Note

Go to the end to download the full example code.

Homog2#

Perform homogenization on several RVE.

  • RVE
  • Periodic boundary conditions
  • E [GPa]
  • v
  • $ux$
  • $uy$
  • $\sigma_{xx}$
  • $\sigma_{yy}$
  • $\sigma_{xy}$
c1111 = 35530825959.242096
c1122 = 16497629783.81898
c1212 = 9516598087.711594


 13 import matplotlib.pyplot as plt
 14 import numpy as np
 15
 16 from EasyFEA import Display, Models, ElemType, Simulations
 17 from EasyFEA.Geoms import Point, Points, Line
 18 from EasyFEA.FEM import FeArray
 19
 20 from Homog1 import Compute_ukl
 21
 22 if __name__ == "__main__":
 23     Display.Clear()
 24
 25     # ----------------------------------------------
 26     # Configuration
 27     # ----------------------------------------------
 28
 29     # use Periodic boundary conditions ?
 30     usePBC = True
 31
 32     geom = "D666"  # hexagon
 33     # geom = "D2"  # rectangle
 34     # geom = "D6"
 35
 36     hollowInclusion = True
 37
 38     # ----------------------------------------------
 39     # Mesh
 40     # ----------------------------------------------
 41     N = 5
 42     elemType = ElemType.TRI6
 43
 44     if geom == "D666":
 45         a = 1
 46         R = 2 * a / np.sqrt(3)
 47         r = R / np.sqrt(2) / 2
 48         phi = np.pi / 6
 49
 50         cos_phi = np.cos(phi)
 51         sin_phi = np.sin(phi)
 52
 53         # Create the contour geometrie
 54         p0 = Point(0, R)
 55         p1 = Point(-cos_phi * R, sin_phi * R)
 56         p2 = Point(-cos_phi * R, -sin_phi * R)
 57         p3 = Point(0, -R)
 58         p4 = Point(cos_phi * R, -sin_phi * R)
 59         p5 = Point(cos_phi * R, sin_phi * R)
 60         # edge length and area
 61         s = Line(p0, p1).length
 62         area = 3 * np.sqrt(3) / 2 * s**2
 63
 64         contour = Points([p0, p1, p2, p3, p4, p5], s / N)
 65         corners = contour.points
 66
 67         # Create the inclusion
 68         p6 = Point(0, (R - r))
 69         p7 = Point(-cos_phi * (R - r), sin_phi * (R - r))
 70         p8 = Point(-cos_phi * (R - r), -sin_phi * (R - r))
 71         p9 = Point(0, -(R - r))
 72         p10 = Point(cos_phi * (R - r), -sin_phi * (R - r))
 73         p11 = Point(cos_phi * (R - r), sin_phi * (R - r))
 74         inclusions = [Points([p6, p7, p8, p9, p10, p11], s / N, not hollowInclusion)]
 75
 76     elif geom == "D2":
 77         a = 1  # width
 78         b = 1.4  # height
 79         e = 1 / 10  # thickness
 80         area = a * b
 81         meshSize = e / N * 2
 82
 83         # Create the contour geometry
 84         p0 = Point(-a / 2, b / 2)
 85         p1 = Point(-a / 2, -b / 2)
 86         p2 = Point(a / 2, -b / 2)
 87         p3 = Point(a / 2, b / 2)
 88         contour = Points([p0, p1, p2, p3], meshSize)
 89         corners = contour.points
 90
 91         # Create the inclusion geometry
 92         p4 = p0 + [e, -e]
 93         p5 = p1 + [e, e]
 94         p6 = p2 + [-e, e]
 95         p7 = p3 + [-e, -e]
 96         inclusions = [Points([p4, p5, p6, p7], meshSize, hollowInclusion)]
 97
 98     elif geom == "D6":
 99         a = 1  # height
