.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/PhaseField/PlateWithHole.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_PhaseField_PlateWithHole.py: PlateWithHole ============= Damage simulation for a plate with a hole subjected to compression. .. GENERATED FROM PYTHON SOURCE LINES 12-271 .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_001.gif :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_001.gif :class: sphx-glr-single-img .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_002.png :alt: $\phi$ :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_003.png :alt: TRI3: Ne = 12578, Nn = 6444 :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_004.png :alt: Boundary conditions :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_005.png :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_006.png :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_007.png :alt: Summary :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_007.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.0/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh 2: 0.000μm, [0.00e+00;0.00e+00], 1:1.094s, tol=0.00e+00 3: 0.800μm, [0.00e+00;0.00e+00], 1:1.072s, tol=1.00e+00 3.20% -> 1.08m 4: 1.600μm, [0.00e+00;0.00e+00], 1:1.071s, tol=7.50e-01 6.40% -> 46.99s 5: 2.400μm, [0.00e+00;0.00e+00], 1:1.051s, tol=5.56e-01 9.60% -> 39.59s 6: 3.200μm, [0.00e+00;0.00e+00], 1:1.073s, tol=4.37e-01 12.80% -> 36.56s 7: 4.000μm, [0.00e+00;0.00e+00], 1:1.056s, tol=3.60e-01 16.00% -> 33.26s 8: 4.800μm, [0.00e+00;0.00e+00], 1:1.085s, tol=3.06e-01 19.20% -> 31.98s 9: 5.600μm, [0.00e+00;0.00e+00], 1:1.051s, tol=2.65e-01 22.40% -> 29.13s 10: 6.400μm, [0.00e+00;0.00e+00], 1:1.076s, tol=2.34e-01 25.60% -> 28.14s 11: 7.200μm, [0.00e+00;0.00e+00], 1:1.059s, tol=2.10e-01 28.80% -> 26.19s 12: 8.000μm, [0.00e+00;0.00e+00], 1:1.073s, tol=1.90e-01 32.00% -> 25.08s 13: 8.800μm, [0.00e+00;0.00e+00], 1:1.054s, tol=1.74e-01 35.20% -> 23.27s 14: 9.600μm, [0.00e+00;0.00e+00], 1:1.073s, tol=1.60e-01 38.40% -> 22.38s 15: 10.400μm, [0.00e+00;0.00e+00], 1:1.054s, tol=1.48e-01 41.60% -> 20.72s 16: 11.200μm, [0.00e+00;0.00e+00], 1:1.074s, tol=1.38e-01 44.80% -> 19.85s 17: 12.000μm, [0.00e+00;0.00e+00], 1:1.068s, tol=1.29e-01 48.00% -> 18.51s 18: 12.800μm, [0.00e+00;0.00e+00], 1:1.078s, tol=1.21e-01 51.20% -> 17.47s 19: 13.600μm, [0.00e+00;0.00e+00], 1:1.053s, tol=1.14e-01 54.40% -> 15.88s 20: 14.400μm, [0.00e+00;0.00e+00], 1:1.090s, tol=1.08e-01 57.60% -> 15.24s 21: 15.200μm, [0.00e+00;0.00e+00], 1:1.056s, tol=1.02e-01 60.80% -> 13.61s 22: 16.000μm, [0.00e+00;0.00e+00], 1:1.072s, tol=9.75e-02 64.00% -> 12.66s 23: 16.800μm, [9.79e-13;4.17e-03], 1:1.068s, tol=9.30e-02 67.20% -> 11.46s 24: 17.600μm, [1.02e-11;4.36e-02], 1:1.101s, tol=8.88e-02 70.40% -> 10.64s 25: 18.400μm, [1.09e-11;1.18e-01], 1:1.082s, tol=8.49e-02 73.60% -> 9.32s 26: 19.200μm, [7.24e-12;2.14e-01], 1:1.074s, tol=8.14e-02 76.80% -> 8.11s 27: 20.000μm, [4.04e-12;3.31e-01], 1:1.076s, tol=7.82e-02 80.00% -> 6.99s 28: 20.800μm, [2.70e-12;5.03e-01], 1:1.084s, tol=7.50e-02 83.20% -> 5.91s 29: 21.600μm, [1.48e-12;7.38e-01], 1:1.046s, tol=7.24e-02 86.40% -> 4.61s 30: 21.800μm, [1.01e-12;9.34e-01], 1:1.079s, tol=2.12e-02 87.20% -> 4.59s 31: 22.000μm, [1.11e-12;9.98e-01], 1:1.066s, tol=3.65e-02 88.00% -> 4.36s 32: 22.200μm, [1.22e-12;1.01e+00], 1:1.084s, tol=6.20e-02 88.80% -> 4.24s 33: 22.400μm, [1.36e-12;1.00e+00], 1:1.070s, tol=8.53e-02 89.60% -> 3.97s 34: 22.600μm, [1.35e-12;1.01e+00], 1:1.077s, tol=9.21e-02 90.40% -> 3.77s 35: 22.800μm, [1.44e-12;1.00e+00], 1:1.055s, tol=1.06e-01 91.20% -> 3.46s 36: 23.000μm, [1.49e-12;1.00e+00], 1:1.073s, tol=1.02e-01 92.00% -> 3.27s 37: 23.200μm, [1.63e-12;1.00e+00], 1:1.059s, tol=9.02e-02 92.80% -> 2.96s 38: 23.400μm, [1.67e-12;1.00e+00], 1:1.083s, tol=9.63e-02 93.60% -> 2.74s 39: 23.600μm, [1.75e-12;1.00e+00], 1:1.055s, tol=1.02e-01 94.40% -> 2.38s 40: 23.800μm, [1.79e-12;1.00e+00], 1:1.072s, tol=1.03e-01 95.20% -> 2.11s 41: 24.000μm, [1.68e-12;1.00e+00], 1:1.059s, tol=1.01e-01 96.00% -> 1.76s 42: 24.200μm, [1.43e-12;1.00e+00], 1:1.073s, tol=9.89e-02 96.80% -> 1.45s 43: 24.400μm, [1.26e-12;1.00e+00], 1:1.070s, tol=1.01e-01 97.60% -> 1.10s 44: 24.600μm, [1.10e-12;1.00e+00], 1:1.077s, tol=9.77e-02 98.40% -> 753.31ms 45: 24.800μm, [1.02e-12;1.00e+00], 1:1.071s, tol=9.67e-02 99.20% -> 380.15ms 46: 25.000μm, [9.35e-13;1.00e+00], 1:1.079s, tol=9.35e-02 100.00% -> -0.00µs Saved mesh data in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.0/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh/Meshes Saved simulation and summary in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.0/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh Generate movie 01/23 (4.35 %) 3.54 s Generate movie 02/23 (8.70 %) 3.31 s Generate movie 03/23 (13.04 %) 3.14 s Generate movie 04/23 (17.39 %) 2.98 s Generate movie 05/23 (21.74 %) 2.82 s Generate movie 06/23 (26.09 %) 2.69 s Generate movie 07/23 (30.43 %) 2.51 s Generate movie 08/23 (34.78 %) 2.41 s Generate movie 09/23 (39.13 %) 2.20 s Generate movie 10/23 (43.48 %) 2.04 s Generate movie 11/23 (47.83 %) 1.97 s Generate movie 12/23 (52.17 %) 1.74 s Generate movie 13/23 (56.52 %) 1.57 s Generate movie 14/23 (60.87 %) 1.43 s Generate movie 15/23 (65.22 %) 1.26 s Generate movie 16/23 (69.57 %) 1.12 s Generate movie 17/23 (73.91 %) 954.07 ms Generate movie 18/23 (78.26 %) 795.82 ms Generate movie 19/23 (82.61 %) 635.25 ms Generate movie 20/23 (86.96 %) 492.82 ms Generate movie 21/23 (91.30 %) 319.25 ms Generate movie 22/23 (95.65 %) 159.98 ms Generate movie 23/23 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 13 import matplotlib.pyplot as plt import numpy as np from EasyFEA import ( Display, Folder, Models, Tic, ElemType, Simulations, PyVista, Paraview, ) from EasyFEA.Geoms import Domain, Circle if __name__ == "__main__": # ---------------------------------------------- # Configuration # ---------------------------------------------- # simu options doSimu = True meshTest = True optimMesh = True # outputs folder = Folder.Results_Dir() plotMesh = False plotIter = False plotEnergy = False makeParaview = False makeMovie = True # Available splits: Bourdin, Amor, Miehe, Stress (isotropic) # He, AnisotStrain, AnisotStress, Zhang (anisotropic) split = Models.PhaseField.SplitType.Miehe # Available regus: AT1, AT2 regu = Models.PhaseField.ReguType.AT1 l0 = 0.12e-3 # m # convergence solver = ( Models.PhaseField.SolverType.History ) # History, HistoryDamage, BoundConstrain maxIter = 1000 tolConv = 1e-0 # load units unitU = "μm" unitF = "kN/mm" unit = 1e6 # ---------------------------------------------- # Geometry # ---------------------------------------------- L = 15e-3 # m h = 30e-3 # m diam = 6e-3 # m thickness = 1 # ---------------------------------------------- # Material # ---------------------------------------------- E = 12e9 # Pa v = 0.3 Gc = 1.4 # J/m2 folder_save = Simulations.PhaseField.Folder( folder, "", split, regu, "", tolConv, solver, meshTest, optimMesh ) Display.MyPrint(folder_save, "green", end="\n") if doSimu: # ---------------------------------------------- # Mesh # ---------------------------------------------- clC = l0 * 2 if meshTest else l0 / 2 if optimMesh: clD = l0 * 4 refineZone = diam * 1.5 / 2 if split in ( Models.PhaseField.SplitType.Bourdin, Models.PhaseField.SplitType.Amor, ): refineGeom = Domain( (0, h / 2 - refineZone), (L, h / 2 + refineZone), clC ) else: refineGeom = Domain( (L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC ) else: clD = l0 if meshTest else l0 / 2 refineGeom = None domain = Domain((0, 0), (L, h), clD) circle = Circle((L / 2, h / 2), diam, clD) mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom]) # Nodes nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0) nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h) nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0)) nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"]) # ---------------------------------------------- # Material # ---------------------------------------------- material = Models.Elastic.Isotropic( 2, E, v, planeStress=False, thickness=thickness ) pfm = Models.PhaseField(material, split, regu, Gc, l0, solver=solver) # ---------------------------------------------- # Boundary conditions # ---------------------------------------------- threshold = 0.6 u_max = 25e-6 uinc0 = 8e-7 if meshTest else 8e-8 uinc1 = 2e-7 if meshTest else 2e-8 config = f""" uinc0 = {uinc0:.1e}; uinc1 = {uinc1:.1e} while ud <= u_max = {u_max:.2e}: ud += uinc0 if simu.damage.max() < {threshold} else uinc1 simu.add_dirichlet(nodes_lower, [0], ["y"]) simu.add_dirichlet(nodes_x0y0, [0], ["x"]) simu.add_dirichlet(nodes_upper, [-ud], ["y"]) """ # ---------------------------------------------- # Simulation # ---------------------------------------------- simu = Simulations.PhaseField(mesh, pfm, folder=folder_save) simu.Results_Set_Bc_Summary(config) dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"]) def Apply_BC(ud: float): simu.Bc_Init() simu.add_dirichlet(nodes_lower, [0], ["y"]) simu.add_dirichlet(nodes_x0y0, [0], ["x"]) simu.add_dirichlet(nodes_upper, [-ud], ["y"]) list_dep = [] list_f = [] ud = -uinc0 iter = 0 nDetect = 0 if plotIter: axIter = Display.Plot_Result(simu, "damage", nodeValues=True) _, axLoad = plt.subplots() axLoad.set_xlabel(f"ud [{unitU}]") axLoad.set_ylabel(f"f [{unitF}]") axLoad.grid() while ud <= u_max: iter += 1 ud += uinc0 if simu.damage.max() < threshold else uinc1 Apply_BC(ud) u, d, convergence = simu.Solve(tolConv, maxIter) simu.Save_Iter() if not convergence: break f = np.sum(simu.Calc_Reaction(dofsY_upper, "elastic")) simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True) if simu.Detect_Damage(nodes_edges, 1): nDetect += 1 if nDetect == 10: break list_dep.append(ud) list_f.append(f) if plotIter: Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter) plt.figure(axIter.figure) plt.pause(1e-12) axLoad.clear() axLoad.plot(np.abs(list_dep * unit), np.abs(list_f / unit), c="blue") axLoad.set_xlabel(f"ud [{unitU}]") axLoad.set_ylabel(f"f [{unitF}]") axLoad.grid() plt.figure(axLoad.figure) plt.pause(1e-12) # ---------------------------------------------- # Saving # ---------------------------------------------- print() Simulations.Save_pickle((list_f, list_dep), folder_save, "force-displacement") simu.Save(folder_save) else: simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save) list_f, list_dep = Simulations.Load_pickle(folder_save, "force-displacement") # ---------------------------------------------- # Results # ---------------------------------------------- if plotEnergy: Display.Plot_Energy(simu, list_f, list_dep, N=400, folder=folder_save) Display.Plot_Result( simu, "damage", nodeValues=True, colorbarIsClose=True, folder=folder_save, filename="damage", ) Display.Plot_Mesh(simu) Display.Plot_BoundaryConditions(simu) Display.Plot_Iter_Summary(simu, folder_save, None, None) ax = Display.Init_Axes() ax.plot(np.abs(list_dep) * unit, np.abs(list_f) / unit, c="blue") ax.set_xlabel(f"ud [{unitU}]") ax.set_ylabel(f"f [{unitF}]") ax.grid() Display.Save_fig(folder_save, "force-displacement") if plotMesh: Display.Plot_Mesh(simu.mesh) if makeParaview: Paraview.Save_simu(simu, folder_save) if makeMovie: simu.Set_Iter(-1) deformFactor = L * 0.05 / simu.Result("displacement_norm").max() iterations = np.arange(0, simu.Niter, simu.Niter // 20) def Func(plotter, iter): simu.Set_Iter(iterations[iter]) thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8)) PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1)) PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif") if doSimu: Tic.Plot_History(folder_save, details=False) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 55.665 seconds) .. _sphx_glr_download_examples_PhaseField_PlateWithHole.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: PlateWithHole.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: PlateWithHole.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: PlateWithHole.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_