Mathematical analysis of successive linear approximation for Mooney-Rivlin material model in finite elasticity

For calculating large deformations in finite elasticity, we have proposed a method of successive linear approximation, by considering the relative descriptional formulation. In this article we briefly describe this method and we prove the existence and uniqueness of weak solutions for boundary value problems for nearly incompressible Mooney-Rivlin materials, that arise in each step of the method.


Introduction
The constitutive equation of a solid body is usually expressed relative to a preferred reference configuration which exhibits specific material symmetries such as isotropy. The constitutive functions are in general nonlinear and linearizations can be used as valid approximation only for small deformations. Therefore, the problem for large deformations leads to boundary value problems involving systems of nonlinear partial differential equations.
In order to circumvent the difficulties due to the nonlinearities, we have proposed a new method for solving numerically the boundary value problem for large deformations. It is based on a successive linear approximation by considering the relative descriptional formulation. Roughly speaking, the constitutive equations are calculated at each state which will be regarded as the reference configuration for the next state, and assuming that the deformation to the next state is small, the updated constitutive equations can be linearized.
As examples for the proposed method, numerical simulations were done (see [1], [2]) for two classical problems concerning Mooney-Rivlin materials, for which the exact solutions are known, namely, the pure shear of a square and the bending of a retangular block into a circular section. The comparison of the numerical results with the exact solutions of these two examples confirms the efficiency of our method.
In the present paper we consider the mathematical analysis of the boundary value problem obtained by linearizing the constitutive equations of nearly incompressible Mooney-Rivlin materials relative to the present configuration and prove the existence and uniqueness of weak solutions.
We organize this paper as follows. In Section 2 we introduce briefly the notion of relative description and we formally deduce in Section 3 the linearization of the constitutive function of a nearly incompressible Mooney-Rivlin material. In Section 4 we consider a boundary value problem involving a system of partial differential equations, that are obtained by linearizing the constitutive equations of a nearly incompressible Mooney-Rivlin material, and which corresponds to one of the steps of the successive linear approximation method. The main result of this paper is contained in section 5, where we prove the existence and uniqueness of weak solutions of this boundary value problem, by considering its variational formulation. For simplicity, we restrict our analysis to the two-dimensional case, but the arguments presented can be extended to three dimensions.

Relative description and successive linear approximation
In this section we introduce the notion of relative description, we formally obtain the linearization of a general constitutive equation of a solid body respect to this configuration and we describe the successive linear approximation method.
Let κ 0 be a reference configuration of a solid body B, B 0 = κ 0 (B), and let be the parametrization of its deformation. Let κ t be the deformed configuration at time t (which we shall always refer as the present time), B t = κ t (B), and be the deformation gradient with respect to the configuration κ 0 .
Let κ τ be the deformed configuration at time τ > t. We define the relative deformation from κ t to κ τ as the function χ t : B t → B τ given by and the corresponding relative displacement as Taking the gradient relative to x in both sides of (2.2), we obtain where I is the identity tensor and are called the displacement gradient and the deformation gradient in the relative description, relative to the present configuration.
On the other hand, taking the gradient relative to X in both sides of (2.2), we obtain from (2.1) and the chain rule, from which we get easily F (X, τ ) = I + H t (x , τ ) F (X, t).
We can represent this situation by the following picture: where and hereafter, for simplicity, sometimes only the time dependence is indicated. Position dependence is usually self-evident and will be indicated for clarity only if necessary.
By considering the time τ = t + ∆t for small enough ∆t, we can assume that the relative displacement gradient is small, Let T be the Cauchy stress tensor given by the constitutive equation Assuming that the operator F κ0 is differentiable, we can calculate the linearization of the constitutive equation (2.5) relative to the current configuration κ t , and assuming that ∆t is small enough, we have formally where dF κ0 (F ) denotes the Fréchet-differential of F κ0 calculated at F . For convenience, we shall write (2.6) as where L(F )[H] := dF κ0 (F )[HF ] defines a fourth order elasticity tensor L(F ) relative to the current configuration κ t .
The successive linear approximation method is the discrete construction of the parametrization χ (X, t) based on the previous arguments. More precisely, let t 0 < · · · < t n−1 < t n < t n+1 < · · · be a sequence of steps with small enough constant spacing ∆t, where at which step we set t = t n and τ = t n+1 . Let the deformation gradient F (X, t n ) and the elastic Cauchy stress tensor T (X, t n ), relative to the preferred configuration κ t0 , assumed to be known. If in any way we calculate the relative displacement u tn (x , t n+1 ), x ∈ κ tn (B), it allows us the update the new reference configuration κ tn+1 relative to the next step by using (2.1) and (2.2), i.e., while the deformation gradient (2.4) and the Cauchy stress (2.7), relative to the preferred configuration κ t0 , can be determined at instant t n+1 respectively by Therefore, after updating the boundary data and the eventual body forces acting on the body, we repeat the cycle from the updated reference configuration κ tn+1 .
We remark that this method can easily be extended to constitutive equation T = F (F,Ḟ ) for viscoelastic solid bodies in general [3].

