.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/LinearizedElasticity/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_LinearizedElasticity_Contact1.py: Contact1 ======== Performing a 'Hertz contact problem' with the assumption of frictionless contact The master mesh is considered non-deformable. WARNING: The assumption of small displacements is highly questionable for this simulation. .. GENERATED FROM PYTHON SOURCE LINES 13-196 .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Contact1_001.gif :alt: Contact1 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Contact1_001.gif :class: sphx-glr-single-img .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Contact1_002.png :alt: Contact1 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Contact1_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.1/docs/examples/LinearizedElasticity/images/sphx_glr_Contact1_002.vtksz .. rst-class:: sphx-glr-script-out .. code-block:: none Eps max = 0.38 % Eps max = 0.76 % Eps max = 1.15 % Eps max = 1.76 % Eps max = 2.20 % Eps max = 2.63 % Eps max = 3.58 % Eps max = 3.79 % Eps max = 3.99 % Eps max = 4.19 % Eps max = 4.38 % Eps max = 4.58 % Eps max = 5.03 % Eps max = 5.53 % Eps max = 6.02 % Eps max = 7.10 % Eps max = 7.49 % Eps max = 7.92 % Eps max = 8.36 % Eps max = 8.98 % Eps max = 9.61 % Eps max = 10.24 % Eps max = 10.87 % Eps max = 11.50 % Eps max = 12.12 % Eps max = 12.75 % Eps max = 13.38 % Eps max = 14.01 % Eps max = 14.64 % Eps max = 15.26 % ==================== Mesh ==================== Element type: QUAD4 Ne = 400, Nn = 441 ==================== Model ==================== ElasIsot: E = 2.10e+05, v = 0.3 planeStress = True thickness = 3.33e+00 solver : pypardiso ============= Boundary Conditions ============= Unspecified. =================== Results =================== W def = 202448.77 Svm max = 72054.31 Evm max = 44.60 % Ux max = 1.39e-01 Ux min = -1.56e-01 Uy max = 0.00e+00 Uy min = -9.83e-01 =================== TicTac =================== Mesh : 600.763 ms Boundary Conditions : 1.952 ms Matrix : 4.039 s Solver : 2.329 s Display : 869.805 ms PostProcessing : 365.494 ms Resolution hyperelastic : 6.224 s PyVista_Interface : 8.941 s Generate movie 0/29 Generate movie 1/29 (3.45 %) 16.46 s Generate movie 2/29 (6.90 %) 12.59 s Generate movie 3/29 (10.34 %) 12.07 s Generate movie 4/29 (13.79 %) 11.52 s Generate movie 5/29 (17.24 %) 11.07 s Generate movie 6/29 (20.69 %) 10.81 s Generate movie 7/29 (24.14 %) 10.41 s Generate movie 8/29 (27.59 %) 9.31 s Generate movie 9/29 (31.03 %) 9.52 s Generate movie 10/29 (34.48 %) 8.99 s Generate movie 11/29 (37.93 %) 8.47 s Generate movie 12/29 (41.38 %) 7.81 s Generate movie 13/29 (44.83 %) 7.45 s Generate movie 14/29 (48.28 %) 6.85 s Generate movie 15/29 (51.72 %) 6.51 s Generate movie 16/29 (55.17 %) 6.06 s Generate movie 17/29 (58.62 %) 5.62 s Generate movie 18/29 (62.07 %) 5.19 s Generate movie 19/29 (65.52 %) 4.70 s Generate movie 20/29 (68.97 %) 4.24 s Generate movie 21/29 (72.41 %) 3.71 s Generate movie 22/29 (75.86 %) 3.26 s Generate movie 23/29 (79.31 %) 2.79 s Generate movie 24/29 (82.76 %) 2.33 s Generate movie 25/29 (86.21 %) 1.87 s Generate movie 26/29 (89.66 %) 1.39 s Generate movie 27/29 (93.10 %) 937.71 ms Generate movie 28/29 (96.55 %) 460.00 ms Generate movie 29/29 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 13 # TODO: Compare results with analytical values ? from EasyFEA import Display, Folder, Models, plt, np, ElemType, Simulations, PyVista from EasyFEA.Geoms import Point, Domain, Points if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- dim = 2 # outputs folder = Folder.Join(Folder.RESULTS_DIR, "Contact", "Contact1") pltIter = False makeMovie = True result = "uy" # geom R = 10 height = R thickness = R / 3 # load N = 30 inc = 1e-0 / N cx, cy = 0, -1 # ---------------------------------------------- # Meshes # ---------------------------------------------- meshSize = R / 20 # slave mesh contour_slave = Domain(Point(-R / 2, 0), Point(R / 2, height), meshSize) if dim == 2: mesh_slave = contour_slave.Mesh_2D([], ElemType.QUAD4, isOrganised=True) else: mesh_slave = contour_slave.Mesh_Extrude( [], [0, 0, -thickness], [4], ElemType.HEXA8, isOrganised=True ) # nodes_slave = mesh_slave.Get_list_groupElem(dim-1)[0].nodes nodes_slave = mesh_slave.Nodes_Conditions(lambda x, y, z: y == height) nodes_y0 = mesh_slave.Nodes_Conditions(lambda x, y, z: y == 0) # master mesh r = R / 2 p0 = Point(-R / 2, height, r=r) p1 = Point(R / 2, height, r=r) p2 = Point(R / 2, height + R) p3 = Point(-R / 2, height + R) contour_master = Points([p0, p1, p2, p3]) yMax = height + np.abs(r) if dim == 2: master_mesh = contour_master.Mesh_2D([], ElemType.TRI3) else: master_mesh = contour_master.Mesh_Extrude( [], [0, 0, -thickness - 2], [4], ElemType.TETRA4 ) groupMaster = master_mesh.Get_list_groupElem(dim - 1)[0] if len(master_mesh.Get_list_groupElem(dim - 1)) > 1: Display.MyPrintError( f"The {groupMaster.elemType.name} element group is used. In 3D, TETRA AND HEXA elements are recommended." ) master_mesh.Translate(dz=-(master_mesh.center[2] - mesh_slave.center[2])) # Display.Plot_Tags(mesh_master, alpha=0.1, showId=True) # get master nodes # nodes_master = mesh_master.Get_list_groupElem(dim-1)[0].nodes if dim == 2: nodes_master = master_mesh.Nodes_Tags(["L0", "L1"]) else: nodes_master = master_mesh.Nodes_Tags(["S1", "S2"]) # # plot meshes # ax = Display.Plot_Mesh(master_mesh, alpha=0) # Display.Plot_Mesh(mesh_slave, ax=ax, alpha=0) # # add nodes interface # ax.scatter(*mesh_slave.coord[nodes_slave, :dim].T, label="slave nodes") # ax.scatter(*master_mesh.coord[nodes_master, :dim].T, label="master nodes") # ax.legend() # ax.set_title("Contact nodes") # ---------------------------------------------- # Simulation # ---------------------------------------------- material = Models.ElasIsot( dim, E=210000, v=0.3, planeStress=True, thickness=thickness ) simu = Simulations.ElasticSimu(mesh_slave, material) list_master_mesh = [master_mesh] if pltIter: ax = Display.Plot_Result(simu, result, deformFactor=1) for i in range(N): master_mesh = master_mesh.copy() master_mesh.Translate(cx * inc, cy * inc) list_master_mesh.append(master_mesh) convergence = False coordo_old = simu.Results_displacement_matrix() + simu.mesh.coord while not convergence: # apply new boundary conditions simu.Bc_Init() simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns()) nodes, newU = simu.Get_contact(master_mesh, nodes_slave, nodes_master) if nodes.size > 0: simu.add_dirichlet(nodes, [newU[:, 0], newU[:, 1]], ["x", "y"]) simu.Solve() # check if there is no new nodes in the master mesh oldSize = nodes.size nodes, __ = simu.Get_contact(master_mesh, nodes_slave, nodes_master) convergence = oldSize == nodes.size simu.Save_Iter() print(f"Eps max = {simu.Result('Strain').max() * 100:3.2f} %") if pltIter: Display.Plot_Result(simu, result, plotMesh=True, deformFactor=1, ax=ax) Display.Plot_Mesh(master_mesh, alpha=0, ax=ax) ax.set_title(result) if dim == 3: Display._Axis_equal_3D( ax, np.concatenate((master_mesh.coord, mesh_slave.coord), 0) ) # # Plot arrows # if nodes.size >0: # # get the nodes coordinates on the interface # coordinates = groupMaster.Get_GaussCoordinates_e_p('mass').reshape(-1,3) # ax.scatter(*coordinates[:,:dim].T) # coordo_new = simu.Results_displacement_matrix() + simu.mesh.coordo # ax.scatter(*coordo_old[nodes,:dim].T) # incU = coordo_new - coordo_oldq # [ax.arrow(*coordo_old[node, :dim], *incU[node,:dim],length_includes_head=True) for node in nodes] plt.pause(1e-12) print(simu) # ---------------------------------------------- # Results # ---------------------------------------------- if makeMovie: def DoAnim(plotter, n): simu.Set_Iter(n) PyVista.Plot( simu, "Svm", 1, style="surface", color="k", plotter=plotter, nColors=10, show_grid=True, ) PyVista.Plot(list_master_mesh[n], plotter=plotter, plotMesh=True, alpha=0.2) PyVista.Movie_func(DoAnim, N, folder=folder, filename=f"{result}.gif") if not pltIter: plotter = PyVista.Plot(simu, result, plotMesh=True, deformFactor=1) PyVista.Plot_Mesh(master_mesh, alpha=0.4, plotter=plotter) plotter.show() Display.plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 19.606 seconds) .. _sphx_glr_download_examples_LinearizedElasticity_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 `_