Note

Go to the end to download the full example code.

CT#

Performs damage simulation on a CT specimen.

CT
  • $\phi$
  • $uy$
  • elemType: TRI3, Ne: 8982, Nn: 4661
  • CT
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v3.0.0/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
   1: 0.000mm, [0.00e+00;0.00e+00], 1:0.112s, tol=0.00e+00
   2: 0.031mm, [0.00e+00;0.00e+00], 1:0.090s, tol=1.00e+00  2.00% -> 4.42s
   3: 0.061mm, [0.00e+00;0.00e+00], 1:0.091s, tol=7.50e-01  4.00% -> 4.36s
   4: 0.092mm, [0.00e+00;0.00e+00], 1:0.091s, tol=5.56e-01  6.00% -> 4.28s
   5: 0.122mm, [0.00e+00;0.00e+00], 1:0.090s, tol=4.37e-01  8.00% -> 4.13s
   6: 0.153mm, [0.00e+00;0.00e+00], 1:0.090s, tol=3.60e-01  10.00% -> 4.03s
   7: 0.184mm, [0.00e+00;0.00e+00], 1:0.091s, tol=3.06e-01  12.00% -> 4.02s
   8: 0.214mm, [0.00e+00;0.00e+00], 1:0.090s, tol=2.65e-01  14.00% -> 3.86s
   9: 0.245mm, [1.63e-04;2.46e-02], 1:0.091s, tol=2.32e-01  16.00% -> 3.80s
  10: 0.276mm, [2.76e-04;8.92e-02], 1:0.092s, tol=2.04e-01  18.00% -> 3.75s
  11: 0.306mm, [3.39e-04;1.99e-01], 1:0.092s, tol=1.81e-01  20.00% -> 3.66s
  12: 0.337mm, [3.21e-04;3.42e-01], 1:0.093s, tol=1.65e-01  22.00% -> 3.62s
  13: 0.367mm, [2.96e-04;5.27e-01], 1:0.091s, tol=1.52e-01  24.00% -> 3.46s
  14: 0.398mm, [2.68e-04;7.44e-01], 1:0.091s, tol=1.59e-01  26.00% -> 3.38s
  15: 0.429mm, [2.44e-04;9.33e-01], 1:0.091s, tol=2.32e-01  28.00% -> 3.28s
  16: 0.459mm, [2.15e-04;1.02e+00], 1:0.091s, tol=3.72e-01  30.00% -> 3.19s
  17: 0.490mm, [1.82e-04;1.02e+00], 1:0.092s, tol=4.13e-01  32.00% -> 3.12s
  18: 0.520mm, [1.52e-04;1.01e+00], 1:0.091s, tol=3.98e-01  34.00% -> 3.01s
  19: 0.551mm, [1.24e-04;1.01e+00], 1:0.092s, tol=2.99e-01  36.00% -> 2.95s
  20: 0.582mm, [1.03e-04;1.01e+00], 1:0.091s, tol=2.73e-01  38.00% -> 2.82s
  21: 0.612mm, [8.73e-05;1.01e+00], 1:0.092s, tol=2.61e-01  40.00% -> 2.75s
  22: 0.643mm, [7.52e-05;1.01e+00], 1:0.094s, tol=2.47e-01  42.00% -> 2.71s
  23: 0.673mm, [6.58e-05;1.01e+00], 1:0.092s, tol=2.42e-01  44.00% -> 2.57s
  24: 0.704mm, [5.74e-05;1.00e+00], 1:0.093s, tol=2.07e-01  46.00% -> 2.51s
  25: 0.735mm, [5.02e-05;1.00e+00], 1:0.091s, tol=2.11e-01  48.00% -> 2.37s
  26: 0.765mm, [4.43e-05;1.00e+00], 1:0.093s, tol=1.77e-01  50.00% -> 2.31s
  27: 0.796mm, [3.96e-05;1.01e+00], 1:0.091s, tol=1.67e-01  52.00% -> 2.20s
  28: 0.827mm, [3.54e-05;1.01e+00], 1:0.092s, tol=1.74e-01  54.00% -> 2.11s
  29: 0.857mm, [3.16e-05;1.00e+00], 1:0.091s, tol=1.49e-01  56.00% -> 2.01s
  30: 0.888mm, [2.87e-05;1.00e+00], 1:0.092s, tol=1.39e-01  58.00% -> 1.92s
  31: 0.918mm, [2.61e-05;1.00e+00], 1:0.092s, tol=1.29e-01  60.00% -> 1.83s
  32: 0.949mm, [2.38e-05;1.00e+00], 1:0.091s, tol=1.31e-01  62.00% -> 1.73s
  33: 0.980mm, [2.20e-05;1.00e+00], 1:0.092s, tol=1.26e-01  64.00% -> 1.66s
  34: 1.010mm, [2.02e-05;1.00e+00], 1:0.092s, tol=1.12e-01  66.00% -> 1.56s
  35: 1.041mm, [1.86e-05;1.00e+00], 1:0.092s, tol=1.11e-01  68.00% -> 1.47s
  36: 1.071mm, [1.75e-05;1.00e+00], 1:0.093s, tol=1.09e-01  70.00% -> 1.39s
  37: 1.102mm, [1.64e-05;1.01e+00], 1:0.091s, tol=1.08e-01  72.00% -> 1.27s
  38: 1.133mm, [1.54e-05;1.01e+00], 1:0.094s, tol=1.09e-01  74.00% -> 1.22s
  39: 1.163mm, [1.48e-05;1.01e+00], 1:0.092s, tol=1.13e-01  76.00% -> 1.10s
  40: 1.194mm, [1.45e-05;1.01e+00], 1:0.092s, tol=1.24e-01  78.00% -> 1.02s
  41: 1.224mm, [1.46e-05;1.01e+00], 1:0.092s, tol=1.55e-01  80.00% -> 916.04ms
  42: 1.255mm, [1.48e-05;1.00e+00], 1:0.093s, tol=4.20e-01  82.00% -> 835.96ms
  43: 1.286mm, [1.54e-05;1.00e+00], 1:0.092s, tol=7.43e-01  84.00% -> 737.98ms
  44: 1.316mm, [1.59e-05;1.00e+00], 1:0.094s, tol=1.25e-01  86.00% -> 658.32ms
  45: 1.347mm, [1.64e-05;1.01e+00], 1:0.094s, tol=1.33e-01  88.00% -> 561.81ms
  46: 1.378mm, [1.71e-05;1.00e+00], 1:0.093s, tol=1.92e-02  90.00% -> 465.44ms
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v3.0.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/v3.0.0/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100

