Note

Go to the end to download the full example code.

LShape#

Damage simulation for a L-part.

LShape
  • $\phi$
  • TRI3: Ne = 4234, Nn = 2178
  • Boundary conditions
  • LShape
  • LShape
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.0/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50
   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.085s, tol=0.00e+00
   2 : 0.040 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=1.00e+00  4.00% -> 1.79s
   3 : 0.080 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=9.95e-01  8.00% -> 1.73s
   4 : 0.120 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=5.56e-01  12.00% -> 1.65s
   5 : 0.160 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=4.38e-01  16.00% -> 1.58s
   6 : 0.200 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=3.60e-01  20.00% -> 1.50s
   7 : 0.240 mm, [0.00e+00; 0.00e+00], 1:0.075s, tol=3.06e-01  24.00% -> 1.42s
   8 : 0.280 mm, [0.00e+00; 8.98e-03], 1:0.075s, tol=2.65e-01  28.00% -> 1.35s
   9 : 0.320 mm, [0.00e+00; 4.57e-02], 1:0.075s, tol=2.32e-01  32.00% -> 1.28s
  10 : 0.360 mm, [0.00e+00; 1.10e-01], 1:0.075s, tol=2.07e-01  36.00% -> 1.20s
  11 : 0.400 mm, [0.00e+00; 2.07e-01], 1:0.076s, tol=1.86e-01  40.00% -> 1.14s
  12 : 0.440 mm, [0.00e+00; 3.48e-01], 1:0.076s, tol=1.69e-01  44.00% -> 1.06s
  13 : 0.480 mm, [0.00e+00; 5.46e-01], 1:0.075s, tol=1.55e-01  48.00% -> 977.04ms
  14 : 0.520 mm, [0.00e+00; 7.78e-01], 1:0.075s, tol=1.62e-01  52.00% -> 900.21ms
  15 : 0.540 mm, [0.00e+00; 9.48e-01], 1:0.075s, tol=2.22e-01  54.00% -> 895.63ms
  16 : 0.560 mm, [0.00e+00; 9.93e-01], 1:0.075s, tol=3.80e-01  56.00% -> 883.93ms
  17 : 0.580 mm, [0.00e+00; 1.00e+00], 1:0.075s, tol=4.44e-01  58.00% -> 870.12ms
  18 : 0.600 mm, [0.00e+00; 1.00e+00], 1:0.075s, tol=3.82e-01  60.00% -> 852.34ms
  19 : 0.620 mm, [0.00e+00; 1.00e+00], 1:0.076s, tol=3.28e-01  62.00% -> 835.47ms
  20 : 0.640 mm, [0.00e+00; 1.00e+00], 1:0.075s, tol=2.91e-01  64.00% -> 806.00ms
  21 : 0.660 mm, [0.00e+00; 1.00e+00], 1:0.076s, tol=2.78e-01  66.00% -> 783.28ms
  22 : 0.680 mm, [0.00e+00; 1.00e+00], 1:0.079s, tol=2.54e-01  68.00% -> 776.42ms
  23 : 0.700 mm, [0.00e+00; 1.00e+00], 1:0.076s, tol=2.53e-01  70.00% -> 715.46ms
  24 : 0.720 mm, [0.00e+00; 1.00e+00], 1:0.078s, tol=2.02e-01  72.00% -> 698.52ms
  25 : 0.740 mm, [0.00e+00; 1.00e+00], 1:0.079s, tol=1.76e-01  74.00% -> 666.90ms
  26 : 0.760 mm, [0.00e+00; 1.00e+00], 1:0.079s, tol=1.87e-01  76.00% -> 620.17ms
  27 : 0.780 mm, [0.00e+00; 1.00e+00], 1:0.084s, tol=1.72e-01  78.00% -> 618.84ms
  28 : 0.800 mm, [0.00e+00; 1.00e+00], 1:0.080s, tol=1.65e-01  80.00% -> 538.23ms
  29 : 0.820 mm, [0.00e+00; 1.00e+00], 1:0.081s, tol=1.90e-01  82.00% -> 500.06ms
  30 : 0.840 mm, [0.00e+00; 1.00e+00], 1:0.081s, tol=1.65e-01  84.00% -> 445.09ms
  31 : 0.860 mm, [0.00e+00; 1.00e+00], 1:0.088s, tol=1.51e-01  86.00% -> 427.44ms
  32 : 0.880 mm, [0.00e+00; 1.00e+00], 1:0.083s, tol=1.15e-01  88.00% -> 352.51ms
  33 : 0.900 mm, [0.00e+00; 1.00e+00], 1:0.079s, tol=1.10e-01  90.00% -> 280.04ms
  34 : 0.920 mm, [0.00e+00; 1.00e+00], 1:0.083s, tol=9.93e-02  92.00% -> 239.08ms
  35 : 0.940 mm, [0.00e+00; 1.00e+00], 1:0.083s, tol=9.94e-02  94.00% -> 180.11ms
  36 : 0.960 mm, [0.00e+00; 1.00e+00], 1:0.083s, tol=9.12e-02  96.00% -> 121.28ms
  37 : 0.980 mm, [0.00e+00; 1.00e+00], 1:0.084s, tol=1.02e-01  98.00% -> 62.05ms
  38 : 1.000 mm, [0.00e+00; 1.00e+00], 1:0.087s, 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.0/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50

