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/envs/v1.6.0/lib/python3.11/site-packages/results/PhaseField/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.279 s, tol=0.00e+00
   2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.275 s, tol=1.00e+00  2.00 % -> 13.49 s
   3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.258 s, tol=1.00e+00  4.00 % -> 12.39 s
   4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.275 s, tol=5.56e-01  6.00 % -> 12.92 s
   5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.292 s, tol=4.37e-01  8.00 % -> 13.41 s
   6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.272 s, tol=3.60e-01  10.00 % -> 12.22 s
   7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.299 s, tol=3.06e-01  12.00 % -> 13.14 s
   8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.289 s, tol=2.65e-01  14.00 % -> 12.42 s
   9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.304 s, tol=2.32e-01  16.00 % -> 12.78 s
  10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.305 s, tol=2.04e-01  18.00 % -> 12.51 s
  11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.323 s, tol=1.81e-01  20.00 % -> 12.92 s
  12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.327 s, tol=1.65e-01  22.00 % -> 12.73 s
  13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.306 s, tol=1.52e-01  24.00 % -> 11.64 s
  14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.278 s, tol=1.59e-01  26.00 % -> 10.28 s
  15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.317 s, tol=2.32e-01  28.00 % -> 11.41 s
  16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.306 s, tol=3.72e-01  30.00 % -> 10.72 s
  17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.317 s, tol=4.13e-01  32.00 % -> 10.78 s
  18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.310 s, tol=3.98e-01  34.00 % -> 10.24 s
  19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.307 s, tol=2.98e-01  36.00 % -> 9.81 s
  20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.305 s, tol=2.73e-01  38.00 % -> 9.46 s
  21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.312 s, tol=2.61e-01  40.00 % -> 9.37 s
  22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.303 s, tol=2.47e-01  42.00 % -> 8.78 s
  23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.307 s, tol=2.42e-01  44.00 % -> 8.58 s
  24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.321 s, tol=2.07e-01  46.00 % -> 8.67 s
  25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.310 s, tol=2.11e-01  48.00 % -> 8.05 s
  26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.315 s, tol=1.77e-01  50.00 % -> 7.88 s
  27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.311 s, tol=1.67e-01  52.00 % -> 7.46 s
  28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.299 s, tol=1.74e-01  54.00 % -> 6.88 s
  29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.288 s, tol=1.49e-01  56.00 % -> 6.33 s
  30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.179 s, tol=1.39e-01  58.00 % -> 3.75 s
  31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.307 s, tol=1.29e-01  60.00 % -> 6.13 s
  32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.313 s, tol=1.31e-01  62.00 % -> 5.95 s
  33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:0.336 s, tol=1.26e-01  64.00 % -> 6.04 s
  34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.316 s, tol=1.12e-01  66.00 % -> 5.37 s
  35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:0.322 s, tol=1.11e-01  68.00 % -> 5.15 s
  36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:0.313 s, tol=1.09e-01  70.00 % -> 4.69 s
  37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:0.311 s, tol=1.08e-01  72.00 % -> 4.36 s
  38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:0.318 s, tol=1.09e-01  74.00 % -> 4.13 s
  39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:0.310 s, tol=1.13e-01  76.00 % -> 3.72 s
  40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:0.314 s, tol=1.24e-01  78.00 % -> 3.45 s
  41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:0.320 s, tol=1.55e-01  80.00 % -> 3.20 s
  42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:0.317 s, tol=4.20e-01  82.00 % -> 2.85 s
  43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:0.311 s, tol=7.43e-01  84.00 % -> 2.49 s
  44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:0.318 s, tol=1.25e-01  86.00 % -> 2.22 s
  45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:0.310 s, tol=1.33e-01  88.00 % -> 1.86 s
  46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:0.313 s, tol=1.92e-02  90.00 % -> 1.57 s
Saved:
/results/PhaseField/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100/simulation.pickle
Saved:
/results/PhaseField/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100/summary.txt

Generate movie 0/22
Generate movie 1/22 (4.55 %) 3.54 s
Generate movie 2/22 (9.09 %) 3.33 s
Generate movie 3/22 (13.64 %) 3.23 s
Generate movie 4/22 (18.18 %) 3.00 s
Generate movie 5/22 (22.73 %) 2.84 s
Generate movie 6/22 (27.27 %) 2.66 s
Generate movie 7/22 (31.82 %) 2.50 s
Generate movie 8/22 (36.36 %) 2.37 s
Generate movie 9/22 (40.91 %) 2.20 s
Generate movie 10/22 (45.45 %) 2.00 s
Generate movie 11/22 (50.00 %) 1.84 s
Generate movie 12/22 (54.55 %) 1.66 s
Generate movie 13/22 (59.09 %) 1.52 s
Generate movie 14/22 (63.64 %) 1.33 s
Generate movie 15/22 (68.18 %) 1.16 s
Generate movie 16/22 (72.73 %) 996.88 ms
Generate movie 17/22 (77.27 %) 837.18 ms
Generate movie 18/22 (81.82 %) 660.25 ms
Generate movie 19/22 (86.36 %) 496.19 ms
Generate movie 20/22 (90.91 %) 339.74 ms
Generate movie 21/22 (95.45 %) 167.13 ms
Generate movie 22/22 (100.00 %) 0.00 µs


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

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

Gallery generated by Sphinx-Gallery