Application to nearly incompressible Mooney-Rivlin materials
From now on we consider a Mooney-Rivlin material whose constitutive equation relative to the preferred reference configuration κ 0 is given by where B = F F T is the left Cauchy-Green strain tensor and the material parameters s 1 and s 2 are constant satisfying s 1 > 0 and s 2 < s 1 . (3.1) Remark 3.1: It is usually assumed s 2 ≤ 0 < s 1 , the so-called E-inequalities (see [4] and [5]), based on the assumption that the free energy function is positive definite for any deformation. Liu [6] has pointed out that "any" deformation is unrealistic from physical point of view, and a thermodynamical stability anaysis only requires s 2 < s 1 . Therefore, we shall include the case 0 < s 2 < s 1 in our analysis.
A direct calculation of the Fréchet-differential of F at F gives For compressible body in general, the pressure p may depend on the deformation gradient F . However, for compressible elastic bodies, we shall assume that the pressure depends only on the determinant of the deformation gradient or, by the mass balance, depends only on the mass density ρ, where ρ 0 denotes the mass density in the referential configuration κ 0 .
For time τ = t + ∆t and from (2.4), we have where tr H means the trace of H and o(2) denotes higher order terms in the small displacement gradient H(τ ). Therefore, assuming that p is differentiable as function of ρ, we have where β(t) = ρ(t)(dp/dρ) t is a material parameter depending on the mass density ρ.
A body is called nearly incompressible if its density is nearly insensitive to change of pressure. Hence, if we regard the density as a function of pressure, ρ = ρ(p), then its derivative with respect to the pressure is nearly zero. This means that, for nearly incompressible materials, the parameter β must be large, Therefore, the Cauchy stress tensor relative to the current configuration κ t is given by and the first Piola-Kirchhoff stress tensor at time τ relative to the current configuration κ t is given by

Linearized boundary value problem and its variational formulation
For simplicity, we denote by κ the current configuration κ t , Ω = κ(B) be the bounded domain of R 3 representing the interior of the region occupied by the body at current configuration κ at the present time t, T 0 = T (t) and B 0 = B(t). Let ∂Ω = Γ 1 ∪ Γ 2 ∪ Γ 3 , n κ be the exterior unit normal to ∂Ω and g be the gravitational force (per unit mass).
We consider the following boundary value problem for the relative displacement u = u(x , τ ), where div is the divergence operator with respect to x , T κ (τ ) = T κ (x , τ ) is the Piolla-Kirchhoff stress tensor at time τ relative to configuration κ at the present time t, which, up to linear terms in relative displacement gradient H = H(τ ) = ∇ x u(τ ), is given by (see (3.2)) and f (τ ) is the surface traction (per unit surface area).
Remark 4.1: At every time step, the idea of formulating the boundary value problem in the form (4.1) is similar to the theory of small deformations superposed on finite deformations (see [7], [5]). In this manner, either we are interested in the evolution of solutions with gradually changing boundary conditions resulting in large deformation, or, we can treat the boundary values of finite elasticity as the final value of a successive small incremental boundary values at each time step (see [8]).
The boundary value problem (4.1) can be formulated as a variational problem. Indeed, let Ω be a smooth enough bounded domain in R 3 and define the space Taking formally the inner product of both sides of the equation in (4.1) by w ∈ V and integrating over Ω, we obtain after integration by parts, Therefore, for u , w ∈ V we consider respectively the bilinear and the linear forms: We notice that the forms L and N can be written in terms of coordinates by where in the above formulas we have used the standard summation convention for repeated indices.
Then, the variational problem is to find the solution u ∈ V such that In order to prove that the solutions of (4.3) is a weak solution of (4.1), the following result concerning existence of a normal trace is useful.
It is well-known that, for a given ϕ ∈ H 1/2 (∂Ω) we may choose ψ ∈ H 1 (Ω) such that γ 0 (ψ) = ϕ and such that ψ H 1 ≤ C 1 ϕ H 1/2 , where the constant C 1 depends only on Ω. Hence, This means that When F is no longer in C 1 (Ω) ∩ C 0 (Ω), using a density argument (see [9]), we can find a sequence {F n } n∈N in C 1 (Ω) ∩ C 0 (Ω) such that Inequality (4.4) shows that {F n · n κ } n∈N is a Cauchy sequence in H −1/2 (∂Ω), whose limit, which is independent of the particular choice of the sequence {F n } n∈N , will be denoted by F · n κ . This finishes the proof.