Generate movie 01/38 (2.63 %) 3.90 s
Generate movie 02/38 (5.26 %) 3.37 s
Generate movie 03/38 (7.89 %) 3.05 s
Generate movie 04/38 (10.53 %) 2.99 s
Generate movie 05/38 (13.16 %) 2.88 s
Generate movie 06/38 (15.79 %) 2.91 s
Generate movie 07/38 (18.42 %) 2.70 s
Generate movie 08/38 (21.05 %) 2.66 s
Generate movie 09/38 (23.68 %) 2.65 s
Generate movie 10/38 (26.32 %) 2.55 s
Generate movie 11/38 (28.95 %) 2.35 s
Generate movie 12/38 (31.58 %) 2.32 s
Generate movie 13/38 (34.21 %) 2.26 s
Generate movie 14/38 (36.84 %) 2.18 s
Generate movie 15/38 (39.47 %) 2.02 s
Generate movie 16/38 (42.11 %) 1.92 s
Generate movie 17/38 (44.74 %) 1.91 s
Generate movie 18/38 (47.37 %) 1.75 s
Generate movie 19/38 (50.00 %) 1.71 s
Generate movie 20/38 (52.63 %) 1.60 s
Generate movie 21/38 (55.26 %) 1.47 s
Generate movie 22/38 (57.89 %) 1.41 s
Generate movie 23/38 (60.53 %) 1.32 s
Generate movie 24/38 (63.16 %) 1.23 s
Generate movie 25/38 (65.79 %) 1.19 s
Generate movie 26/38 (68.42 %) 1.08 s
Generate movie 27/38 (71.05 %) 988.13 ms
Generate movie 28/38 (73.68 %) 889.48 ms
Generate movie 29/38 (76.32 %) 795.45 ms
Generate movie 30/38 (78.95 %) 717.02 ms
Generate movie 31/38 (81.58 %) 642.88 ms
Generate movie 32/38 (84.21 %) 531.94 ms
Generate movie 33/38 (86.84 %) 442.18 ms
Generate movie 34/38 (89.47 %) 354.35 ms
Generate movie 35/38 (92.11 %) 274.31 ms
Generate movie 36/38 (94.74 %) 177.50 ms
Generate movie 37/38 (97.37 %) 91.74 ms
Generate movie 38/38 (100.00 %) 0.00 µs


 13 import matplotlib.pyplot as plt
 14 import numpy as np
 15
 16 from EasyFEA import Display, Folder, Models, ElemType, Simulations, Paraview, PyVista
 17 from EasyFEA.Geoms import Point, Points, Domain, Circle
 18
 19 if __name__ == "__main__":
 20     Display.Clear()
 21
 22     # ----------------------------------------------
 23     # Configuration
 24     # ----------------------------------------------
 25
 26     # simu options
 27     doSimu = True
 28     meshTest = True
 29     optimMesh = True
 30
 31     # outputs
 32     folder = Folder.Results_Dir()
 33     pltIter = False
 34     pltLoad = False
 35     makeMovie = True
 36     makeParaview = False
 37
 38     # geom
 39     dim = 2
 40     L = 250  # mm
 41     nL = 50
 42     ep = 100
 43
 44     # model
 45     E = 2e4  # MPa
 46     v = 0.18
 47
 48     split = Models.PhaseField.SplitType.Miehe
 49     regu = Models.PhaseField.ReguType.AT1
 50     Gc = 130  # J/m2
 51     Gc *= 1000 / 1e6  # mJ/mm2
 52
 53     # convergence
 54     tolConv = 1e-0
 55     convOption = 2
 56
 57     # load
 58     uMax = 1  # mm
 59     inc0 = uMax / 25
 60     inc1 = inc0 / 2
 61
 62     config = f"""
 63     uMax = {uMax}
 64
 65     inc0 = {inc0}
 66     inc1 = {inc1}
 67
 68     while ud <= uMax:
 69
 70     ud += inc0 if simu.damage.max() < 0.6 else inc1
 71
 72     simu.add_dirichlet(nodes_circle, [0], ['d'], "damage")
 73     simu.add_dirichlet(nodes_y0, [0]*dim, simu.Get_dofs())
 74     simu.add_dirichlet(nodes_load, [ud], ['y'])
 75     """
 76
 77     # ----------------------------------------------
 78     # Mesh
 79     # ----------------------------------------------
 80     l0 = L / nL
 81
 82     if meshTest:
 83         hC = l0
 84     else:
 85         hC = 0.5
 86         # hC = 0.25
 87
 88     p1 = Point()
 89     p2 = Point(L, 0)
 90     p3 = Point(L, L)
 91     p4 = Point(2 * L - 30, L)
 92     p5 = Point(2 * L, L)
 93     p6 = Point(2 * L, 2 * L)
 94     p7 = Point(0, 2 * L)
 95     if optimMesh:
 96         h = 100
 97         refineDomain = Domain(Point(0, L - h / 3), Point(L + h / 3, L + h), hC)
 98         hD = hC * 5
 99     else:
