.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/Contact/Contact3.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_Contact3.py: Contact3 ======== Frictionless contact between a thin elastic arch strip and a rigid block, solved with the penalty method and Newton-Raphson. The rigid block 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 (arch strip) and 3D (extruded arch strip); set ``dim`` below. The non-linear simulation (``ElasticContact``) lives in ``_utils.py``. .. GENERATED FROM PYTHON SOURCE LINES 16-125 .. image-sg:: /examples/Contact/images/sphx_glr_Contact3_001.gif :alt: Contact3 :srcset: /examples/Contact/images/sphx_glr_Contact3_001.gif :class: sphx-glr-single-img .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/Contact/images/sphx_glr_Contact3_002.png :alt: Contact3 :srcset: /examples/Contact/images/sphx_glr_Contact3_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_Contact3_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 3.317591015983e+05 At Newton iteration 2 norm is 8.823302971212e+05 At Newton iteration 3 norm is 4.817575525634e+04 At Newton iteration 4 norm is 5.654339121181e+03 At Newton iteration 5 norm is 1.309684810071e+03 At Newton iteration 6 norm is 1.076949440762e+02 At Newton iteration 7 norm is 9.287557323790e+01 At Newton iteration 8 norm is 8.100208244980e+01 At Newton iteration 9 norm is 8.255348762812e+01 At Newton iteration 10 norm is 9.830625477488e+01 At Newton iteration 11 norm is 1.035701793339e+02 At Newton iteration 12 norm is 4.170797681642e+01 At Newton iteration 13 norm is 3.374784836984e+01 At Newton iteration 14 norm is 3.550637878922e+01 At Newton iteration 15 norm is 3.129502005272e+01 At Newton iteration 16 norm is 1.878599545516e+01 At Newton iteration 17 norm is 1.073427713889e-10 ======= elastic problem at iteration 1 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 5.115550059265e-10 ======= elastic problem at iteration 2 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 5.516151581912e-10 ======= elastic problem at iteration 3 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 6.235578438922e-10 ======= elastic problem at iteration 4 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.499225455071e+03 At Newton iteration 3 norm is 2.880391331552e+02 At Newton iteration 4 norm is 5.116356946139e-10 ======= elastic problem at iteration 5 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 7.986856982269e-10 ======= elastic problem at iteration 6 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 8.355332756958e-10 ======= elastic problem at iteration 7 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 9.131250719347e-10 ======= elastic problem at iteration 8 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.033042535117e-09 ======= elastic problem at iteration 9 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 2.127042476046e+02 At Newton iteration 3 norm is 2.297943455031e+00 At Newton iteration 4 norm is 3.068947162576e-01 At Newton iteration 5 norm is 1.046315291892e-09 ======= elastic problem at iteration 10 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 3.083780453789e+03 At Newton iteration 3 norm is 2.526106928419e+02 At Newton iteration 4 norm is 1.189880401945e+02 At Newton iteration 5 norm is 1.183178684717e-09 ======= elastic problem at iteration 11 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.320396084363e-09 ======= elastic problem at iteration 12 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.441451878807e-09 ======= elastic problem at iteration 13 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.538779091218e-09 ======= elastic problem at iteration 14 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.678832810242e-09 ======= elastic problem at iteration 15 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.833354552873e-09 ======= elastic problem at iteration 16 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.846337922863e-09 ======= elastic problem at iteration 17 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.895384037304e-09 ======= elastic problem at iteration 18 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 1.874177123996e-09 ======= elastic problem at iteration 19 ======= At Newton iteration 1 norm is 3.317591015983e+05 At Newton iteration 2 norm is 2.823474082624e+02 At Newton iteration 3 norm is 2.143841630050e-09 ==================== Mesh ==================== elemType: TETRA4 Ne: 2880 Nn: 960 ==================== Model ==================== Isotropic: E = 2.10e+05, v = 0.3 solver:scipy ============= Boundary Conditions ============= Unspecified. =================== Results =================== W def = 692.92 Svm max = 3977.16 Evm max = 2.46 % Ux max = 1.72e+00 Ux min = -1.70e+00 Uy max = 2.28e+00 Uy min = -2.00e+00 Uz max = 1.31e-01 Uz min = -1.82e-03 =================== TicTac =================== Mesh: 116.684 ms Boundary Conditions: 140.462 µs Matrix: 1.367 s Solver: 1.004 s Resolution elastic: 2.478 s PostProcessing: 688.898 µs Generate movie 01/20 (5.00 %) 4.02 s Generate movie 02/20 (10.00 %) 5.55 s Generate movie 03/20 (15.00 %) 5.30 s Generate movie 04/20 (20.00 %) 5.01 s Generate movie 05/20 (25.00 %) 4.68 s Generate movie 06/20 (30.00 %) 4.31 s Generate movie 07/20 (35.00 %) 4.08 s Generate movie 08/20 (40.00 %) 3.71 s Generate movie 09/20 (45.00 %) 3.40 s Generate movie 10/20 (50.00 %) 3.13 s Generate movie 11/20 (55.00 %) 2.82 s Generate movie 12/20 (60.00 %) 2.49 s Generate movie 13/20 (65.00 %) 2.21 s Generate movie 14/20 (70.00 %) 1.88 s Generate movie 15/20 (75.00 %) 1.57 s Generate movie 16/20 (80.00 %) 1.26 s Generate movie 17/20 (85.00 %) 946.56 ms Generate movie 18/20 (90.00 %) 643.66 ms Generate movie 19/20 (95.00 %) 317.98 ms Generate movie 20/20 (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" e = 10 L = 3 * e t = 1 # strip thickness h = 10 # strip height r = 3 # fillet radius thickness = e / 2 # out-of-plane extent (3D) mS = t / 5 if dim == 2 else t N = 20 # load steps delta = 2 * t penalty = 1e6 # contact stiffness folder = Folder.Results_Dir() # ---------------------------------------------- # Mesh # ---------------------------------------------- # body # ---- p1 = Point(-L / 2 - e) p2 = Point(-L / 2, r=r) p3 = Point(-e / 2, h - t, r=r) p4 = Point(e / 2, h - t, r=r) p5 = Point(L / 2, r=r) p6 = Point(L / 2 + e) lower = Points([p1, p2, p3, p4, p5, p6]) upper = lower.copy() upper.Translate(dy=t) contour = Points(list(lower.points) + list(upper.points[::-1]), mS) if dim == 2: mesh = contour.Mesh_2D([], ElemType.TRI3) else: nz = max(1, round(thickness / mS)) mesh = contour.Mesh_Extrude([], [0, 0, thickness], [nz], ElemType.TETRA4) # block # ----- domain = Domain((-L / 2 - 2 * e, -5 * t), (L / 2 + 2 * e, 0)) if dim == 2: block = domain.Mesh_2D([], ElemType.QUAD4, isOrganised=True) else: block = domain.Mesh_Extrude( [], [0, 0, thickness * 2], [1], ElemType.HEXA8, isOrganised=True ) block.Translate(dz=-block.center[2] / 2) nodes_top = mesh.Nodes_Conditions(lambda x, y, z: y == h) nodes_contact = block.Nodes_Conditions(lambda x, y, z: y == 0) block.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) simu._contactMesh = block print(f"Penalty contact solve in {dim}D (Newton per step):") for i in range(N): # udpate load dep = [0.0] * simu.dim dep[1] = -(i + 1) / N * delta # solve contact simu.Bc_Init() simu.add_dirichlet(nodes_top, dep, 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(block, color="gray", alpha=0.4, plotter=plotter) PyVista.Plot_Elements(block, color="k", dimElem=1, linewidth=2, plotter=plotter) PyVista.Movie_func(Plot_Iter, N, folder=folder, filename="contact.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 10.076 seconds) .. _sphx_glr_download_examples_Contact_Contact3.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Contact3.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Contact3.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Contact3.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_