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.5/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
   1: 0.000mm, [0.00e+00;0.00e+00], 1:0.438s, tol=0.00e+00
   2: 0.031mm, [0.00e+00;0.00e+00], 1:0.442s, tol=1.00e+00  2.00% -> 21.68s
   3: 0.061mm, [0.00e+00;0.00e+00], 1:0.450s, tol=7.50e-01  4.00% -> 21.61s
   4: 0.092mm, [0.00e+00;0.00e+00], 1:0.441s, tol=5.56e-01  6.00% -> 20.74s
   5: 0.122mm, [0.00e+00;0.00e+00], 1:0.445s, tol=4.37e-01  8.00% -> 20.48s
   6: 0.153mm, [0.00e+00;0.00e+00], 1:0.442s, tol=3.60e-01  10.00% -> 19.89s
   7: 0.184mm, [0.00e+00;0.00e+00], 1:0.449s, tol=3.06e-01  12.00% -> 19.74s
   8: 0.214mm, [0.00e+00;0.00e+00], 1:0.442s, tol=2.65e-01  14.00% -> 19.02s
   9: 0.245mm, [1.63e-04;2.46e-02], 1:0.444s, tol=2.32e-01  16.00% -> 18.64s
  10: 0.276mm, [2.76e-04;8.92e-02], 1:0.441s, tol=2.04e-01  18.00% -> 18.09s
  11: 0.306mm, [3.39e-04;1.99e-01], 1:0.450s, tol=1.81e-01  20.00% -> 18.01s
  12: 0.337mm, [3.21e-04;3.42e-01], 1:0.447s, tol=1.65e-01  22.00% -> 17.42s
  13: 0.367mm, [2.96e-04;5.27e-01], 1:0.444s, tol=1.52e-01  24.00% -> 16.88s
  14: 0.398mm, [2.68e-04;7.44e-01], 1:0.440s, tol=1.59e-01  26.00% -> 16.29s
  15: 0.429mm, [2.44e-04;9.33e-01], 1:0.450s, tol=2.32e-01  28.00% -> 16.20s
  16: 0.459mm, [2.15e-04;1.02e+00], 1:0.449s, tol=3.72e-01  30.00% -> 15.73s
  17: 0.490mm, [1.82e-04;1.02e+00], 1:0.446s, tol=4.13e-01  32.00% -> 15.15s
  18: 0.520mm, [1.52e-04;1.01e+00], 1:0.442s, tol=3.98e-01  34.00% -> 14.59s
  19: 0.551mm, [1.24e-04;1.01e+00], 1:0.454s, tol=2.99e-01  36.00% -> 14.53s
  20: 0.582mm, [1.03e-04;1.01e+00], 1:0.444s, tol=2.73e-01  38.00% -> 13.76s
  21: 0.612mm, [8.73e-05;1.01e+00], 1:0.445s, tol=2.61e-01  40.00% -> 13.36s
  22: 0.643mm, [7.52e-05;1.01e+00], 1:0.445s, tol=2.47e-01  42.00% -> 12.90s
  23: 0.673mm, [6.58e-05;1.01e+00], 1:0.448s, tol=2.42e-01  44.00% -> 12.54s
  24: 0.704mm, [5.74e-05;1.00e+00], 1:0.447s, tol=2.07e-01  46.00% -> 12.06s
  25: 0.735mm, [5.02e-05;1.00e+00], 1:0.471s, tol=2.11e-01  48.00% -> 12.25s
  26: 0.765mm, [4.43e-05;1.00e+00], 1:0.467s, tol=1.77e-01  50.00% -> 11.66s
  27: 0.796mm, [3.96e-05;1.01e+00], 1:0.475s, tol=1.67e-01  52.00% -> 11.40s
  28: 0.827mm, [3.54e-05;1.01e+00], 1:0.474s, tol=1.74e-01  54.00% -> 10.91s
  29: 0.857mm, [3.16e-05;1.00e+00], 1:0.475s, tol=1.49e-01  56.00% -> 10.44s
  30: 0.888mm, [2.87e-05;1.00e+00], 1:0.465s, tol=1.39e-01  58.00% -> 9.77s
  31: 0.918mm, [2.61e-05;1.00e+00], 1:0.478s, tol=1.29e-01  60.00% -> 9.55s
  32: 0.949mm, [2.38e-05;1.00e+00], 1:0.475s, tol=1.31e-01  62.00% -> 9.03s
  33: 0.980mm, [2.20e-05;1.00e+00], 1:0.500s, tol=1.26e-01  64.00% -> 9.01s
  34: 1.010mm, [2.02e-05;1.00e+00], 1:0.497s, tol=1.12e-01  66.00% -> 8.45s
  35: 1.041mm, [1.86e-05;1.00e+00], 1:0.504s, tol=1.11e-01  68.00% -> 8.06s
  36: 1.071mm, [1.75e-05;1.00e+00], 1:0.591s, tol=1.09e-01  70.00% -> 8.86s
  37: 1.102mm, [1.64e-05;1.01e+00], 1:0.594s, tol=1.08e-01  72.00% -> 8.31s
  38: 1.133mm, [1.54e-05;1.01e+00], 1:0.589s, tol=1.09e-01  74.00% -> 7.66s
  39: 1.163mm, [1.48e-05;1.01e+00], 1:0.598s, tol=1.13e-01  76.00% -> 7.18s
  40: 1.194mm, [1.45e-05;1.01e+00], 1:0.597s, tol=1.24e-01  78.00% -> 6.56s
  41: 1.224mm, [1.46e-05;1.01e+00], 1:0.596s, tol=1.55e-01  80.00% -> 5.96s
  42: 1.255mm, [1.48e-05;1.00e+00], 1:0.598s, tol=4.20e-01  82.00% -> 5.38s
  43: 1.286mm, [1.54e-05;1.00e+00], 1:0.620s, tol=7.43e-01  84.00% -> 4.96s
  44: 1.316mm, [1.59e-05;1.00e+00], 1:0.630s, tol=1.25e-01  86.00% -> 4.41s
  45: 1.347mm, [1.64e-05;1.01e+00], 1:0.625s, tol=1.33e-01  88.00% -> 3.75s
  46: 1.378mm, [1.71e-05;1.00e+00], 1:0.627s, tol=1.92e-02  90.00% -> 3.13s
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/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.8.5/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100

