.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/Contact/Contact1.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_Contact1.py: Contact1 ======== A rigid indenter is pressed into an elastic half-space (2D, **plane strain**) and the finite-element contact pressure is compared to the closed-form solutions, with ``E* = E/(1-ν²)``: - **parabola** ``z = x²/2R`` (Hertz line contact): ``p(x) = p₀√(1-x²/a²)``, ``p₀ = a·E*/(2R)`` — a semi-ellipse. - **wedge** ``z = |x|·tanθ``: ``p(x) = (E*·tanθ/π)·arccosh(a/|x|)`` — a log peak. The contact half-width ``a`` is read from the FE solution; the analytical curve is plotted with that same ``a`` (so the comparison is shape + magnitude). The contact pressure is obtained directly from the penalty contact as ``p = εₙ⟨-g⟩`` on the surface. .. GENERATED FROM PYTHON SOURCE LINES 20-241 .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/Contact/images/sphx_glr_Contact1_001.png :alt: Contact1 :srcset: /examples/Contact/images/sphx_glr_Contact1_001.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_Contact1_001.vtksz .. image-sg:: /examples/Contact/images/sphx_glr_Contact1_002.gif :alt: Contact1 :srcset: /examples/Contact/images/sphx_glr_Contact1_002.gif :class: sphx-glr-single-img .. image-sg:: /examples/Contact/images/sphx_glr_Contact1_003.gif :alt: Contact1 :srcset: /examples/Contact/images/sphx_glr_Contact1_003.gif :class: sphx-glr-single-img .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/Contact/images/sphx_glr_Contact1_004.png :alt: Contact1 :srcset: /examples/Contact/images/sphx_glr_Contact1_004.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_Contact1_004.vtksz .. image-sg:: /examples/Contact/images/sphx_glr_Contact1_005.png :alt: parabola (a = 0.419), wedge (a = 0.473) :srcset: /examples/Contact/images/sphx_glr_Contact1_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [parabola] pressing the rigid indenter (Newton per step): ======= elastic problem at iteration 0 ======= At Newton iteration 1 norm is 5.536805892087e+03 At Newton iteration 2 norm is 1.336298576388e+02 At Newton iteration 3 norm is 4.858189561363e+01 At Newton iteration 4 norm is 1.965806277962e+01 At Newton iteration 5 norm is 6.146824670310e+00 At Newton iteration 6 norm is 3.769446750457e-12 ======= elastic problem at iteration 1 ======= At Newton iteration 1 norm is 6.924212846549e+03 At Newton iteration 2 norm is 1.507335478405e+02 At Newton iteration 3 norm is 7.419047099757e+01 At Newton iteration 4 norm is 3.478426175903e+01 At Newton iteration 5 norm is 1.885129836358e+01 At Newton iteration 6 norm is 2.824259659696e+00 At Newton iteration 7 norm is 8.257509454179e-12 ======= elastic problem at iteration 2 ======= At Newton iteration 1 norm is 7.711403046518e+03 At Newton iteration 2 norm is 1.774220663138e+02 At Newton iteration 3 norm is 6.682079946032e+01 At Newton iteration 4 norm is 2.735547015955e+01 At Newton iteration 5 norm is 9.302345940108e+00 At Newton iteration 6 norm is 1.340417369513e-11 ======= elastic problem at iteration 3 ======= At Newton iteration 1 norm is 8.313675151848e+03 At Newton iteration 2 norm is 1.803032438628e+02 At Newton iteration 3 norm is 6.533342138905e+01 At Newton iteration 4 norm is 2.405611452238e+01 At Newton iteration 5 norm is 3.948698407378e+00 At Newton iteration 6 norm is 1.821760791086e-11 ======= elastic problem at iteration 4 ======= At Newton iteration 1 norm is 8.818907274369e+03 At Newton iteration 2 norm is 1.594557717343e+02 At Newton iteration 3 norm is 1.004618205320e+02 At Newton iteration 4 norm is 4.640474059448e+01 At Newton iteration 5 norm is 2.459138270659e+01 At Newton iteration 6 norm is 3.643228916354e+00 At Newton iteration 7 norm is 2.353048230135e-11 ======= elastic problem at iteration 5 ======= At Newton iteration 1 norm is 9.254979631209e+03 At Newton iteration 2 norm is 1.693663162754e+02 At Newton iteration 3 norm is 1.093546985568e+02 At Newton iteration 4 norm is 5.030820094599e+01 At Newton iteration 5 norm is 5.338180308909e+00 At Newton iteration 6 norm is 2.800552478699e-11 ======= elastic problem at iteration 6 ======= At Newton iteration 1 norm is 9.644403746802e+03 At Newton iteration 2 norm is 1.837055602393e+02 At Newton iteration 3 norm is 1.199733597807e+02 At Newton iteration 4 norm is 5.977487725063e+01 At Newton iteration 5 norm is 9.775216534900e+00 At Newton iteration 6 norm is 3.320535179354e-11 ======= elastic problem at iteration 7 ======= At Newton iteration 1 norm is 9.994142186081e+03 At Newton iteration 2 norm is 1.643649443689e+02 At Newton iteration 3 norm is 7.100589654474e+01 At Newton iteration 4 norm is 1.885584229551e+01 At Newton iteration 5 norm is 4.040741247376e-11 ======= elastic problem at iteration 8 ======= At Newton iteration 1 norm is 1.031556506142e+04 At Newton iteration 2 norm is 1.846666676678e+02 At Newton iteration 3 norm is 8.634216166540e+01 At Newton iteration 4 norm is 2.843868240077e+01 At Newton iteration 5 norm is 4.901650786398e+00 At Newton iteration 6 norm is 4.610704407389e-11 ======= elastic problem at iteration 9 ======= At Newton iteration 1 norm is 1.061520518188e+04 At Newton iteration 2 norm is 1.735699098822e+02 At Newton iteration 3 norm is 1.051266955796e+02 At Newton iteration 4 norm is 4.206082216874e+01 At Newton iteration 5 norm is 1.602554912949e+01 At Newton iteration 6 norm is 5.186825359889e-11 ======= elastic problem at iteration 10 ======= At Newton iteration 1 norm is 1.089468953721e+04 At Newton iteration 2 norm is 1.967415724883e+02 At Newton iteration 3 norm is 1.258063637510e+02 At Newton iteration 4 norm is 5.776602967803e+01 At Newton iteration 5 norm is 3.032635689800e+01 At Newton iteration 6 norm is 3.437445123131e+00 At Newton iteration 7 norm is 5.970618529298e-11 ======= elastic problem at iteration 11 ======= At Newton iteration 1 norm is 1.116302764789e+04 At Newton iteration 2 norm is 1.817153901348e+02 At Newton iteration 3 norm is 7.847938864843e+01 At Newton iteration 4 norm is 1.951932688858e+01 At Newton iteration 5 norm is 6.314203785870e-11 [parabola] a=0.419 error vs analytical over 0.1a-0.85a: mean 0% max 1% Generate movie 01/12 (8.33 %) 2.66 s Generate movie 02/12 (16.67 %) 3.29 s Generate movie 03/12 (25.00 %) 2.78 s Generate movie 04/12 (33.33 %) 2.50 s Generate movie 05/12 (41.67 %) 2.18 s Generate movie 06/12 (50.00 %) 1.89 s Generate movie 07/12 (58.33 %) 1.56 s Generate movie 08/12 (66.67 %) 1.28 s Generate movie 09/12 (75.00 %) 965.52 ms Generate movie 10/12 (83.33 %) 628.89 ms Generate movie 11/12 (91.67 %) 313.71 ms Generate movie 12/12 (100.00 %) 0.00 µs [wedge] pressing the rigid indenter (Newton per step): ======= elastic problem at iteration 0 ======= At Newton iteration 1 norm is 6.736266687598e+03 At Newton iteration 2 norm is 7.238066934825e+02 At Newton iteration 3 norm is 1.455169662229e+02 At Newton iteration 4 norm is 8.527654385270e+01 At Newton iteration 5 norm is 6.179022523118e-12 ======= elastic problem at iteration 1 ======= At Newton iteration 1 norm is 1.044830094297e+04 At Newton iteration 2 norm is 7.218267191013e+02 At Newton iteration 3 norm is 1.666526524070e+02 At Newton iteration 4 norm is 5.905000527297e+01 At Newton iteration 5 norm is 1.391936726019e-11 ======= elastic problem at iteration 2 ======= At Newton iteration 1 norm is 1.310754115797e+04 At Newton iteration 2 norm is 7.024703210645e+02 At Newton iteration 3 norm is 1.293800126937e+02 At Newton iteration 4 norm is 8.889876310227e+01 At Newton iteration 5 norm is 4.066926771745e+01 At Newton iteration 6 norm is 2.356942823635e-11 ======= elastic problem at iteration 3 ======= At Newton iteration 1 norm is 1.533010007974e+04 At Newton iteration 2 norm is 6.850037859907e+02 At Newton iteration 3 norm is 1.535887444124e+02 At Newton iteration 4 norm is 6.162716758653e+01 At Newton iteration 5 norm is 1.626376027190e+01 At Newton iteration 6 norm is 3.372581795418e-11 ======= elastic problem at iteration 4 ======= At Newton iteration 1 norm is 1.735857879453e+04 At Newton iteration 2 norm is 6.684958294746e+02 At Newton iteration 3 norm is 1.491381961254e+02 At Newton iteration 4 norm is 6.807775395264e+01 At Newton iteration 5 norm is 3.113727608740e+01 At Newton iteration 6 norm is 4.471076356498e-11 ======= elastic problem at iteration 5 ======= At Newton iteration 1 norm is 1.926247104997e+04 At Newton iteration 2 norm is 6.706878251850e+02 At Newton iteration 3 norm is 1.570830389009e+02 At Newton iteration 4 norm is 7.411624608189e+01 At Newton iteration 5 norm is 3.763981102320e+01 At Newton iteration 6 norm is 1.547822566224e+00 At Newton iteration 7 norm is 5.602731477211e-11 ======= elastic problem at iteration 6 ======= At Newton iteration 1 norm is 2.106920123204e+04 At Newton iteration 2 norm is 6.495395617115e+02 At Newton iteration 3 norm is 1.428840484679e+02 At Newton iteration 4 norm is 7.049824480479e+01 At Newton iteration 5 norm is 4.051662604050e+01 At Newton iteration 6 norm is 6.864918385469e+00 At Newton iteration 7 norm is 6.818971394499e-11 ======= elastic problem at iteration 7 ======= At Newton iteration 1 norm is 2.276869018696e+04 At Newton iteration 2 norm is 6.489914597175e+02 At Newton iteration 3 norm is 1.388657797508e+02 At Newton iteration 4 norm is 6.516506777748e+01 At Newton iteration 5 norm is 5.443723725811e+00 At Newton iteration 6 norm is 7.965796497607e-11 ======= elastic problem at iteration 8 ======= At Newton iteration 1 norm is 2.444117929985e+04 At Newton iteration 2 norm is 6.277568072516e+02 At Newton iteration 3 norm is 1.165282274664e+02 At Newton iteration 4 norm is 5.552777913147e+01 At Newton iteration 5 norm is 2.950391665648e+01 At Newton iteration 6 norm is 9.563669248353e-11 ======= elastic problem at iteration 9 ======= At Newton iteration 1 norm is 2.604666450243e+04 At Newton iteration 2 norm is 6.234931195458e+02 At Newton iteration 3 norm is 1.335466139163e+02 At Newton iteration 4 norm is 4.520648797953e+01 At Newton iteration 5 norm is 1.786100859015e+01 At Newton iteration 6 norm is 1.096941931236e-10 ======= elastic problem at iteration 10 ======= At Newton iteration 1 norm is 2.761454055922e+04 At Newton iteration 2 norm is 6.169702866905e+02 At Newton iteration 3 norm is 1.441537783234e+02 At Newton iteration 4 norm is 5.703103385829e+01 At Newton iteration 5 norm is 6.072912967074e+00 At Newton iteration 6 norm is 1.222048763634e-10 ======= elastic problem at iteration 11 ======= At Newton iteration 1 norm is 2.915203671188e+04 At Newton iteration 2 norm is 6.096796538288e+02 At Newton iteration 3 norm is 1.495307684225e+02 At Newton iteration 4 norm is 6.623428440921e+01 At Newton iteration 5 norm is 1.522177670189e+01 At Newton iteration 6 norm is 1.395931823722e-10 [wedge] a=0.473 error vs analytical over 0.1a-0.85a: mean 0% max 0% Generate movie 01/12 (8.33 %) 1.87 s Generate movie 02/12 (16.67 %) 3.18 s Generate movie 03/12 (25.00 %) 2.83 s Generate movie 04/12 (33.33 %) 2.56 s Generate movie 05/12 (41.67 %) 2.20 s Generate movie 06/12 (50.00 %) 1.91 s Generate movie 07/12 (58.33 %) 1.60 s Generate movie 08/12 (66.67 %) 1.28 s Generate movie 09/12 (75.00 %) 945.10 ms Generate movie 10/12 (83.33 %) 638.28 ms Generate movie 11/12 (91.67 %) 316.75 ms Generate movie 12/12 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 22 import numpy as np import matplotlib.pyplot as plt from EasyFEA import Terminal, Folder, Models, ElemType, Mesh, PyVista from EasyFEA.Geoms import Domain, Points, Circle from EasyFEA.FEM import MatrixType from _utils import RigidContact # ---------------------------------------------- # Rigid indenter profiles z(x), spanning [0, W] and rising away from x=0 # ---------------------------------------------- def parabola_indenter() -> Mesh: # cluster samples near x=0: a curved indenter is gauged by nearest-sample # projection, so it must be resolved finely over the small contact (a << W), # otherwise that sampling — not the body mesh — dominates the pressure error xs = W * np.linspace(0, 1, 200) ** 2 geom = Points( [ *[(x, x**2 / (2 * Rc)) for x in xs], # lower bound *[(W, W**2 / (2 * Rc) + 2), (0, 2)], ] ) mesh = geom.Mesh_2D([], ElemType.TRI3) nodes = mesh.Nodes_Conditions(lambda x, y, z: y <= x**2 / (2 * Rc) + 1e-9) mesh.Set_Tag(nodes, "contact") return mesh def wedge_indenter() -> Mesh: geom = Points( [ (0, 0), (W, W * np.tan(theta)), (W, W * np.tan(theta) + 2), (0, 2), ] ) mesh = geom.Mesh_2D([], ElemType.TRI3) nodes = mesh.Nodes_Conditions(lambda x, y, z: y <= x * np.tan(theta) + 1e-9) mesh.Set_Tag(nodes, "contact") return mesh def parabola_analytic(x, a): return a * Es / (2 * Rc) * np.sqrt(np.clip(1 - (x / a) ** 2, 0, None)) def wedge_analytic(x, a): return Es * np.tan(theta) / np.pi * np.arccosh(np.clip(a / np.abs(x), 1.0, None)) indenter_cases = { "parabola": (parabola_indenter, parabola_analytic), "wedge": (wedge_indenter, wedge_analytic), } def contact_pressure(simu: RigidContact, indenter: Mesh): """Penalty contact pressure ``εₙ⟨-g⟩`` and x-coordinate on the surface Gauss points.""" matrixType = MatrixType.mass list_p = [] list_x = [] for group1d in simu.mesh.Get_list_groupElem(1): N_pg = group1d.Get_N_pg(matrixType)[:, 0, :] # x = X + u elements = group1d.Get_Elements_Tag("top") X_e_pg = group1d.Get_GaussCoordinates_e_pg(matrixType)[elements] u_e_pg = simu.displacement.reshape(simu.mesh.Nn, 2) x_e_pg = X_e_pg.copy() x_e_pg[..., :2] += np.einsum( "pn,enc->epc", N_pg, u_e_pg[group1d.connect[elements]] ) list_x.extend(x_e_pg[..., 0].ravel()) # get pressure for groupIndent in indenter.Get_list_groupElem(1): elements = ( groupIndent.Get_Elements_Tag("contact") if "contact" in groupIndent.elementTags else None ) gap_e_pg, _ = groupIndent._Get_gap_and_normal( x_e_pg=x_e_pg, elements=elements, coord=indenter.center, matrixType=matrixType, ) p_e_pg = simu.penalty * np.maximum(-gap_e_pg, 0.0) list_p.extend(p_e_pg.ravel()) return np.asarray(list_p), np.asarray(list_x) if __name__ == "__main__": Terminal.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- folder = Folder.Results_Dir() result = "Svm" E, v = 210000.0, 0.3 Es = E / (1 - v**2) W, D = 6.0, 6.0 # elastic half-space (symmetry about x=0): wide & deep vs contact meshSize = W / 10 N = 12 # indentation steps penalty = 100 * Es Rc = 6.0 # parabola radius of curvature theta = np.deg2rad(6) # wedge half-angle indeter_delta = { "parabola": 0.05, "wedge": 0.12, } # indentation depth per case # ---------------------------------------------- # Elastic half-space (surface at y=0) # ---------------------------------------------- body = Domain((0, -D), (W, 0), meshSize) refineGeoms = [ Circle((0, 0), W * coef * 2, W / N) for coef, N in zip( [0.2, 0.4], [200, 100], ) ] mesh = body.Mesh_2D([], ElemType.TRI3, refineGeoms=refineGeoms) nodes_bottom = mesh.Nodes_Conditions(lambda x, y, z: y == -D) nodes_sym = mesh.Nodes_Conditions(lambda x, y, z: x == 0) nodes_top = mesh.Nodes_Conditions(lambda x, y, z: y == 0) mesh.Set_Tag(nodes_top, "top") PyVista.Plot_Mesh(mesh).show() material = Models.Elastic.Isotropic(2, E=E, v=v, planeStress=False) # ---------------------------------------------- # Solve each case and compare FE pressure to theory # ---------------------------------------------- fig, axes = plt.subplots(1, len(indenter_cases), figsize=(11, 4)) for ax, (name, (build, analytic)) in zip(axes, indenter_cases.items()): simu = RigidContact(mesh, material, penalty) indenter = build() list_indeter = [indenter] delta = indeter_delta[name] print(f"\n[{name}] pressing the rigid indenter (Newton per step):") for i in range(N): # update indeter indenter = list_indeter[0] indenter = indenter.copy() indenter.Translate(dy=-(i + 1) / N * delta) list_indeter.append(indenter) simu._contactMesh = indenter # solve contact simu.Bc_Init() simu.add_dirichlet(nodes_bottom, [0, 0], simu.Get_unknowns()) simu.add_dirichlet(nodes_sym, [0], ["x"]) simu.Solve() simu.Save_Iter() pg, xg = contact_pressure(simu, indenter) order = np.argsort(xg) xg, pg = xg[order], pg[order] a = xg[pg > 1e-6 * pg.max()].max() # FE contact half-width # relative error over the inner patch (skip the singular center & edge), # normalised by the FE peak pressure band = (xg > 0.1 * a) & (xg < 0.85 * a) rel = np.abs(pg[band] - analytic(xg[band], a)) / pg.max() print( f"\n[{name}] a={a:.3f} error vs analytical over 0.1a-0.85a: " f"mean {rel.mean() * 100:.0f}% max {rel.max() * 100:.0f}%" ) # ---------------------------------------------- # Results # ---------------------------------------------- ax.plot(xg, pg, "o", ms=3, label="FE εₙ⟨-g⟩") xa = np.linspace(1e-3, a, 200) ax.plot(xa, analytic(xa, a), "k-", label="analytical") ax.set_title(f"{name} (a = {a:.3f})") ax.set_xlabel("x") ax.set_ylabel("contact pressure") ax.set_xlim(0, 1.6 * a) ax.legend() ax.grid(True) def Plot_Iter(plotter: PyVista.pv.Plotter, n): simu.Set_Iter(n) PyVista.Plot( simu, result, 1, color="k", nColors=10, show_grid=True, plotter=plotter ) PyVista.Plot(list_indeter[n], color="gray", plotter=plotter) PyVista.Plot_Elements( list_indeter[n], dimElem=1, color="k", linewidth=1.5, plotter=plotter ) PyVista._setCameraPosition(plotter, 2) PyVista.Movie_func(Plot_Iter, N, folder=folder, filename=f"{name}.gif") fig.tight_layout() plotter = PyVista._Plotter() Plot_Iter(plotter, -1) plotter.show() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 16.557 seconds) .. _sphx_glr_download_examples_Contact_Contact1.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Contact1.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Contact1.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Contact1.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_