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.8.1/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.431s, tol=0.00e+00
   2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.411s, tol=1.00e+00  2.00% -> 20.12s
   3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.411s, tol=1.00e+00  4.00% -> 19.75s
   4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.413s, tol=5.56e-01  6.00% -> 19.42s
   5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.411s, tol=4.37e-01  8.00% -> 18.89s
   6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.411s, tol=3.60e-01  10.00% -> 18.50s
   7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.411s, tol=3.06e-01  12.00% -> 18.10s
   8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.410s, tol=2.65e-01  14.00% -> 17.63s
   9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.412s, tol=2.32e-01  16.00% -> 17.29s
  10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.411s, tol=2.04e-01  18.00% -> 16.87s
  11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.411s, tol=1.81e-01  20.00% -> 16.43s
  12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.413s, tol=1.65e-01  22.00% -> 16.09s
  13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.411s, tol=1.52e-01  24.00% -> 15.63s
  14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.412s, tol=1.59e-01  26.00% -> 15.23s
  15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.412s, tol=2.32e-01  28.00% -> 14.83s
  16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.417s, tol=3.72e-01  30.00% -> 14.61s
  17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.411s, tol=4.13e-01  32.00% -> 13.99s
  18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.411s, tol=3.98e-01  34.00% -> 13.56s
  19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.410s, tol=2.98e-01  36.00% -> 13.14s
  20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.410s, tol=2.73e-01  38.00% -> 12.70s
  21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.414s, tol=2.61e-01  40.00% -> 12.43s
  22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.411s, tol=2.47e-01  42.00% -> 11.93s
  23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.412s, tol=2.42e-01  44.00% -> 11.53s
  24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.414s, tol=2.07e-01  46.00% -> 11.18s
  25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.436s, tol=2.11e-01  48.00% -> 11.34s
  26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.438s, tol=1.77e-01  50.00% -> 10.95s
  27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.437s, tol=1.67e-01  52.00% -> 10.49s
  28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.442s, tol=1.74e-01  54.00% -> 10.16s
  29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.443s, tol=1.49e-01  56.00% -> 9.74s
  30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.436s, tol=1.39e-01  58.00% -> 9.15s
  31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.439s, tol=1.29e-01  60.00% -> 8.78s
  32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.440s, tol=1.31e-01  62.00% -> 8.36s
  33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:0.466s, tol=1.26e-01  64.00% -> 8.39s
  34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.464s, tol=1.12e-01  66.00% -> 7.88s
  35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:0.466s, tol=1.11e-01  68.00% -> 7.46s
  36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:0.548s, tol=1.09e-01  70.00% -> 8.22s
  37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:0.551s, tol=1.08e-01  72.00% -> 7.71s
  38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:0.548s, tol=1.09e-01  74.00% -> 7.12s
  39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:0.554s, tol=1.13e-01  76.00% -> 6.65s
  40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:0.553s, tol=1.24e-01  78.00% -> 6.08s
  41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:0.559s, tol=1.55e-01  80.00% -> 5.59s
  42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:0.561s, tol=4.20e-01  82.00% -> 5.05s
  43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:0.570s, tol=7.43e-01  84.00% -> 4.56s
  44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:0.588s, tol=1.25e-01  86.00% -> 4.11s
  45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:0.578s, tol=1.33e-01  88.00% -> 3.47s
  46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:0.586s, tol=1.92e-02  90.00% -> 2.93s
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.1/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100

Generate movie 01/23 (4.35 %) 2.56 s
Generate movie 02/23 (8.70 %) 2.18 s
Generate movie 03/23 (13.04 %) 2.03 s
Generate movie 04/23 (17.39 %) 2.11 s
Generate movie 05/23 (21.74 %) 1.96 s
Generate movie 06/23 (26.09 %) 1.74 s
Generate movie 07/23 (30.43 %) 1.65 s
Generate movie 08/23 (34.78 %) 1.55 s
Generate movie 09/23 (39.13 %) 1.56 s
Generate movie 10/23 (43.48 %) 1.30 s
Generate movie 11/23 (47.83 %) 1.20 s
Generate movie 12/23 (52.17 %) 1.13 s
Generate movie 13/23 (56.52 %) 1.02 s
Generate movie 14/23 (60.87 %) 950.68 ms
Generate movie 15/23 (65.22 %) 820.37 ms
Generate movie 16/23 (69.57 %) 719.31 ms
Generate movie 17/23 (73.91 %) 611.70 ms
Generate movie 18/23 (78.26 %) 531.85 ms
Generate movie 19/23 (82.61 %) 418.13 ms
Generate movie 20/23 (86.96 %) 299.84 ms
Generate movie 21/23 (91.30 %) 206.40 ms
Generate movie 22/23 (95.65 %) 99.44 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 = Models.PhaseField.SplitType.Miehe
 54     regu = Models.PhaseField.ReguType.AT1
 55
 56     folder_save = Simulations.PhaseField.Folder(
 57         f"{folder}{dim}D", "", split, regu, "", 1, "", meshTest, optimMesh, nL=nL
 58     )
 59
 60     Display.MyPrint(folder_save, "green", end="\n")
 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, 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 simu.Detect_Damage(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         deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
223
224         iterations = np.arange(0, simu.Niter, simu.Niter // 20)
225
226         def Func(plotter, iter):
227             simu.Set_Iter(iterations[iter])
228             thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
229             PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
230
231         PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
232
233     Display.Plot_Result(simu, "damage", folder=folder_save)
234     Display.Plot_Result(simu, "uy", deformFactor=1)
235     Display.Plot_Mesh(mesh)
236     Display.Plot_Tags(mesh, folder=folder_save)
237
238     plt.show()

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

Gallery generated by Sphinx-Gallery