Existence and uniqueness of solution in two-dimensions
Let Ω be a bounded Lipschitz domain of R 2 whose boundary ∂Ω = Γ 1 ∪ Γ 2 ∪ Γ 3 , with meas(Γ i ) = 0 for i = 1, 2, 3, and consider the space For u, w ∈ V, we introduce where L 2 is the usual L 2 -norm. It is well-known that the Poincaré inequality holds if meas(Γ 3 ) = 0, i.e., there exists a constante C such that In this case, ·; · and V define an inner product and a norm in V, respectively. From now on we assume that where by S + 2 (R) we denote the set of all symmetric and positive definite 2 × 2 matrix, and we set It is clear that the forms L and N defined in (4.2) are continuous in V.
Recalling that H and W are 2 × 2 matrix whose entries are given by the bilinear form L(u , w ) defined in (4.2) can be written as In particular, for W = H we have  Hence, to prove that L is coercive, it is sufficient to show that there exists α > 0 such that i.e., it suffices to show that the bilinear form A(x ; H, W ) is uniformly coercive as function of 2 × 2 matrices. Furthermore, a direct calculation (see the Appendix) gives that the coercivity of A(x , H, W ) is equivalent to the semipositivity of the matrix A(x ) − αI, for all x ∈ Ω, with A(x ) given by where γ 1 e γ 2 are the eigenvalues of B 0 and I is the 4 × 4 identity matrix. Therefore, we can also write where, following the notation introduced in the Appendix, We notice that if A(x ) − αI is uniformly semipositive in Ω, then and consequently, In order to analyze the matrix (5.5) and in view of the conditions (3.1), we must distinguish two cases: s 2 < 0 < s 1 and 0 ≤ s 2 < s 1 . In both cases, we fix a constant k > max{0, s 2 s −1 1 } and take ε := s 1 − s 2 k −1 . Now, let a 0 = a 0 (x ) and b 0 = b 0 (x ) be the functions defined by Assuming that det B 0 ≥ k, we have the following inequalities Notice that the above inequalities indicate that the interval [a 0 (x ), b 0 (x )] has nonempty interior for all x ∈ Ω if det B 0 ≥ k, and this will be essential in proving the next theorem.
A necessary and sufficient condition for the matrix A − αI be positive semidefinite is It is clear that the condition (5.11) implies, in particular, A 22 − α ≥ 0 and A 44 − α ≥ 0, since A is symmetric.
Step 1: Analysis of the first block of A: In the case s 2 < 0, we have f (γ) ≥ −s 1 s 2 for all γ > 0. So, In the case s 2 ≥ 0, we can assume without loss of generality that γ 1 ≥ γ 2 . Then, as and Therefore, in the two cases, On the other hand, if we denote f i = f (γ i ) and g i = g(γ i ), i = 1, 2, we get . Since we have Step 2: Analysis of the second block of A: By the definition of T 0 , we have Hence, A 33 − α ≥ 0 if, and only if, 14) The conditions (5.13) and (5.14) can be expressed by It is noteworthy that (5.12) implies (5.13). Moreover, the interval is not empty, because if we denote and hence, from (5.10): which implies that the interval defined by (5.15) is not empty. Now, according to the notation introduced above, we have where we are denoting Hence, The above inequality holds if the discriminant of the binomial (5.16) is positive. In fact, we have Note that To simplify the notation, consider Then, Note also that, from (5.8) we have: which implies that, in both cases, 0 < αD/2C < 1. So, by calculating the roots a α and b α of the binomial (5.16), we get, and the condition ( Therefore, from (5.8) and (5.6), we have b α ≥ b 0 − αd and a α ≤ a 0 + αd.
Therefore, in the case s 2 < 0, we have from which we conclude that On the other hand, it is easy to see that which implies that, and so, a * ≤ a 0 .
In the case s 2 ≥ 0, we have Since tr B 0 ≥ 2 √ det B 0 and d ≥ 2, it follows that b * − 2α ≥ b 0 − αd holds for all α > 0. Likewise, and a * + α ≤ a 0 + αd holds for all α > 0. This finishes the proof. Remark 5.5: Theorem 5.1 gives a sufficient condition for the existence of a unique weak solution of the boundary value problem (4.1) corresponding to each time step of the successive approximation. As we are supposing that the material is nearly incompressible, it is reasonable to expect that det B 0 ≈ 1. This implies that, if γ 1 and γ 2 are the eigenvalue of B 0 , γ 1 ≈ 1/γ 2 and tr B −1 0 ≈ γ 1 + 1/γ 1 . So, the hypothesis (5.9) does not means that we are assuming that p 0 is small. On the other hand, numerical experiments show that the hypothesis (5.9) can be very restrictive in the presence of gravitational body forces. In this case, we can incorporate the potential of the gravitational force into the pressure, and analyze the re-formulated problem.
Remark 5.6: The previous results hold if we assume that Γ 3 = ∅. In fact, unlike the space V introduced in (5.1), we must consider V = u ∈ (H 1 (Ω)) 2 ; u · n κ = 0 on Γ 2 . (5.18) However, in this case, it is necessary to assume that the domain Ω satisfies a geometric property to ensure that (5.2) is a norm. This can be done by supposing that Ω has the following property: There is no constant vector c ∈ R 2 such that c · n κ (x ) = 0, ∀x ∈ Γ 2 .

Appendix
Without loss of generality, we can assume that B 0 is a diagonal matrix, given by