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.5/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50
   1: 0.000mm, [0.00e+00;0.00e+00], 1:0.084s, tol=0.00e+00
   2: 0.040mm, [0.00e+00;0.00e+00], 1:0.074s, tol=1.00e+00  4.00% -> 1.79s
   3: 0.080mm, [0.00e+00;0.00e+00], 1:0.075s, tol=7.50e-01  8.00% -> 1.72s
   4: 0.120mm, [0.00e+00;0.00e+00], 1:0.074s, tol=5.56e-01  12.00% -> 1.64s
   5: 0.160mm, [0.00e+00;0.00e+00], 1:0.074s, tol=4.37e-01  16.00% -> 1.56s
   6: 0.200mm, [0.00e+00;0.00e+00], 1:0.075s, tol=3.60e-01  20.00% -> 1.49s
   7: 0.240mm, [0.00e+00;0.00e+00], 1:0.075s, tol=3.06e-01  24.00% -> 1.42s
   8: 0.280mm, [0.00e+00;8.98e-03], 1:0.075s, tol=2.65e-01  28.00% -> 1.34s
   9: 0.320mm, [0.00e+00;4.57e-02], 1:0.075s, tol=2.32e-01  32.00% -> 1.27s
  10: 0.360mm, [0.00e+00;1.10e-01], 1:0.075s, tol=2.07e-01  36.00% -> 1.20s
  11: 0.400mm, [0.00e+00;2.07e-01], 1:0.075s, tol=1.86e-01  40.00% -> 1.12s
  12: 0.440mm, [0.00e+00;3.48e-01], 1:0.075s, tol=1.69e-01  44.00% -> 1.05s
  13: 0.480mm, [0.00e+00;5.46e-01], 1:0.075s, tol=1.55e-01  48.00% -> 972.70ms
  14: 0.520mm, [0.00e+00;7.78e-01], 1:0.075s, tol=1.62e-01  52.00% -> 895.84ms
  15: 0.540mm, [0.00e+00;9.48e-01], 1:0.075s, tol=2.22e-01  54.00% -> 891.78ms
  16: 0.560mm, [0.00e+00;9.93e-01], 1:0.075s, tol=3.80e-01  56.00% -> 882.24ms
  17: 0.580mm, [0.00e+00;1.00e+00], 1:0.077s, tol=4.44e-01  58.00% -> 896.45ms
  18: 0.600mm, [0.00e+00;1.00e+00], 1:0.075s, tol=3.82e-01  60.00% -> 847.73ms
  19: 0.620mm, [0.00e+00;1.00e+00], 1:0.075s, tol=3.28e-01  62.00% -> 825.26ms
  20: 0.640mm, [0.00e+00;1.00e+00], 1:0.075s, tol=2.91e-01  64.00% -> 804.07ms
  21: 0.660mm, [0.00e+00;1.00e+00], 1:0.075s, tol=2.78e-01  66.00% -> 774.51ms
  22: 0.680mm, [0.00e+00;1.00e+00], 1:0.078s, tol=2.54e-01  68.00% -> 772.57ms
  23: 0.700mm, [0.00e+00;1.00e+00], 1:0.076s, tol=2.53e-01  70.00% -> 716.52ms
  24: 0.720mm, [0.00e+00;1.00e+00], 1:0.078s, tol=2.03e-01  72.00% -> 697.20ms
  25: 0.740mm, [0.00e+00;1.00e+00], 1:0.079s, tol=1.76e-01  74.00% -> 665.82ms
  26: 0.760mm, [0.00e+00;1.00e+00], 1:0.078s, tol=1.87e-01  76.00% -> 618.67ms
  27: 0.780mm, [0.00e+00;1.00e+00], 1:0.084s, tol=1.72e-01  78.00% -> 617.04ms
  28: 0.800mm, [0.00e+00;1.00e+00], 1:0.080s, tol=1.65e-01  80.00% -> 537.92ms
  29: 0.820mm, [0.00e+00;1.00e+00], 1:0.096s, tol=1.90e-01  82.00% -> 591.18ms
  30: 0.840mm, [0.00e+00;1.00e+00], 1:0.080s, tol=1.64e-01  84.00% -> 441.18ms
  31: 0.860mm, [0.00e+00;1.00e+00], 1:0.087s, tol=1.51e-01  86.00% -> 423.48ms
  32: 0.880mm, [0.00e+00;1.00e+00], 1:0.083s, tol=1.15e-01  88.00% -> 350.52ms
  33: 0.900mm, [0.00e+00;1.00e+00], 1:0.079s, tol=1.10e-01  90.00% -> 279.12ms
  34: 0.920mm, [0.00e+00;1.00e+00], 1:0.084s, tol=9.94e-02  92.00% -> 241.01ms
  35: 0.940mm, [0.00e+00;1.00e+00], 1:0.083s, tol=9.95e-02  94.00% -> 180.00ms
  36: 0.960mm, [0.00e+00;1.00e+00], 1:0.083s, tol=9.12e-02  96.00% -> 120.54ms
  37: 0.980mm, [0.00e+00;1.00e+00], 1:0.084s, tol=1.02e-01  98.00% -> 61.84ms
  38: 1.000mm, [0.00e+00;1.00e+00], 1:0.087s, tol=9.35e-02  100.00% -> -0.00µs
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50

