Note
Go to the end to download the full example code.
Contact1#
A rigid indenter is pressed into an elastic half-space (2D, plane strain) and the finite-element contact pressure is compared to the closed-form solutions, with E* = E/(1-ν²):
parabola
z = x²/2R(Hertz line contact):p(x) = p₀√(1-x²/a²),p₀ = a·E*/(2R)— a semi-ellipse.wedge
z = |x|·tanθ:p(x) = (E*·tanθ/π)·arccosh(a/|x|)— a log peak.
The contact half-width a is read from the FE solution; the analytical curve is plotted with that same a (so the comparison is shape + magnitude).
The contact pressure is obtained directly from the penalty contact as p = εₙ⟨-g⟩ on the surface.
[parabola] pressing the rigid indenter (Newton per step):
======= elastic problem at iteration 0 =======
At Newton iteration 1 norm is 5.536805892087e+03
At Newton iteration 2 norm is 1.336298576388e+02
At Newton iteration 3 norm is 4.858189561363e+01
At Newton iteration 4 norm is 1.965806277962e+01
At Newton iteration 5 norm is 6.146824670310e+00
At Newton iteration 6 norm is 3.769446750457e-12
======= elastic problem at iteration 1 =======
At Newton iteration 1 norm is 6.924212846549e+03
At Newton iteration 2 norm is 1.507335478405e+02
At Newton iteration 3 norm is 7.419047099757e+01
At Newton iteration 4 norm is 3.478426175903e+01
At Newton iteration 5 norm is 1.885129836358e+01
At Newton iteration 6 norm is 2.824259659696e+00
At Newton iteration 7 norm is 8.257509454179e-12
======= elastic problem at iteration 2 =======
At Newton iteration 1 norm is 7.711403046518e+03
At Newton iteration 2 norm is 1.774220663138e+02
At Newton iteration 3 norm is 6.682079946032e+01
At Newton iteration 4 norm is 2.735547015955e+01
At Newton iteration 5 norm is 9.302345940108e+00
At Newton iteration 6 norm is 1.340417369513e-11
======= elastic problem at iteration 3 =======
At Newton iteration 1 norm is 8.313675151848e+03
At Newton iteration 2 norm is 1.803032438628e+02
At Newton iteration 3 norm is 6.533342138905e+01
At Newton iteration 4 norm is 2.405611452238e+01
At Newton iteration 5 norm is 3.948698407378e+00
At Newton iteration 6 norm is 1.821760791086e-11
======= elastic problem at iteration 4 =======
At Newton iteration 1 norm is 8.818907274369e+03
At Newton iteration 2 norm is 1.594557717343e+02
At Newton iteration 3 norm is 1.004618205320e+02
At Newton iteration 4 norm is 4.640474059448e+01
At Newton iteration 5 norm is 2.459138270659e+01
At Newton iteration 6 norm is 3.643228916354e+00
At Newton iteration 7 norm is 2.353048230135e-11
======= elastic problem at iteration 5 =======
At Newton iteration 1 norm is 9.254979631209e+03
At Newton iteration 2 norm is 1.693663162754e+02
At Newton iteration 3 norm is 1.093546985568e+02
At Newton iteration 4 norm is 5.030820094599e+01
At Newton iteration 5 norm is 5.338180308909e+00
At Newton iteration 6 norm is 2.800552478699e-11
======= elastic problem at iteration 6 =======
At Newton iteration 1 norm is 9.644403746802e+03
At Newton iteration 2 norm is 1.837055602393e+02
At Newton iteration 3 norm is 1.199733597807e+02
At Newton iteration 4 norm is 5.977487725063e+01
At Newton iteration 5 norm is 9.775216534900e+00
At Newton iteration 6 norm is 3.320535179354e-11
======= elastic problem at iteration 7 =======
At Newton iteration 1 norm is 9.994142186081e+03
At Newton iteration 2 norm is 1.643649443689e+02
At Newton iteration 3 norm is 7.100589654474e+01
At Newton iteration 4 norm is 1.885584229551e+01
At Newton iteration 5 norm is 4.040741247376e-11
======= elastic problem at iteration 8 =======
At Newton iteration 1 norm is 1.031556506142e+04
At Newton iteration 2 norm is 1.846666676678e+02
At Newton iteration 3 norm is 8.634216166540e+01
At Newton iteration 4 norm is 2.843868240077e+01
At Newton iteration 5 norm is 4.901650786398e+00
At Newton iteration 6 norm is 4.610704407389e-11
======= elastic problem at iteration 9 =======
At Newton iteration 1 norm is 1.061520518188e+04
At Newton iteration 2 norm is 1.735699098822e+02
At Newton iteration 3 norm is 1.051266955796e+02
At Newton iteration 4 norm is 4.206082216874e+01
At Newton iteration 5 norm is 1.602554912949e+01
At Newton iteration 6 norm is 5.186825359889e-11
======= elastic problem at iteration 10 =======
At Newton iteration 1 norm is 1.089468953721e+04
At Newton iteration 2 norm is 1.967415724883e+02
At Newton iteration 3 norm is 1.258063637510e+02
At Newton iteration 4 norm is 5.776602967803e+01
At Newton iteration 5 norm is 3.032635689800e+01
At Newton iteration 6 norm is 3.437445123131e+00
At Newton iteration 7 norm is 5.970618529298e-11
======= elastic problem at iteration 11 =======
At Newton iteration 1 norm is 1.116302764789e+04
At Newton iteration 2 norm is 1.817153901348e+02
At Newton iteration 3 norm is 7.847938864843e+01
At Newton iteration 4 norm is 1.951932688858e+01
At Newton iteration 5 norm is 6.314203785870e-11
[parabola] a=0.419 error vs analytical over 0.1a-0.85a: mean 0% max 1%
Generate movie 01/12 (8.33 %) 2.66 s
Generate movie 02/12 (16.67 %) 3.29 s
Generate movie 03/12 (25.00 %) 2.78 s
Generate movie 04/12 (33.33 %) 2.50 s
Generate movie 05/12 (41.67 %) 2.18 s
Generate movie 06/12 (50.00 %) 1.89 s
Generate movie 07/12 (58.33 %) 1.56 s
Generate movie 08/12 (66.67 %) 1.28 s
Generate movie 09/12 (75.00 %) 965.52 ms
Generate movie 10/12 (83.33 %) 628.89 ms
Generate movie 11/12 (91.67 %) 313.71 ms
Generate movie 12/12 (100.00 %) 0.00 µs
[wedge] pressing the rigid indenter (Newton per step):
======= elastic problem at iteration 0 =======
At Newton iteration 1 norm is 6.736266687598e+03
At Newton iteration 2 norm is 7.238066934825e+02
At Newton iteration 3 norm is 1.455169662229e+02
At Newton iteration 4 norm is 8.527654385270e+01
At Newton iteration 5 norm is 6.179022523118e-12
======= elastic problem at iteration 1 =======
At Newton iteration 1 norm is 1.044830094297e+04
At Newton iteration 2 norm is 7.218267191013e+02
At Newton iteration 3 norm is 1.666526524070e+02
At Newton iteration 4 norm is 5.905000527297e+01
At Newton iteration 5 norm is 1.391936726019e-11
======= elastic problem at iteration 2 =======
At Newton iteration 1 norm is 1.310754115797e+04
At Newton iteration 2 norm is 7.024703210645e+02
At Newton iteration 3 norm is 1.293800126937e+02
At Newton iteration 4 norm is 8.889876310227e+01
At Newton iteration 5 norm is 4.066926771745e+01
At Newton iteration 6 norm is 2.356942823635e-11
======= elastic problem at iteration 3 =======
At Newton iteration 1 norm is 1.533010007974e+04
At Newton iteration 2 norm is 6.850037859907e+02
At Newton iteration 3 norm is 1.535887444124e+02
At Newton iteration 4 norm is 6.162716758653e+01
At Newton iteration 5 norm is 1.626376027190e+01
At Newton iteration 6 norm is 3.372581795418e-11
======= elastic problem at iteration 4 =======
At Newton iteration 1 norm is 1.735857879453e+04
At Newton iteration 2 norm is 6.684958294746e+02
At Newton iteration 3 norm is 1.491381961254e+02
At Newton iteration 4 norm is 6.807775395264e+01
At Newton iteration 5 norm is 3.113727608740e+01
At Newton iteration 6 norm is 4.471076356498e-11
======= elastic problem at iteration 5 =======
At Newton iteration 1 norm is 1.926247104997e+04
At Newton iteration 2 norm is 6.706878251850e+02
At Newton iteration 3 norm is 1.570830389009e+02
At Newton iteration 4 norm is 7.411624608189e+01
At Newton iteration 5 norm is 3.763981102320e+01
At Newton iteration 6 norm is 1.547822566224e+00
At Newton iteration 7 norm is 5.602731477211e-11
======= elastic problem at iteration 6 =======
At Newton iteration 1 norm is 2.106920123204e+04
At Newton iteration 2 norm is 6.495395617115e+02
At Newton iteration 3 norm is 1.428840484679e+02
At Newton iteration 4 norm is 7.049824480479e+01
At Newton iteration 5 norm is 4.051662604050e+01
At Newton iteration 6 norm is 6.864918385469e+00
At Newton iteration 7 norm is 6.818971394499e-11
======= elastic problem at iteration 7 =======
At Newton iteration 1 norm is 2.276869018696e+04
At Newton iteration 2 norm is 6.489914597175e+02
At Newton iteration 3 norm is 1.388657797508e+02
At Newton iteration 4 norm is 6.516506777748e+01
At Newton iteration 5 norm is 5.443723725811e+00
At Newton iteration 6 norm is 7.965796497607e-11
======= elastic problem at iteration 8 =======
At Newton iteration 1 norm is 2.444117929985e+04
At Newton iteration 2 norm is 6.277568072516e+02
At Newton iteration 3 norm is 1.165282274664e+02
At Newton iteration 4 norm is 5.552777913147e+01
At Newton iteration 5 norm is 2.950391665648e+01
At Newton iteration 6 norm is 9.563669248353e-11
======= elastic problem at iteration 9 =======
At Newton iteration 1 norm is 2.604666450243e+04
At Newton iteration 2 norm is 6.234931195458e+02
At Newton iteration 3 norm is 1.335466139163e+02
At Newton iteration 4 norm is 4.520648797953e+01
At Newton iteration 5 norm is 1.786100859015e+01
At Newton iteration 6 norm is 1.096941931236e-10
======= elastic problem at iteration 10 =======
At Newton iteration 1 norm is 2.761454055922e+04
At Newton iteration 2 norm is 6.169702866905e+02
At Newton iteration 3 norm is 1.441537783234e+02
At Newton iteration 4 norm is 5.703103385829e+01
At Newton iteration 5 norm is 6.072912967074e+00
At Newton iteration 6 norm is 1.222048763634e-10
======= elastic problem at iteration 11 =======
At Newton iteration 1 norm is 2.915203671188e+04
At Newton iteration 2 norm is 6.096796538288e+02
At Newton iteration 3 norm is 1.495307684225e+02
At Newton iteration 4 norm is 6.623428440921e+01
At Newton iteration 5 norm is 1.522177670189e+01
At Newton iteration 6 norm is 1.395931823722e-10
[wedge] a=0.473 error vs analytical over 0.1a-0.85a: mean 0% max 0%
Generate movie 01/12 (8.33 %) 1.87 s
Generate movie 02/12 (16.67 %) 3.18 s
Generate movie 03/12 (25.00 %) 2.83 s
Generate movie 04/12 (33.33 %) 2.56 s
Generate movie 05/12 (41.67 %) 2.20 s
Generate movie 06/12 (50.00 %) 1.91 s
Generate movie 07/12 (58.33 %) 1.60 s
Generate movie 08/12 (66.67 %) 1.28 s
Generate movie 09/12 (75.00 %) 945.10 ms
Generate movie 10/12 (83.33 %) 638.28 ms
Generate movie 11/12 (91.67 %) 316.75 ms
Generate movie 12/12 (100.00 %) 0.00 µs
22 import numpy as np
23 import matplotlib.pyplot as plt
24
25 from EasyFEA import Terminal, Folder, Models, ElemType, Mesh, PyVista
26 from EasyFEA.Geoms import Domain, Points, Circle
27 from EasyFEA.FEM import MatrixType
28
29 from _utils import RigidContact
30
31
32 # ----------------------------------------------
33 # Rigid indenter profiles z(x), spanning [0, W] and rising away from x=0
34 # ----------------------------------------------
35 def parabola_indenter() -> Mesh:
36 # cluster samples near x=0: a curved indenter is gauged by nearest-sample
37 # projection, so it must be resolved finely over the small contact (a << W),
38 # otherwise that sampling — not the body mesh — dominates the pressure error
39 xs = W * np.linspace(0, 1, 200) ** 2
40 geom = Points(
41 [
42 *[(x, x**2 / (2 * Rc)) for x in xs], # lower bound
43 *[(W, W**2 / (2 * Rc) + 2), (0, 2)],
44 ]
45 )
46 mesh = geom.Mesh_2D([], ElemType.TRI3)
47 nodes = mesh.Nodes_Conditions(lambda x, y, z: y <= x**2 / (2 * Rc) + 1e-9)
48 mesh.Set_Tag(nodes, "contact")
49 return mesh
50
51
52 def wedge_indenter() -> Mesh:
53 geom = Points(
54 [
55 (0, 0),
56 (W, W * np.tan(theta)),
57 (W, W * np.tan(theta) + 2),
58 (0, 2),
59 ]
60 )
61 mesh = geom.Mesh_2D([], ElemType.TRI3)
62 nodes = mesh.Nodes_Conditions(lambda x, y, z: y <= x * np.tan(theta) + 1e-9)
63 mesh.Set_Tag(nodes, "contact")
64 return mesh
65
66
67 def parabola_analytic(x, a):
68 return a * Es / (2 * Rc) * np.sqrt(np.clip(1 - (x / a) ** 2, 0, None))
69
70
71 def wedge_analytic(x, a):
72 return Es * np.tan(theta) / np.pi * np.arccosh(np.clip(a / np.abs(x), 1.0, None))
73
74
75 indenter_cases = {
76 "parabola": (parabola_indenter, parabola_analytic),
77 "wedge": (wedge_indenter, wedge_analytic),
78 }
79
80
81 def contact_pressure(simu: RigidContact, indenter: Mesh):
82 """Penalty contact pressure ``εₙ⟨-g⟩`` and x-coordinate on the surface Gauss points."""
83
84 matrixType = MatrixType.mass
85
86 list_p = []
87 list_x = []
88
89 for group1d in simu.mesh.Get_list_groupElem(1):
90 N_pg = group1d.Get_N_pg(matrixType)[:, 0, :]
91 # x = X + u
92 elements = group1d.Get_Elements_Tag("top")
93 X_e_pg = group1d.Get_GaussCoordinates_e_pg(matrixType)[elements]
94 u_e_pg = simu.displacement.reshape(simu.mesh.Nn, 2)
95 x_e_pg = X_e_pg.copy()
96 x_e_pg[..., :2] += np.einsum(
97 "pn,enc->epc", N_pg, u_e_pg[group1d.connect[elements]]
98 )
99 list_x.extend(x_e_pg[..., 0].ravel())
100 # get pressure
101 for groupIndent in indenter.Get_list_groupElem(1):
102 elements = (
103 groupIndent.Get_Elements_Tag("contact")
104 if "contact" in groupIndent.elementTags
105 else None
106 )
107 gap_e_pg, _ = groupIndent._Get_gap_and_normal(
108 x_e_pg=x_e_pg,
109 elements=elements,
110 coord=indenter.center,
111 matrixType=matrixType,
112 )
113 p_e_pg = simu.penalty * np.maximum(-gap_e_pg, 0.0)
114 list_p.extend(p_e_pg.ravel())
115
116 return np.asarray(list_p), np.asarray(list_x)
117
118
119 if __name__ == "__main__":
120 Terminal.Clear()
121
122 # ----------------------------------------------
123 # Configuration
124 # ----------------------------------------------
125
126 folder = Folder.Results_Dir()
127 result = "Svm"
128
129 E, v = 210000.0, 0.3
130 Es = E / (1 - v**2)
131
132 W, D = 6.0, 6.0 # elastic half-space (symmetry about x=0): wide & deep vs contact
133 meshSize = W / 10
134 N = 12 # indentation steps
135 penalty = 100 * Es
136
137 Rc = 6.0 # parabola radius of curvature
138 theta = np.deg2rad(6) # wedge half-angle
139 indeter_delta = {
140 "parabola": 0.05,
141 "wedge": 0.12,
142 } # indentation depth per case
143
144 # ----------------------------------------------
145 # Elastic half-space (surface at y=0)
146 # ----------------------------------------------
147 body = Domain((0, -D), (W, 0), meshSize)
148
149 refineGeoms = [
150 Circle((0, 0), W * coef * 2, W / N)
151 for coef, N in zip(
152 [0.2, 0.4],
153 [200, 100],
154 )
155 ]
156
157 mesh = body.Mesh_2D([], ElemType.TRI3, refineGeoms=refineGeoms)
158 nodes_bottom = mesh.Nodes_Conditions(lambda x, y, z: y == -D)
159 nodes_sym = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
160 nodes_top = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
161 mesh.Set_Tag(nodes_top, "top")
162
163 PyVista.Plot_Mesh(mesh).show()
164
165 material = Models.Elastic.Isotropic(2, E=E, v=v, planeStress=False)
166
167 # ----------------------------------------------
168 # Solve each case and compare FE pressure to theory
169 # ----------------------------------------------
170
171 fig, axes = plt.subplots(1, len(indenter_cases), figsize=(11, 4))
172 for ax, (name, (build, analytic)) in zip(axes, indenter_cases.items()):
173 simu = RigidContact(mesh, material, penalty)
174 indenter = build()
175
176 list_indeter = [indenter]
177 delta = indeter_delta[name]
178
179 print(f"\n[{name}] pressing the rigid indenter (Newton per step):")
180 for i in range(N):
181 # update indeter
182 indenter = list_indeter[0]
183 indenter = indenter.copy()
184 indenter.Translate(dy=-(i + 1) / N * delta)
185 list_indeter.append(indenter)
186 simu._contactMesh = indenter
187 # solve contact
188 simu.Bc_Init()
189 simu.add_dirichlet(nodes_bottom, [0, 0], simu.Get_unknowns())
190 simu.add_dirichlet(nodes_sym, [0], ["x"])
191 simu.Solve()
192 simu.Save_Iter()
193
194 pg, xg = contact_pressure(simu, indenter)
195 order = np.argsort(xg)
196 xg, pg = xg[order], pg[order]
197 a = xg[pg > 1e-6 * pg.max()].max() # FE contact half-width
198 # relative error over the inner patch (skip the singular center & edge),
199 # normalised by the FE peak pressure
200 band = (xg > 0.1 * a) & (xg < 0.85 * a)
201 rel = np.abs(pg[band] - analytic(xg[band], a)) / pg.max()
202 print(
203 f"\n[{name}] a={a:.3f} error vs analytical over 0.1a-0.85a: "
204 f"mean {rel.mean() * 100:.0f}% max {rel.max() * 100:.0f}%"
205 )
206
207 # ----------------------------------------------
208 # Results
209 # ----------------------------------------------
210
211 ax.plot(xg, pg, "o", ms=3, label="FE εₙ⟨-g⟩")
212 xa = np.linspace(1e-3, a, 200)
213 ax.plot(xa, analytic(xa, a), "k-", label="analytical")
214 ax.set_title(f"{name} (a = {a:.3f})")
215 ax.set_xlabel("x")
216 ax.set_ylabel("contact pressure")
217 ax.set_xlim(0, 1.6 * a)
218 ax.legend()
219 ax.grid(True)
220
221 def Plot_Iter(plotter: PyVista.pv.Plotter, n):
222 simu.Set_Iter(n)
223 PyVista.Plot(
224 simu, result, 1, color="k", nColors=10, show_grid=True, plotter=plotter
225 )
226 PyVista.Plot(list_indeter[n], color="gray", plotter=plotter)
227 PyVista.Plot_Elements(
228 list_indeter[n], dimElem=1, color="k", linewidth=1.5, plotter=plotter
229 )
230 PyVista._setCameraPosition(plotter, 2)
231
232 PyVista.Movie_func(Plot_Iter, N, folder=folder, filename=f"{name}.gif")
233
234 fig.tight_layout()
235
236 plotter = PyVista._Plotter()
237 Plot_Iter(plotter, -1)
238 plotter.show()
239
240 plt.show()
Total running time of the script: (0 minutes 16.557 seconds)