Note
Go to the end to download the full example code.
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.
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 : 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 : 229.037 ms
Boundary Conditions : 737.429 µs
Matrix : 1.928 s
Solver : 1.995 s
Display : 336.066 ms
PostProcessing : 126.329 ms
Resolution hyperelastic : 3.863 s
PyVista_Interface : 1.650 s
Generate movie 01/30 (3.33 %) 13.52 s
Generate movie 02/30 (6.67 %) 6.61 s
Generate movie 03/30 (10.00 %) 4.94 s
Generate movie 04/30 (13.33 %) 4.95 s
Generate movie 05/30 (16.67 %) 4.76 s
Generate movie 06/30 (20.00 %) 4.66 s
Generate movie 07/30 (23.33 %) 4.22 s
Generate movie 08/30 (26.67 %) 4.06 s
Generate movie 09/30 (30.00 %) 4.41 s
Generate movie 10/30 (33.33 %) 3.75 s
Generate movie 11/30 (36.67 %) 3.57 s
Generate movie 12/30 (40.00 %) 3.35 s
Generate movie 13/30 (43.33 %) 3.23 s
Generate movie 14/30 (46.67 %) 3.08 s
Generate movie 15/30 (50.00 %) 2.82 s
Generate movie 16/30 (53.33 %) 2.66 s
Generate movie 17/30 (56.67 %) 2.43 s
Generate movie 18/30 (60.00 %) 2.26 s
Generate movie 19/30 (63.33 %) 2.06 s
Generate movie 20/30 (66.67 %) 1.90 s
Generate movie 21/30 (70.00 %) 1.71 s
Generate movie 22/30 (73.33 %) 1.51 s
Generate movie 23/30 (76.67 %) 1.31 s
Generate movie 24/30 (80.00 %) 1.12 s
Generate movie 25/30 (83.33 %) 959.23 ms
Generate movie 26/30 (86.67 %) 745.57 ms
Generate movie 27/30 (90.00 %) 607.61 ms
Generate movie 28/30 (93.33 %) 385.19 ms
Generate movie 29/30 (96.67 %) 188.72 ms
Generate movie 30/30 (100.00 %) 0.00 µs
13 # TODO: Compare results with analytical values ?
14
15 from EasyFEA import Display, Folder, Models, plt, np, ElemType, Simulations, PyVista
16 from EasyFEA.Geoms import Point, Domain, Points
17
18 if __name__ == "__main__":
19 Display.Clear()
20
21 # ----------------------------------------------
22 # Configuration
23 # ----------------------------------------------
24 dim = 2
25
26 # outputs
27 folder = Folder.Join(Folder.RESULTS_DIR, "Contact", "Contact1")
28 pltIter = False
29 makeMovie = True
30 result = "uy"
31
32 # geom
33 R = 10
34 height = R
35 thickness = R / 3
36
37 # load
38 N = 30
39 inc = 1e-0 / N
40 cx, cy = 0, -1
41
42 # ----------------------------------------------
43 # Meshes
44 # ----------------------------------------------
45
46 meshSize = R / 20
47
48 # slave mesh
49 contour_slave = Domain(Point(-R / 2, 0), Point(R / 2, height), meshSize)
50 if dim == 2:
51 mesh_slave = contour_slave.Mesh_2D([], ElemType.QUAD4, isOrganised=True)
52 else:
53 mesh_slave = contour_slave.Mesh_Extrude(
54 [], [0, 0, -thickness], [4], ElemType.HEXA8, isOrganised=True
55 )
56
57 # nodes_slave = mesh_slave.Get_list_groupElem(dim-1)[0].nodes
58 nodes_slave = mesh_slave.Nodes_Conditions(lambda x, y, z: y == height)
59 nodes_y0 = mesh_slave.Nodes_Conditions(lambda x, y, z: y == 0)
60
61 # master mesh
62 r = R / 2
63 p0 = Point(-R / 2, height, r=r)
64 p1 = Point(R / 2, height, r=r)
65 p2 = Point(R / 2, height + R)
66 p3 = Point(-R / 2, height + R)
67 contour_master = Points([p0, p1, p2, p3])
68
69 yMax = height + np.abs(r)
70 if dim == 2:
71 master_mesh = contour_master.Mesh_2D([], ElemType.TRI3)
72 else:
73 master_mesh = contour_master.Mesh_Extrude(
74 [], [0, 0, -thickness - 2], [4], ElemType.TETRA4
75 )
76 groupMaster = master_mesh.Get_list_groupElem(dim - 1)[0]
77 if len(master_mesh.Get_list_groupElem(dim - 1)) > 1:
78 Display.MyPrintError(
79 f"The {groupMaster.elemType.name} element group is used. In 3D, TETRA AND HEXA elements are recommended."
80 )
81 master_mesh.Translate(dz=-(master_mesh.center[2] - mesh_slave.center[2]))
82
83 # Display.Plot_Tags(mesh_master, alpha=0.1, showId=True)
84
85 # get master nodes
86 # nodes_master = mesh_master.Get_list_groupElem(dim-1)[0].nodes
87 if dim == 2:
88 nodes_master = master_mesh.Nodes_Tags(["L0", "L1"])
89 else:
90 nodes_master = master_mesh.Nodes_Tags(["S1", "S2"])
91
92 # # plot meshes
93 # ax = Display.Plot_Mesh(master_mesh, alpha=0)
94 # Display.Plot_Mesh(mesh_slave, ax=ax, alpha=0)
95 # # add nodes interface
96 # ax.scatter(*mesh_slave.coord[nodes_slave, :dim].T, label="slave nodes")
97 # ax.scatter(*master_mesh.coord[nodes_master, :dim].T, label="master nodes")
98 # ax.legend()
99 # ax.set_title("Contact nodes")
100
101 # ----------------------------------------------
102 # Simulation
103 # ----------------------------------------------
104 material = Models.ElasIsot(
105 dim, E=210000, v=0.3, planeStress=True, thickness=thickness
106 )
107 simu = Simulations.ElasticSimu(mesh_slave, material)
108
109 list_master_mesh = [master_mesh]
110
111 if pltIter:
112 ax = Display.Plot_Result(simu, result, deformFactor=1)
113
114 for i in range(N):
115 master_mesh = master_mesh.copy()
116 master_mesh.Translate(cx * inc, cy * inc)
117
118 list_master_mesh.append(master_mesh)
119
120 convergence = False
121
122 coordo_old = simu.Results_displacement_matrix() + simu.mesh.coord
123
124 while not convergence:
125 # apply new boundary conditions
126 simu.Bc_Init()
127 simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
128
129 nodes, newU = simu.Get_contact(master_mesh, nodes_slave, nodes_master)
130
131 if nodes.size > 0:
132 simu.add_dirichlet(nodes, [newU[:, 0], newU[:, 1]], ["x", "y"])
133
134 simu.Solve()
135
136 # check if there is no new nodes in the master mesh
137 oldSize = nodes.size
138 nodes, __ = simu.Get_contact(master_mesh, nodes_slave, nodes_master)
139 convergence = oldSize == nodes.size
140
141 simu.Save_Iter()
142
143 print(f"Eps max = {simu.Result('Strain').max() * 100:3.2f} %")
144
145 if pltIter:
146 Display.Plot_Result(simu, result, plotMesh=True, deformFactor=1, ax=ax)
147 Display.Plot_Mesh(master_mesh, alpha=0, ax=ax)
148 ax.set_title(result)
149 if dim == 3:
150 Display._Axis_equal_3D(
151 ax, np.concatenate((master_mesh.coord, mesh_slave.coord), 0)
152 )
153
154 # # Plot arrows
155 # if nodes.size >0:
156 # # get the nodes coordinates on the interface
157 # coordinates = groupMaster.Get_GaussCoordinates_e_p('mass').reshape(-1,3)
158 # ax.scatter(*coordinates[:,:dim].T)
159
160 # coordo_new = simu.Results_displacement_matrix() + simu.mesh.coord
161 # ax.scatter(*coordo_old[nodes,:dim].T)
162 # incU = coordo_new - coordo_oldq
163 # [ax.arrow(*coordo_old[node, :dim], *incU[node,:dim],length_includes_head=True) for node in nodes]
164
165 plt.pause(1e-12)
166
167 print(simu)
168
169 # ----------------------------------------------
170 # Results
171 # ----------------------------------------------
172 if makeMovie:
173
174 def DoAnim(plotter, n):
175 simu.Set_Iter(n)
176 PyVista.Plot(
177 simu,
178 "Svm",
179 1,
180 style="surface",
181 color="k",
182 plotter=plotter,
183 nColors=10,
184 show_grid=True,
185 )
186 PyVista.Plot(list_master_mesh[n], plotter=plotter, plotMesh=True, alpha=0.2)
187
188 PyVista.Movie_func(DoAnim, N, folder=folder, filename=f"{result}.gif")
189
190 if not pltIter:
191 plotter = PyVista.Plot(simu, result, plotMesh=True, deformFactor=1)
192 PyVista.Plot_Mesh(master_mesh, alpha=0.4, plotter=plotter)
193 plotter.show()
194
195 Display.plt.show()
Total running time of the script: (0 minutes 8.471 seconds)