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.5/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.369 s, tol=0.00e+00
3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.360 s, tol=1.00e+00 3.20 % -> 21.77 s
4 : 1.600 μm, [3.98e-09; 3.63e-02], 1:0.392 s, tol=8.04e-01 6.40 % -> 17.18 s
5 : 2.400 μm, [3.74e-05; 3.63e-02], 1:0.408 s, tol=5.54e-01 9.60 % -> 15.36 s
6 : 3.200 μm, [8.41e-05; 3.63e-02], 1:0.415 s, tol=4.37e-01 12.80 % -> 14.13 s
7 : 4.000 μm, [1.50e-04; 3.63e-02], 1:0.416 s, tol=3.60e-01 16.00 % -> 13.10 s
8 : 4.800 μm, [2.34e-04; 3.63e-02], 1:0.292 s, tol=3.05e-01 19.20 % -> 8.60 s
9 : 5.600 μm, [3.37e-04; 3.63e-02], 1:0.387 s, tol=2.65e-01 22.40 % -> 10.71 s
10 : 6.400 μm, [4.60e-04; 3.64e-02], 1:0.418 s, tol=2.34e-01 25.60 % -> 10.93 s
11 : 7.200 μm, [6.01e-04; 4.57e-02], 1:0.415 s, tol=2.10e-01 28.80 % -> 10.25 s
12 : 8.000 μm, [7.63e-04; 5.81e-02], 1:0.354 s, tol=1.90e-01 32.00 % -> 8.27 s
13 : 8.800 μm, [9.46e-04; 7.21e-02], 1:0.383 s, tol=1.73e-01 35.20 % -> 8.45 s
14 : 9.600 μm, [1.15e-03; 8.78e-02], 1:0.270 s, tol=1.59e-01 38.40 % -> 5.63 s
15 : 10.400 μm, [1.38e-03; 1.05e-01], 1:0.317 s, tol=1.47e-01 41.60 % -> 6.22 s
16 : 11.200 μm, [1.62e-03; 1.25e-01], 1:0.423 s, tol=1.37e-01 44.80 % -> 7.82 s
17 : 12.000 μm, [1.89e-03; 1.46e-01], 1:0.431 s, tol=1.28e-01 48.00 % -> 7.47 s
18 : 12.800 μm, [2.19e-03; 1.70e-01], 1:0.410 s, tol=1.21e-01 51.20 % -> 6.64 s
19 : 13.600 μm, [2.50e-03; 1.97e-01], 1:0.353 s, tol=1.14e-01 54.40 % -> 5.33 s
20 : 14.400 μm, [2.85e-03; 2.29e-01], 1:0.416 s, tol=1.07e-01 57.60 % -> 5.81 s
21 : 15.200 μm, [3.21e-03; 2.65e-01], 1:0.407 s, tol=1.02e-01 60.80 % -> 5.24 s
22 : 16.000 μm, [3.61e-03; 3.04e-01], 1:0.416 s, tol=9.68e-02 64.00 % -> 4.91 s
23 : 16.800 μm, [4.03e-03; 3.48e-01], 1:0.403 s, tol=9.23e-02 67.20 % -> 4.33 s
24 : 17.600 μm, [4.47e-03; 3.96e-01], 1:0.310 s, tol=8.82e-02 70.40 % -> 3.00 s
25 : 18.400 μm, [4.95e-03; 4.49e-01], 1:0.410 s, tol=8.44e-02 73.60 % -> 3.53 s
26 : 19.200 μm, [5.45e-03; 5.05e-01], 1:0.417 s, tol=8.10e-02 76.80 % -> 3.15 s
27 : 20.000 μm, [5.98e-03; 5.61e-01], 1:0.411 s, tol=7.80e-02 80.00 % -> 2.67 s
28 : 20.800 μm, [6.53e-03; 6.38e-01], 1:0.413 s, tol=7.53e-02 83.20 % -> 2.25 s
29 : 21.000 μm, [7.11e-03; 7.36e-01], 1:0.321 s, tol=2.04e-02 84.00 % -> 1.71 s
30 : 21.200 μm, [7.30e-03; 8.35e-01], 1:0.385 s, tol=2.10e-02 84.80 % -> 2.00 s
31 : 21.400 μm, [7.46e-03; 9.20e-01], 1:0.411 s, tol=2.34e-02 85.60 % -> 2.07 s
32 : 21.600 μm, [7.64e-03; 9.75e-01], 1:0.392 s, tol=2.99e-02 86.40 % -> 1.91 s
33 : 21.800 μm, [7.87e-03; 9.98e-01], 1:0.355 s, tol=4.11e-02 87.20 % -> 1.67 s
34 : 22.000 μm, [8.18e-03; 1.00e+00], 1:0.405 s, tol=5.57e-02 88.00 % -> 1.82 s
35 : 22.200 μm, [8.57e-03; 1.00e+00], 1:0.414 s, tol=6.35e-02 88.80 % -> 1.77 s
36 : 22.400 μm, [8.97e-03; 1.00e+00], 1:0.408 s, tol=8.48e-02 89.60 % -> 1.66 s
37 : 22.600 μm, [9.40e-03; 1.00e+00], 1:0.413 s, tol=8.07e-02 90.40 % -> 1.58 s
38 : 22.800 μm, [9.82e-03; 1.00e+00], 1:0.393 s, tol=7.54e-02 91.20 % -> 1.40 s
39 : 23.000 μm, [1.01e-02; 1.00e+00], 1:0.249 s, tol=7.34e-02 92.00 % -> 821.96 ms
40 : 23.200 μm, [1.04e-02; 1.00e+00], 1:0.314 s, tol=7.26e-02 92.80 % -> 949.36 ms
41 : 23.400 μm, [1.06e-02; 1.00e+00], 1:0.420 s, tol=7.46e-02 93.60 % -> 1.15 s
42 : 23.600 μm, [1.08e-02; 1.00e+00], 1:0.419 s, tol=7.51e-02 94.40 % -> 1.02 s
43 : 23.800 μm, [1.10e-02; 1.00e+00], 1:0.411 s, tol=7.56e-02 95.20 % -> 869.76 ms
44 : 24.000 μm, [1.11e-02; 1.00e+00], 1:0.364 s, tol=7.53e-02 96.00 % -> 651.72 ms
45 : 24.200 μm, [1.12e-02; 1.00e+00], 1:0.421 s, tol=7.52e-02 96.80 % -> 611.98 ms
46 : 24.400 μm, [1.14e-02; 1.00e+00], 1:0.412 s, tol=7.46e-02 97.60 % -> 455.61 ms
47 : 24.600 μm, [1.15e-02; 1.00e+00], 1:0.417 s, tol=7.42e-02 98.40 % -> 311.72 ms
48 : 24.800 μm, [1.16e-02; 1.00e+00], 1:0.409 s, tol=7.27e-02 99.20 % -> 154.92 ms
49 : 25.000 μm, [1.17e-02; 1.00e+00], 1:0.337 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
Generate movie 0/23
Generate movie 1/23 (4.35 %) 3.69 s
Generate movie 2/23 (8.70 %) 3.54 s
Generate movie 3/23 (13.04 %) 3.38 s
Generate movie 4/23 (17.39 %) 3.19 s
Generate movie 5/23 (21.74 %) 3.01 s
Generate movie 6/23 (26.09 %) 2.81 s
Generate movie 7/23 (30.43 %) 2.67 s
Generate movie 8/23 (34.78 %) 2.52 s
Generate movie 9/23 (39.13 %) 2.35 s
Generate movie 10/23 (43.48 %) 2.19 s
Generate movie 11/23 (47.83 %) 2.01 s
Generate movie 12/23 (52.17 %) 1.83 s
Generate movie 13/23 (56.52 %) 1.67 s
Generate movie 14/23 (60.87 %) 1.52 s
Generate movie 15/23 (65.22 %) 1.35 s
Generate movie 16/23 (69.57 %) 1.17 s
Generate movie 17/23 (73.91 %) 1.01 s
Generate movie 18/23 (78.26 %) 807.42 ms
Generate movie 19/23 (82.61 %) 672.81 ms
Generate movie 20/23 (86.96 %) 509.05 ms
Generate movie 21/23 (91.30 %) 335.59 ms
Generate movie 22/23 (95.65 %) 168.95 ms
Generate movie 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._pvMesh(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 25.682 seconds)