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