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.

Static Scene

Interactive Scene

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

# 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()