.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/Contact/Contact2.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_Contact_Contact2.py: Contact2 ======== Frictionless contact between a deformable block and a rigid indenter, solved as a genuine non-linear problem with Newton-Raphson and a **penalty** regularisation of the unilateral (Signorini) condition. The rigid indenter is treated as an obstacle: at every contact-surface Gauss point the normal gap ``gₙ`` to the obstacle surface is measured and, where it is negative (penetration), a penalty traction ``εₙ⟨-gₙ⟩ n`` resists it.These contributions are added to the elastic residual/tangent through ``Operators.NonLinear.PenaltyContact`` and the non-linear system ``A(u) Δu = -R(u)`` is solved with Newton at each load step. Runs in 2D (block + half-disc punch) and 3D (box + half-cylinder punch); set ``dim`` below. The non-linear simulation (``ElasticContact``) lives in ``_utils.py``. .. GENERATED FROM PYTHON SOURCE LINES 16-132 .. image-sg:: /examples/Contact/images/sphx_glr_Contact2_001.gif :alt: Contact2 :srcset: /examples/Contact/images/sphx_glr_Contact2_001.gif :class: sphx-glr-single-img .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/Contact/images/sphx_glr_Contact2_002.png :alt: Contact2 :srcset: /examples/Contact/images/sphx_glr_Contact2_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v3.1.0/docs/examples/Contact/images/sphx_glr_Contact2_002.vtksz .. rst-class:: sphx-glr-script-out .. code-block:: none Penalty contact solve in 3D (Newton per step): ======= elastic problem at iteration 0 ======= At Newton iteration 1 norm is 1.185313683882e+05 At Newton iteration 2 norm is 5.812774982115e+02 At Newton iteration 3 norm is 9.901128618810e-09 ======= elastic problem at iteration 1 ======= At Newton iteration 1 norm is 1.166673631212e+05 At Newton iteration 2 norm is 1.162554996434e+03 At Newton iteration 3 norm is 7.749237811253e-09 ======= elastic problem at iteration 2 ======= At Newton iteration 1 norm is 1.148067412636e+05 At Newton iteration 2 norm is 1.743832494649e+03 At Newton iteration 3 norm is 5.100563749518e-09 ======= elastic problem at iteration 3 ======= At Newton iteration 1 norm is 1.129496700202e+05 At Newton iteration 2 norm is 2.325109992869e+03 At Newton iteration 3 norm is 1.008345954691e-08 ======= elastic problem at iteration 4 ======= At Newton iteration 1 norm is 1.110963274459e+05 At Newton iteration 2 norm is 2.906387491082e+03 At Newton iteration 3 norm is 1.210719790661e-08 ======= elastic problem at iteration 5 ======= At Newton iteration 1 norm is 1.092469033083e+05 At Newton iteration 2 norm is 3.487664989292e+03 At Newton iteration 3 norm is 9.346021680056e-09 ======= elastic problem at iteration 6 ======= At Newton iteration 1 norm is 1.074016000308e+05 At Newton iteration 2 norm is 4.068942487516e+03 At Newton iteration 3 norm is 8.768863505757e-09 ======= elastic problem at iteration 7 ======= At Newton iteration 1 norm is 1.055606337234e+05 At Newton iteration 2 norm is 4.650219985722e+03 At Newton iteration 3 norm is 7.563893860032e-09 ======= elastic problem at iteration 8 ======= At Newton iteration 1 norm is 1.037242353127e+05 At Newton iteration 2 norm is 5.231497483937e+03 At Newton iteration 3 norm is 8.001390881827e-09 ======= elastic problem at iteration 9 ======= At Newton iteration 1 norm is 1.018926517788e+05 At Newton iteration 2 norm is 5.812774982156e+03 At Newton iteration 3 norm is 1.478495779020e-08 ======= elastic problem at iteration 10 ======= At Newton iteration 1 norm is 1.000661475126e+05 At Newton iteration 2 norm is 6.394052480372e+03 At Newton iteration 3 norm is 1.041112153236e-08 ======= elastic problem at iteration 11 ======= At Newton iteration 1 norm is 9.824500580546e+04 At Newton iteration 2 norm is 6.975329978588e+03 At Newton iteration 3 norm is 8.464284619152e-09 ======= elastic problem at iteration 12 ======= At Newton iteration 1 norm is 9.642953048542e+04 At Newton iteration 2 norm is 7.556607476799e+03 At Newton iteration 3 norm is 8.758018624992e-09 ======= elastic problem at iteration 13 ======= At Newton iteration 1 norm is 9.462004771586e+04 At Newton iteration 2 norm is 8.137884975020e+03 At Newton iteration 3 norm is 1.209765291054e-08 ======= elastic problem at iteration 14 ======= At Newton iteration 1 norm is 9.281690797513e+04 At Newton iteration 2 norm is 8.719162473235e+03 At Newton iteration 3 norm is 9.616481315885e-09 ======= elastic problem at iteration 15 ======= At Newton iteration 1 norm is 9.102048823516e+04 At Newton iteration 2 norm is 9.300439971447e+03 At Newton iteration 3 norm is 9.611662141835e-09 ======= elastic problem at iteration 16 ======= At Newton iteration 1 norm is 8.923119436215e+04 At Newton iteration 2 norm is 9.881717469665e+03 At Newton iteration 3 norm is 7.262590546939e-09 ======= elastic problem at iteration 17 ======= At Newton iteration 1 norm is 8.744946376172e+04 At Newton iteration 2 norm is 1.046299496788e+04 At Newton iteration 3 norm is 1.431849773446e-08 ======= elastic problem at iteration 18 ======= At Newton iteration 1 norm is 8.567576829723e+04 At Newton iteration 2 norm is 1.104427246609e+04 At Newton iteration 3 norm is 9.918191748510e-09 ======= elastic problem at iteration 19 ======= At Newton iteration 1 norm is 8.391061750845e+04 At Newton iteration 2 norm is 1.162554996431e+04 At Newton iteration 3 norm is 9.706379453461e-09 ======= elastic problem at iteration 20 ======= At Newton iteration 1 norm is 8.215456216286e+04 At Newton iteration 2 norm is 1.220682746252e+04 At Newton iteration 3 norm is 4.965093719761e-09 ======= elastic problem at iteration 21 ======= At Newton iteration 1 norm is 8.040819817587e+04 At Newton iteration 2 norm is 1.278810496073e+04 At Newton iteration 3 norm is 1.387281033992e-08 ======= elastic problem at iteration 22 ======= At Newton iteration 1 norm is 7.871095643409e+04 At Newton iteration 2 norm is 1.180866389997e+04 At Newton iteration 3 norm is 1.076353168632e-08 ======= elastic problem at iteration 23 ======= At Newton iteration 1 norm is 7.830569184190e+04 At Newton iteration 2 norm is 6.483218661311e-09 ======= elastic problem at iteration 24 ======= At Newton iteration 1 norm is 7.830569184189e+04 At Newton iteration 2 norm is 8.257101205397e-09 ======= elastic problem at iteration 25 ======= At Newton iteration 1 norm is 7.830569184183e+04 At Newton iteration 2 norm is 1.016022596422e-08 ======= elastic problem at iteration 26 ======= At Newton iteration 1 norm is 7.830569184190e+04 At Newton iteration 2 norm is 7.109110221733e-09 ======= elastic problem at iteration 27 ======= At Newton iteration 1 norm is 7.830569184190e+04 At Newton iteration 2 norm is 8.274836961939e-09 ======= elastic problem at iteration 28 ======= At Newton iteration 1 norm is 7.830569184190e+04 At Newton iteration 2 norm is 7.358984570085e-09 ======= elastic problem at iteration 29 ======= At Newton iteration 1 norm is 7.830569184189e+04 At Newton iteration 2 norm is 4.735839961635e-09 ==================== Mesh ==================== elemType: HEXA8 Ne: 192 Nn: 324 ==================== Model ==================== Isotropic: E = 2.10e+05, v = 0.3 solver:scipy ============= Boundary Conditions ============= Unspecified. =================== Results =================== W def = 2618.76 Svm max = 2910.49 Evm max = 1.80 % Ux max = 1.62e-02 Ux min = -1.62e-02 Uy max = 0.00e+00 Uy min = -1.06e-01 Uz max = 1.26e-02 Uz min = -1.26e-02 =================== TicTac =================== Mesh: 56.295 ms Boundary Conditions: 208.735 µs Matrix: 651.949 ms Solver: 623.681 ms Resolution elastic: 1.380 s PostProcessing: 433.372 µs Generate movie 01/30 (3.33 %) 6.05 s Generate movie 02/30 (6.67 %) 10.47 s Generate movie 03/30 (10.00 %) 9.94 s Generate movie 04/30 (13.33 %) 9.59 s Generate movie 05/30 (16.67 %) 9.50 s Generate movie 06/30 (20.00 %) 8.89 s Generate movie 07/30 (23.33 %) 8.64 s Generate movie 08/30 (26.67 %) 8.30 s Generate movie 09/30 (30.00 %) 7.81 s Generate movie 10/30 (33.33 %) 7.41 s Generate movie 11/30 (36.67 %) 7.14 s Generate movie 12/30 (40.00 %) 6.76 s Generate movie 13/30 (43.33 %) 6.38 s Generate movie 14/30 (46.67 %) 6.05 s Generate movie 15/30 (50.00 %) 5.66 s Generate movie 16/30 (53.33 %) 5.31 s Generate movie 17/30 (56.67 %) 4.89 s Generate movie 18/30 (60.00 %) 4.60 s Generate movie 19/30 (63.33 %) 4.21 s Generate movie 20/30 (66.67 %) 3.80 s Generate movie 21/30 (70.00 %) 3.41 s Generate movie 22/30 (73.33 %) 2.99 s Generate movie 23/30 (76.67 %) 2.70 s Generate movie 24/30 (80.00 %) 2.30 s Generate movie 25/30 (83.33 %) 1.90 s Generate movie 26/30 (86.67 %) 1.52 s Generate movie 27/30 (90.00 %) 1.14 s Generate movie 28/30 (93.33 %) 765.82 ms Generate movie 29/30 (96.67 %) 381.38 ms Generate movie 30/30 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 17 from EasyFEA import Terminal, Folder, Models, ElemType, PyVista from EasyFEA.Geoms import Point, Domain, Points from _utils import RigidContact if __name__ == "__main__": Terminal.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- dim = 3 # 2 or 3 result = "Svm" R = 10 # block size thickness = R / 3 # out-of-plane extent (3D) meshSize = R / 20 if dim == 2 else R / 8 N = 30 # load steps delta = 0.1 penalty = 1e7 # contact stiffness (larger -> less penetration) folder = Folder.Results_Dir() # ---------------------------------------------- # Body # ---------------------------------------------- body = Domain( (-R / 2, 0), (R / 2, R), meshSize, ) if dim == 2: mesh = body.Mesh_2D( [], ElemType.QUAD4, isOrganised=True, ) else: nz = max(1, round(thickness / meshSize)) mesh = body.Mesh_Extrude( [], [0, 0, thickness], [nz], ElemType.HEXA8, isOrganised=True, ) nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0) # ---------------------------------------------- # Indenter # ---------------------------------------------- r = R / 3 contour = Points( [ Point(-R / 2, R, r=r), Point(R / 2, R, r=r), Point(R / 2, 2 * R), Point(-R / 2, 2 * R), ] ) if dim == 2: indenter = contour.Mesh_2D([], ElemType.TRI3) else: indenter = contour.Mesh_Extrude([], [0, 0, thickness], [nz], ElemType.TETRA4) # lower (contact) surface of the punch: below the end of the fillets nodes_contact = indenter.Nodes_Conditions(lambda x, y, z: y <= R + r + 1e-6) indenter.Set_Tag(nodes_contact, "contact") # ---------------------------------------------- # Simulation # ---------------------------------------------- material = Models.Elastic.Isotropic( dim, E=210000, v=0.3, planeStress=True, thickness=thickness ) simu = RigidContact(mesh, material, penalty) list_indenter = [indenter] print(f"Penalty contact solve in {dim}D (Newton per step):") for i in range(N): # update indeter indenter = list_indenter[0].copy() indenter.Translate(dy=-(i + 1) / N * delta) # lower the rigid indenter list_indenter.append(indenter) simu._contactMesh = indenter # solve contact simu.Bc_Init() simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns()) simu.Solve() simu.Save_Iter() print(simu) # ---------------------------------------------- # Results # ---------------------------------------------- def Plot_Iter(plotter, n): simu.Set_Iter(n) PyVista.Plot( simu, result, 1, color="k", nColors=21, show_grid=True, plotter=plotter ) PyVista.Plot(list_indenter[n], color="gray", alpha=0.4, plotter=plotter) PyVista.Plot_Elements( list_indenter[n], color="k", dimElem=1, linewidth=2, plotter=plotter ) PyVista._setCameraPosition(plotter, mesh.inDim) PyVista.Movie_func(Plot_Iter, N, folder=folder, filename=f"{result}.gif") plotter = PyVista._Plotter() result = "uy" Plot_Iter(plotter, -1) plotter.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.068 seconds) .. _sphx_glr_download_examples_Contact_Contact2.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Contact2.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Contact2.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Contact2.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_