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