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 ====================
elemType: 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.506 ms
Boundary Conditions: 310.898 µs
PostProcessing: 193.908 ms
Matrix: 9.120 ms
Solver: 82.753 ms
Generate movie 01/30 (3.33 %) 8.15 s
Generate movie 02/30 (6.67 %) 8.87 s
Generate movie 03/30 (10.00 %) 6.88 s
Generate movie 04/30 (13.33 %) 6.51 s
Generate movie 05/30 (16.67 %) 6.28 s
Generate movie 06/30 (20.00 %) 6.04 s
Generate movie 07/30 (23.33 %) 5.86 s
Generate movie 08/30 (26.67 %) 5.58 s
Generate movie 09/30 (30.00 %) 5.39 s
Generate movie 10/30 (33.33 %) 5.14 s
Generate movie 11/30 (36.67 %) 4.92 s
Generate movie 12/30 (40.00 %) 4.60 s
Generate movie 13/30 (43.33 %) 4.37 s
Generate movie 14/30 (46.67 %) 4.07 s
Generate movie 15/30 (50.00 %) 3.87 s
Generate movie 16/30 (53.33 %) 3.62 s
Generate movie 17/30 (56.67 %) 3.36 s
Generate movie 18/30 (60.00 %) 3.08 s
Generate movie 19/30 (63.33 %) 2.82 s
Generate movie 20/30 (66.67 %) 2.59 s
Generate movie 21/30 (70.00 %) 2.34 s
Generate movie 22/30 (73.33 %) 2.04 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.29 s
Generate movie 26/30 (86.67 %) 1.03 s
Generate movie 27/30 (90.00 %) 763.33 ms
Generate movie 28/30 (93.33 %) 512.28 ms
Generate movie 29/30 (96.67 %) 256.85 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 Terminal, Matplotlib, Folder, Models, ElemType, Simulations, PyVista
23 from EasyFEA.Geoms import Point, Domain, Points
24
25 if __name__ == "__main__":
26 Terminal.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 Terminal.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 # Matplotlib.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 = Matplotlib.Plot_Mesh(master_mesh, alpha=0)
101 # Matplotlib.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 = Matplotlib.Plot(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 Matplotlib.Plot(simu, result, plotMesh=True, deformFactor=1, ax=ax)
154 Matplotlib.Plot_Mesh(master_mesh, alpha=0, ax=ax)
155 ax.set_title(result)
156 if dim == 3:
157 Matplotlib._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.410 seconds)