Updated Lagrangian for Compressible Hyperelastic Material with Frictionless Contact

The Updated Lagrangian method for nonlinear elasticity with contact is presented. First, we describe the Total Lagrangian for a compressible Neo-Hookean material. Next, we introduce the Updated Lagrangian formulation for Neo-Hookean and Ogden compressible materials with contact. An advantage of this approach is that at each iteration only a linear system is solved. The linear problem to be solved is written in the updated domain. Numerical results are presented: compression of a Hertz half ball and of a hyperelastic ring against a flat rigid foundation, and contact of an elastic cube and a ball.

The incremental method (see [17], Sections 6.10-6.12) varies the forces and the imposed displacement by small increments from zero to desired values to successively solve linearized problems written in the undeformed domain. The Updated Lagrangian method is similar, but the linear problem to be solved is written in the updated domain (see [18] Sections 2.6-2.8 or [19] Section 14.8). This method was developed initially for nonlinear elasticity, but recently it was successfully employed for dynamic fluid-structure interaction [20]. An advantage of this approach is that at each time step, only a linear system is solved. A stability result is obtained in [21].
Nonlinear elasticity equations with frictionless contact can be formulated in term of a constrained nonlinear optimization problem: the nonlinear cost function is the deformation elastic energy, and the constraints are the non-penetration condition of the elastic structure into the obstacle and the imposed displacement on some boundary. The Lagrange multiplier and augmented Lagrangian methods, as well as the mortar, penalty and active set methods come from the constrained nonlinear programming algorithms, where the cost function is written in the undeformed domain.
By introducing a positive function defined on the contact zone, in Nitsche's method the problem is formulated as a system of equations which can be solved by generalized Newton's method. The semismooth Newton method can be considered related: the problem is reformulated using non-differentiable approximating equations.
The purpose of this paper is to present the Updated Lagrangian method for nonlinear elasticity with contact. The novelty is to use this method for contact problem. We can also highlight that the linearized problem written in the undeformed domain for Neo-Hookean and Ogden compressible materials are derived.
In the second section we describe the Total Lagrangian for compressible Neo-Hookean material. In Sections 3 and 4, we introduce the Updated Lagrangian formulation for Neo-Hookean and Ogden compressible materials with contact. The last section is devoted to the numerical results.

Contact without Friction in Non-Linear Elasticity Using Total Lagrangian Framework
We consider Ω 0 ⊂ R 2 an undeformed structure domain, and we assume that it is an open, bounded and connected subset. Its boundary is Lipschitz and admits the decomposition ∂Ω 0 = Γ D ∪ Γ N ∪ Γ C , where Γ D , Γ N and Γ C are relatively open subsets, mutually disjointed. On Γ D we impose a given displacement U D , in Ω 0 volume forces f are applied, and it is subjected to surface loads h on Γ N . A portion of Γ C will be in contact with a rigid foundation after deformation.
A particle of the structure whose initial position was the point X = (X 1 , X 2 ) will occupy the position If A is a square matrix, we denote by det A, tr(A), A −1 , A T the determinant, the trace, the inverse and the transpose matrix of A, respectively. We write A −T = A −1 T , and We denote by F(X) = I + ∇ X U(X) the gradient of the deformation, where I is the unity matrix, and we write C = F T F, J(X) = det F(X). The first and the second Piola-Kirchhoff stress tensors are denoted by Π and Σ, respectively, and the following equality holds: Π = FΣ. For the hyperelastic material, there exists a strain energy function W such that ∂W ∂F = Π and 2 ∂W ∂C = Σ (see [22], Chapter 6). The Cauchy stress tensor σ is computed by σ(x) = 1 J FΣF T (X), where x = X + U(X). The rigid foundation is modeled as the graph of a function ψ ∈ C 1 (R), and we denote its graph by graph(ψ) = (X 1 , X 2 ) ∈ R 2 , X 2 = ψ(X 1 ) and its epigraph by We assume that the undeformed structure domain Ω 0 is into epi(ψ). The problem to solve is

