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 = 61038330264.452805
c1122 = 28341243079.475285
c1212 = 16348543592.488798


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

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

Gallery generated by Sphinx-Gallery