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.7.1/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh 2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:1.752 s, tol=0.00e+00
3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:1.686 s, tol=1.00e+00 3.20 % -> 1.70 m
4 : 1.600 μm, [1.03e-07; 3.70e-02], 1:1.686 s, tol=8.04e-01 6.40 % -> 1.23 m
5 : 2.400 μm, [3.74e-05; 3.70e-02], 1:1.700 s, tol=5.54e-01 9.60 % -> 1.07 m
6 : 3.200 μm, [8.41e-05; 3.70e-02], 1:1.724 s, tol=4.37e-01 12.80 % -> 58.74 s
7 : 4.000 μm, [1.50e-04; 3.70e-02], 1:1.697 s, tol=3.60e-01 16.00 % -> 53.45 s
8 : 4.800 μm, [2.34e-04; 3.70e-02], 1:1.684 s, tol=3.05e-01 19.20 % -> 49.60 s
9 : 5.600 μm, [3.37e-04; 3.70e-02], 1:1.688 s, tol=2.65e-01 22.40 % -> 46.77 s
10 : 6.400 μm, [4.60e-04; 3.71e-02], 1:1.704 s, tol=2.34e-01 25.60 % -> 44.56 s
11 : 7.200 μm, [6.01e-04; 4.57e-02], 1:1.627 s, tol=2.10e-01 28.80 % -> 40.22 s
12 : 8.000 μm, [7.63e-04; 5.81e-02], 1:1.684 s, tol=1.90e-01 32.00 % -> 39.36 s
13 : 8.800 μm, [9.46e-04; 7.21e-02], 1:1.682 s, tol=1.73e-01 35.20 % -> 37.15 s
14 : 9.600 μm, [1.15e-03; 8.78e-02], 1:1.695 s, tol=1.59e-01 38.40 % -> 35.34 s
15 : 10.400 μm, [1.38e-03; 1.05e-01], 1:1.695 s, tol=1.47e-01 41.60 % -> 33.31 s
16 : 11.200 μm, [1.62e-03; 1.25e-01], 1:1.684 s, tol=1.37e-01 44.80 % -> 31.12 s
17 : 12.000 μm, [1.89e-03; 1.46e-01], 1:1.700 s, tol=1.28e-01 48.00 % -> 29.46 s
18 : 12.800 μm, [2.19e-03; 1.70e-01], 1:1.671 s, tol=1.21e-01 51.20 % -> 27.07 s
19 : 13.600 μm, [2.50e-03; 1.97e-01], 1:1.690 s, tol=1.14e-01 54.40 % -> 25.50 s
20 : 14.400 μm, [2.85e-03; 2.29e-01], 1:1.708 s, tol=1.07e-01 57.60 % -> 23.88 s
21 : 15.200 μm, [3.21e-03; 2.65e-01], 1:1.700 s, tol=1.02e-01 60.80 % -> 21.93 s
22 : 16.000 μm, [3.61e-03; 3.04e-01], 1:1.694 s, tol=9.68e-02 64.00 % -> 20.01 s
23 : 16.800 μm, [4.03e-03; 3.48e-01], 1:1.698 s, tol=9.23e-02 67.20 % -> 18.23 s
24 : 17.600 μm, [4.47e-03; 3.96e-01], 1:1.597 s, tol=8.82e-02 70.40 % -> 15.44 s
25 : 18.400 μm, [4.95e-03; 4.49e-01], 1:1.702 s, tol=8.44e-02 73.60 % -> 14.65 s
26 : 19.200 μm, [5.45e-03; 5.05e-01], 1:1.687 s, tol=8.10e-02 76.80 % -> 12.74 s
27 : 20.000 μm, [5.98e-03; 5.61e-01], 1:1.675 s, tol=7.80e-02 80.00 % -> 10.89 s
28 : 20.800 μm, [6.53e-03; 6.38e-01], 1:1.678 s, tol=7.53e-02 83.20 % -> 9.15 s
29 : 21.000 μm, [7.11e-03; 7.36e-01], 1:1.679 s, tol=2.04e-02 84.00 % -> 8.95 s
30 : 21.200 μm, [7.30e-03; 8.35e-01], 1:1.551 s, tol=2.10e-02 84.80 % -> 8.06 s
31 : 21.400 μm, [7.46e-03; 9.20e-01], 1:1.670 s, tol=2.34e-02 85.60 % -> 8.43 s
32 : 21.600 μm, [7.64e-03; 9.75e-01], 1:1.685 s, tol=2.99e-02 86.40 % -> 8.22 s
33 : 21.800 μm, [7.87e-03; 9.98e-01], 1:1.681 s, tol=4.11e-02 87.20 % -> 7.89 s
34 : 22.000 μm, [8.18e-03; 1.00e+00], 1:1.670 s, tol=5.57e-02 88.00 % -> 7.51 s
35 : 22.200 μm, [8.57e-03; 1.00e+00], 1:1.682 s, tol=6.35e-02 88.80 % -> 7.21 s
36 : 22.400 μm, [8.97e-03; 1.00e+00], 1:1.680 s, tol=8.48e-02 89.60 % -> 6.82 s
37 : 22.600 μm, [9.40e-03; 1.00e+00], 1:1.667 s, tol=8.07e-02 90.40 % -> 6.37 s
38 : 22.800 μm, [9.82e-03; 1.00e+00], 1:1.670 s, tol=7.54e-02 91.20 % -> 5.96 s
39 : 23.000 μm, [1.01e-02; 1.00e+00], 1:1.663 s, tol=7.34e-02 92.00 % -> 5.49 s
40 : 23.200 μm, [1.04e-02; 1.00e+00], 1:1.686 s, tol=7.26e-02 92.80 % -> 5.10 s
41 : 23.400 μm, [1.06e-02; 1.00e+00], 1:1.678 s, tol=7.46e-02 93.60 % -> 4.59 s
42 : 23.600 μm, [1.08e-02; 1.00e+00], 1:1.656 s, tol=7.51e-02 94.40 % -> 4.03 s
43 : 23.800 μm, [1.10e-02; 1.00e+00], 1:1.668 s, tol=7.56e-02 95.20 % -> 3.53 s
44 : 24.000 μm, [1.11e-02; 1.00e+00], 1:1.693 s, tol=7.53e-02 96.00 % -> 3.03 s
45 : 24.200 μm, [1.12e-02; 1.00e+00], 1:1.694 s, tol=7.52e-02 96.80 % -> 2.46 s
46 : 24.400 μm, [1.14e-02; 1.00e+00], 1:1.692 s, tol=7.46e-02 97.60 % -> 1.87 s
47 : 24.600 μm, [1.15e-02; 1.00e+00], 1:1.695 s, tol=7.42e-02 98.40 % -> 1.27 s
48 : 24.800 μm, [1.16e-02; 1.00e+00], 1:1.688 s, tol=7.27e-02 99.20 % -> 639.69 ms
49 : 25.000 μm, [1.17e-02; 1.00e+00], 1:1.704 s, tol=7.13e-02 100.00 % -> -0.00 µs
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/force-displacement.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/simulation.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/summary.txt
Generate movie 01/24 (4.17 %) 6.10 s
Generate movie 02/24 (8.33 %) 5.58 s
Generate movie 03/24 (12.50 %) 5.28 s
Generate movie 04/24 (16.67 %) 5.02 s
Generate movie 05/24 (20.83 %) 4.74 s
Generate movie 06/24 (25.00 %) 4.69 s
Generate movie 07/24 (29.17 %) 4.27 s
Generate movie 08/24 (33.33 %) 4.02 s
Generate movie 09/24 (37.50 %) 3.75 s
Generate movie 10/24 (41.67 %) 3.52 s
Generate movie 11/24 (45.83 %) 3.27 s
Generate movie 12/24 (50.00 %) 3.03 s
Generate movie 13/24 (54.17 %) 2.67 s
Generate movie 14/24 (58.33 %) 2.49 s
Generate movie 15/24 (62.50 %) 2.27 s
Generate movie 16/24 (66.67 %) 2.00 s
Generate movie 17/24 (70.83 %) 1.75 s
Generate movie 18/24 (75.00 %) 1.45 s
Generate movie 19/24 (79.17 %) 1.25 s
Generate movie 20/24 (83.33 %) 1.01 s
Generate movie 21/24 (87.50 %) 741.85 ms
Generate movie 22/24 (91.67 %) 499.95 ms
Generate movie 23/24 (95.83 %) 251.19 ms
Generate movie 24/24 (100.00 %) 0.00 µs
12 import matplotlib.pyplot as plt
13 import numpy as np
14
15 from EasyFEA import (
16 Display,
17 Folder,
18 Models,
19 Tic,
20 ElemType,
21 Mesh,
22 Simulations,
23 PyVista,
24 Paraview,
25 )
26 from EasyFEA.Geoms import Domain, Circle
27
28 import multiprocessing
29
30 # Display.Clear()
31
32
33 # ----------------------------------------------
34 # Configuration
35 # ----------------------------------------------
36
37 # simu options
38 doSimu = True
39 meshTest = True
40 optimMesh = True
41 useParallel = False
42 nProcs = 4 # number of processes in parallel
43
44 # outputs
45 folder = Folder.Results_Dir()
46 plotMesh = False
47 plotIter = False
48 plotResult = True
49 plotEnergy = False
50 showFig = True
51
52 saveParaview = False
53 makeMovie = True
54
55 # models
56
57 # splits = ["Bourdin","Amor","Miehe","Stress"] # Splits Isotropes
58 # splits = ["He","AnisotStrain","AnisotStress","Zhang"] # Splits Anisotropes
59 # splits = ["Bourdin","Amor","Miehe","Stress","He","AnisotStrain","AnisotStress","Zhang"]
60 # splits = ["Zhang"]
61 # splits = ["AnisotStrain","AnisotStress","Zhang"]
62 splits = ["Miehe"]
63
64 regus = ["AT2"] # ["AT1", "AT2"]
65 # regus = ["AT1", "AT2"]
66
67 l0 = 0.12e-3
68
69 # convergence
70 solver = (
71 Models.PhaseField.SolverType.History
72 ) # ["History", "HistoryDamage", "BoundConstrain"]
73 maxIter = 1000
74 tolConv = 1e-0
75
76 # ----------------------------------------------
77 # Mesh
78 # ----------------------------------------------
79
80
81 def DoMesh(
82 L: float, h: float, diam: float, thickness: float, l0: float, split: str
83 ) -> Mesh:
84 clC = l0 * 2 if meshTest else l0 / 2
85 if optimMesh:
86 clD = l0 * 4
87 refineZone = diam * 1.5 / 2
88 if split in ["Bourdin", "Amor"]:
89 refineGeom = Domain((0, h / 2 - refineZone), (L, h / 2 + refineZone), clC)
90 else:
91 refineGeom = Domain((L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC)
92 else:
93 clD = l0 if meshTest else l0 / 2
94 refineGeom = None
95
96 domain = Domain((0, 0), (L, h), clD)
97 circle = Circle((L / 2, h / 2), diam, clD, isHollow=True)
98
99 # ax = Display.Init_Axes()
100 # domain.Plot(ax, color="k", plotPoints=False)
101 # circle.Plot(ax, color="k", plotPoints=False)
102 # # if refineGeom != None:
103 # # refineGeom.Plot(ax, color='k', plotPoints=False)
104 # # ax.scatter(((L+diam)/2, L/2), (h/2, (h+diam)/2), c='k')
105 # ax.axis("off")
106 # Display.Save_fig(folder, "sample", True)
107
108 mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom])
109
110 # ax = Display.Plot_Mesh(mesh, lw=0.3, facecolors="white")
111 # ax.axis("off")
112 # ax.set_title("")
113 # Display.Save_fig(folder, "mesh", transparent=True)
114
115 return mesh
116
117
118 # ----------------------------------------------
119 # Do Simu
120 # ----------------------------------------------
121 def DoSimu(split: str, regu: str):
122
123 folder_save = Simulations.PhaseField.Folder(
124 folder, "Elas_Isot", split, regu, "DP", tolConv, solver, meshTest, optimMesh
125 )
126
127 Display.MyPrint(folder_save, "green")
128
129 # ----------------------------------------------
130 # Geom
131 # ----------------------------------------------
132
133 L = 15e-3
134 h = 30e-3
135 thickness = 1
136 diam = 6e-3
137
138 # load units
139 unitU = "μm"
140 unitF = "kN/mm"
141 unit = 1e6
142
143 # ----------------------------------------------
144 # Mesh
145 # ----------------------------------------------
146
147 if doSimu:
148 mesh = DoMesh(L, h, diam, thickness, l0, split)
149
150 # Get Nodes
151 nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
152 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
153 nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0))
154 nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"])
155 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
156
157 # ----------------------------------------------
158 # Boundary conditions
159 # ----------------------------------------------
160
161 threshold = 0.6
162 u_max = 25e-6
163 # u_max = 35e-6
164
165 uinc0 = 8e-7 if meshTest else 8e-8
166 uinc1 = 2e-7 if meshTest else 2e-8
167
168 config = f"""
169 while ud <= u_max:
170
171 ud += uinc0 if simu.damage.max() < threshold else uinc1
172
173 u_max = {u_max}
174 uinc0 = {uinc0:.1e} (simu.damage.max() < {threshold})
175 uinc1 = {uinc1:.1e}
176
177 simu.add_dirichlet(nodes_lower, [0], ["y"])
178 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
179 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
180 if dim == 3:
181 simu.add_dirichlet(nodes_y0z0, [0], ["z"])
182 """
183
184 # ----------------------------------------------
185 # Material
186 # ----------------------------------------------
187 E = 12e9
188 v = 0.3
189 planeStress = False
190 material = Models.Elastic.Isotropic(2, E, v, planeStress, thickness)
191
192 gc = 1.4
193 pfm = Models.PhaseField(material, split, regu, gc, l0, solver=solver)
194
195 # ----------------------------------------------
196 # Simulation
197 # ----------------------------------------------
198 simu = Simulations.PhaseField(mesh, pfm, verbosity=False)
199
200 simu.Results_Set_Bc_Summary(config)
201
202 dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"])
203
204 def Apply_BC(ud: float):
205 simu.Bc_Init()
206 simu.add_dirichlet(nodes_lower, [0], ["y"])
207 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
208 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
209
210 # INIT
211 displacement = []
212 force = []
213 ud = -uinc0
214 iter = 0
215 nDetect = 0
216
217 if plotIter:
218 axIter = Display.Plot_Result(simu, "damage", nodeValues=True)
219
220 force = np.asarray(force)
221 displacement = np.asarray(displacement)
222 _, axLoad = Display.Plot_Force_Displacement(
223 force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]"
224 )
225
226 while ud <= u_max:
227 iter += 1
228 ud += uinc0 if simu.damage.max() < threshold else uinc1
229
230 Apply_BC(ud)
231
232 u, d, Ku, convergence = simu.Solve(tolConv, maxIter)
233 simu.Save_Iter()
234
235 # stop if the simulation does not converge
236 if not convergence:
237 break
238
239 f = np.sum(Ku[dofsY_upper, :] @ u)
240
241 simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True)
242
243 # Detect damaged edges
244 if np.any(d[nodes_edges] >= 1):
245 nDetect += 1
246 if nDetect == 10:
247 break
248
249 displacement = np.concatenate((displacement, [ud]))
250 force = np.concatenate((force, [f]))
251
252 if plotIter:
253 Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter)
254 plt.figure(axIter.figure)
255 plt.pause(1e-12)
256
257 force = np.asarray(force)
258 displacement = np.asarray(displacement)
259 Display.Plot_Force_Displacement(
260 force / unit,
261 displacement * unit,
262 f"ud [{unitU}]",
263 f"f [{unitF}]",
264 ax=axLoad,
265 )[1]
266 plt.figure(axLoad.figure)
267 plt.pause(1e-12)
268
269 # ----------------------------------------------
270 # Saving
271 # ----------------------------------------------
272 force = np.asarray(force)
273 displacement = np.asarray(displacement)
274 print()
275 Simulations.Save_pickle(
276 (force, displacement), folder_save, "force-displacement"
277 )
278 simu.Save(folder_save)
279
280 else:
281 # ----------------------------------------------
282 # Load
283 # ----------------------------------------------
284 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
285 force, displacement = Simulations.Load_pickle(folder_save, "force-displacement")
286
287 # ----------------------------------------------
288 # Results
289 # ---------------------------------------------
290 if plotEnergy:
291 Display.Plot_Energy(simu, force, displacement, N=400, folder=folder_save)
292
293 if plotResult:
294 Display.Plot_Result(
295 simu,
296 "damage",
297 nodeValues=True,
298 colorbarIsClose=True,
299 folder=folder_save,
300 filename="damage",
301 )
302 Display.Plot_Mesh(simu)
303 Display.Plot_BoundaryConditions(simu)
304 Display.Plot_Iter_Summary(simu, folder_save, None, None)
305 Display.Plot_Force_Displacement(
306 force / unit,
307 displacement * unit,
308 f"ud [{unitU}]",
309 f"f [{unitF}]",
310 folder_save,
311 )
312
313 if plotMesh:
314 Display.Plot_Mesh(mesh)
315
316 if saveParaview:
317 Paraview.Save_simu(simu, folder_save)
318
319 if makeMovie:
320 simu.Set_Iter(-1)
321 depMax = simu.Result("displacement_norm").max()
322 deformFactor = L * 0.05 / depMax
323
324 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
325
326 def Func(plotter, iter):
327 simu.Set_Iter(iterations[iter])
328
329 grid = PyVista._pvMesh(simu, "damage", deformFactor)
330
331 tresh = grid.threshold((0, 0.8))
332
333 PyVista.Plot(
334 tresh,
335 "damage",
336 deformFactor,
337 plotMesh=True,
338 plotter=plotter,
339 clim=(0, 1),
340 )
341
342 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
343
344 if doSimu:
345 Tic.Plot_History(folder_save, details=False)
346
347 if showFig:
348 plt.show()
349
350 # Display.plt.close("all")
351 Tic.Clear()
352
353
354 if __name__ == "__main__":
355 # generates configs
356 Splits = []
357 Regus = []
358 for split in splits.copy():
359 for regu in regus.copy():
360 Splits.append(split)
361 Regus.append(regu)
362
363 if useParallel:
364 items = [(split, regu) for split, regu in zip(Splits, Regus)]
365 with multiprocessing.Pool(nProcs) as pool:
366 for result in pool.starmap(DoSimu, items):
367 pass
368 else:
369 [DoSimu(split, regu) for split, regu in zip(Splits, Regus)]
Total running time of the script: (1 minutes 31.841 seconds)