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