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.9.0/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
   1: 0.000mm, [0.00e+00;0.00e+00], 1:0.624s, tol=0.00e+00
   2: 0.031mm, [0.00e+00;0.00e+00], 1:0.597s, tol=1.00e+00  2.00% -> 29.25s
   3: 0.061mm, [0.00e+00;0.00e+00], 1:0.599s, tol=7.50e-01  4.00% -> 28.78s
   4: 0.092mm, [0.00e+00;0.00e+00], 1:0.606s, tol=5.56e-01  6.00% -> 28.50s
   5: 0.122mm, [0.00e+00;0.00e+00], 1:0.598s, tol=4.38e-01  8.00% -> 27.50s
   6: 0.153mm, [0.00e+00;0.00e+00], 1:0.598s, tol=3.60e-01  10.00% -> 26.93s
   7: 0.184mm, [0.00e+00;0.00e+00], 1:0.600s, tol=3.06e-01  12.00% -> 26.41s
   8: 0.214mm, [0.00e+00;0.00e+00], 1:0.597s, tol=2.65e-01  14.00% -> 25.66s
   9: 0.245mm, [1.63e-04;2.46e-02], 1:0.599s, tol=2.32e-01  16.00% -> 25.15s
  10: 0.276mm, [2.76e-04;8.92e-02], 1:0.600s, tol=2.04e-01  18.00% -> 24.59s
  11: 0.306mm, [3.39e-04;1.99e-01], 1:0.595s, tol=1.81e-01  20.00% -> 23.81s
  12: 0.337mm, [3.21e-04;3.42e-01], 1:0.599s, tol=1.65e-01  22.00% -> 23.36s
  13: 0.367mm, [2.96e-04;5.27e-01], 1:0.601s, tol=1.52e-01  24.00% -> 22.83s
  14: 0.398mm, [2.68e-04;7.44e-01], 1:0.596s, tol=1.59e-01  26.00% -> 22.05s
  15: 0.429mm, [2.44e-04;9.33e-01], 1:0.597s, tol=2.32e-01  28.00% -> 21.51s
  16: 0.459mm, [2.15e-04;1.02e+00], 1:0.596s, tol=3.72e-01  30.00% -> 20.87s
  17: 0.490mm, [1.82e-04;1.02e+00], 1:0.607s, tol=4.13e-01  32.00% -> 20.65s
  18: 0.520mm, [1.52e-04;1.01e+00], 1:0.599s, tol=3.98e-01  34.00% -> 19.76s
  19: 0.551mm, [1.24e-04;1.01e+00], 1:0.602s, tol=2.99e-01  36.00% -> 19.26s
  20: 0.582mm, [1.03e-04;1.01e+00], 1:0.602s, tol=2.73e-01  38.00% -> 18.67s
  21: 0.612mm, [8.73e-05;1.01e+00], 1:0.604s, tol=2.61e-01  40.00% -> 18.12s
  22: 0.643mm, [7.52e-05;1.01e+00], 1:0.595s, tol=2.47e-01  42.00% -> 17.27s
  23: 0.673mm, [6.58e-05;1.01e+00], 1:0.605s, tol=2.42e-01  44.00% -> 16.95s
  24: 0.704mm, [5.74e-05;1.00e+00], 1:0.601s, tol=2.07e-01  46.00% -> 16.23s
  25: 0.735mm, [5.02e-05;1.00e+00], 1:0.630s, tol=2.11e-01  48.00% -> 16.37s
  26: 0.765mm, [4.43e-05;1.00e+00], 1:0.629s, tol=1.77e-01  50.00% -> 15.73s
  27: 0.796mm, [3.96e-05;1.01e+00], 1:0.632s, tol=1.67e-01  52.00% -> 15.16s
  28: 0.827mm, [3.54e-05;1.01e+00], 1:0.635s, tol=1.74e-01  54.00% -> 14.60s
  29: 0.857mm, [3.16e-05;1.00e+00], 1:0.646s, tol=1.49e-01  56.00% -> 14.20s
  30: 0.888mm, [2.87e-05;1.00e+00], 1:0.656s, tol=1.39e-01  58.00% -> 13.77s
  31: 0.918mm, [2.61e-05;1.00e+00], 1:0.642s, tol=1.29e-01  60.00% -> 12.84s
  32: 0.949mm, [2.38e-05;1.00e+00], 1:0.640s, tol=1.31e-01  62.00% -> 12.16s
  33: 0.980mm, [2.20e-05;1.00e+00], 1:0.668s, tol=1.26e-01  64.00% -> 12.03s
  34: 1.010mm, [2.02e-05;1.00e+00], 1:0.672s, tol=1.12e-01  66.00% -> 11.43s
  35: 1.041mm, [1.86e-05;1.00e+00], 1:0.667s, tol=1.11e-01  68.00% -> 10.68s
  36: 1.071mm, [1.75e-05;1.00e+00], 1:0.771s, tol=1.09e-01  70.00% -> 11.56s
  37: 1.102mm, [1.64e-05;1.01e+00], 1:0.765s, tol=1.08e-01  72.00% -> 10.72s
  38: 1.133mm, [1.54e-05;1.01e+00], 1:0.771s, tol=1.09e-01  74.00% -> 10.03s
  39: 1.163mm, [1.48e-05;1.01e+00], 1:0.772s, tol=1.13e-01  76.00% -> 9.26s
  40: 1.194mm, [1.45e-05;1.01e+00], 1:0.783s, tol=1.24e-01  78.00% -> 8.62s
  41: 1.224mm, [1.46e-05;1.01e+00], 1:0.779s, tol=1.55e-01  80.00% -> 7.79s
  42: 1.255mm, [1.48e-05;1.00e+00], 1:0.782s, tol=4.20e-01  82.00% -> 7.04s
  43: 1.286mm, [1.54e-05;1.00e+00], 1:0.808s, tol=7.43e-01  84.00% -> 6.47s
  44: 1.316mm, [1.59e-05;1.00e+00], 1:0.828s, tol=1.25e-01  86.00% -> 5.80s
  45: 1.347mm, [1.64e-05;1.01e+00], 1:0.809s, tol=1.33e-01  88.00% -> 4.85s
  46: 1.378mm, [1.71e-05;1.00e+00], 1:0.817s, tol=1.92e-02  90.00% -> 4.09s
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.0/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.0/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100