Updated Lagrangian for Compressible Neo-Hookean Material with Contact
We suppose that the material is homogeneous, isotropic and can be described by the compressible Neo-Hookean model ( [18], p. 239); the strain energy function is and the second Piola-Kirchhoff stress is where λ, µ are the Lamé constants of the linearized theory (see [23], Chapter 5).
We denote by Ω n the image of Ω 0 via the map X → X + U n (X), and we set Ω = Ω n the computational domain at the increment n. The map from Ω 0 to Ω n+1 defined by X → x = X + U n+1 (X); the composition of the map from Ω 0 to Ω is defined by X → x = X + U n (X), and the map from Ω to Ω n+1 defined by With the notations F = I + ∇ x u and J = det F, J n = det F n , we obtain For the Neo-Hoohean material, we have σ(x) = 1 J λ(ln J)I + 1 J µ FF T − I , and we set [17], Section 2.6), and taking into account (5), we get

Using (4), it follows that
If A is a square matrix, by linearization, we have: (see [23], Chapter 3.2). We can linearize the map u → F Σ by We have For simplicity, we assume that the displacement U D on Γ D , the volume forces f in Ω 0 and the surface loads h on Γ N are constant. Let N ∈ N * be the number of increments. We will solve successively N linearized problems written in the updated configuration Ω n , for 0 ≤ n ≤ N − 1.
The boundaries Γ n D , Γ n N , Γ n C of Ω n are obtained, respectively, from Γ D , Γ N , Γ C via the map X → X + U n (X). Let us introduce the linear application of the increments of the forces In the case of small displacements u, the constraint of non-penetration of the elastic structure into the rigid obstacle can be approached by see [1]. We introduce the convex set The linearized problem written in the updated configuration to be solved is the variational inequality: find u ∈ K such that Proposition 1. The bi-linear application ( u, w) → a nh ( u, w) is symmetric. Let as introduce the quadratic optimization problem with affine constraints

Proof. If
A solution of (11) is also a solution of (10). If a nh is coercive, the variational inequality (10) has a unique solution which is also the unique solution of optimization problem (11) (see [5]).
Problem (11) will be solved numerically by the Interior Point algorithm implemented in the software FreeFem++ (see [24]). The novelty of this approach is that Problem (11) is written in the updated configuration.

Updated Lagrangian for Compressible Ogden Material with Contact
We suppose that the compressible material is of type Ogden [25] with the strain energy function with I 1 = tr(C), I 2 = 1 2 (tr(C)) 2 − tr C 2 , I 3 = det(C) and c 1 , c 2 , a > 0. The first two terms correspond to the Mooney-Rivlin material, and the volumetric part of strain-energy functions used here a(I 3 − 1) − (c 1 + c 2 + a) ln I 3 was proposed in [26] in order to obtain polyconvexity and the coerciveness of the strain energy function.
We have From the Cayley-Hamilton theorem in 2D, we have and using that C is symmetric, we get As in the preceding section, using (5) and (4), we obtain We have We have a similar result as in Section 3.
As previously, the linearized problem written in the updated configuration to be solved is the variational inequality: find u ∈ K such that The associated optimization problem is

Numerical Results
Let T h be a triangulation of Ω of size h, with nv vertices. We set φ i : T h → R the shape function associated with vertex A i , which is a piecewise linear function and is globally continuous. For the two-dimension displacements, we introduce the basis 0), for i = 1, . . . , nv and φ nv+i = (0, φ i ), for i = 1, . . . , nv.
We define the matrix A ∈ R 2 nv×2 nv and the vector b ∈ R 2 nv by The constraint w ∈ K will be treated weakly. Thus, we introduce matrix C ∈ R ngC×2 nv , where ngC is the number of vertex A i ∈ Γ n C and the vector d ∈ R ngC by The discrete version of (11) is For the numerical tests, we employed the software FreeFem++ (see [24]). The optimization problem (15)-(17) is solved by the library IPOPT "Interior Point OPTimizer", which has an interface in FreeFem++.

