.. 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.1/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh 2: 0.000μm, [0.00e+00;0.00e+00], 1:0.862s, tol=0.00e+00 3: 0.800μm, [0.00e+00;0.00e+00], 1:0.832s, tol=1.00e+00 3.20% -> 50.31s 4: 1.600μm, [0.00e+00;0.00e+00], 1:0.835s, tol=7.50e-01 6.40% -> 36.63s 5: 2.400μm, [0.00e+00;0.00e+00], 1:0.846s, tol=5.56e-01 9.60% -> 31.87s 6: 3.200μm, [0.00e+00;0.00e+00], 1:0.840s, tol=4.38e-01 12.80% -> 28.62s 7: 4.000μm, [0.00e+00;0.00e+00], 1:0.838s, tol=3.60e-01 16.00% -> 26.41s 8: 4.800μm, [0.00e+00;0.00e+00], 1:0.840s, tol=3.06e-01 19.20% -> 24.75s 9: 5.600μm, [0.00e+00;0.00e+00], 1:0.839s, tol=2.65e-01 22.40% -> 23.26s 10: 6.400μm, [0.00e+00;0.00e+00], 1:0.835s, tol=2.34e-01 25.60% -> 21.83s 11: 7.200μm, [0.00e+00;0.00e+00], 1:0.837s, tol=2.10e-01 28.80% -> 20.68s 12: 8.000μm, [0.00e+00;0.00e+00], 1:0.835s, tol=1.90e-01 32.00% -> 19.51s 13: 8.800μm, [0.00e+00;0.00e+00], 1:0.831s, tol=1.74e-01 35.20% -> 18.36s 14: 9.600μm, [0.00e+00;0.00e+00], 1:0.832s, tol=1.60e-01 38.40% -> 17.35s 15: 10.400μm, [0.00e+00;0.00e+00], 1:0.831s, tol=1.48e-01 41.60% -> 16.33s 16: 11.200μm, [0.00e+00;0.00e+00], 1:0.834s, tol=1.38e-01 44.80% -> 15.41s 17: 12.000μm, [0.00e+00;0.00e+00], 1:0.835s, tol=1.29e-01 48.00% -> 14.47s 18: 12.800μm, [0.00e+00;0.00e+00], 1:0.831s, tol=1.21e-01 51.20% -> 13.46s 19: 13.600μm, [0.00e+00;0.00e+00], 1:0.833s, tol=1.14e-01 54.40% -> 12.57s 20: 14.400μm, [0.00e+00;0.00e+00], 1:0.836s, tol=1.08e-01 57.60% -> 11.69s 21: 15.200μm, [0.00e+00;0.00e+00], 1:0.834s, tol=1.02e-01 60.80% -> 10.76s 22: 16.000μm, [0.00e+00;0.00e+00], 1:0.835s, tol=9.75e-02 64.00% -> 9.86s 23: 16.800μm, [9.79e-13;4.17e-03], 1:0.832s, tol=9.30e-02 67.20% -> 8.94s 24: 17.600μm, [1.02e-11;4.36e-02], 1:0.836s, tol=8.88e-02 70.40% -> 8.08s 25: 18.400μm, [1.09e-11;1.18e-01], 1:0.846s, tol=8.49e-02 73.60% -> 7.28s 26: 19.200μm, [7.24e-12;2.14e-01], 1:0.836s, tol=8.14e-02 76.80% -> 6.32s 27: 20.000μm, [4.04e-12;3.31e-01], 1:0.833s, tol=7.82e-02 80.00% -> 5.41s 28: 20.800μm, [2.70e-12;5.03e-01], 1:0.833s, tol=7.50e-02 83.20% -> 4.54s 29: 21.600μm, [1.48e-12;7.38e-01], 1:0.832s, tol=7.24e-02 86.40% -> 3.67s 30: 21.800μm, [1.01e-12;9.34e-01], 1:0.836s, tol=2.12e-02 87.20% -> 3.56s 31: 22.000μm, [1.11e-12;9.98e-01], 1:0.834s, tol=3.65e-02 88.00% -> 3.41s 32: 22.200μm, [1.22e-12;1.01e+00], 1:0.834s, tol=6.20e-02 88.80% -> 3.26s 33: 22.400μm, [1.36e-12;1.00e+00], 1:0.839s, tol=8.53e-02 89.60% -> 3.12s 34: 22.600μm, [1.35e-12;1.01e+00], 1:0.835s, tol=9.21e-02 90.40% -> 2.93s 35: 22.800μm, [1.44e-12;1.00e+00], 1:0.838s, tol=1.06e-01 91.20% -> 2.75s 36: 23.000μm, [1.49e-12;1.00e+00], 1:0.834s, tol=1.02e-01 92.00% -> 2.54s 37: 23.200μm, [1.63e-12;1.00e+00], 1:0.833s, tol=9.02e-02 92.80% -> 2.33s 38: 23.400μm, [1.67e-12;1.00e+00], 1:0.835s, tol=9.63e-02 93.60% -> 2.11s 39: 23.600μm, [1.75e-12;1.00e+00], 1:0.835s, tol=1.02e-01 94.40% -> 1.88s 40: 23.800μm, [1.79e-12;1.00e+00], 1:0.832s, tol=1.03e-01 95.20% -> 1.64s 41: 24.000μm, [1.68e-12;1.00e+00], 1:0.838s, tol=1.01e-01 96.00% -> 1.40s 42: 24.200μm, [1.43e-12;1.00e+00], 1:0.833s, tol=9.89e-02 96.80% -> 1.13s 43: 24.400μm, [1.26e-12;1.00e+00], 1:0.845s, tol=1.01e-01 97.60% -> 872.98ms 44: 24.600μm, [1.10e-12;1.00e+00], 1:0.833s, tol=9.77e-02 98.40% -> 582.57ms 45: 24.800μm, [1.02e-12;1.00e+00], 1:0.839s, tol=9.67e-02 99.20% -> 297.64ms 46: 25.000μm, [9.35e-13;1.00e+00], 1:0.837s, tol=9.35e-02 100.00% -> -0.00µs Saved mesh data in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.1/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.1/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh Generate movie 01/23 (4.35 %) 2.67 s Generate movie 02/23 (8.70 %) 2.24 s Generate movie 03/23 (13.04 %) 2.12 s Generate movie 04/23 (17.39 %) 2.06 s Generate movie 05/23 (21.74 %) 1.93 s Generate movie 06/23 (26.09 %) 1.81 s Generate movie 07/23 (30.43 %) 1.71 s Generate movie 08/23 (34.78 %) 1.60 s Generate movie 09/23 (39.13 %) 1.52 s Generate movie 10/23 (43.48 %) 1.37 s Generate movie 11/23 (47.83 %) 1.28 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 %) 951.06 ms Generate movie 15/23 (65.22 %) 860.32 ms Generate movie 16/23 (69.57 %) 841.49 ms Generate movie 17/23 (73.91 %) 651.04 ms Generate movie 18/23 (78.26 %) 537.90 ms Generate movie 19/23 (82.61 %) 427.68 ms Generate movie 20/23 (86.96 %) 329.83 ms Generate movie 21/23 (91.30 %) 217.79 ms Generate movie 22/23 (95.65 %) 106.49 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 43.486 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 `_