100         refineDomain = None
101         hD = hC
102
103     contour = Points([p1, p2, p3, p4, p5, p6, p7], hD)
104
105     circle = Circle(p5, 100)
106
107     if dim == 2:
108         mesh = contour.Mesh_2D([], ElemType.TRI3, refineGeoms=[refineDomain])
109     else:
110         mesh = contour.Mesh_Extrude(
111             [], [0, 0, -ep], [3], ElemType.HEXA8, refineGeoms=[refineDomain]
112         )
113
114     # Display.Plot_Mesh(mesh)
115     # Display.Plot_Tags(mesh)
116     # from EasyFEA import PyVista
117     # PyVista.Plot_Mesh(mesh).show()
118
119     nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
120     nodes_load = mesh.Nodes_Conditions(lambda x, y, z: (y == L) & (x >= 2 * L - 30))
121     node3 = mesh.Nodes_Point(p3)
122     node4 = mesh.Nodes_Point(p4)
123     nodes_circle = mesh.Nodes_Cylinder(circle, direction=[0, 0, ep])
124     nodes_edges = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) | (y == 0))
125
126     # ----------------------------------------------
127     # Simulation
128     # ----------------------------------------------
129     material = Models.Elastic.Isotropic(dim, E, v, True, ep)
130     pfm = Models.PhaseField(material, split, regu, Gc, l0)
131
132     folder_save = Simulations.PhaseField.Folder(
133         f"{folder}{dim}D",
134         "",
135         pfm.split,
136         pfm.regularization,
137         "",
138         tolConv,
139         "",
140         meshTest,
141         optimMesh,
142         nL=nL,
143     )
144
145     Display.MyPrint(folder_save, "green", end="\n")
146
147     if doSimu:
148         simu = Simulations.PhaseField(mesh, pfm)
149         simu.Results_Set_Bc_Summary(config)
150
151         dofsY_load = simu.Bc_dofs_nodes(nodes_load, ["y"])
152
153         if pltIter:
154             axIter = Display.Plot_Result(simu, "damage")
155
156             axLoad = Display.Init_Axes()
157             axLoad.set_xlabel("displacement [mm]")
158             axLoad.set_ylabel("load [kN]")
159
160         list_ud = []
161         list_f = []
162         ud = -inc0
163         iter = -1
164
165         while ud <= uMax:
166             # update displacement
167             iter += 1
168             ud += inc0 if simu.damage.max() < 0.6 else inc1
169
170             # update boundary conditions
171             simu.Bc_Init()
172             simu.add_dirichlet(nodes_circle, [0], ["d"], "damage")
173             simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
174             simu.add_dirichlet(nodes_load, [ud], ["y"])
175
176             # solve
177             u, d, convergence = simu.Solve(tolConv, 500, convOption)
178             # calc load
179             f = np.sum(simu.Calc_Reaction(dofsY_load, "elastic"))
180
181             # saves load and displacement
182             list_ud.append(ud)
183             list_f.append(f)
184
185             # print iter
186             simu.Results_Set_Iteration_Summary(iter, ud, "mm", ud / uMax, True)
187
188             # saves iteration
189             simu.Save_Iter()
190
191             if pltIter:
192                 plt.figure(axIter.figure)
193                 Display.Plot_Result(simu, "damage", ax=axIter)
194                 plt.pause(1e-12)
195
196                 plt.figure(axLoad.figure)
197                 axLoad.scatter(ud, f / 1000, c="black")
198                 plt.pause(1e-12)
199
200             if not convergence or simu.Detect_Damage(nodes_edges):
201                 # stops simulation if damage occurs on edges or convergence has not been reached
202                 break
203
204         # saves load and displacement
205         Simulations.Save_pickle((list_f, list_ud), folder_save, "force-displacement")
206
207         # saves the simulation
208         simu.Save(folder_save)
209
210     else:
211         simu = Simulations.Load_Simu(folder_save)
212         mesh = simu.mesh
213
214     list_f, list_ud = Simulations.Load_pickle(folder_save, "force-displacement")
215
216     # ----------------------------------------------
217     # Results
218     # ----------------------------------------------
219     Display.Plot_Result(simu, "damage", folder=folder_save)
220     Display.Plot_Mesh(simu)
221     Display.Plot_BoundaryConditions(simu)
222
223     axLoad = Display.Init_Axes()
224     axLoad.set_xlabel("displacement [mm]")
225     axLoad.set_ylabel("load [kN]")
226     axLoad.plot(list_ud, np.abs(list_f) / 1000, c="blue")
227     Display.Save_fig(folder_save, "forcedep")
228
229     Display.Plot_Iter_Summary(simu, folder_save)
230
231     if makeMovie:
232         simu.Set_Iter(-1)
233         deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
234
235         iterations = np.arange(0, simu.Niter, simu.Niter // 20)
236
237         def Func(plotter, iter):
238             simu.Set_Iter(iterations[iter])
239             thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
240             PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
241
242         PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
243
244     if makeParaview:
245         Paraview.Save_simu(simu, folder_save)
246
247     plt.show()

Total running time of the script: (0 minutes 10.006 seconds)

Gallery generated by Sphinx-Gallery