Note

Go to the end to download the full example code.

CT#

Performs damage simulation on a CT specimen.

CT
  • $\phi$
  • $uy$
  • TRI3: Ne = 8982, Nn = 4661
  • CT
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.466 s, tol=0.00e+00
   2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.470 s, tol=1.00e+00  2.00 % -> 23.03 s
   3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.445 s, tol=1.00e+00  4.00 % -> 21.36 s
   4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.457 s, tol=5.56e-01  6.00 % -> 21.49 s
   5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.456 s, tol=4.37e-01  8.00 % -> 20.97 s
   6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.458 s, tol=3.60e-01  10.00 % -> 20.60 s
   7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.460 s, tol=3.06e-01  12.00 % -> 20.23 s
   8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.462 s, tol=2.65e-01  14.00 % -> 19.87 s
   9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.457 s, tol=2.32e-01  16.00 % -> 19.18 s
  10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.465 s, tol=2.04e-01  18.00 % -> 19.05 s
  11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.459 s, tol=1.81e-01  20.00 % -> 18.34 s
  12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.455 s, tol=1.65e-01  22.00 % -> 17.75 s
  13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.456 s, tol=1.52e-01  24.00 % -> 17.33 s
  14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.464 s, tol=1.59e-01  26.00 % -> 17.16 s
  15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.451 s, tol=2.32e-01  28.00 % -> 16.22 s
  16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.452 s, tol=3.72e-01  30.00 % -> 15.82 s
  17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.452 s, tol=4.13e-01  32.00 % -> 15.37 s
  18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.446 s, tol=3.98e-01  34.00 % -> 14.71 s
  19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.455 s, tol=2.98e-01  36.00 % -> 14.56 s
  20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.462 s, tol=2.73e-01  38.00 % -> 14.32 s
  21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.451 s, tol=2.61e-01  40.00 % -> 13.53 s
  22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.465 s, tol=2.47e-01  42.00 % -> 13.47 s
  23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.470 s, tol=2.42e-01  44.00 % -> 13.17 s
  24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.464 s, tol=2.07e-01  46.00 % -> 12.52 s
  25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.484 s, tol=2.11e-01  48.00 % -> 12.60 s
  26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.492 s, tol=1.77e-01  50.00 % -> 12.29 s
  27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.473 s, tol=1.67e-01  52.00 % -> 11.35 s
  28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.492 s, tol=1.74e-01  54.00 % -> 11.32 s
  29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.483 s, tol=1.49e-01  56.00 % -> 10.63 s
  30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.476 s, tol=1.39e-01  58.00 % -> 10.01 s
  31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.483 s, tol=1.29e-01  60.00 % -> 9.66 s
  32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.483 s, tol=1.31e-01  62.00 % -> 9.19 s
  33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:0.512 s, tol=1.26e-01  64.00 % -> 9.21 s
  34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.518 s, tol=1.12e-01  66.00 % -> 8.81 s
  35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:0.522 s, tol=1.11e-01  68.00 % -> 8.35 s
  36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:0.608 s, tol=1.09e-01  70.00 % -> 9.11 s
  37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:0.614 s, tol=1.08e-01  72.00 % -> 8.60 s
  38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:0.602 s, tol=1.09e-01  74.00 % -> 7.82 s
  39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:0.594 s, tol=1.13e-01  76.00 % -> 7.12 s
  40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:0.610 s, tol=1.24e-01  78.00 % -> 6.70 s
  41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:0.613 s, tol=1.55e-01  80.00 % -> 6.13 s
  42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:0.612 s, tol=4.20e-01  82.00 % -> 5.51 s
  43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:0.640 s, tol=7.43e-01  84.00 % -> 5.12 s
  44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:0.663 s, tol=1.25e-01  86.00 % -> 4.64 s
  45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:0.646 s, tol=1.33e-01  88.00 % -> 3.88 s
  46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:0.652 s, tol=1.92e-02  90.00 % -> 3.26 s
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100/simulation.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100/summary.txt

