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, dof = 882
==================== 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 : 123.352 ms
Boundary Conditions : 664.234 µs
Matrix : 39.871 ms
Solver : 62.654 ms
Display : 335.538 ms
PostProcessing : 126.215 ms
Movie_func 0/29
Movie_func 1/29 (3.45 %) 8.55 s
Movie_func 2/29 (6.90 %) 6.45 s
Movie_func 3/29 (10.34 %) 6.31 s
Movie_func 4/29 (13.79 %) 5.96 s
Movie_func 5/29 (17.24 %) 5.68 s
Movie_func 6/29 (20.69 %) 5.56 s
Movie_func 7/29 (24.14 %) 5.29 s
Movie_func 8/29 (27.59 %) 5.17 s
Movie_func 9/29 (31.03 %) 4.79 s
Movie_func 10/29 (34.48 %) 4.58 s
Movie_func 11/29 (37.93 %) 4.35 s
Movie_func 12/29 (41.38 %) 4.09 s
Movie_func 13/29 (44.83 %) 3.83 s
Movie_func 14/29 (48.28 %) 3.60 s
Movie_func 15/29 (51.72 %) 3.33 s
Movie_func 16/29 (55.17 %) 3.19 s
Movie_func 17/29 (58.62 %) 2.89 s
Movie_func 18/29 (62.07 %) 2.64 s
Movie_func 19/29 (65.52 %) 2.42 s
Movie_func 20/29 (68.97 %) 2.22 s
Movie_func 21/29 (72.41 %) 1.91 s
Movie_func 22/29 (75.86 %) 1.69 s
Movie_func 23/29 (79.31 %) 1.45 s
Movie_func 24/29 (82.76 %) 1.22 s
Movie_func 25/29 (86.21 %) 980.14 ms
Movie_func 26/29 (89.66 %) 736.16 ms
Movie_func 27/29 (93.10 %) 489.93 ms
Movie_func 28/29 (96.55 %) 244.87 ms
Movie_func 29/29 (100.00 %) 0.00 µs
13 # TODO: Compare results with analytical values ?
14
15 from EasyFEA import (
16 Display,
17 Folder,
18 Models,
19 plt,
20 np,
21 ElemType,
22 Simulations,
23 PyVista,
24 )
25 from EasyFEA.Geoms import Point, Domain, Points
26
27 if __name__ == "__main__":
28 Display.Clear()
29
30 # ----------------------------------------------
31 # Configuration
32 # ----------------------------------------------
33 dim = 2
34
35 # outputs
36 folder = Folder.Join(Folder.RESULTS_DIR, "Contact", "Contact1")
37 pltIter = False
38 makeMovie = True
39 result = "uy"
40
41 # geom
42 R = 10
43 height = R
44 thickness = R / 3
45
46 # load
47 N = 30
48 inc = 1e-0 / N
49 cx, cy = 0, -1
50
51 # ----------------------------------------------
52 # Meshes
53 # ----------------------------------------------
54
55 meshSize = R / 20
56
57 # slave mesh
58 contour_slave = Domain(Point(-R / 2, 0), Point(R / 2, height), meshSize)
59 if dim == 2:
60 mesh_slave = contour_slave.Mesh_2D([], ElemType.QUAD4, isOrganised=True)
61 else:
62 mesh_slave = contour_slave.Mesh_Extrude(
63 [], [0, 0, -thickness], [4], ElemType.HEXA8, isOrganised=True
64 )
65
66 # nodes_slave = mesh_slave.Get_list_groupElem(dim-1)[0].nodes
67 nodes_slave = mesh_slave.Nodes_Conditions(lambda x, y, z: y == height)
68 nodes_y0 = mesh_slave.Nodes_Conditions(lambda x, y, z: y == 0)
69
70 # master mesh
71 r = R / 2
72 p0 = Point(-R / 2, height, r=r)
73 p1 = Point(R / 2, height, r=r)
74 p2 = Point(R / 2, height + R)
75 p3 = Point(-R / 2, height + R)
76 contour_master = Points([p0, p1, p2, p3])
77
78 yMax = height + np.abs(r)
79 if dim == 2:
80 master_mesh = contour_master.Mesh_2D([], ElemType.TRI3)
81 else:
82 master_mesh = contour_master.Mesh_Extrude(
83 [], [0, 0, -thickness - 2], [4], ElemType.TETRA4
84 )
85 groupMaster = master_mesh.Get_list_groupElem(dim - 1)[0]
86 if len(master_mesh.Get_list_groupElem(dim - 1)) > 1:
87 Display.MyPrintError(
88 f"The {groupMaster.elemType.name} element group is used. In 3D, TETRA AND HEXA elements are recommended."
89 )
90 master_mesh.Translate(dz=-(master_mesh.center[2] - mesh_slave.center[2]))
91
92 # Display.Plot_Tags(mesh_master, alpha=0.1, showId=True)
93
94 # get master nodes
95 # nodes_master = mesh_master.Get_list_groupElem(dim-1)[0].nodes
96 if dim == 2:
97 nodes_master = master_mesh.Nodes_Tags(["L0", "L1"])
98 else:
99 nodes_master = master_mesh.Nodes_Tags(["S1", "S2"])
100
101 # # plot meshes
102 # ax = Display.Plot_Mesh(master_mesh, alpha=0)
103 # Display.Plot_Mesh(mesh_slave, ax=ax, alpha=0)
104 # # add nodes interface
105 # ax.scatter(*mesh_slave.coord[nodes_slave, :dim].T, label="slave nodes")
106 # ax.scatter(*master_mesh.coord[nodes_master, :dim].T, label="master nodes")
107 # ax.legend()
108 # ax.set_title("Contact nodes")
109
110 # ----------------------------------------------
111 # Simulation
112 # ----------------------------------------------
113 material = Models.ElasIsot(
114 dim, E=210000, v=0.3, planeStress=True, thickness=thickness
115 )
116 simu = Simulations.ElasticSimu(mesh_slave, material)
117
118 list_master_mesh = [master_mesh]
119
120 if pltIter:
121 ax = Display.Plot_Result(simu, result, deformFactor=1)
122
123 for i in range(N):
124 master_mesh = master_mesh.copy()
125 master_mesh.Translate(cx * inc, cy * inc)
126
127 list_master_mesh.append(master_mesh)
128
129 convergence = False
130
131 coordo_old = simu.Results_displacement_matrix() + simu.mesh.coord
132
133 while not convergence:
134 # apply new boundary conditions
135 simu.Bc_Init()
136 simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
137
138 nodes, newU = simu.Get_contact(master_mesh, nodes_slave, nodes_master)
139
140 if nodes.size > 0:
141 simu.add_dirichlet(nodes, [newU[:, 0], newU[:, 1]], ["x", "y"])
142
143 simu.Solve()
144
145 # check if there is no new nodes in the master mesh
146 oldSize = nodes.size
147 nodes, __ = simu.Get_contact(master_mesh, nodes_slave, nodes_master)
148 convergence = oldSize == nodes.size
149
150 simu.Save_Iter()
151
152 print(f"Eps max = {simu.Result('Strain').max() * 100:3.2f} %")
153
154 if pltIter:
155 Display.Plot_Result(simu, result, plotMesh=True, deformFactor=1, ax=ax)
156 Display.Plot_Mesh(master_mesh, alpha=0, ax=ax)
157 ax.set_title(result)
158 if dim == 3:
159 Display._Axis_equal_3D(
160 ax, np.concatenate((master_mesh.coord, mesh_slave.coord), 0)
161 )
162
163 # # Plot arrows
164 # if nodes.size >0:
165 # # get the nodes coordinates on the interface
166 # coordinates = groupMaster.Get_GaussCoordinates_e_p('mass').reshape(-1,3)
167 # ax.scatter(*coordinates[:,:dim].T)
168
169 # coordo_new = simu.Results_displacement_matrix() + simu.mesh.coordo
170 # ax.scatter(*coordo_old[nodes,:dim].T)
171 # incU = coordo_new - coordo_oldq
172 # [ax.arrow(*coordo_old[node, :dim], *incU[node,:dim],length_includes_head=True) for node in nodes]
173
174 plt.pause(1e-12)
175
176 print(simu)
177
178 # ----------------------------------------------
179 # Results
180 # ----------------------------------------------
181 if makeMovie:
182
183 def DoAnim(plotter, n):
184 simu.Set_Iter(n)
185 PyVista.Plot(
186 simu,
187 "Svm",
188 1,
189 style="surface",
190 color="k",
191 plotter=plotter,
192 nColors=10,
193 show_grid=True,
194 )
195 PyVista.Plot(list_master_mesh[n], plotter=plotter, plotMesh=True, alpha=0.2)
196
197 PyVista.Movie_func(DoAnim, N, folder=folder, filename=f"{result}.gif")
198
199 if not pltIter:
200 plotter = PyVista.Plot(simu, result, plotMesh=True, deformFactor=1)
201 PyVista.Plot_Mesh(master_mesh, alpha=0.4, plotter=plotter)
202 plotter.show()
203
204 Display.plt.show()
Total running time of the script: (0 minutes 10.550 seconds)