Generate movie 01/38 (2.63 %) 4.00 s
Generate movie 02/38 (5.26 %) 3.41 s
Generate movie 03/38 (7.89 %) 3.13 s
Generate movie 04/38 (10.53 %) 3.08 s
Generate movie 05/38 (13.16 %) 2.96 s
Generate movie 06/38 (15.79 %) 2.88 s
Generate movie 07/38 (18.42 %) 2.77 s
Generate movie 08/38 (21.05 %) 2.64 s
Generate movie 09/38 (23.68 %) 2.63 s
Generate movie 10/38 (26.32 %) 2.48 s
Generate movie 11/38 (28.95 %) 2.40 s
Generate movie 12/38 (31.58 %) 2.37 s
Generate movie 13/38 (34.21 %) 2.24 s
Generate movie 14/38 (36.84 %) 2.18 s
Generate movie 15/38 (39.47 %) 2.06 s
Generate movie 16/38 (42.11 %) 2.01 s
Generate movie 17/38 (44.74 %) 1.88 s
Generate movie 18/38 (47.37 %) 1.78 s
Generate movie 19/38 (50.00 %) 1.73 s
Generate movie 20/38 (52.63 %) 1.62 s
Generate movie 21/38 (55.26 %) 1.53 s
Generate movie 22/38 (57.89 %) 1.51 s
Generate movie 23/38 (60.53 %) 1.38 s
Generate movie 24/38 (63.16 %) 1.31 s
Generate movie 25/38 (65.79 %) 1.21 s
Generate movie 26/38 (68.42 %) 1.11 s
Generate movie 27/38 (71.05 %) 996.33 ms
Generate movie 28/38 (73.68 %) 926.50 ms
Generate movie 29/38 (76.32 %) 821.65 ms
Generate movie 30/38 (78.95 %) 738.82 ms
Generate movie 31/38 (81.58 %) 633.96 ms
Generate movie 32/38 (84.21 %) 548.93 ms
Generate movie 33/38 (86.84 %) 461.24 ms
Generate movie 34/38 (89.47 %) 367.90 ms
Generate movie 35/38 (92.11 %) 282.41 ms
Generate movie 36/38 (94.74 %) 179.90 ms
Generate movie 37/38 (97.37 %) 90.92 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     inc0 = {inc0}
 65     inc1 = {inc1}
 66
 67     while ud <= uMax:
 68
 69         ud += inc0 if simu.damage.max() < 0.6 else inc1
 70
 71         simu.add_dirichlet(nodes_circle, [0], ['d'], "damage")
 72         simu.add_dirichlet(nodes_y0, [0]*dim, simu.Get_dofs())
 73         simu.add_dirichlet(nodes_load, [ud], ['y'])
 74     """
 75
 76     # ----------------------------------------------
 77     # Mesh
 78     # ----------------------------------------------
 79     l0 = L / nL
 80
 81     if meshTest:
 82         hC = l0
 83     else:
 84         hC = 0.5
 85         # hC = 0.25
 86
 87     p1 = Point()
 88     p2 = Point(L, 0)
 89     p3 = Point(L, L)
 90     p4 = Point(2 * L - 30, L)
 91     p5 = Point(2 * L, L)
 92     p6 = Point(2 * L, 2 * L)
 93     p7 = Point(0, 2 * L)
 94     if optimMesh:
 95         h = 100
 96         refineDomain = Domain(Point(0, L - h / 3), Point(L + h / 3, L + h), hC)
 97         hD = hC * 5
 98     else:
 99         refineDomain = None
100         hD = hC
101
102     contour = Points([p1, p2, p3, p4, p5, p6, p7], hD)
103
104     circle = Circle(p5, 100)
105
106     if dim == 2:
107         mesh = contour.Mesh_2D([], ElemType.TRI3, refineGeoms=[refineDomain])
108     else:
109         mesh = contour.Mesh_Extrude(
110             [], [0, 0, -ep], [3], ElemType.HEXA8, refineGeoms=[refineDomain]
111         )
112
113     # Display.Plot_Mesh(mesh)
114     # Display.Plot_Tags(mesh)
115     # from EasyFEA import PyVista
116     # PyVista.Plot_Mesh(mesh).show()
117
118     nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
119     nodes_load = mesh.Nodes_Conditions(lambda x, y, z: (y == L) & (x >= 2 * L - 30))
120     node3 = mesh.Nodes_Point(p3)
121     node4 = mesh.Nodes_Point(p4)
122     nodes_circle = mesh.Nodes_Cylinder(circle, direction=[0, 0, ep])
123     nodes_edges = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) | (y == 0))
124
125     # ----------------------------------------------
126     # Simulation
127     # ----------------------------------------------
128     material = Models.Elastic.Isotropic(dim, E, v, True, ep)
129     pfm = Models.PhaseField(material, split, regu, Gc, l0)
130
131     folder_save = Simulations.PhaseField.Folder(
132         f"{folder}{dim}D",
133         "",
134         pfm.split,
135         pfm.regularization,
136         "",
137         tolConv,
138         "",
139         meshTest,
140         optimMesh,
141         nL=nL,
142     )
143
144     Display.MyPrint(folder_save, "green", end="\n")
145
146     if doSimu:
147         simu = Simulations.PhaseField(mesh, pfm, folder=folder_save)
148         simu.Results_Set_Bc_Summary(config)
149
150         dofsY_load = simu.Bc_dofs_nodes(nodes_load, ["y"])
151
152         if pltIter:
153             axIter = Display.Plot_Result(simu, "damage")
154
155             axLoad = Display.Init_Axes()
156             axLoad.set_xlabel("displacement [mm]")
157             axLoad.set_ylabel("load [kN]")
158
159         list_ud = []
160         list_f = []
161         ud = -inc0
162         iter = -1
163
164         while ud <= uMax:
165             # update displacement
166             iter += 1
167             ud += inc0 if simu.damage.max() < 0.6 else inc1
168
169             # update boundary conditions
170             simu.Bc_Init()
171             simu.add_dirichlet(nodes_circle, [0], ["d"], "damage")
172             simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
173             simu.add_dirichlet(nodes_load, [ud], ["y"])
174
175             # solve
176             u, d, convergence = simu.Solve(tolConv, 500, convOption)
177             # calc load
178             f = np.sum(simu.Calc_Reaction(dofsY_load, "elastic"))
179
180             # saves load and displacement
181             list_ud.append(ud)
182             list_f.append(f)
183
184             # print iter
185             simu.Results_Set_Iteration_Summary(iter, ud, "mm", ud / uMax, True)
186
187             # saves iteration
188             simu.Save_Iter()
189
190             if pltIter:
191                 plt.figure(axIter.figure)
192                 Display.Plot_Result(simu, "damage", ax=axIter)
193                 plt.pause(1e-12)
194
195                 plt.figure(axLoad.figure)
196                 axLoad.scatter(ud, f / 1000, c="black")
197                 plt.pause(1e-12)
198
199             if not convergence or simu.Detect_Damage(nodes_edges):
200                 # stops simulation if damage occurs on edges or convergence has not been reached
201                 break
202
203         # saves load and displacement
204         Simulations.Save_pickle((list_f, list_ud), folder_save, "force-displacement")
205
206         # saves the simulation
207         simu.Save(folder_save)
208
209     else:
210         simu = Simulations.Load_Simu(folder_save)
211         mesh = simu.mesh
212
213     list_f, list_ud = Simulations.Load_pickle(folder_save, "force-displacement")
214
215     # ----------------------------------------------
216     # Results
217     # ----------------------------------------------
218     Display.Plot_Result(simu, "damage", folder=folder_save)
219     Display.Plot_Mesh(simu)
220     Display.Plot_BoundaryConditions(simu)
221
222     axLoad = Display.Init_Axes()
223     axLoad.set_xlabel("displacement [mm]")
224     axLoad.set_ylabel("load [kN]")
225     axLoad.plot(list_ud, np.abs(list_f) / 1000, c="blue")
226     Display.Save_fig(folder_save, "forcedep")
227
228     Display.Plot_Iter_Summary(simu, folder_save)
229
230     if makeMovie:
231         simu.Set_Iter(-1)
232         deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
233
234         iterations = np.arange(0, simu.Niter, simu.Niter // 20)
235
236         def Func(plotter, iter):
237             simu.Set_Iter(iterations[iter])
238             thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
239             PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
240
241         PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
242
243     if makeParaview:
244         Paraview.Save_simu(simu, folder_save)
245
246     plt.show()

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

Gallery generated by Sphinx-Gallery