Generate movie 01/23 (4.35 %) 2.90 s
Generate movie 02/23 (8.70 %) 2.21 s
Generate movie 03/23 (13.04 %) 2.16 s
Generate movie 04/23 (17.39 %) 2.03 s
Generate movie 05/23 (21.74 %) 1.87 s
Generate movie 06/23 (26.09 %) 1.78 s
Generate movie 07/23 (30.43 %) 1.63 s
Generate movie 08/23 (34.78 %) 1.59 s
Generate movie 09/23 (39.13 %) 1.56 s
Generate movie 10/23 (43.48 %) 1.34 s
Generate movie 11/23 (47.83 %) 1.23 s
Generate movie 12/23 (52.17 %) 1.21 s
Generate movie 13/23 (56.52 %) 1.08 s
Generate movie 14/23 (60.87 %) 973.31 ms
Generate movie 15/23 (65.22 %) 857.36 ms
Generate movie 16/23 (69.57 %) 706.91 ms
Generate movie 17/23 (73.91 %) 595.96 ms
Generate movie 18/23 (78.26 %) 520.66 ms
Generate movie 19/23 (82.61 %) 434.85 ms
Generate movie 20/23 (86.96 %) 316.16 ms
Generate movie 21/23 (91.30 %) 206.36 ms
Generate movie 22/23 (95.65 %) 105.62 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 Display, Folder, Models, ElemType, Simulations, PyVista, Paraview
 17 from EasyFEA.Geoms import Point, Points, Circle, Line, Contour, Domain
 18
 19 if __name__ == "__main__":
 20     Display.Clear()
 21
 22     # ----------------------------------------------
 23     # Config
 24     # ----------------------------------------------
 25
 26     dim = 2
 27
 28     # simu options
 29     doSimu = True
 30     meshTest = True
 31     optimMesh = True
 32
 33     # outputs
 34     folder = Folder.Results_Dir()
 35     plotGeom = False
 36     plotIter = False
 37     makeParaview = False
 38     makeMovie = True
 39
 40     # geom
 41     L = 60  # mm
 42     e = 4
 43     t = 30
 44     r = 2
 45     t2 = 15
 46     diam = 8
 47     thickness = 8
 48
 49     # model
 50     Gc = 100
 51     nL = 100
 52     l0 = L / nL
 53     split = "Miehe"
 54     regu = "AT1"
 55
 56     folder_save = Simulations.PhaseField.Folder(
 57         f"{folder}{dim}D", "Isot", split, regu, "", 1, "", meshTest, optimMesh, nL=nL
 58     )
 59
 60     Display.MyPrint(folder_save, "green")
 61
 62     # ----------------------------------------------
 63     # Geom
 64     # ----------------------------------------------
 65     clC = l0 if meshTest else l0 / 2  # meshSize on the crack
 66     clD = l0 * 2 if optimMesh else clC
 67     mS = l0 / 2
 68
 69     pt1 = Point(0, -L / 2)
 70     pt2 = Point(L, -L / 2)
 71     pt3 = Point(L, L / 2)
 72     pt4 = Point(0, L / 2)
 73     pt5 = Point(0, e / 2)
 74     pt6 = Point(t + r, e / 2, r=r)
 75     pt7 = Point(t + r, -e / 2, r=r)
 76     pt8 = Point(0, -e / 2)
 77     points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], clD)
 78
 79     contour = points.Get_Contour()
 80
 81     circle1 = Circle(Point(t2, -L / 2 + t2), diam, clD)
 82     circle2 = Circle(Point(t2, L / 2 - t2), diam, clD)
 83
 84     if plotGeom:
 85         ax = contour.Plot()
 86         contour.Plot_Geoms([circle1, circle2], ax=ax)
 87
 88     # ----------------------------------------------
 89     # Mesh
 90     # ----------------------------------------------
 91
 92     refineGeom = (
 93         Domain(Point(t, -e * 1.5), Point(L, e * 1.5, thickness), clC)
 94         if optimMesh
 95         else None
 96     )
 97
 98     if plotGeom:
 99         refineGeom.Plot(ax=ax)
