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/envs/v1.5.4/lib/python3.11/site-packages/results/PhaseField/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh 2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:0.160 s, tol=0.00e+00
3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.131 s, tol=1.00e+00 3.20 % -> 7.95 s
4 : 1.600 μm, [3.98e-09; 3.63e-02], 1:0.118 s, tol=8.04e-01 6.40 % -> 5.20 s
5 : 2.400 μm, [3.74e-05; 3.63e-02], 1:0.131 s, tol=5.54e-01 9.60 % -> 4.95 s
6 : 3.200 μm, [8.41e-05; 3.63e-02], 1:0.132 s, tol=4.37e-01 12.80 % -> 4.51 s
7 : 4.000 μm, [1.50e-04; 3.63e-02], 1:0.140 s, tol=3.60e-01 16.00 % -> 4.42 s
8 : 4.800 μm, [2.34e-04; 3.63e-02], 1:0.140 s, tol=3.05e-01 19.20 % -> 4.14 s
9 : 5.600 μm, [3.37e-04; 3.63e-02], 1:0.132 s, tol=2.65e-01 22.40 % -> 3.66 s
10 : 6.400 μm, [4.60e-04; 3.64e-02], 1:0.132 s, tol=2.34e-01 25.60 % -> 3.46 s
11 : 7.200 μm, [6.01e-04; 4.57e-02], 1:0.131 s, tol=2.10e-01 28.80 % -> 3.24 s
12 : 8.000 μm, [7.63e-04; 5.81e-02], 1:0.132 s, tol=1.90e-01 32.00 % -> 3.07 s
13 : 8.800 μm, [9.46e-04; 7.21e-02], 1:0.131 s, tol=1.73e-01 35.20 % -> 2.90 s
14 : 9.600 μm, [1.15e-03; 8.78e-02], 1:0.132 s, tol=1.59e-01 38.40 % -> 2.74 s
15 : 10.400 μm, [1.38e-03; 1.05e-01], 1:0.120 s, tol=1.47e-01 41.60 % -> 2.36 s
16 : 11.200 μm, [1.62e-03; 1.25e-01], 1:0.131 s, tol=1.37e-01 44.80 % -> 2.43 s
17 : 12.000 μm, [1.89e-03; 1.46e-01], 1:0.132 s, tol=1.28e-01 48.00 % -> 2.28 s
18 : 12.800 μm, [2.19e-03; 1.70e-01], 1:0.133 s, tol=1.21e-01 51.20 % -> 2.15 s
19 : 13.600 μm, [2.50e-03; 1.97e-01], 1:0.132 s, tol=1.14e-01 54.40 % -> 1.98 s
20 : 14.400 μm, [2.85e-03; 2.29e-01], 1:0.132 s, tol=1.07e-01 57.60 % -> 1.84 s
21 : 15.200 μm, [3.21e-03; 2.65e-01], 1:0.142 s, tol=1.02e-01 60.80 % -> 1.83 s
22 : 16.000 μm, [3.61e-03; 3.04e-01], 1:0.132 s, tol=9.68e-02 64.00 % -> 1.56 s
23 : 16.800 μm, [4.03e-03; 3.48e-01], 1:0.131 s, tol=9.23e-02 67.20 % -> 1.41 s
24 : 17.600 μm, [4.47e-03; 3.96e-01], 1:0.133 s, tol=8.82e-02 70.40 % -> 1.29 s
25 : 18.400 μm, [4.95e-03; 4.49e-01], 1:0.132 s, tol=8.44e-02 73.60 % -> 1.14 s
26 : 19.200 μm, [5.45e-03; 5.05e-01], 1:0.131 s, tol=8.10e-02 76.80 % -> 992.88 ms
27 : 20.000 μm, [5.98e-03; 5.61e-01], 1:0.132 s, tol=7.80e-02 80.00 % -> 854.91 ms
28 : 20.800 μm, [6.53e-03; 6.38e-01], 1:0.132 s, tol=7.53e-02 83.20 % -> 716.97 ms
29 : 21.000 μm, [7.11e-03; 7.36e-01], 1:0.132 s, tol=2.04e-02 84.00 % -> 701.39 ms
30 : 21.200 μm, [7.30e-03; 8.35e-01], 1:0.132 s, tol=2.10e-02 84.80 % -> 684.91 ms
31 : 21.400 μm, [7.46e-03; 9.20e-01], 1:0.145 s, tol=2.34e-02 85.60 % -> 733.54 ms
32 : 21.600 μm, [7.64e-03; 9.75e-01], 1:0.132 s, tol=2.99e-02 86.40 % -> 642.21 ms
33 : 21.800 μm, [7.87e-03; 9.98e-01], 1:0.133 s, tol=4.11e-02 87.20 % -> 625.53 ms
34 : 22.000 μm, [8.18e-03; 1.00e+00], 1:0.132 s, tol=5.57e-02 88.00 % -> 593.28 ms
35 : 22.200 μm, [8.57e-03; 1.00e+00], 1:0.133 s, tol=6.35e-02 88.80 % -> 570.60 ms
36 : 22.400 μm, [8.97e-03; 1.00e+00], 1:0.133 s, tol=8.48e-02 89.60 % -> 540.98 ms
37 : 22.600 μm, [9.40e-03; 1.00e+00], 1:0.132 s, tol=8.07e-02 90.40 % -> 503.66 ms
38 : 22.800 μm, [9.82e-03; 1.00e+00], 1:0.132 s, tol=7.54e-02 91.20 % -> 470.65 ms
39 : 23.000 μm, [1.01e-02; 1.00e+00], 1:0.132 s, tol=7.34e-02 92.00 % -> 437.56 ms
40 : 23.200 μm, [1.04e-02; 1.00e+00], 1:0.121 s, tol=7.26e-02 92.80 % -> 367.54 ms
41 : 23.400 μm, [1.06e-02; 1.00e+00], 1:0.140 s, tol=7.46e-02 93.60 % -> 383.77 ms
42 : 23.600 μm, [1.08e-02; 1.00e+00], 1:0.132 s, tol=7.51e-02 94.40 % -> 320.31 ms
43 : 23.800 μm, [1.10e-02; 1.00e+00], 1:0.132 s, tol=7.56e-02 95.20 % -> 279.29 ms
44 : 24.000 μm, [1.11e-02; 1.00e+00], 1:0.135 s, tol=7.53e-02 96.00 % -> 242.61 ms
45 : 24.200 μm, [1.12e-02; 1.00e+00], 1:0.132 s, tol=7.52e-02 96.80 % -> 191.71 ms
46 : 24.400 μm, [1.14e-02; 1.00e+00], 1:0.121 s, tol=7.46e-02 97.60 % -> 133.52 ms
47 : 24.600 μm, [1.15e-02; 1.00e+00], 1:0.132 s, tol=7.42e-02 98.40 % -> 98.63 ms
48 : 24.800 μm, [1.16e-02; 1.00e+00], 1:0.135 s, tol=7.27e-02 99.20 % -> 51.28 ms
49 : 25.000 μm, [1.17e-02; 1.00e+00], 1:0.121 s, tol=7.13e-02 100.00 % -> -0.00 µs
Saved:
/results/PhaseField/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/force-displacement.pickle
Saved:
/results/PhaseField/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/simulation.pickle
Saved:
/results/PhaseField/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/summary.txt
Movie_func 0/23
Movie_func 1/23 (4.35 %) 3.11 s
Movie_func 2/23 (8.70 %) 2.95 s
Movie_func 3/23 (13.04 %) 2.77 s
Movie_func 4/23 (17.39 %) 2.63 s
Movie_func 5/23 (21.74 %) 2.50 s
Movie_func 6/23 (26.09 %) 2.35 s
Movie_func 7/23 (30.43 %) 2.24 s
Movie_func 8/23 (34.78 %) 2.09 s
Movie_func 9/23 (39.13 %) 1.97 s
Movie_func 10/23 (43.48 %) 1.86 s
Movie_func 11/23 (47.83 %) 1.72 s
Movie_func 12/23 (52.17 %) 1.54 s
Movie_func 13/23 (56.52 %) 1.42 s
Movie_func 14/23 (60.87 %) 1.28 s
Movie_func 15/23 (65.22 %) 1.12 s
Movie_func 16/23 (69.57 %) 986.97 ms
Movie_func 17/23 (73.91 %) 845.75 ms
Movie_func 18/23 (78.26 %) 706.71 ms
Movie_func 19/23 (82.61 %) 568.21 ms
Movie_func 20/23 (86.96 %) 427.83 ms
Movie_func 21/23 (91.30 %) 282.40 ms
Movie_func 22/23 (95.65 %) 145.32 ms
Movie_func 23/23 (100.00 %) 0.00 µs
12 from EasyFEA import (
13 Display,
14 Folder,
15 Models,
16 plt,
17 np,
18 Tic,
19 ElemType,
20 Mesh,
21 Simulations,
22 PyVista,
23 Paraview,
24 )
25 from EasyFEA.Geoms import Domain, Circle
26
27 import multiprocessing
28
29 # Display.Clear()
30
31
32 # ----------------------------------------------
33 # Configuration
34 # ----------------------------------------------
35
36 # simu options
37 doSimu = True
38 meshTest = True
39 optimMesh = True
40 useParallel = False
41 nProcs = 4 # number of processes in parallel
42
43 # outputs
44 plotMesh = False
45 plotIter = False
46 plotResult = True
47 plotEnergy = False
48 showFig = True
49
50 saveParaview = False
51 makeMovie = True
52
53 # models
54
55 # splits = ["Bourdin","Amor","Miehe","Stress"] # Splits Isotropes
56 # splits = ["He","AnisotStrain","AnisotStress","Zhang"] # Splits Anisotropes
57 # splits = ["Bourdin","Amor","Miehe","Stress","He","AnisotStrain","AnisotStress","Zhang"]
58 # splits = ["Zhang"]
59 # splits = ["AnisotStrain","AnisotStress","Zhang"]
60 splits = ["Miehe"]
61
62 regus = ["AT2"] # ["AT1", "AT2"]
63 # regus = ["AT1", "AT2"]
64
65 l0 = 0.12e-3
66
67 # convergence
68 solver = (
69 Models.PhaseField.SolverType.History
70 ) # ["History", "HistoryDamage", "BoundConstrain"]
71 maxIter = 1000
72 tolConv = 1e-0
73
74 # ----------------------------------------------
75 # Mesh
76 # ----------------------------------------------
77
78
79 def DoMesh(
80 L: float, h: float, diam: float, thickness: float, l0: float, split: str
81 ) -> Mesh:
82 clC = l0 * 2 if meshTest else l0 / 2
83 if optimMesh:
84 clD = l0 * 4
85 refineZone = diam * 1.5 / 2
86 if split in ["Bourdin", "Amor"]:
87 refineGeom = Domain((0, h / 2 - refineZone), (L, h / 2 + refineZone), clC)
88 else:
89 refineGeom = Domain((L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC)
90 else:
91 clD = l0 if meshTest else l0 / 2
92 refineGeom = None
93
94 domain = Domain((0, 0), (L, h), clD)
95 circle = Circle((L / 2, h / 2), diam, clD, isHollow=True)
96
97 # ax = Display.Init_Axes()
98 # domain.Plot(ax, color="k", plotPoints=False)
99 # circle.Plot(ax, color="k", plotPoints=False)
100 # # if refineGeom != None:
101 # # refineGeom.Plot(ax, color='k', plotPoints=False)
102 # # ax.scatter(((L+diam)/2, L/2), (h/2, (h+diam)/2), c='k')
103 # ax.axis("off")
104 # Display.Save_fig(Folder.RESULTS_DIR, "sample", True)
105
106 mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom])
107
108 # ax = Display.Plot_Mesh(mesh, lw=0.3, facecolors="white")
109 # ax.axis("off")
110 # ax.set_title("")
111 # Display.Save_fig(Folder.RESULTS_DIR, "mesh", transparent=True)
112
113 return mesh
114
115
116 # ----------------------------------------------
117 # Do Simu
118 # ----------------------------------------------
119 def DoSimu(split: str, regu: str):
120 folder_save = Folder.PhaseField_Folder(
121 "PlateWithHole",
122 "Elas_Isot",
123 split,
124 regu,
125 "DP",
126 tolConv,
127 solver,
128 meshTest,
129 optimMesh,
130 )
131
132 Display.MyPrint(folder_save, "green")
133
134 # ----------------------------------------------
135 # Geom
136 # ----------------------------------------------
137
138 L = 15e-3
139 h = 30e-3
140 thickness = 1
141 diam = 6e-3
142
143 # load units
144 unitU = "μm"
145 unitF = "kN/mm"
146 unit = 1e6
147
148 # ----------------------------------------------
149 # Mesh
150 # ----------------------------------------------
151
152 if doSimu:
153 mesh = DoMesh(L, h, diam, thickness, l0, split)
154
155 # Get Nodes
156 nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
157 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
158 nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0))
159 nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"])
160 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
161
162 # ----------------------------------------------
163 # Boundary conditions
164 # ----------------------------------------------
165
166 threshold = 0.6
167 u_max = 25e-6
168 # u_max = 35e-6
169
170 uinc0 = 8e-7 if meshTest else 8e-8
171 uinc1 = 2e-7 if meshTest else 2e-8
172
173 config = f"""
174 while ud <= u_max:
175
176 ud += uinc0 if simu.damage.max() < threshold else uinc1
177
178 u_max = {u_max}
179 uinc0 = {uinc0:.1e} (simu.damage.max() < {threshold})
180 uinc1 = {uinc1:.1e}
181
182 simu.add_dirichlet(nodes_lower, [0], ["y"])
183 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
184 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
185 if dim == 3:
186 simu.add_dirichlet(nodes_y0z0, [0], ["z"])
187 """
188
189 # ----------------------------------------------
190 # Material
191 # ----------------------------------------------
192 E = 12e9
193 v = 0.3
194 planeStress = False
195 material = Models.ElasIsot(2, E, v, planeStress, thickness)
196
197 gc = 1.4
198 pfm = Models.PhaseField(material, split, regu, gc, l0, solver=solver)
199
200 # ----------------------------------------------
201 # Simulation
202 # ----------------------------------------------
203 simu = Simulations.PhaseFieldSimu(mesh, pfm, verbosity=False)
204
205 simu.Results_Set_Bc_Summary(config)
206
207 dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"])
208
209 def Apply_BC(ud: float):
210 simu.Bc_Init()
211 simu.add_dirichlet(nodes_lower, [0], ["y"])
212 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
213 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
214
215 # INIT
216 displacement = []
217 force = []
218 ud = -uinc0
219 iter = 0
220 nDetect = 0
221
222 if plotIter:
223 axIter = Display.Plot_Result(simu, "damage", nodeValues=True)
224
225 force = np.asarray(force)
226 displacement = np.asarray(displacement)
227 _, axLoad = Display.Plot_Force_Displacement(
228 force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]"
229 )
230
231 while ud <= u_max:
232 iter += 1
233 ud += uinc0 if simu.damage.max() < threshold else uinc1
234
235 Apply_BC(ud)
236
237 u, d, Kglob, convergence = simu.Solve(tolConv, maxIter)
238 simu.Save_Iter()
239
240 # stop if the simulation does not converge
241 if not convergence:
242 break
243
244 f = np.sum(Kglob[dofsY_upper, :] @ u)
245
246 simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True)
247
248 # Detect damaged edges
249 if np.any(d[nodes_edges] >= 1):
250 nDetect += 1
251 if nDetect == 10:
252 break
253
254 displacement = np.concatenate((displacement, [ud]))
255 force = np.concatenate((force, [f]))
256
257 if plotIter:
258 Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter)
259 plt.figure(axIter.figure)
260 plt.pause(1e-12)
261
262 force = np.asarray(force)
263 displacement = np.asarray(displacement)
264 Display.Plot_Force_Displacement(
265 force / unit,
266 displacement * unit,
267 f"ud [{unitU}]",
268 f"f [{unitF}]",
269 ax=axLoad,
270 )[1]
271 plt.figure(axLoad.figure)
272 plt.pause(1e-12)
273
274 # ----------------------------------------------
275 # Saving
276 # ----------------------------------------------
277 force = np.asarray(force)
278 displacement = np.asarray(displacement)
279 print()
280 Simulations.Save_pickle(
281 (force, displacement), folder_save, "force-displacement"
282 )
283 simu.Save(folder_save)
284
285 else:
286 # ----------------------------------------------
287 # Load
288 # ----------------------------------------------
289 simu: Simulations.PhaseFieldSimu = Simulations.Load_Simu(folder_save)
290 force, displacement = Simulations.Load_pickle(folder_save, "force-displacement")
291
292 # ----------------------------------------------
293 # Results
294 # ---------------------------------------------
295 if plotEnergy:
296 Display.Plot_Energy(simu, force, displacement, N=400, folder=folder_save)
297
298 if plotResult:
299 Display.Plot_Result(
300 simu,
301 "damage",
302 nodeValues=True,
303 colorbarIsClose=True,
304 folder=folder_save,
305 filename="damage",
306 )
307 Display.Plot_Mesh(simu)
308 Display.Plot_BoundaryConditions(simu)
309 Display.Plot_Iter_Summary(simu, folder_save, None, None)
310 Display.Plot_Force_Displacement(
311 force / unit,
312 displacement * unit,
313 f"ud [{unitU}]",
314 f"f [{unitF}]",
315 folder_save,
316 )
317
318 if plotMesh:
319 Display.Plot_Mesh(mesh)
320
321 if saveParaview:
322 Paraview.Save_simu(simu, folder_save)
323
324 if makeMovie:
325 simu.Set_Iter(-1)
326 depMax = simu.Result("displacement_norm").max()
327 deformFactor = L * 0.05 / depMax
328
329 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
330
331 def Func(plotter, iter):
332 simu.Set_Iter(iterations[iter])
333
334 grid = PyVista._pyVistaMesh(simu, "damage", deformFactor)
335
336 tresh = grid.threshold((0, 0.8))
337
338 PyVista.Plot(
339 tresh,
340 "damage",
341 deformFactor,
342 plotMesh=True,
343 plotter=plotter,
344 clim=(0, 1),
345 )
346
347 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
348
349 if doSimu:
350 Tic.Plot_History(folder_save, details=False)
351
352 if showFig:
353 plt.show()
354
355 # plt.close("all")
356 Tic.Clear()
357
358
359 if __name__ == "__main__":
360 # generates configs
361 Splits = []
362 Regus = []
363 for split in splits.copy():
364 for regu in regus.copy():
365 Splits.append(split)
366 Regus.append(regu)
367
368 if useParallel:
369 items = [(split, regu) for split, regu in zip(Splits, Regus)]
370 with multiprocessing.Pool(nProcs) as pool:
371 for result in pool.starmap(DoSimu, items):
372 pass
373 else:
374 [DoSimu(split, regu) for split, regu in zip(Splits, Regus)]
Total running time of the script: (0 minutes 12.863 seconds)