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