100
101     if dim == 2:
102         crack = Line(
103             Point(t + r, 0, isOpen=True), Point(t + r + 6, 0), clC, isOpen=True
104         )
105
106         mesh = contour.Mesh_2D(
107             [circle1, circle2],
108             ElemType.TRI3,
109             cracks=[crack],
110             refineGeoms=[refineGeom],
111         )
112     else:
113         elemType = ElemType.PRISM6
114
115         pc1 = Point(t + r, 0, 0, True)
116         pc2 = Point(t + r + 6, 0, 0)
117         pc3 = pc2 + [0, 0, thickness]
118         pc4 = Point(t + r, 0, thickness, True)
119         #
120         line1 = Line(pc1, pc2, clC, True)
121         line2 = Line(pc2, pc3, clC, False)
122         line3 = Line(pc3, pc4, clC, True)
123         line4 = Line(pc4, pc1, clC, True)
124         #
125         crack = Contour([line1, line2, line3, line4], isOpen=True)
126         cracks = [crack]
127
128         if plotGeom:
129             crack.Plot(ax=ax)
130
131         mesh = contour.Mesh_Extrude(
132             [circle1, circle2],
133             [0, 0, thickness],
134             [4],
135             elemType,
136             cracks=cracks,
137             additionalLines=[line1],
138             refineGeoms=[refineGeom],
139         )
140
141     # PyVista.Plot_Mesh(mesh).show()
142     # PyVista.Plot_Nodes(mesh, mesh.orphanNodes).show()
143
144     nodes_1 = mesh.Nodes_Cylinder(circle1)
145     nodes_2 = mesh.Nodes_Cylinder(circle2)
146
147     nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
148
149     # ----------------------------------------------
150     # Simu
151     # ----------------------------------------------
152
153     mat = Models.Elastic.Isotropic(dim, thickness=thickness, planeStress=True)
154
155     pfm = Models.PhaseField(mat, split, regu, Gc, l0)
156
157     if doSimu:
158         displacements = np.linspace(0, L / 40, 50)
159
160         config = """
161         displacements = np.linspace(0, L/40, 50)
162
163         for i, dep in enumerate(displacements):
164
165         if dim == 2:
166             simu.add_dirichlet(nodes_1, [0,-dep], ["x","y"])
167             simu.add_dirichlet(nodes_2, [0,dep], ["x","y"])
168         else:
169             simu.add_dirichlet(nodes_1, [0,-dep, 0], ["x", "y", "z"])
170             simu.add_dirichlet(nodes_2, [0,dep, 0], ["x", "y", "z"])
171         """
172
173         simu = Simulations.PhaseField(mesh, pfm)
174         simu.Results_Set_Bc_Summary(config)
175
176         if plotIter:
177             ax = Display.Plot_Result(simu, "damage", 1, plotMesh=True)
178
179         for i, dep in enumerate(displacements):
180             simu.Bc_Init()
181
182             if dim == 2:
183                 simu.add_dirichlet(nodes_1, [0, -dep], ["x", "y"])
184                 simu.add_dirichlet(nodes_2, [0, dep], ["x", "y"])
185             else:
186                 simu.add_dirichlet(nodes_1, [0, -dep, 0], ["x", "y", "z"])
187                 simu.add_dirichlet(nodes_2, [0, dep, 0], ["x", "y", "z"])
188
189             # PyVista.Plot_BoundaryConditions(simu).show()
190
191             u, d, K, converg = simu.Solve()
192
193             simu.Results_Set_Iteration_Summary(
194                 i, dep, "mm", i / displacements.size, remove=True
195             )
196
197             assert converg
198
199             simu.Save_Iter()
200
201             if np.any(d[nodes_xL] >= 1):
202                 break
203
204             if plotIter:
205                 Display.Plot_Result(simu, "damage", 1, plotMesh=True, ax=ax)
206                 plt.pause(1e-12)
207
208         simu.Save(folder_save)
209
210     else:
211         simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
212
213     # ----------------------------------------------
214     # Results
215     # ----------------------------------------------
216
217     if makeParaview:
218         Paraview.Save_simu(simu, folder_save)
219
220     if makeMovie:
221         simu.Set_Iter(-1)
222         depMax = simu.Result("displacement_norm").max()
223         deformFactor = L * 0.05 / depMax
224
225         iterations = np.arange(0, simu.Niter, simu.Niter // 20)
226
227         def Func(plotter, iter):
228             simu.Set_Iter(iterations[iter])
229
230             grid = PyVista._pvMesh(simu, "damage", deformFactor)
231
232             tresh = grid.threshold((0, 0.8))
233
234             PyVista.Plot(
235                 tresh,
236                 "damage",
237                 deformFactor,
238                 plotMesh=True,
239                 plotter=plotter,
240                 clim=(0, 1),
241             )
242
243         PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
244
245     Display.Plot_Result(simu, "damage", folder=folder_save)
246     Display.Plot_Result(simu, "uy", deformFactor=1)
247     Display.Plot_Mesh(mesh)
248     Display.Plot_Tags(mesh, folder=folder_save)
249
250     plt.show()

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

Gallery generated by Sphinx-Gallery