.. 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 16-199 .. 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.7.0/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 ==================== Isotropic: E = 2.10e+05, v = 0.3 planeStress = True thickness = 3.33e+00 solver : scipy ============= Boundary Conditions ============= Unspecified. =================== Results =================== W def = 202448.77 Svm max = 75667.12 Evm max = 46.01 % Ux max = 1.39e-01 Ux min = -1.56e-01 Uy max = 0.00e+00 Uy min = -9.83e-01 =================== TicTac =================== Mesh : 307.982 ms Boundary Conditions : 1.088 ms Matrix : 3.160 s Solver : 2.807 s Display : 492.430 ms PostProcessing : 243.154 ms Resolution hyperelastic : 5.848 s PyVista_Interface : 2.394 s Generate movie 01/30 (3.33 %) 17.33 s Generate movie 02/30 (6.67 %) 9.05 s Generate movie 03/30 (10.00 %) 7.08 s Generate movie 04/30 (13.33 %) 6.85 s Generate movie 05/30 (16.67 %) 6.58 s Generate movie 06/30 (20.00 %) 6.39 s Generate movie 07/30 (23.33 %) 5.98 s Generate movie 08/30 (26.67 %) 5.77 s Generate movie 09/30 (30.00 %) 5.56 s Generate movie 10/30 (33.33 %) 5.29 s Generate movie 11/30 (36.67 %) 5.15 s Generate movie 12/30 (40.00 %) 5.10 s Generate movie 13/30 (43.33 %) 5.05 s Generate movie 14/30 (46.67 %) 4.35 s Generate movie 15/30 (50.00 %) 3.99 s Generate movie 16/30 (53.33 %) 3.87 s Generate movie 17/30 (56.67 %) 3.42 s Generate movie 18/30 (60.00 %) 3.18 s Generate movie 19/30 (63.33 %) 2.94 s Generate movie 20/30 (66.67 %) 2.72 s Generate movie 21/30 (70.00 %) 2.45 s Generate movie 22/30 (73.33 %) 2.12 s Generate movie 23/30 (76.67 %) 1.86 s Generate movie 24/30 (80.00 %) 1.57 s Generate movie 25/30 (83.33 %) 1.31 s Generate movie 26/30 (86.67 %) 1.07 s Generate movie 27/30 (90.00 %) 809.84 ms Generate movie 28/30 (93.33 %) 542.88 ms Generate movie 29/30 (96.67 %) 266.37 ms Generate movie 30/30 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 16 # 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.Results_Dir() 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.Elastic.Isotropic( dim, E=210000, v=0.3, planeStress=True, thickness=thickness ) simu = Simulations.Elastic(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.coord # 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 11.557 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 `_