Generate movie 01/23 (4.35 %) 2.84 s
Generate movie 02/23 (8.70 %) 2.21 s
Generate movie 03/23 (13.04 %) 2.09 s
Generate movie 04/23 (17.39 %) 1.95 s
Generate movie 05/23 (21.74 %) 1.83 s
Generate movie 06/23 (26.09 %) 1.74 s
Generate movie 07/23 (30.43 %) 1.64 s
Generate movie 08/23 (34.78 %) 1.56 s
Generate movie 09/23 (39.13 %) 1.41 s
Generate movie 10/23 (43.48 %) 1.32 s
Generate movie 11/23 (47.83 %) 1.22 s
Generate movie 12/23 (52.17 %) 1.19 s
Generate movie 13/23 (56.52 %) 1.05 s
Generate movie 14/23 (60.87 %) 964.45 ms
Generate movie 15/23 (65.22 %) 839.84 ms
Generate movie 16/23 (69.57 %) 735.61 ms
Generate movie 17/23 (73.91 %) 646.79 ms
Generate movie 18/23 (78.26 %) 529.15 ms
Generate movie 19/23 (82.61 %) 440.93 ms
Generate movie 20/23 (86.96 %) 319.21 ms
Generate movie 21/23 (91.30 %) 211.63 ms
Generate movie 22/23 (95.65 %) 107.92 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 28.864 seconds)

Gallery generated by Sphinx-Gallery