Note
Go to the end to download the full example code.
Homog5#
Conduct 3d homogenization on a periodic mesh generated with microgen.
Warning: The file contains tag data that couldn't be processed.
Successfully imported the mesh in EasyFEA.
Element type: TETRA4
Ne = 11112, Nn = 3826
TETRA4 -> Ne = 11112
C_hom =
7.977e-02 3.394e-02 3.390e-02 2.593e-05 5.758e-06 1.544e-05
3.394e-02 7.973e-02 3.389e-02 4.211e-05 3.200e-05 7.147e-05
3.390e-02 3.389e-02 7.977e-02 -8.281e-06 2.238e-05 1.388e-05
2.593e-05 4.211e-05 -8.281e-06 6.687e-02 -2.169e-05 -2.921e-05
5.758e-06 3.200e-05 2.238e-05 -2.169e-05 6.698e-02 -6.466e-05
1.544e-05 7.147e-05 1.388e-05 -2.921e-05 -6.466e-05 6.694e-02
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Folder, Models, Simulations, MeshIO, PyVista
17 from EasyFEA.FEM import FeArray
18
19 from Homog4 import Compute_ukl, Get_nodes, Get_pairedNodes
20
21 if __name__ == "__main__":
22 Display.Clear()
23
24 # ----------------------------------------------
25 # Configuration
26 # ----------------------------------------------
27
28 # use Periodic boundary conditions ?
29 usePBC = True
30 plotPBC = False
31 plotSurfaces = False
32
33 folderResults = Folder.Results_Dir()
34 meshes_dir = Folder.Join(Folder.Dir(n=2), "_meshes")
35
36 # ----------------------------------------------
37 # Mesh
38 # ----------------------------------------------
39
40 gmshFile = Folder.Join(meshes_dir, "octet_truss.msh")
41 mesh = MeshIO.Gmsh_to_EasyFEA(gmshFile)
42 mesh.Translate(*-mesh.center) # center mesh on 0,0,0
43
44 plotter = PyVista.Plot_Mesh(mesh)
45 plotter.show_grid()
46 plotter.add_title("RVE")
47 plotter.show()
48
49 # ----------------------------------------------
50 # Get paired nodes
51 # ----------------------------------------------
52
53 tuple_nodes = Get_nodes(mesh, plotSurfaces=plotSurfaces)
54 if usePBC:
55 nodesKUBC = None
56 pairedNodes = Get_pairedNodes(mesh, *tuple_nodes, plotPBC=plotPBC)
57 else:
58 nodesKUBC = set(np.concatenate(tuple_nodes))
59 nodesKUBC = list(nodesKUBC)
60 pairedNodes = None
61
62 # ----------------------------------------------
63 # Material and Simulation
64 # ----------------------------------------------
65 material = Models.Elastic.Isotropic(3, E=1, v=0.3)
66
67 simu = Simulations.Elastic(mesh, material)
68
69 # ----------------------------------------------
70 # Homogenization
71 # ----------------------------------------------
72 r2 = np.sqrt(2)
73 E1 = np.array(
74 [
75 [1, 0, 0],
76 [0, 0, 0],
77 [0, 0, 0],
78 ]
79 )
80 E2 = np.array(
81 [
82 [0, 0, 0],
83 [0, 1, 0],
84 [0, 0, 0],
85 ]
86 )
87 E3 = np.array(
88 [
89 [0, 0, 0],
90 [0, 0, 0],
91 [0, 0, 1],
92 ]
93 )
94 E12 = np.array(
95 [
96 [0, 1 / r2, 0],
97 [1 / r2, 0, 0],
98 [0, 0, 0],
99 ]
100 )
101 E13 = np.array(
102 [
103 [0, 0, 1 / r2],
104 [0, 0, 0],
105 [1 / r2, 0, 0],
106 ]
107 )
108 E23 = np.array(
109 [
110 [0, 0, 0],
111 [0, 0, 1 / r2],
112 [0, 1 / r2, 0],
113 ]
114 )
115
116 u11 = Compute_ukl(simu, E1, nodesKUBC, pairedNodes, True)
117 u22 = Compute_ukl(simu, E2, nodesKUBC, pairedNodes)
118 u33 = Compute_ukl(simu, E3, nodesKUBC, pairedNodes)
119 u12 = Compute_ukl(simu, E12, nodesKUBC, pairedNodes, True)
120 u13 = Compute_ukl(simu, E13, nodesKUBC, pairedNodes)
121 u23 = Compute_ukl(simu, E23, nodesKUBC, pairedNodes)
122
123 u11_e = mesh.Locates_sol_e(u11, asFeArray=True)
124 u22_e = mesh.Locates_sol_e(u22, asFeArray=True)
125 u33_e = mesh.Locates_sol_e(u33, asFeArray=True)
126 u12_e = mesh.Locates_sol_e(u12, asFeArray=True)
127 u13_e = mesh.Locates_sol_e(u13, asFeArray=True)
128 u23_e = mesh.Locates_sol_e(u23, asFeArray=True)
129
130 # ----------------------------------------------
131 # Effective elasticity tensor (C_hom)
132 # ----------------------------------------------
133 U_e = FeArray.zeros(*u11_e.shape, 6)
134
135 U_e[..., 0] = u11_e
136 U_e[..., 1] = u22_e
137 U_e[..., 2] = u33_e
138 U_e[..., 3] = u23_e
139 U_e[..., 4] = u13_e
140 U_e[..., 5] = u12_e
141
142 matrixType = "mass"
143 wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType)
144 B_e_pg = mesh.Get_B_e_pg(matrixType)
145
146 C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2])
147
148 xMin, yMin, zMin = mesh.coord.min(axis=0)
149 xMax, yMax, zMax = mesh.coord.max(axis=0)
150 volume = (xMax - xMin) * (yMax - yMin) * (zMax - zMin)
151
152 C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / volume
153
154 formatted_array = ""
155 for i in range(6):
156 formatted_array += "\n"
157 for j in range(6):
158 formatted_array += f"{C_hom[i,j]:10.3e} "
159
160 print("C_hom =", formatted_array)
161
162 plt.show()
Total running time of the script: (0 minutes 9.578 seconds)