Note
Go to the end to download the full example code.
Contact2#
Frictionless contact between a deformable block and a rigid indenter, solved as a genuine non-linear problem with Newton-Raphson and a penalty regularisation of the unilateral (Signorini) condition.
The rigid indenter is treated as an obstacle: at every contact-surface Gauss point the normal gap gₙ to the obstacle surface is measured and, where it is negative (penetration), a penalty traction εₙ⟨-gₙ⟩ n resists it.These contributions are added to the elastic residual/tangent through Operators.NonLinear.PenaltyContact and the non-linear system A(u) Δu = -R(u) is solved with Newton at each load step.
Runs in 2D (block + half-disc punch) and 3D (box + half-cylinder punch); set dim below. The non-linear simulation (ElasticContact) lives in _utils.py.
Penalty contact solve in 3D (Newton per step):
======= elastic problem at iteration 0 =======
At Newton iteration 1 norm is 1.185313683882e+05
At Newton iteration 2 norm is 5.812774982115e+02
At Newton iteration 3 norm is 9.901128618810e-09
======= elastic problem at iteration 1 =======
At Newton iteration 1 norm is 1.166673631212e+05
At Newton iteration 2 norm is 1.162554996434e+03
At Newton iteration 3 norm is 7.749237811253e-09
======= elastic problem at iteration 2 =======
At Newton iteration 1 norm is 1.148067412636e+05
At Newton iteration 2 norm is 1.743832494649e+03
At Newton iteration 3 norm is 5.100563749518e-09
======= elastic problem at iteration 3 =======
At Newton iteration 1 norm is 1.129496700202e+05
At Newton iteration 2 norm is 2.325109992869e+03
At Newton iteration 3 norm is 1.008345954691e-08
======= elastic problem at iteration 4 =======
At Newton iteration 1 norm is 1.110963274459e+05
At Newton iteration 2 norm is 2.906387491082e+03
At Newton iteration 3 norm is 1.210719790661e-08
======= elastic problem at iteration 5 =======
At Newton iteration 1 norm is 1.092469033083e+05
At Newton iteration 2 norm is 3.487664989292e+03
At Newton iteration 3 norm is 9.346021680056e-09
======= elastic problem at iteration 6 =======
At Newton iteration 1 norm is 1.074016000308e+05
At Newton iteration 2 norm is 4.068942487516e+03
At Newton iteration 3 norm is 8.768863505757e-09
======= elastic problem at iteration 7 =======
At Newton iteration 1 norm is 1.055606337234e+05
At Newton iteration 2 norm is 4.650219985722e+03
At Newton iteration 3 norm is 7.563893860032e-09
======= elastic problem at iteration 8 =======
At Newton iteration 1 norm is 1.037242353127e+05
At Newton iteration 2 norm is 5.231497483937e+03
At Newton iteration 3 norm is 8.001390881827e-09
======= elastic problem at iteration 9 =======
At Newton iteration 1 norm is 1.018926517788e+05
At Newton iteration 2 norm is 5.812774982156e+03
At Newton iteration 3 norm is 1.478495779020e-08
======= elastic problem at iteration 10 =======
At Newton iteration 1 norm is 1.000661475126e+05
At Newton iteration 2 norm is 6.394052480372e+03
At Newton iteration 3 norm is 1.041112153236e-08
======= elastic problem at iteration 11 =======
At Newton iteration 1 norm is 9.824500580546e+04
At Newton iteration 2 norm is 6.975329978588e+03
At Newton iteration 3 norm is 8.464284619152e-09
======= elastic problem at iteration 12 =======
At Newton iteration 1 norm is 9.642953048542e+04
At Newton iteration 2 norm is 7.556607476799e+03
At Newton iteration 3 norm is 8.758018624992e-09
======= elastic problem at iteration 13 =======
At Newton iteration 1 norm is 9.462004771586e+04
At Newton iteration 2 norm is 8.137884975020e+03
At Newton iteration 3 norm is 1.209765291054e-08
======= elastic problem at iteration 14 =======
At Newton iteration 1 norm is 9.281690797513e+04
At Newton iteration 2 norm is 8.719162473235e+03
At Newton iteration 3 norm is 9.616481315885e-09
======= elastic problem at iteration 15 =======
At Newton iteration 1 norm is 9.102048823516e+04
At Newton iteration 2 norm is 9.300439971447e+03
At Newton iteration 3 norm is 9.611662141835e-09
======= elastic problem at iteration 16 =======
At Newton iteration 1 norm is 8.923119436215e+04
At Newton iteration 2 norm is 9.881717469665e+03
At Newton iteration 3 norm is 7.262590546939e-09
======= elastic problem at iteration 17 =======
At Newton iteration 1 norm is 8.744946376172e+04
At Newton iteration 2 norm is 1.046299496788e+04
At Newton iteration 3 norm is 1.431849773446e-08
======= elastic problem at iteration 18 =======
At Newton iteration 1 norm is 8.567576829723e+04
At Newton iteration 2 norm is 1.104427246609e+04
At Newton iteration 3 norm is 9.918191748510e-09
======= elastic problem at iteration 19 =======
At Newton iteration 1 norm is 8.391061750845e+04
At Newton iteration 2 norm is 1.162554996431e+04
At Newton iteration 3 norm is 9.706379453461e-09
======= elastic problem at iteration 20 =======
At Newton iteration 1 norm is 8.215456216286e+04
At Newton iteration 2 norm is 1.220682746252e+04
At Newton iteration 3 norm is 4.965093719761e-09
======= elastic problem at iteration 21 =======
At Newton iteration 1 norm is 8.040819817587e+04
At Newton iteration 2 norm is 1.278810496073e+04
At Newton iteration 3 norm is 1.387281033992e-08
======= elastic problem at iteration 22 =======
At Newton iteration 1 norm is 7.871095643409e+04
At Newton iteration 2 norm is 1.180866389997e+04
At Newton iteration 3 norm is 1.076353168632e-08
======= elastic problem at iteration 23 =======
At Newton iteration 1 norm is 7.830569184190e+04
At Newton iteration 2 norm is 6.483218661311e-09
======= elastic problem at iteration 24 =======
At Newton iteration 1 norm is 7.830569184189e+04
At Newton iteration 2 norm is 8.257101205397e-09
======= elastic problem at iteration 25 =======
At Newton iteration 1 norm is 7.830569184183e+04
At Newton iteration 2 norm is 1.016022596422e-08
======= elastic problem at iteration 26 =======
At Newton iteration 1 norm is 7.830569184190e+04
At Newton iteration 2 norm is 7.109110221733e-09
======= elastic problem at iteration 27 =======
At Newton iteration 1 norm is 7.830569184190e+04
At Newton iteration 2 norm is 8.274836961939e-09
======= elastic problem at iteration 28 =======
At Newton iteration 1 norm is 7.830569184190e+04
At Newton iteration 2 norm is 7.358984570085e-09
======= elastic problem at iteration 29 =======
At Newton iteration 1 norm is 7.830569184189e+04
At Newton iteration 2 norm is 4.735839961635e-09
==================== Mesh ====================
elemType: HEXA8
Ne: 192
Nn: 324
==================== Model ====================
Isotropic:
E = 2.10e+05, v = 0.3
solver:scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
W def = 2618.76
Svm max = 2910.49
Evm max = 1.80 %
Ux max = 1.62e-02
Ux min = -1.62e-02
Uy max = 0.00e+00
Uy min = -1.06e-01
Uz max = 1.26e-02
Uz min = -1.26e-02
=================== TicTac ===================
Mesh: 56.295 ms
Boundary Conditions: 208.735 µs
Matrix: 651.949 ms
Solver: 623.681 ms
Resolution elastic: 1.380 s
PostProcessing: 433.372 µs
Generate movie 01/30 (3.33 %) 6.05 s
Generate movie 02/30 (6.67 %) 10.47 s
Generate movie 03/30 (10.00 %) 9.94 s
Generate movie 04/30 (13.33 %) 9.59 s
Generate movie 05/30 (16.67 %) 9.50 s
Generate movie 06/30 (20.00 %) 8.89 s
Generate movie 07/30 (23.33 %) 8.64 s
Generate movie 08/30 (26.67 %) 8.30 s
Generate movie 09/30 (30.00 %) 7.81 s
Generate movie 10/30 (33.33 %) 7.41 s
Generate movie 11/30 (36.67 %) 7.14 s
Generate movie 12/30 (40.00 %) 6.76 s
Generate movie 13/30 (43.33 %) 6.38 s
Generate movie 14/30 (46.67 %) 6.05 s
Generate movie 15/30 (50.00 %) 5.66 s
Generate movie 16/30 (53.33 %) 5.31 s
Generate movie 17/30 (56.67 %) 4.89 s
Generate movie 18/30 (60.00 %) 4.60 s
Generate movie 19/30 (63.33 %) 4.21 s
Generate movie 20/30 (66.67 %) 3.80 s
Generate movie 21/30 (70.00 %) 3.41 s
Generate movie 22/30 (73.33 %) 2.99 s
Generate movie 23/30 (76.67 %) 2.70 s
Generate movie 24/30 (80.00 %) 2.30 s
Generate movie 25/30 (83.33 %) 1.90 s
Generate movie 26/30 (86.67 %) 1.52 s
Generate movie 27/30 (90.00 %) 1.14 s
Generate movie 28/30 (93.33 %) 765.82 ms
Generate movie 29/30 (96.67 %) 381.38 ms
Generate movie 30/30 (100.00 %) 0.00 µs
17 from EasyFEA import Terminal, Folder, Models, ElemType, PyVista
18 from EasyFEA.Geoms import Point, Domain, Points
19
20 from _utils import RigidContact
21
22 if __name__ == "__main__":
23 Terminal.Clear()
24
25 # ----------------------------------------------
26 # Configuration
27 # ----------------------------------------------
28
29 dim = 3 # 2 or 3
30 result = "Svm"
31
32 R = 10 # block size
33 thickness = R / 3 # out-of-plane extent (3D)
34 meshSize = R / 20 if dim == 2 else R / 8
35
36 N = 30 # load steps
37 delta = 0.1
38 penalty = 1e7 # contact stiffness (larger -> less penetration)
39
40 folder = Folder.Results_Dir()
41
42 # ----------------------------------------------
43 # Body
44 # ----------------------------------------------
45 body = Domain(
46 (-R / 2, 0),
47 (R / 2, R),
48 meshSize,
49 )
50 if dim == 2:
51 mesh = body.Mesh_2D(
52 [],
53 ElemType.QUAD4,
54 isOrganised=True,
55 )
56 else:
57 nz = max(1, round(thickness / meshSize))
58 mesh = body.Mesh_Extrude(
59 [],
60 [0, 0, thickness],
61 [nz],
62 ElemType.HEXA8,
63 isOrganised=True,
64 )
65 nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
66
67 # ----------------------------------------------
68 # Indenter
69 # ----------------------------------------------
70 r = R / 3
71 contour = Points(
72 [
73 Point(-R / 2, R, r=r),
74 Point(R / 2, R, r=r),
75 Point(R / 2, 2 * R),
76 Point(-R / 2, 2 * R),
77 ]
78 )
79 if dim == 2:
80 indenter = contour.Mesh_2D([], ElemType.TRI3)
81 else:
82 indenter = contour.Mesh_Extrude([], [0, 0, thickness], [nz], ElemType.TETRA4)
83 # lower (contact) surface of the punch: below the end of the fillets
84 nodes_contact = indenter.Nodes_Conditions(lambda x, y, z: y <= R + r + 1e-6)
85
86 indenter.Set_Tag(nodes_contact, "contact")
87
88 # ----------------------------------------------
89 # Simulation
90 # ----------------------------------------------
91 material = Models.Elastic.Isotropic(
92 dim, E=210000, v=0.3, planeStress=True, thickness=thickness
93 )
94 simu = RigidContact(mesh, material, penalty)
95
96 list_indenter = [indenter]
97 print(f"Penalty contact solve in {dim}D (Newton per step):")
98 for i in range(N):
99 # update indeter
100 indenter = list_indenter[0].copy()
101 indenter.Translate(dy=-(i + 1) / N * delta) # lower the rigid indenter
102 list_indenter.append(indenter)
103 simu._contactMesh = indenter
104 # solve contact
105 simu.Bc_Init()
106 simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
107 simu.Solve()
108 simu.Save_Iter()
109
110 print(simu)
111
112 # ----------------------------------------------
113 # Results
114 # ----------------------------------------------
115 def Plot_Iter(plotter, n):
116 simu.Set_Iter(n)
117 PyVista.Plot(
118 simu, result, 1, color="k", nColors=21, show_grid=True, plotter=plotter
119 )
120 PyVista.Plot(list_indenter[n], color="gray", alpha=0.4, plotter=plotter)
121 PyVista.Plot_Elements(
122 list_indenter[n], color="k", dimElem=1, linewidth=2, plotter=plotter
123 )
124 PyVista._setCameraPosition(plotter, mesh.inDim)
125
126 PyVista.Movie_func(Plot_Iter, N, folder=folder, filename=f"{result}.gif")
127
128 plotter = PyVista._Plotter()
129 result = "uy"
130 Plot_Iter(plotter, -1)
131 plotter.show()
Total running time of the script: (0 minutes 15.068 seconds)