100         b = 2  # width
101         c = np.sqrt(a**2 + b**2)
102
103         e = b / 10  # thickness
104         l1 = b / 2
105
106         area = a * b
107
108         theta = np.arctan(a / b)
109         alpha = (np.pi / 2 - theta) / 2
110         cos_alpha = np.cos(alpha)
111         sin_alpha = np.sin(alpha)
112         phi = np.pi / 3
113         cos_phi = np.cos(phi)
114         sin_phi = np.sin(phi)
115
116         l2 = (b - l1 * sin_alpha) / 2
117         hx = e / cos_phi / 4
118         hy = e / sin_phi / 4
119
120         # symmetry functions
121         def Sym_x(point: Point) -> Point:
122             return Point(-point.x, point.y)
123
124         def Sym_y(point: Point) -> Point:
125             return Point(point.x, -point.y)
126
127         # points in the non-rotated base
128         p0 = Point(
129             l1 / 2 + l2 * cos_phi + e / 2 * cos_alpha, l2 * sin_phi - e / 2 * sin_alpha
130         )
131         p1 = p0 + [-e * cos_alpha, e * sin_alpha]
132         p2 = Point(l1 / 2 - hy, hx)
133         p3 = Sym_x(p2)
134         p4 = Sym_x(p1)
135         p5 = Sym_x(p0)
136         p6 = Point(-l1 / 2 - np.sqrt(hx**2 + hy**2))
137         p7 = Sym_y(p5)
138         p8 = Sym_y(p4)
139         p9 = Sym_y(p3)
140         p10 = Sym_y(p2)
141         p11 = Sym_y(p1)
142         p12 = Sym_y(p0)
143         p13 = Sym_x(p6)
144
145         # do some tests to check if the geometry has been created correctly
146         t1 = Line(p2, p10).length
147         t2 = Line(p2, p13).length
148         t3 = Line(p10, p13).length
149         assert (
150             np.abs(e - (t1 + t2 + t3) / 3) / e <= 1e-12
151         )  # check that t1 = t2 = t3 = e
152         t4 = Line(p0, p1).length
153         assert np.abs(t4 - e) / e <= 1e-12  # check that t4 = e
154
155         alpha = -alpha
156         rot = np.array(
157             [
158                 [np.cos(alpha), -np.sin(alpha), 0],
159                 [np.sin(alpha), np.cos(alpha), 0],
160                 [0, 0, 1],
161             ]
162         )
163
164         rotate_points = []
165         ax = Display.Init_Axes()
166         for p, point in enumerate(
167             [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13]
168         ):
169             assert isinstance(point, Point)
170
171             newCoord = rot @ point.coord
172
173             ax.scatter(*newCoord[:2], c="black")
174             ax.text(*newCoord[:2], f"p{p}", c="black")
175
176             rotate_points.append(Point(*newCoord))
177
178         corners = [rotate_points[p] for p in [0, 1, 4, 5, 7, 8, 11, 12]]
179
180         hollowInclusion = True  # dont change
181
182         contour = Points(rotate_points, e / N * 2)
183
184         inclusions = []
185
186     else:
187         raise Exception("Unknown geom")
188
189     mesh = contour.Mesh_2D(inclusions, elemType)
190
191     Display.Plot_Mesh(mesh, title="RVE")
192     # Display.Plot_Tags(mesh)
193
194     nodes_matrix = mesh.Nodes_Tags(["S0"])
195     elements_matrix = mesh.Elements_Nodes(nodes_matrix)
196
197     if not hollowInclusion:
198         nodes_inclusion = mesh.Nodes_Tags(["S1"])
199         elements_inclusion = mesh.Elements_Nodes(nodes_inclusion)
200
201     nCorners = len(corners)
202     nEdges = nCorners // 2
203
204     if usePBC:
205         nodes_kubc = mesh.Nodes_Point(corners)
206         paired_nodes = mesh.Get_Paired_Nodes(nodes_kubc, True)
207     else:
208         nodes_kubc = mesh.Nodes_Tags([f"L{i}" for i in range(6)])
209         paired_nodes = None
210
211     # ----------------------------------------------
212     # Material and simulation
213     # ----------------------------------------------
214     E = np.ones(mesh.Ne) * 70 * 1e9
215     v = np.ones(mesh.Ne) * 0.45
216
217     if not hollowInclusion:
218         E[elements_inclusion] = 200 * 1e9
219         v[elements_inclusion] = 0.3
220
221     Display.Plot_Result(mesh, E * 1e-9, nodeValues=False, title="E [GPa]")
222     Display.Plot_Result(mesh, v, nodeValues=False, title="v")
223
224     material = Models.Elastic.Isotropic(2, E, v, planeStress=False)
225
226     simu = Simulations.Elastic(mesh, material)
227
228     # ----------------------------------------------
229     # Homogenization
230     # ----------------------------------------------
231     r2 = np.sqrt(2)
232     E11 = np.array([[1, 0], [0, 0]])
233     E22 = np.array([[0, 0], [0, 1]])
234     E12 = np.array([[0, 1 / r2], [1 / r2, 0]])
235
236     u11 = Compute_ukl(simu, nodes_kubc, E11, paired_nodes)
237     u22 = Compute_ukl(simu, nodes_kubc, E22, paired_nodes)
238     u12 = Compute_ukl(simu, nodes_kubc, E12, paired_nodes, True)
239
240     u11_e = mesh.Locates_sol_e(u11, asFeArray=True)
241     u22_e = mesh.Locates_sol_e(u22, asFeArray=True)
242     u12_e = mesh.Locates_sol_e(u12, asFeArray=True)
243
244     # ----------------------------------------------
245     # Effective elasticity tensor (C_hom)
246     # ----------------------------------------------
247     U_e = FeArray.zeros(*u11_e.shape, 3)
248
249     U_e[..., 0] = u11_e
250     U_e[..., 1] = u22_e
251     U_e[..., 2] = u12_e
252
253     matrixType = "mass"
254     wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType)
255     B_e_pg = mesh.Get_B_e_pg(matrixType)
256
257     C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2])
258
259     C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / area
260
261     print(f"c1111 = {C_hom[0, 0]}")
262     print(f"c1122 = {C_hom[0, 1]}")
263     print(f"c1212 = {C_hom[2, 2] / 2}")
264
265     plt.show()

Total running time of the script: (0 minutes 0.876 seconds)

Gallery generated by Sphinx-Gallery