.. 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-272 .. 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.8.0/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh 2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:0.887s, tol=0.00e+00 3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.840s, tol=1.00e+00 3.20% -> 50.79s 4 : 1.600 μm, [0.00e+00; 0.00e+00], 1:0.858s, tol=8.04e-01 6.40% -> 37.66s 5 : 2.400 μm, [0.00e+00; 0.00e+00], 1:0.838s, tol=5.56e-01 9.60% -> 31.56s 6 : 3.200 μm, [0.00e+00; 0.00e+00], 1:0.862s, tol=4.38e-01 12.80% -> 29.38s 7 : 4.000 μm, [0.00e+00; 0.00e+00], 1:0.838s, tol=3.60e-01 16.00% -> 26.39s 8 : 4.800 μm, [0.00e+00; 0.00e+00], 1:0.859s, tol=3.06e-01 19.20% -> 25.31s 9 : 5.600 μm, [0.00e+00; 0.00e+00], 1:0.836s, tol=2.65e-01 22.40% -> 23.17s 10 : 6.400 μm, [0.00e+00; 0.00e+00], 1:0.858s, tol=2.34e-01 25.60% -> 22.45s 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.863s, tol=1.90e-01 32.00% -> 20.18s 13 : 8.800 μm, [0.00e+00; 0.00e+00], 1:0.841s, tol=1.74e-01 35.20% -> 18.58s 14 : 9.600 μm, [0.00e+00; 0.00e+00], 1:0.858s, tol=1.60e-01 38.40% -> 17.90s 15 : 10.400 μm, [0.00e+00; 0.00e+00], 1:0.838s, tol=1.48e-01 41.60% -> 16.48s 16 : 11.200 μm, [0.00e+00; 0.00e+00], 1:0.858s, tol=1.38e-01 44.80% -> 15.87s 17 : 12.000 μm, [0.00e+00; 0.00e+00], 1:0.845s, tol=1.29e-01 48.00% -> 14.64s 18 : 12.800 μm, [0.00e+00; 0.00e+00], 1:0.855s, tol=1.21e-01 51.20% -> 13.85s 19 : 13.600 μm, [0.00e+00; 0.00e+00], 1:0.836s, tol=1.14e-01 54.40% -> 12.62s 20 : 14.400 μm, [0.00e+00; 0.00e+00], 1:0.859s, tol=1.08e-01 57.60% -> 12.02s 21 : 15.200 μm, [0.00e+00; 0.00e+00], 1:0.851s, tol=1.02e-01 60.80% -> 10.98s 22 : 16.000 μm, [0.00e+00; 0.00e+00], 1:0.864s, tol=9.75e-02 64.00% -> 10.20s 23 : 16.800 μm, [9.79e-13; 4.17e-03], 1:0.837s, tol=9.30e-02 67.20% -> 8.98s 24 : 17.600 μm, [1.02e-11; 4.36e-02], 1:0.860s, tol=8.88e-02 70.40% -> 8.32s 25 : 18.400 μm, [1.09e-11; 1.18e-01], 1:0.838s, tol=8.49e-02 73.60% -> 7.22s 26 : 19.200 μm, [7.24e-12; 2.14e-01], 1:0.859s, tol=8.14e-02 76.80% -> 6.49s 27 : 20.000 μm, [4.04e-12; 3.31e-01], 1:0.841s, tol=7.82e-02 80.00% -> 5.47s 28 : 20.800 μm, [2.70e-12; 5.03e-01], 1:0.861s, tol=7.50e-02 83.20% -> 4.69s 29 : 21.600 μm, [1.48e-12; 7.38e-01], 1:0.846s, tol=7.24e-02 86.40% -> 3.73s 30 : 21.800 μm, [1.01e-12; 9.34e-01], 1:0.863s, tol=2.12e-02 87.20% -> 3.67s 31 : 22.000 μm, [1.11e-12; 9.98e-01], 1:0.840s, tol=3.65e-02 88.00% -> 3.43s 32 : 22.200 μm, [1.22e-12; 1.01e+00], 1:0.857s, tol=6.20e-02 88.80% -> 3.35s 33 : 22.400 μm, [1.36e-12; 1.00e+00], 1:0.840s, tol=8.53e-02 89.60% -> 3.12s 34 : 22.600 μm, [1.35e-12; 1.01e+00], 1:0.858s, tol=9.21e-02 90.40% -> 3.01s 35 : 22.800 μm, [1.44e-12; 1.00e+00], 1:0.842s, tol=1.06e-01 91.20% -> 2.76s 36 : 23.000 μm, [1.49e-12; 1.00e+00], 1:0.858s, tol=1.02e-01 92.00% -> 2.61s 37 : 23.200 μm, [1.63e-12; 1.00e+00], 1:0.843s, tol=9.02e-02 92.80% -> 2.35s 38 : 23.400 μm, [1.67e-12; 1.00e+00], 1:0.866s, tol=9.63e-02 93.60% -> 2.19s 39 : 23.600 μm, [1.75e-12; 1.00e+00], 1:0.839s, tol=1.02e-01 94.40% -> 1.89s 40 : 23.800 μm, [1.79e-12; 1.00e+00], 1:0.861s, tol=1.03e-01 95.20% -> 1.69s 41 : 24.000 μm, [1.68e-12; 1.00e+00], 1:0.841s, tol=1.01e-01 96.00% -> 1.40s 42 : 24.200 μm, [1.43e-12; 1.00e+00], 1:0.862s, tol=9.89e-02 96.80% -> 1.17s 43 : 24.400 μm, [1.26e-12; 1.00e+00], 1:0.839s, tol=1.01e-01 97.60% -> 866.16ms 44 : 24.600 μm, [1.10e-12; 1.00e+00], 1:0.862s, tol=9.77e-02 98.40% -> 603.04ms 45 : 24.800 μm, [1.02e-12; 1.00e+00], 1:0.844s, tol=9.67e-02 99.20% -> 299.57ms 46 : 25.000 μm, [9.35e-13; 1.00e+00], 1:0.867s, tol=9.35e-02 100.00% -> -0.00µs Saved simulation and summary in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.0/examples/PhaseField/results/PlateWithHole/Test/Miehe_AT1_optimMesh Generate movie 01/23 (4.35 %) 2.54 s Generate movie 02/23 (8.70 %) 2.21 s Generate movie 03/23 (13.04 %) 2.15 s Generate movie 04/23 (17.39 %) 2.05 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.70 s Generate movie 08/23 (34.78 %) 1.58 s Generate movie 09/23 (39.13 %) 1.45 s Generate movie 10/23 (43.48 %) 1.35 s Generate movie 11/23 (47.83 %) 1.28 s Generate movie 12/23 (52.17 %) 1.16 s Generate movie 13/23 (56.52 %) 1.06 s Generate movie 14/23 (60.87 %) 933.13 ms Generate movie 15/23 (65.22 %) 867.55 ms Generate movie 16/23 (69.57 %) 744.38 ms Generate movie 17/23 (73.91 %) 627.48 ms Generate movie 18/23 (78.26 %) 536.81 ms Generate movie 19/23 (82.61 %) 426.74 ms Generate movie 20/23 (86.96 %) 318.88 ms Generate movie 21/23 (91.30 %) 212.87 ms Generate movie 22/23 (95.65 %) 107.35 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, isHollow=True) 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""" E = {E:.2e} Pa; v = {v}; Gc = {Gc} J/m2; l0 = {l0:.2e} m while ud <= u_max = {u_max:.2e}: ud += uinc0 if simu.damage.max() < {threshold} else uinc1 uinc0 = {uinc0:.1e}; uinc1 = {uinc1:.1e} 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, verbosity=False) 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 44.053 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 `_