Note
Go to the end to download the full example code.
PlateWithHole#
Damage simulation for a plate with a hole subjected to compression.
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh
2: 0.000μm, [0.00e+00;0.00e+00], 1:1.178s, tol=0.00e+00
3: 0.800μm, [0.00e+00;0.00e+00], 1:1.120s, tol=1.00e+00 3.20% -> 1.13m
4: 1.600μm, [0.00e+00;0.00e+00], 1:1.145s, tol=7.50e-01 6.40% -> 50.26s
5: 2.400μm, [0.00e+00;0.00e+00], 1:1.121s, tol=5.56e-01 9.60% -> 42.22s
6: 3.200μm, [0.00e+00;0.00e+00], 1:1.122s, tol=4.37e-01 12.80% -> 38.21s
7: 4.000μm, [0.00e+00;0.00e+00], 1:1.119s, tol=3.60e-01 16.00% -> 35.24s
8: 4.800μm, [0.00e+00;0.00e+00], 1:1.112s, tol=3.06e-01 19.20% -> 32.74s
9: 5.600μm, [0.00e+00;0.00e+00], 1:1.129s, tol=2.65e-01 22.40% -> 31.28s
10: 6.400μm, [0.00e+00;0.00e+00], 1:1.117s, tol=2.34e-01 25.60% -> 29.21s
11: 7.200μm, [0.00e+00;0.00e+00], 1:1.100s, tol=2.10e-01 28.80% -> 27.20s
12: 8.000μm, [0.00e+00;0.00e+00], 1:1.098s, tol=1.90e-01 32.00% -> 25.66s
13: 8.800μm, [0.00e+00;0.00e+00], 1:1.090s, tol=1.74e-01 35.20% -> 24.07s
14: 9.600μm, [0.00e+00;0.00e+00], 1:1.093s, tol=1.60e-01 38.40% -> 22.79s
15: 10.400μm, [0.00e+00;0.00e+00], 1:1.100s, tol=1.48e-01 41.60% -> 21.62s
16: 11.200μm, [0.00e+00;0.00e+00], 1:1.104s, tol=1.38e-01 44.80% -> 20.40s
17: 12.000μm, [0.00e+00;0.00e+00], 1:1.132s, tol=1.29e-01 48.00% -> 19.63s
18: 12.800μm, [0.00e+00;0.00e+00], 1:1.112s, tol=1.21e-01 51.20% -> 18.02s
19: 13.600μm, [0.00e+00;0.00e+00], 1:1.113s, tol=1.14e-01 54.40% -> 16.79s
20: 14.400μm, [0.00e+00;0.00e+00], 1:1.113s, tol=1.08e-01 57.60% -> 15.56s
21: 15.200μm, [0.00e+00;0.00e+00], 1:1.121s, tol=1.02e-01 60.80% -> 14.46s
22: 16.000μm, [0.00e+00;0.00e+00], 1:1.114s, tol=9.75e-02 64.00% -> 13.16s
23: 16.800μm, [9.79e-13;4.17e-03], 1:1.145s, tol=9.30e-02 67.20% -> 12.30s
24: 17.600μm, [1.02e-11;4.36e-02], 1:1.139s, tol=8.88e-02 70.40% -> 11.02s
25: 18.400μm, [1.09e-11;1.18e-01], 1:1.117s, tol=8.49e-02 73.60% -> 9.62s
26: 19.200μm, [7.24e-12;2.14e-01], 1:1.137s, tol=8.14e-02 76.80% -> 8.58s
27: 20.000μm, [4.04e-12;3.31e-01], 1:1.134s, tol=7.82e-02 80.00% -> 7.37s
28: 20.800μm, [2.70e-12;5.03e-01], 1:1.133s, tol=7.50e-02 83.20% -> 6.18s
29: 21.600μm, [1.48e-12;7.38e-01], 1:1.129s, tol=7.24e-02 86.40% -> 4.97s
30: 21.800μm, [1.01e-12;9.34e-01], 1:1.109s, tol=2.12e-02 87.20% -> 4.72s
31: 22.000μm, [1.11e-12;9.98e-01], 1:1.109s, tol=3.65e-02 88.00% -> 4.54s
32: 22.200μm, [1.22e-12;1.01e+00], 1:1.104s, tol=6.20e-02 88.80% -> 4.32s
33: 22.400μm, [1.36e-12;1.00e+00], 1:1.105s, tol=8.53e-02 89.60% -> 4.10s
34: 22.600μm, [1.35e-12;1.01e+00], 1:1.096s, tol=9.21e-02 90.40% -> 3.84s
35: 22.800μm, [1.44e-12;1.00e+00], 1:1.101s, tol=1.06e-01 91.20% -> 3.61s
36: 23.000μm, [1.49e-12;1.00e+00], 1:1.102s, tol=1.02e-01 92.00% -> 3.35s
37: 23.200μm, [1.63e-12;1.00e+00], 1:1.104s, tol=9.02e-02 92.80% -> 3.08s
38: 23.400μm, [1.67e-12;1.00e+00], 1:1.122s, tol=9.63e-02 93.60% -> 2.84s
39: 23.600μm, [1.75e-12;1.00e+00], 1:1.122s, tol=1.02e-01 94.40% -> 2.53s
40: 23.800μm, [1.79e-12;1.00e+00], 1:1.130s, tol=1.03e-01 95.20% -> 2.22s
41: 24.000μm, [1.68e-12;1.00e+00], 1:1.130s, tol=1.01e-01 96.00% -> 1.88s
42: 24.200μm, [1.43e-12;1.00e+00], 1:1.144s, tol=9.89e-02 96.80% -> 1.55s
43: 24.400μm, [1.26e-12;1.00e+00], 1:1.139s, tol=1.01e-01 97.60% -> 1.18s
44: 24.600μm, [1.10e-12;1.00e+00], 1:1.133s, tol=9.77e-02 98.40% -> 792.50ms
45: 24.800μm, [1.02e-12;1.00e+00], 1:1.125s, tol=9.67e-02 99.20% -> 399.23ms
46: 25.000μm, [9.35e-13;1.00e+00], 1:1.112s, tol=9.35e-02 100.00% -> -0.00µs
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh
Generate movie 01/23 (4.35 %) 3.62 s
Generate movie 02/23 (8.70 %) 3.34 s
Generate movie 03/23 (13.04 %) 3.16 s
Generate movie 04/23 (17.39 %) 3.01 s
Generate movie 05/23 (21.74 %) 2.84 s
Generate movie 06/23 (26.09 %) 2.69 s
Generate movie 07/23 (30.43 %) 2.53 s
Generate movie 08/23 (34.78 %) 2.36 s
Generate movie 09/23 (39.13 %) 2.22 s
Generate movie 10/23 (43.48 %) 2.05 s
Generate movie 11/23 (47.83 %) 1.90 s
Generate movie 12/23 (52.17 %) 1.76 s
Generate movie 13/23 (56.52 %) 1.61 s
Generate movie 14/23 (60.87 %) 1.44 s
Generate movie 15/23 (65.22 %) 1.28 s
Generate movie 16/23 (69.57 %) 1.11 s
Generate movie 17/23 (73.91 %) 950.40 ms
Generate movie 18/23 (78.26 %) 786.54 ms
Generate movie 19/23 (82.61 %) 639.74 ms
Generate movie 20/23 (86.96 %) 479.39 ms
Generate movie 21/23 (91.30 %) 316.77 ms
Generate movie 22/23 (95.65 %) 157.97 ms
Generate movie 23/23 (100.00 %) 0.00 µs
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import (
17 Display,
18 Folder,
19 Models,
20 Tic,
21 ElemType,
22 Simulations,
23 PyVista,
24 Paraview,
25 )
26 from EasyFEA.Geoms import Domain, Circle
27
28 if __name__ == "__main__":
29
30 # ----------------------------------------------
31 # Configuration
32 # ----------------------------------------------
33
34 # simu options
35 doSimu = True
36 meshTest = True
37 optimMesh = True
38
39 # outputs
40 folder = Folder.Results_Dir()
41 plotMesh = False
42 plotIter = False
43 plotEnergy = False
44
45 makeParaview = False
46 makeMovie = True
47
48 # Available splits: Bourdin, Amor, Miehe, Stress (isotropic)
49 # He, AnisotStrain, AnisotStress, Zhang (anisotropic)
50 split = Models.PhaseField.SplitType.Miehe
51
52 # Available regus: AT1, AT2
53 regu = Models.PhaseField.ReguType.AT1
54
55 l0 = 0.12e-3 # m
56
57 # convergence
58 solver = (
59 Models.PhaseField.SolverType.History
60 ) # History, HistoryDamage, BoundConstrain
61 maxIter = 1000
62 tolConv = 1e-0
63
64 # load units
65 unitU = "μm"
66 unitF = "kN/mm"
67 unit = 1e6
68
69 # ----------------------------------------------
70 # Geometry
71 # ----------------------------------------------
72 L = 15e-3 # m
73 h = 30e-3 # m
74 diam = 6e-3 # m
75 thickness = 1
76
77 # ----------------------------------------------
78 # Material
79 # ----------------------------------------------
80 E = 12e9 # Pa
81 v = 0.3
82 Gc = 1.4 # J/m2
83
84 folder_save = Simulations.PhaseField.Folder(
85 folder, "", split, regu, "", tolConv, solver, meshTest, optimMesh
86 )
87 Display.MyPrint(folder_save, "green", end="\n")
88
89 if doSimu:
90 # ----------------------------------------------
91 # Mesh
92 # ----------------------------------------------
93 clC = l0 * 2 if meshTest else l0 / 2
94 if optimMesh:
95 clD = l0 * 4
96 refineZone = diam * 1.5 / 2
97 if split in (
98 Models.PhaseField.SplitType.Bourdin,
99 Models.PhaseField.SplitType.Amor,
100 ):
101 refineGeom = Domain(
102 (0, h / 2 - refineZone), (L, h / 2 + refineZone), clC
103 )
104 else:
105 refineGeom = Domain(
106 (L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC
107 )
108 else:
109 clD = l0 if meshTest else l0 / 2
110 refineGeom = None
111
112 domain = Domain((0, 0), (L, h), clD)
113 circle = Circle((L / 2, h / 2), diam, clD)
114 mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom])
115
116 # Nodes
117 nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
118 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
119 nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0))
120 nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"])
121
122 # ----------------------------------------------
123 # Material
124 # ----------------------------------------------
125 material = Models.Elastic.Isotropic(
126 2, E, v, planeStress=False, thickness=thickness
127 )
128 pfm = Models.PhaseField(material, split, regu, Gc, l0, solver=solver)
129
130 # ----------------------------------------------
131 # Boundary conditions
132 # ----------------------------------------------
133 threshold = 0.6
134 u_max = 25e-6
135 uinc0 = 8e-7 if meshTest else 8e-8
136 uinc1 = 2e-7 if meshTest else 2e-8
137
138 config = f"""
139 uinc0 = {uinc0:.1e}; uinc1 = {uinc1:.1e}
140
141 while ud <= u_max = {u_max:.2e}:
142
143 ud += uinc0 if simu.damage.max() < {threshold} else uinc1
144
145 simu.add_dirichlet(nodes_lower, [0], ["y"])
146 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
147 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
148 """
149
150 # ----------------------------------------------
151 # Simulation
152 # ----------------------------------------------
153 simu = Simulations.PhaseField(mesh, pfm, folder=folder_save)
154 simu.Results_Set_Bc_Summary(config)
155
156 dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"])
157
158 def Apply_BC(ud: float):
159 simu.Bc_Init()
160 simu.add_dirichlet(nodes_lower, [0], ["y"])
161 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
162 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
163
164 list_dep = []
165 list_f = []
166 ud = -uinc0
167 iter = 0
168 nDetect = 0
169
170 if plotIter:
171 axIter = Display.Plot_Result(simu, "damage", nodeValues=True)
172 _, axLoad = plt.subplots()
173 axLoad.set_xlabel(f"ud [{unitU}]")
174 axLoad.set_ylabel(f"f [{unitF}]")
175 axLoad.grid()
176
177 while ud <= u_max:
178 iter += 1
179 ud += uinc0 if simu.damage.max() < threshold else uinc1
180
181 Apply_BC(ud)
182
183 u, d, convergence = simu.Solve(tolConv, maxIter)
184 simu.Save_Iter()
185
186 if not convergence:
187 break
188
189 f = np.sum(simu.Calc_Reaction(dofsY_upper, "elastic"))
190 simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True)
191
192 if simu.Detect_Damage(nodes_edges, 1):
193 nDetect += 1
194 if nDetect == 10:
195 break
196
197 list_dep.append(ud)
198 list_f.append(f)
199
200 if plotIter:
201 Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter)
202 plt.figure(axIter.figure)
203 plt.pause(1e-12)
204 axLoad.clear()
205 axLoad.plot(np.abs(list_dep * unit), np.abs(list_f / unit), c="blue")
206 axLoad.set_xlabel(f"ud [{unitU}]")
207 axLoad.set_ylabel(f"f [{unitF}]")
208 axLoad.grid()
209 plt.figure(axLoad.figure)
210 plt.pause(1e-12)
211
212 # ----------------------------------------------
213 # Saving
214 # ----------------------------------------------
215 print()
216 Simulations.Save_pickle((list_f, list_dep), folder_save, "force-displacement")
217 simu.Save(folder_save)
218
219 else:
220 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
221 list_f, list_dep = Simulations.Load_pickle(folder_save, "force-displacement")
222
223 # ----------------------------------------------
224 # Results
225 # ----------------------------------------------
226 if plotEnergy:
227 Display.Plot_Energy(simu, list_f, list_dep, N=400, folder=folder_save)
228
229 Display.Plot_Result(
230 simu,
231 "damage",
232 nodeValues=True,
233 colorbarIsClose=True,
234 folder=folder_save,
235 filename="damage",
236 )
237 Display.Plot_Mesh(simu)
238 Display.Plot_BoundaryConditions(simu)
239 Display.Plot_Iter_Summary(simu, folder_save, None, None)
240
241 ax = Display.Init_Axes()
242 ax.plot(np.abs(list_dep) * unit, np.abs(list_f) / unit, c="blue")
243 ax.set_xlabel(f"ud [{unitU}]")
244 ax.set_ylabel(f"f [{unitF}]")
245 ax.grid()
246 Display.Save_fig(folder_save, "force-displacement")
247
248 if plotMesh:
249 Display.Plot_Mesh(simu.mesh)
250
251 if makeParaview:
252 Paraview.Save_simu(simu, folder_save)
253
254 if makeMovie:
255 simu.Set_Iter(-1)
256 deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
257
258 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
259
260 def Func(plotter, iter):
261 simu.Set_Iter(iterations[iter])
262 thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
263 PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
264
265 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
266
267 if doSimu:
268 Tic.Plot_History(folder_save, details=False)
269
270 plt.show()
Total running time of the script: (0 minutes 57.991 seconds)