Generate movie 01/23 (4.35 %) 3.48 s
Generate movie 02/23 (8.70 %) 3.18 s
Generate movie 03/23 (13.04 %) 3.00 s
Generate movie 04/23 (17.39 %) 2.86 s
Generate movie 05/23 (21.74 %) 2.72 s
Generate movie 06/23 (26.09 %) 2.55 s
Generate movie 07/23 (30.43 %) 2.41 s
Generate movie 08/23 (34.78 %) 2.27 s
Generate movie 09/23 (39.13 %) 2.09 s
Generate movie 10/23 (43.48 %) 1.97 s
Generate movie 11/23 (47.83 %) 1.81 s
Generate movie 12/23 (52.17 %) 1.67 s
Generate movie 13/23 (56.52 %) 1.50 s
Generate movie 14/23 (60.87 %) 1.36 s
Generate movie 15/23 (65.22 %) 1.20 s
Generate movie 16/23 (69.57 %) 1.04 s
Generate movie 17/23 (73.91 %) 898.08 ms
Generate movie 18/23 (78.26 %) 745.66 ms
Generate movie 19/23 (82.61 %) 595.48 ms
Generate movie 20/23 (86.96 %) 446.57 ms
Generate movie 21/23 (91.30 %) 298.25 ms
Generate movie 22/23 (95.65 %) 148.40 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 Terminal, Matplotlib, Folder, Models, ElemType, Simulations, PyVista, Paraview
 17 from EasyFEA.Geoms import Point, Points, Circle, Line, Contour, Domain
 18
 19 if __name__ == "__main__":
 20     Terminal.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     Terminal.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 = Matplotlib.Plot(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                 Matplotlib.Plot(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     Matplotlib.Plot(simu, "damage", folder=folder_save)
234     Matplotlib.Plot(simu, "uy", deformFactor=1)
235     Matplotlib.Plot_Mesh(mesh)
236     Matplotlib.Plot_Tags(mesh, folder=folder_save)
237
238     plt.show()

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

Gallery generated by Sphinx-Gallery