Generate movie 01/23 (4.35 %) 3.54 s
Generate movie 02/23 (8.70 %) 3.23 s
Generate movie 03/23 (13.04 %) 3.09 s
Generate movie 04/23 (17.39 %) 2.92 s
Generate movie 05/23 (21.74 %) 2.76 s
Generate movie 06/23 (26.09 %) 2.62 s
Generate movie 07/23 (30.43 %) 2.48 s
Generate movie 08/23 (34.78 %) 2.33 s
Generate movie 09/23 (39.13 %) 2.15 s
Generate movie 10/23 (43.48 %) 2.04 s
Generate movie 11/23 (47.83 %) 1.87 s
Generate movie 12/23 (52.17 %) 1.71 s
Generate movie 13/23 (56.52 %) 1.56 s
Generate movie 14/23 (60.87 %) 1.41 s
Generate movie 15/23 (65.22 %) 1.24 s
Generate movie 16/23 (69.57 %) 1.08 s
Generate movie 17/23 (73.91 %) 935.36 ms
Generate movie 18/23 (78.26 %) 783.46 ms
Generate movie 19/23 (82.61 %) 618.58 ms
Generate movie 20/23 (86.96 %) 464.45 ms
Generate movie 21/23 (91.30 %) 314.56 ms
Generate movie 22/23 (95.65 %) 154.06 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, folder=folder_save)
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 38.604 seconds)

Gallery generated by Sphinx-Gallery