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.8.3/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh
2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:0.826s, tol=0.00e+00
3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.816s, tol=1.00e+00 3.20% -> 49.37s
4 : 1.600 μm, [0.00e+00; 0.00e+00], 1:0.809s, tol=7.50e-01 6.40% -> 35.48s
5 : 2.400 μm, [0.00e+00; 0.00e+00], 1:0.816s, tol=5.56e-01 9.60% -> 30.73s
6 : 3.200 μm, [0.00e+00; 0.00e+00], 1:0.816s, tol=4.38e-01 12.80% -> 27.79s
7 : 4.000 μm, [0.00e+00; 0.00e+00], 1:0.815s, tol=3.60e-01 16.00% -> 25.66s
8 : 4.800 μm, [0.00e+00; 0.00e+00], 1:0.813s, tol=3.06e-01 19.20% -> 23.96s
9 : 5.600 μm, [0.00e+00; 0.00e+00], 1:0.844s, tol=2.65e-01 22.40% -> 23.39s
10 : 6.400 μm, [0.00e+00; 0.00e+00], 1:0.836s, tol=2.34e-01 25.60% -> 21.86s
11 : 7.200 μm, [0.00e+00; 0.00e+00], 1:0.865s, tol=2.10e-01 28.80% -> 21.38s
12 : 8.000 μm, [0.00e+00; 0.00e+00], 1:0.841s, tol=1.90e-01 32.00% -> 19.66s
13 : 8.800 μm, [0.00e+00; 0.00e+00], 1:0.853s, tol=1.74e-01 35.20% -> 18.85s
14 : 9.600 μm, [0.00e+00; 0.00e+00], 1:0.833s, tol=1.60e-01 38.40% -> 17.37s
15 : 10.400 μm, [0.00e+00; 0.00e+00], 1:0.862s, tol=1.48e-01 41.60% -> 16.93s
16 : 11.200 μm, [0.00e+00; 0.00e+00], 1:0.856s, tol=1.38e-01 44.80% -> 15.82s
17 : 12.000 μm, [0.00e+00; 0.00e+00], 1:0.845s, tol=1.29e-01 48.00% -> 14.64s
18 : 12.800 μm, [0.00e+00; 0.00e+00], 1:0.834s, tol=1.21e-01 51.20% -> 13.51s
19 : 13.600 μm, [0.00e+00; 0.00e+00], 1:0.833s, tol=1.14e-01 54.40% -> 12.57s
20 : 14.400 μm, [0.00e+00; 0.00e+00], 1:0.827s, tol=1.08e-01 57.60% -> 11.57s
21 : 15.200 μm, [0.00e+00; 0.00e+00], 1:0.837s, tol=1.02e-01 60.80% -> 10.79s
22 : 16.000 μm, [0.00e+00; 0.00e+00], 1:0.829s, tol=9.75e-02 64.00% -> 9.79s
23 : 16.800 μm, [9.79e-13; 4.17e-03], 1:0.847s, tol=9.30e-02 67.20% -> 9.10s
24 : 17.600 μm, [1.02e-11; 4.36e-02], 1:0.849s, tol=8.88e-02 70.40% -> 8.21s
25 : 18.400 μm, [1.09e-11; 1.18e-01], 1:0.858s, tol=8.49e-02 73.60% -> 7.38s
26 : 19.200 μm, [7.24e-12; 2.14e-01], 1:0.842s, tol=8.14e-02 76.80% -> 6.36s
27 : 20.000 μm, [4.04e-12; 3.31e-01], 1:0.838s, tol=7.82e-02 80.00% -> 5.45s
28 : 20.800 μm, [2.70e-12; 5.03e-01], 1:0.829s, tol=7.50e-02 83.20% -> 4.52s
29 : 21.600 μm, [1.48e-12; 7.38e-01], 1:0.846s, tol=7.24e-02 86.40% -> 3.73s
30 : 21.800 μm, [1.01e-12; 9.34e-01], 1:0.846s, tol=2.12e-02 87.20% -> 3.60s
31 : 22.000 μm, [1.11e-12; 9.98e-01], 1:0.844s, tol=3.65e-02 88.00% -> 3.45s
32 : 22.200 μm, [1.22e-12; 1.01e+00], 1:0.829s, tol=6.20e-02 88.80% -> 3.24s
33 : 22.400 μm, [1.36e-12; 1.00e+00], 1:0.838s, tol=8.53e-02 89.60% -> 3.11s
34 : 22.600 μm, [1.35e-12; 1.01e+00], 1:0.829s, tol=9.21e-02 90.40% -> 2.91s
35 : 22.800 μm, [1.44e-12; 1.00e+00], 1:0.851s, tol=1.06e-01 91.20% -> 2.79s
36 : 23.000 μm, [1.49e-12; 1.00e+00], 1:0.837s, tol=1.02e-01 92.00% -> 2.55s
37 : 23.200 μm, [1.63e-12; 1.00e+00], 1:0.854s, tol=9.02e-02 92.80% -> 2.38s
38 : 23.400 μm, [1.67e-12; 1.00e+00], 1:0.838s, tol=9.63e-02 93.60% -> 2.12s
39 : 23.600 μm, [1.75e-12; 1.00e+00], 1:0.821s, tol=1.02e-01 94.40% -> 1.85s
40 : 23.800 μm, [1.79e-12; 1.00e+00], 1:0.808s, tol=1.03e-01 95.20% -> 1.59s
41 : 24.000 μm, [1.68e-12; 1.00e+00], 1:0.837s, tol=1.01e-01 96.00% -> 1.39s
42 : 24.200 μm, [1.43e-12; 1.00e+00], 1:0.829s, tol=9.89e-02 96.80% -> 1.12s
43 : 24.400 μm, [1.26e-12; 1.00e+00], 1:0.839s, tol=1.01e-01 97.60% -> 866.89ms
44 : 24.600 μm, [1.10e-12; 1.00e+00], 1:0.845s, tol=9.77e-02 98.40% -> 590.54ms
45 : 24.800 μm, [1.02e-12; 1.00e+00], 1:0.844s, tol=9.67e-02 99.20% -> 299.46ms
46 : 25.000 μm, [9.35e-13; 1.00e+00], 1:0.851s, tol=9.35e-02 100.00% -> -0.00µs
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.3/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh
Generate movie 01/23 (4.35 %) 2.65 s
Generate movie 02/23 (8.70 %) 2.30 s
Generate movie 03/23 (13.04 %) 2.19 s
Generate movie 04/23 (17.39 %) 2.07 s
Generate movie 05/23 (21.74 %) 1.97 s
Generate movie 06/23 (26.09 %) 1.83 s
Generate movie 07/23 (30.43 %) 1.89 s
Generate movie 08/23 (34.78 %) 1.67 s
Generate movie 09/23 (39.13 %) 1.57 s
Generate movie 10/23 (43.48 %) 1.47 s
Generate movie 11/23 (47.83 %) 1.38 s
Generate movie 12/23 (52.17 %) 1.24 s
Generate movie 13/23 (56.52 %) 1.13 s
Generate movie 14/23 (60.87 %) 989.28 ms
Generate movie 15/23 (65.22 %) 889.68 ms
Generate movie 16/23 (69.57 %) 774.13 ms
Generate movie 17/23 (73.91 %) 669.00 ms
Generate movie 18/23 (78.26 %) 565.78 ms
Generate movie 19/23 (82.61 %) 455.39 ms
Generate movie 20/23 (86.96 %) 336.75 ms
Generate movie 21/23 (91.30 %) 230.04 ms
Generate movie 22/23 (95.65 %) 113.78 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 E = {E:.2e} Pa; v = {v}; Gc = {Gc} J/m2; l0 = {l0:.2e} m
140
141 while ud <= u_max = {u_max:.2e}:
142
143 ud += uinc0 if simu.damage.max() < {threshold} else uinc1
144 uinc0 = {uinc0:.1e}; uinc1 = {uinc1:.1e}
145
146 simu.add_dirichlet(nodes_lower, [0], ["y"])
147 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
148 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
149 """
150
151 # ----------------------------------------------
152 # Simulation
153 # ----------------------------------------------
154 simu = Simulations.PhaseField(mesh, pfm, verbosity=False)
155 simu.Results_Set_Bc_Summary(config)
156
157 dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"])
158
159 def Apply_BC(ud: float):
160 simu.Bc_Init()
161 simu.add_dirichlet(nodes_lower, [0], ["y"])
162 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
163 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
164
165 list_dep = []
166 list_f = []
167 ud = -uinc0
168 iter = 0
169 nDetect = 0
170
171 if plotIter:
172 axIter = Display.Plot_Result(simu, "damage", nodeValues=True)
173 _, axLoad = plt.subplots()
174 axLoad.set_xlabel(f"ud [{unitU}]")
175 axLoad.set_ylabel(f"f [{unitF}]")
176 axLoad.grid()
177
178 while ud <= u_max:
179 iter += 1
180 ud += uinc0 if simu.damage.max() < threshold else uinc1
181
182 Apply_BC(ud)
183
184 u, d, convergence = simu.Solve(tolConv, maxIter)
185 simu.Save_Iter()
186
187 if not convergence:
188 break
189
190 f = np.sum(simu.Calc_Reaction(dofsY_upper, "elastic"))
191 simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True)
192
193 if simu.Detect_Damage(nodes_edges, 1):
194 nDetect += 1
195 if nDetect == 10:
196 break
197
198 list_dep.append(ud)
199 list_f.append(f)
200
201 if plotIter:
202 Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter)
203 plt.figure(axIter.figure)
204 plt.pause(1e-12)
205 axLoad.clear()
206 axLoad.plot(np.abs(list_dep * unit), np.abs(list_f / unit), c="blue")
207 axLoad.set_xlabel(f"ud [{unitU}]")
208 axLoad.set_ylabel(f"f [{unitF}]")
209 axLoad.grid()
210 plt.figure(axLoad.figure)
211 plt.pause(1e-12)
212
213 # ----------------------------------------------
214 # Saving
215 # ----------------------------------------------
216 print()
217 Simulations.Save_pickle((list_f, list_dep), folder_save, "force-displacement")
218 simu.Save(folder_save)
219
220 else:
221 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
222 list_f, list_dep = Simulations.Load_pickle(folder_save, "force-displacement")
223
224 # ----------------------------------------------
225 # Results
226 # ----------------------------------------------
227 if plotEnergy:
228 Display.Plot_Energy(simu, list_f, list_dep, N=400, folder=folder_save)
229
230 Display.Plot_Result(
231 simu,
232 "damage",
233 nodeValues=True,
234 colorbarIsClose=True,
235 folder=folder_save,
236 filename="damage",
237 )
238 Display.Plot_Mesh(simu)
239 Display.Plot_BoundaryConditions(simu)
240 Display.Plot_Iter_Summary(simu, folder_save, None, None)
241
242 ax = Display.Init_Axes()
243 ax.plot(np.abs(list_dep) * unit, np.abs(list_f) / unit, c="blue")
244 ax.set_xlabel(f"ud [{unitU}]")
245 ax.set_ylabel(f"f [{unitF}]")
246 ax.grid()
247 Display.Save_fig(folder_save, "force-displacement")
248
249 if plotMesh:
250 Display.Plot_Mesh(simu.mesh)
251
252 if makeParaview:
253 Paraview.Save_simu(simu, folder_save)
254
255 if makeMovie:
256 simu.Set_Iter(-1)
257 deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
258
259 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
260
261 def Func(plotter, iter):
262 simu.Set_Iter(iterations[iter])
263 thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
264 PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
265
266 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
267
268 if doSimu:
269 Tic.Plot_History(folder_save, details=False)
270
271 plt.show()
Total running time of the script: (0 minutes 43.580 seconds)