Compression of a Hertz Half Ball against a Foundation
This example is adapted from [11]. The undeformed structure domain Ω 0 is a half ball of radius R = 8 m with center (0, R), and the rigid foundation is given by ψ(X 1 ) = 0. The boundary Γ D is the little segment [AB] = {t ∈ (−0.095, 0.095); x(t) = t, y(t) = R}, Γ C is the half of a circle {t ∈ (π, 2π); x(t) = R cos(t), y(t) = R + R sin(t)}, and Γ N is the rest of ∂Ω 0 .
The quadratic optimization problem with affine constraints is The analytical normal stress in the contact zone given by the Hertz theory is The contact zone is |x 1 | < b and x 2 = 0. The number of nodes on Γ C is 252; the mesh of Ω 0 has 11,759 vertices and 23,096 triangles. For h 2 = −2 and x 1 = 0, the analytical value for b is 1.40621 and for the normal stress is −14.4871, while the computed normal stress is −14.545. In Figure 1, we can see: the initial mesh, the von Mises stress and a zoom of the contact zone. The numerical solution is consistent with the analytical solution. The IPOPT algorithm solves the optimization problem after 10 iterations.

Compression of a Hyperelastic Ring against a Flat Rigid Foundation
This example is adapted from [11]. The undeformed structure domain Ω 0 is a ring of exterior radius R e = 10 m, interior radius R i = 9 m and center (0, R e ), and the rigid foundation is given by ψ(X 1 ) = 0. The boundary Γ D is the arc of a circle AB = {t ∈ ( π 2 − π 48 , π 2 + π 48 ); x(t) = R e cos(t), y(t) = R e + R e sin(t)}, Γ C is the inferior half of a circle {t ∈ (π, 2π); x(t) = R e cos(t), y(t) = R e + R e sin(t)}, and Γ N is the rest of ∂Ω 0 .
The number of nodes on Γ C is 240; the mesh of Ω 0 has 4342 vertices and 7734 triangles. We employ finite element P 1 , and we set N = 14 as the number of increments. The average number of iterations of the IPOPT algorithm in order to solve, at each increment, the optimization problem is 12.
We can see in Figure 2 the final mesh and the von Mises stress for Neo-Hookean-like material, and in Figure 3 the initial, intermediate and final meshes for the Ogden material. We denote by U nh the solution in the case of the Neo-Hookean material and by U og for the Ogden-like material. We have U nh − U og L 2 (Ω 0 ) = 1.88063. Our solution for the Ogden-like material is similar to the one obtained in [11]. (c) (d) Figure 3. Ring, Ogden model: mesh deformation after: 0 (a), 5 (b), 10 (c) and 14 (d) increments.

Contact of an Elastic Cube and a Ball
This example is adapted from [13]. The undeformed structure domain Ω 0 is the cube (0, 1) 3 , and the obstacle is a ball with center (0.5, 0.5, −0.3) and radius 0.3, see Figure 4.
For 3D, we have the same formula for (7), but (9) is replaced by The mesh is controled by the number k of segments on each edge of the cube. The mesh of Ω 0 has: 9261 vertices and 48,000 tetrahedrons for k = 20, 35,937 vertices and 197,608 tetrahedrons for k = 32 and 68,921 vertices and 384,000 tetrahedrons for k = 40. We employ the finite element P 1 , and we set N = 4 as the number of increments.  The average number of iterations of the IPOPT algorithm is 14 for k = 20 and 16.75 for k = 40. The total CPU time on an Intel Sandy Bridge 16 × 3.30 GHz and 64 GB RAM was 6 min for k = 20, 24 min for k = 32 and 50 min for k = 40.

Conclusions
An Updated Lagrangian method for nonlinear elasticity with frictionless contact was presented. The linearized problem written in the updated configuration for Neo-Hookean and Ogden compressible materials were derived. At each iteration, only a linear system was solved. Two-and three-dimensional numerical simulations were performed.

Conclusions
An Updated Lagrangian method for nonlinear elasticity with frictionless contact was presented. The linearized problem written in the updated configuration for Neo-Hookean and Ogden compressible materials were derived. At each iteration, only a linear system was solved. Two-and three-dimensional numerical simulations were performed.