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.45283
c1122 = 28341243079.475273
c1212 = 16348543592.488775


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

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

Gallery generated by Sphinx-Gallery