Note
Go to the end to download the full example code.
Contact3#
Frictionless contact between a thin elastic arch strip and a rigid block, solved with the penalty method and Newton-Raphson.
The rigid block 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 (arch strip) and 3D (extruded arch strip); 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 3.317591015983e+05
At Newton iteration 2 norm is 8.823302971212e+05
At Newton iteration 3 norm is 4.817575525634e+04
At Newton iteration 4 norm is 5.654339121181e+03
At Newton iteration 5 norm is 1.309684810071e+03
At Newton iteration 6 norm is 1.076949440762e+02
At Newton iteration 7 norm is 9.287557323790e+01
At Newton iteration 8 norm is 8.100208244980e+01
At Newton iteration 9 norm is 8.255348762812e+01
At Newton iteration 10 norm is 9.830625477488e+01
At Newton iteration 11 norm is 1.035701793339e+02
At Newton iteration 12 norm is 4.170797681642e+01
At Newton iteration 13 norm is 3.374784836984e+01
At Newton iteration 14 norm is 3.550637878922e+01
At Newton iteration 15 norm is 3.129502005272e+01
At Newton iteration 16 norm is 1.878599545516e+01
At Newton iteration 17 norm is 1.073427713889e-10
======= elastic problem at iteration 1 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 5.115550059265e-10
======= elastic problem at iteration 2 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 5.516151581912e-10
======= elastic problem at iteration 3 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 6.235578438922e-10
======= elastic problem at iteration 4 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.499225455071e+03
At Newton iteration 3 norm is 2.880391331552e+02
At Newton iteration 4 norm is 5.116356946139e-10
======= elastic problem at iteration 5 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 7.986856982269e-10
======= elastic problem at iteration 6 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 8.355332756958e-10
======= elastic problem at iteration 7 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 9.131250719347e-10
======= elastic problem at iteration 8 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.033042535117e-09
======= elastic problem at iteration 9 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 2.127042476046e+02
At Newton iteration 3 norm is 2.297943455031e+00
At Newton iteration 4 norm is 3.068947162576e-01
At Newton iteration 5 norm is 1.046315291892e-09
======= elastic problem at iteration 10 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 3.083780453789e+03
At Newton iteration 3 norm is 2.526106928419e+02
At Newton iteration 4 norm is 1.189880401945e+02
At Newton iteration 5 norm is 1.183178684717e-09
======= elastic problem at iteration 11 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.320396084363e-09
======= elastic problem at iteration 12 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.441451878807e-09
======= elastic problem at iteration 13 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.538779091218e-09
======= elastic problem at iteration 14 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.678832810242e-09
======= elastic problem at iteration 15 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.833354552873e-09
======= elastic problem at iteration 16 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.846337922863e-09
======= elastic problem at iteration 17 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.895384037304e-09
======= elastic problem at iteration 18 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 1.874177123996e-09
======= elastic problem at iteration 19 =======
At Newton iteration 1 norm is 3.317591015983e+05
At Newton iteration 2 norm is 2.823474082624e+02
At Newton iteration 3 norm is 2.143841630050e-09
==================== Mesh ====================
elemType: TETRA4
Ne: 2880
Nn: 960
==================== Model ====================
Isotropic:
E = 2.10e+05, v = 0.3
solver:scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
W def = 692.92
Svm max = 3977.16
Evm max = 2.46 %
Ux max = 1.72e+00
Ux min = -1.70e+00
Uy max = 2.28e+00
Uy min = -2.00e+00
Uz max = 1.31e-01
Uz min = -1.82e-03
=================== TicTac ===================
Mesh: 116.684 ms
Boundary Conditions: 140.462 µs
Matrix: 1.367 s
Solver: 1.004 s
Resolution elastic: 2.478 s
PostProcessing: 688.898 µs
Generate movie 01/20 (5.00 %) 4.02 s
Generate movie 02/20 (10.00 %) 5.55 s
Generate movie 03/20 (15.00 %) 5.30 s
Generate movie 04/20 (20.00 %) 5.01 s
Generate movie 05/20 (25.00 %) 4.68 s
Generate movie 06/20 (30.00 %) 4.31 s
Generate movie 07/20 (35.00 %) 4.08 s
Generate movie 08/20 (40.00 %) 3.71 s
Generate movie 09/20 (45.00 %) 3.40 s
Generate movie 10/20 (50.00 %) 3.13 s
Generate movie 11/20 (55.00 %) 2.82 s
Generate movie 12/20 (60.00 %) 2.49 s
Generate movie 13/20 (65.00 %) 2.21 s
Generate movie 14/20 (70.00 %) 1.88 s
Generate movie 15/20 (75.00 %) 1.57 s
Generate movie 16/20 (80.00 %) 1.26 s
Generate movie 17/20 (85.00 %) 946.56 ms
Generate movie 18/20 (90.00 %) 643.66 ms
Generate movie 19/20 (95.00 %) 317.98 ms
Generate movie 20/20 (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 e = 10
33 L = 3 * e
34 t = 1 # strip thickness
35 h = 10 # strip height
36 r = 3 # fillet radius
37 thickness = e / 2 # out-of-plane extent (3D)
38 mS = t / 5 if dim == 2 else t
39
40 N = 20 # load steps
41 delta = 2 * t
42 penalty = 1e6 # contact stiffness
43
44 folder = Folder.Results_Dir()
45
46 # ----------------------------------------------
47 # Mesh
48 # ----------------------------------------------
49
50 # body
51 # ----
52 p1 = Point(-L / 2 - e)
53 p2 = Point(-L / 2, r=r)
54 p3 = Point(-e / 2, h - t, r=r)
55 p4 = Point(e / 2, h - t, r=r)
56 p5 = Point(L / 2, r=r)
57 p6 = Point(L / 2 + e)
58
59 lower = Points([p1, p2, p3, p4, p5, p6])
60 upper = lower.copy()
61 upper.Translate(dy=t)
62 contour = Points(list(lower.points) + list(upper.points[::-1]), mS)
63
64 if dim == 2:
65 mesh = contour.Mesh_2D([], ElemType.TRI3)
66 else:
67 nz = max(1, round(thickness / mS))
68 mesh = contour.Mesh_Extrude([], [0, 0, thickness], [nz], ElemType.TETRA4)
69
70 # block
71 # -----
72 domain = Domain((-L / 2 - 2 * e, -5 * t), (L / 2 + 2 * e, 0))
73
74 if dim == 2:
75 block = domain.Mesh_2D([], ElemType.QUAD4, isOrganised=True)
76 else:
77 block = domain.Mesh_Extrude(
78 [], [0, 0, thickness * 2], [1], ElemType.HEXA8, isOrganised=True
79 )
80 block.Translate(dz=-block.center[2] / 2)
81
82 nodes_top = mesh.Nodes_Conditions(lambda x, y, z: y == h)
83 nodes_contact = block.Nodes_Conditions(lambda x, y, z: y == 0)
84 block.Set_Tag(nodes_contact, "contact")
85
86 # ----------------------------------------------
87 # Simulation
88 # ----------------------------------------------
89 material = Models.Elastic.Isotropic(
90 dim, E=210000, v=0.3, planeStress=True, thickness=thickness
91 )
92 simu = RigidContact(mesh, material, penalty)
93 simu._contactMesh = block
94
95 print(f"Penalty contact solve in {dim}D (Newton per step):")
96 for i in range(N):
97 # udpate load
98 dep = [0.0] * simu.dim
99 dep[1] = -(i + 1) / N * delta
100 # solve contact
101 simu.Bc_Init()
102 simu.add_dirichlet(nodes_top, dep, simu.Get_unknowns())
103 simu.Solve()
104 simu.Save_Iter()
105
106 print(simu)
107
108 # ----------------------------------------------
109 # Results
110 # ----------------------------------------------
111 def Plot_Iter(plotter, n):
112 simu.Set_Iter(n)
113 PyVista.Plot(
114 simu, result, 1, color="k", nColors=21, show_grid=True, plotter=plotter
115 )
116 PyVista.Plot(block, color="gray", alpha=0.4, plotter=plotter)
117 PyVista.Plot_Elements(block, color="k", dimElem=1, linewidth=2, plotter=plotter)
118
119 PyVista.Movie_func(Plot_Iter, N, folder=folder, filename="contact.gif")
120
121 plotter = PyVista._Plotter()
122 result = "uy"
123 Plot_Iter(plotter, -1)
124 plotter.show()
Total running time of the script: (0 minutes 10.076 seconds)