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.1/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.982 s, tol=0.00e+00
   2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.918 s, tol=1.00e+00  2.00 % -> 44.98 s
   3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.853 s, tol=1.00e+00  4.00 % -> 40.92 s
   4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.875 s, tol=5.56e-01  6.00 % -> 41.11 s
   5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.888 s, tol=4.37e-01  8.00 % -> 40.86 s
   6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.885 s, tol=3.60e-01  10.00 % -> 39.81 s
   7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.872 s, tol=3.06e-01  12.00 % -> 38.39 s
   8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.874 s, tol=2.65e-01  14.00 % -> 37.60 s
   9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.889 s, tol=2.32e-01  16.00 % -> 37.36 s
  10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.816 s, tol=2.04e-01  18.00 % -> 33.46 s
  11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.838 s, tol=1.81e-01  20.00 % -> 33.53 s
  12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.891 s, tol=1.65e-01  22.00 % -> 34.75 s
  13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.890 s, tol=1.52e-01  24.00 % -> 33.82 s
  14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.872 s, tol=1.59e-01  26.00 % -> 32.28 s
  15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.887 s, tol=2.32e-01  28.00 % -> 31.93 s
  16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.922 s, tol=3.72e-01  30.00 % -> 32.28 s
  17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.901 s, tol=4.13e-01  32.00 % -> 30.63 s
  18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.905 s, tol=3.98e-01  34.00 % -> 29.86 s
  19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.876 s, tol=2.98e-01  36.00 % -> 28.03 s
  20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.912 s, tol=2.73e-01  38.00 % -> 28.27 s
  21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.918 s, tol=2.61e-01  40.00 % -> 27.53 s
  22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.913 s, tol=2.47e-01  42.00 % -> 26.49 s
  23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.889 s, tol=2.42e-01  44.00 % -> 24.90 s
  24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.896 s, tol=2.07e-01  46.00 % -> 24.20 s
  25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.948 s, tol=2.11e-01  48.00 % -> 24.65 s
  26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.926 s, tol=1.77e-01  50.00 % -> 23.15 s
  27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.969 s, tol=1.67e-01  52.00 % -> 23.25 s
  28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.919 s, tol=1.74e-01  54.00 % -> 21.14 s
  29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.943 s, tol=1.49e-01  56.00 % -> 20.75 s
  30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.904 s, tol=1.39e-01  58.00 % -> 18.99 s
  31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.941 s, tol=1.29e-01  60.00 % -> 18.81 s
  32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.978 s, tol=1.31e-01  62.00 % -> 18.58 s
  33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:1.002 s, tol=1.26e-01  64.00 % -> 18.03 s
  34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.990 s, tol=1.12e-01  66.00 % -> 16.83 s
  35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:1.008 s, tol=1.11e-01  68.00 % -> 16.12 s
  36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:1.153 s, tol=1.09e-01  70.00 % -> 17.30 s
  37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:1.178 s, tol=1.08e-01  72.00 % -> 16.50 s
  38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:1.162 s, tol=1.09e-01  74.00 % -> 15.11 s
  39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:1.105 s, tol=1.13e-01  76.00 % -> 13.26 s
  40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:1.167 s, tol=1.24e-01  78.00 % -> 12.84 s
  41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:1.163 s, tol=1.55e-01  80.00 % -> 11.63 s
  42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:1.150 s, tol=4.20e-01  82.00 % -> 10.35 s
  43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:1.200 s, tol=7.43e-01  84.00 % -> 9.60 s
  44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:1.208 s, tol=1.25e-01  86.00 % -> 8.45 s
  45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:1.226 s, tol=1.33e-01  88.00 % -> 7.35 s
  46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:1.231 s, tol=1.92e-02  90.00 % -> 6.15 s
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/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.1/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100/summary.txt

Generate movie 01/23 (4.35 %) 5.73 s
Generate movie 02/23 (8.70 %) 5.08 s
Generate movie 03/23 (13.04 %) 4.71 s
Generate movie 04/23 (17.39 %) 4.57 s
Generate movie 05/23 (21.74 %) 4.28 s
Generate movie 06/23 (26.09 %) 4.07 s
Generate movie 07/23 (30.43 %) 3.84 s
Generate movie 08/23 (34.78 %) 3.59 s
Generate movie 09/23 (39.13 %) 3.37 s
Generate movie 10/23 (43.48 %) 3.11 s
Generate movie 11/23 (47.83 %) 2.94 s
Generate movie 12/23 (52.17 %) 2.66 s
Generate movie 13/23 (56.52 %) 2.39 s
Generate movie 14/23 (60.87 %) 2.15 s
Generate movie 15/23 (65.22 %) 1.93 s
Generate movie 16/23 (69.57 %) 1.71 s
Generate movie 17/23 (73.91 %) 1.43 s
Generate movie 18/23 (78.26 %) 1.19 s
Generate movie 19/23 (82.61 %) 950.03 ms
Generate movie 20/23 (86.96 %) 727.02 ms
Generate movie 21/23 (91.30 %) 486.04 ms
Generate movie 22/23 (95.65 %) 241.54 ms
Generate movie 23/23 (100.00 %) 0.00 µs


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

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

Gallery generated by Sphinx-Gallery