F Peyraut

Z.-Q Feng
email: feng@iup.univ-evry.fr

Q.-C He

N Labed

Robust numerical analysis of homogeneous and non-homogeneous deformations

Keywords: Hyperelasticity, Ogden's model, Non-linear analysis, Finite element method

The present paper is devoted to the modeling of finite deformations of hyperelastic bodies described by Ogden's model. For the finite element implementation, the explicit expressions of the tangent operator and the stress tensor are derived. If the eigen-values of the right Cauchy-Green strain tensor are distinct, a simple and compact tensor formula, well adapted for numerical implementation, is proposed. A limiting technique is used to take into account the special case of coalescent eigenvalues in which non-differentiability occurs. Three test examples, including homogeneous and non-homogeneous deformations, are proposed to illustrate the developed formulation.

Introduction

Many industrial applications are concerned by foam-like or rubber-like materials. The tire technology is for example one of the main application fields of rubber-like shells. Such materials are also used to improve the performance of safety glass by using interlayers to join the splinters in the case of a crash [START_REF] Bois | Modelling of safety glass for crash simulation[END_REF]. Another application for the automotive industry is the use of polyurethane foams for a better comfort of car seats [START_REF] Prudent | Simulation numérique du confort de sièges d'automobiles : comportement mécanique 3D de mousses de polyuréthanne ; optimisation 2D d'un profil simplifié de siège[END_REF]. But all these materials are highly non-linear and require specific and robust techniques to perform numerical computation.

In the last decades, many attempts have been made to solve non-linear problems, including hyperelastic materials, with the finite element method. A very detailed review on finite element formulation for non-linear analysis, including two and three dimensional problems and involving isotropic, orthotropic, rubber-like and elasto-plastic materials, has been provided by Sussman and Bathe [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF]. In their paper, one appendix is entirely devoted to the stress and tangent constitutive tensors for the compressible Ogden's material model but only by considering the particular case of two dimensional plane strain or axisymmetric analysis. Simo and Taylor [START_REF] Simo | Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms[END_REF] have provided a closed form expression in the most general case of the spatial and material tangent moduli and applied it to the Ogden's material. Saleeb, Chang and Arnold [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF] have also given a closed form expression but without numerical examples. Lastly, Basar and Itskov [START_REF] Basar | Finite element formulation of the Ogden material model with application to rubber-like shells[END_REF] have obtained a closed form expression for the Ogden's material tangent moduli by using the Cauchy-Green strain tensor invariants to perform derivation and have applied the results for the analysis of rubber-like shells. But all these papers have not intensely studied the algorithmic issue to increase computational efficiency while we have focused on this topic in this work by developing approaches well adapted for numerical implementation.

In isotropic non-linear elasticity, there exist many models to describe the hyperelastic behavior of foam-like or rubber-like materials [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF][START_REF] Ogden | Non-linear Elastic Deformations[END_REF][START_REF] Ciarlet | Three-dimensional Elasticity[END_REF][START_REF] Arruda | A three dimensional constitutive model for the large deformation stretch behavior of rubber elastic materials[END_REF][START_REF] Gent | A new constitutive relation for rubber[END_REF][START_REF] Simo | Computational inelasticity[END_REF][START_REF] Horgan | Finite thermoelasticity with limiting chain extensibility[END_REF]. The model proposed by Ogden [START_REF] Ogden | Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids[END_REF] is one of the most important ones and is expressed in terms of principal stretches of the Cauchy-Green strain tensor. As outlined by Valanis and Landel [START_REF] Valanis | The strain-energy function of a hyperelastic material in terms of the extension ratios[END_REF] this form is more suitable to perform experimental characterization than a strain energy density expressed in terms of the strain invariants. But, for this model, the finite element implementation of the elastic tangent operator usually requires lengthy tedious development because specific computational procedures are needed to deal with the special case of coalescent eigenvalues. For this special case, the perturbation technique [START_REF] Miehe | Computation of isotropic tensor functions[END_REF][START_REF] Miehe | Aspects of the formulation and finite element implementation of large strain isotropic elasticity[END_REF] and the limiting procedure [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF][START_REF] Basar | Finite element formulation of the Ogden material model with application to rubber-like shells[END_REF][START_REF] Duffett | The analysis of incompressible hyperelastic bodies by the finite element method[END_REF][START_REF] Itskov | On the theory of fourth-order tensors and their applications in computational mechanics[END_REF] are two approaches mainly used. The perturbation technique consists in making two identical eigenvalues slightly different by adding a small numerical perturbation to one of them. Formulas available for distinct eigenvalues can then be readily applied. But the error can be of the same order or greater than the perturbation, especially with three dimensional problems. On the other hand, the limiting procedure consists in studying theoretically the limit when the eigenvalues approach one to the others. This procedure provides closed-form formulas for the stress tensor and the tangent operator and can also be applied in a more general context than the problem studied here: for example non-symmetric tensor involved in large strain viscoplasticity, anisotropic and particularly single crystal plasticity [START_REF] Itskov | Computation of the exponential and other isotropic tensor functions and their derivatives[END_REF]. In our work, the limiting procedure is selected because it offers a more stable numerical approach than the perturbation technique.

To calculate the constitutive tangent tensor, one way consists in referring to the principal axes of this tensor. Chadwick and Ogden have drawn attention to the advantages offered by this choice in the field of acceleration wave's theory [START_REF] Chadwick | A theorem of tensor calculus and its application to isotropic elasticity[END_REF][START_REF] Chadwick | On the definition of elastic moduli[END_REF]. Duffett and Reddy have made extensive use of this technique but only by considering plane problems [START_REF] Duffett | The analysis of incompressible hyperelastic bodies by the finite element method[END_REF]. Peng and Chang have also used this technique by focusing on total and updated Lagrangian finite element formulations [START_REF] Peng | A compressible approach in finite element analysis of rubber-elastic materials[END_REF]. But this technique requires stress and material moduli calculation in principal directions first and then transformations in an active coordinate system. Consequently, eigenvectors calculation must be performed previously. To avoid it, a spectral decomposition expressed in terms of eigenvalue-bases instead of eigenvectors can be used [START_REF] Morman | The generalized strain measure with application to nonhomogeneous deformations in rubber-like solids[END_REF][START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF]. This kind of spectral decomposition is very efficient because eigenvalue-bases are directly deduced from the right Cauchy-Green strain tensor by simple algebraic formulas [START_REF] Simo | Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms[END_REF]. We have thus selected it in this work.

To improve numerical results and the computational efficiency, the aim of the present paper is to develop formulae as simple as possible for calculation of the stress and stiffness tangent tensors associated with Ogden's materials. To this end, we use tensor calculus and keep as much as possible eigenvalue-bases expressions. The tensor formula obtained offers two advantages. Firstly, its compact and simple form is suitable for a numerical implementation and leads to less computational times. Secondly, it provides a very convenient way to perform the limiting procedure while theoretical calculations related to this procedure are often considered as lengthy and tedious in the literature. Moreover, we have shown that our formulas are equivalent to those provided by [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF] in the case of coalescent eigenvalues.

Three numerical test examples are performed in this study to show the validity of the model developed. For homogeneous deformations, reliable results are obtained in the case of double or triple coalescence as well as for distinct eigenvalues. Moreover, an excellent agreement with an analytical solution [START_REF] Ogden | Non-linear Elastic Deformations[END_REF] has been found. In the case of a pure torsion problem with a circular cylinder [START_REF] Polignone | Pure torsion of compressible nonlinearly elastic circular cylinders[END_REF], giving non-homogeneous deformations with two identical eigenvalues, the numerical results were successfully compared to the analytical solution. The third example is solved to show non-homogeneous deformations with variable eigenvalues.

Notation.

In order to present the essential structure of subsequent formulae in a clear manner, coordinate-free notation is employed as much as possible. A bold-face Latin lowercase letter, say a, and a bold-face Latin capital letter, say A, will denote a vector and second-order tensor, respectively. An outline Latin capital letter, say A, refers to a fourthorder tensor. If their Cartesian components are noted as a i , A ij and A ij kl , the main coordinate-free symbols used in this paper are related to the corresponding index symbols in the following way:

(a ⊗ b) ij = a i b j ,( A ⊗ B) ij kl = A ij B kl , (A ⊗ B) ij kl = A ik B jl ,( A ⊗ B) ij kl = A il B jk .
For conciseness, the convention of summation over repeated indices is adopted and it will also be convenient to define

A ⊗ B = 1 2 (A ⊗ B + A ⊗ B).

Hyperelastic bodies undergoing large deformations

Rubber-like materials are usually taken to be hyperelastic and often undergo large deformations. To describe the geometrical transformations in R 3 , the deformation gradient tensor is introduced by

F = I +∇u = I + ∂u ∂X , ( 1
)
where I is the identity tensor, X ∈ R 3 the position vector and u ∈ R 3 the displacement vector. Due to large deformations and large rotations, Green-Lagrange strain is adopted for the non-linear relationships between strains and displacements. We denote by C the stretch tensor or the right Cauchy-Green strain tensor: C = F T F. The Green-Lagrange strain tensor E is defined as

E = 1 2 (C -I). (2)
For an hyperelastic law, there exists an elastic potential function W (or strain energy density function) which is a scalar function of one of the strain tensors and whose derivative with respect to one strain component determines the corresponding stress component. This can be expressed by

S = ∂W ∂E = 2 ∂W ∂C , (3)
where S is the second Piola-Kirchoff stress tensor. In particular, for isotropic hyperelasticity [START_REF] Ciarlet | Three-dimensional Elasticity[END_REF], Eq. ( 3) can be written as

S = 2 I 3 ∂W ∂I 3 C -1 + ∂W ∂I 1 + I 1 ∂W ∂I 2 I - ∂W ∂I 2 C , (4)
where I i (i = 1, 2, 3) denote the invariants of the right Cauchy-Green deformation tensor C:

I 1 = tr(C), I 2 = I 2 1 -tr C 2 /2,I 3 = det(C). ( 5
)
or equivalently

I 1 = λ 1 + λ 2 + λ 3 ,I 2 = λ 1 λ 2 + λ 2 λ 3 + λ 1 λ 3 ,I 3 = λ 1 λ 2 λ 3 , (6)
where λ a (a = 1, 2, 3) are the eigenvalues of the right Cauchy-Green tensor C. These eigenvalues are solution of the characteristic polynomial:

λ 3 -I 1 λ 2 + I 2 λ -I 3 = 0, (7)
which can be solved in closed form thanks to the Cardano's formula (see for example Eqs. ( 9), [START_REF] Gent | A new constitutive relation for rubber[END_REF] in [START_REF] Basar | Finite element formulation of the Ogden material model with application to rubber-like shells[END_REF]):

λ a = 1 3 I 1 + 2 I 2 1 -3I 2 cos β + a2π 3 ,a = 1, 2, 3, (8)
with

β = cos -1 2I 3 1 -9I 1 I 2 + 27I 3 2(I 2 1 -3I 2 ) 3/2 . ( 9
)
The Ogden constitutive law is used to model compressible or incompressible rubber-like materials [START_REF] Ogden | Non-linear Elastic Deformations[END_REF]. This model is based on the principal eigenvalues of C. The strain energy density function is given as follows

W = N i=1 μ i α i λ α i /2 1 + λ α i /2 2 + λ α i /2 3 -3 + N i=1 μ i α i β i I -α i β i /2 3 -1 , (10)
where N , μ i , α i and β i are material parameters. The initial shear modulus, G, and the initial bulk modulus, K,are given as

G = 1 2 N i=1 μ i α i ,K = N i=1 μ i α i 1 3 + β i . (11)
The first and second derivatives of W with respect to C are needed for the stresses and tangent modulus calculation (see next sections). It is noted from Eq. ( 10) that the strain energy is split in two parts. The first one depends on eigenvalues of C while the second one depends on the third invariant I 3 . The derivative of the second part can be calculated in a conventional way by the chain rule while the derivative of the first part must require special attention for coalescent eigenvalues.

Calculation of the stress tensor

From Eqs. ( 3), [START_REF] Gent | A new constitutive relation for rubber[END_REF], the stress tensor can be split into two parts:

S = S 1 + S 2 . ( 12
)
The first part S 1 is related to the principal stretches of C and the second part S 2 is related to the third invariant of C. The calculation of S 2 is relatively simple. In fact, by using Eqs. ( 4) and ( 10), S 2 can be easily expressed as

S 2 =- N i=1 μ i I 3 -α i β i /2 C -1 . (13)
However, it is much more delicate to calculate S 1 . Particular attention must be paid when dealing with distinct or coalescent eigenvalues of C.

Distinct eigenvalues

In view of Eqs. ( 3), [START_REF] Gent | A new constitutive relation for rubber[END_REF], S 1 is given by:

S 1 = 3 a=1 S(λ a ) ∂λ a ∂C (14)
where S(λ a ) is defined by:

S(λ a ) = N i=1 μ i λ a α i /2-1 . (15)
We have adopted here the notation S(λ a ) to compare our formulas to those established by Saleeb, Chang and Arnold (see Eq. ( 25) in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF]). The computation of λ a in Eq. ( 14) is performed by using Eqs. ( 8), ( 9) while ∂λ a ∂C is calculated by using the spectral decomposition of C. This decomposition is written as follows (see for example Eq. ( 4) in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF]):

C = 3 a=1 λ a M a ( 16
)
where the eigenvalue-bases M a (a = 1, 2, 3) are related to the eigenvectors n a of C by:

M a = n a ⊗ n a . (17)
A classical fundamental result allows to calculate M a by (see for example Eq. ( 7) in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF] and Eq. (19a) in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF]):

M a = ∂λ a ∂C = C 2 -(I 1 -λ a )C + I 3 λ -1 a I (λ a -λ b )(λ a -λ c ) (18)
where (a,b,c) represents a cyclic permutation of (1, 2, 3). By reporting Eq. ( 18) in Eq. ( 14), the part of the stress tensor related to the principal stretches is thus simplified as follows:

S 1 = 3 a=1 S(λ a )M a . (19)
It is noted that the calculation of eigenvectors n a , used in Eq. ( 17) to define M a , is unnecessary because the later can be directly obtained from Eq. ( 18). To perform the stress evaluation, the best way is to use Eq. ( 18) first and Eq. ( 19) next. But, in order to make comparisons with known results, it is interesting to report Eq. ( 18) in Eq. ( 19) and factorize the terms related to C 2 , C and I:

S 1 = pC 2 + qC + rI (20)
where the parameters p, q and r are defined by

p =-D 3 a=1 (λ a -λ b )S(λ c ), (21)
q = D 3 a=1 λ 2 a -λ 2 b S(λ c ), (22)
r = D 3 a=1 (λ a -λ b )λ a λ b S(λ c ) (23)
with

D = 1 (λ 1 -λ 2 )(λ 2 -λ 3 )(λ 3 -λ 1 ) . (24)
Eqs. ( 20)-( 24) are consistent with Eqs. (26a), [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF] provided in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF]. However, it should be noticed that Eq. ( 18), and consequently Eqs. ( 20)- [START_REF] Schröder | Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions[END_REF], are valid only for distinct eigenvalues. In case of coalescent eigenvalues, a special procedure must be applied.

Double eigenvalues

In this case (say λ 1 = λ 2 = λ 3 ), we recall first a classical matrix relation available for any case (see for example Eq. ( 5) in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF])

M 1 + M 2 + M 3 = I (25)
and a second relation which only holds in case of double eigenvalues (see for example Eq. (19b) in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF]):

M 3 = C -λ 1 I λ 3 -λ 1 . ( 26
)
By restarting from Eq. ( 20) and using Eq. ( 25) it can be easily shown that

S 1 = S(λ 1 )I + S(λ 3 ) -S(λ 1 ) M 3 . (27)
In view of Eqs. ( 26), [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF], we can factorize the terms related to C and I:

S 1 = sC + tI (28)
where the parameters s and t are defined by

s = S(λ 3 ) -S(λ 1 ) λ 3 -λ 1 ,t = λ 3 S(λ 1 ) -λ 1 S(λ 3 ) λ 3 -λ 1 . (29)
Eqs. ( 28), (29) are consistent with Eqs. (26b), [START_REF] Valanis | The strain-energy function of a hyperelastic material in terms of the extension ratios[END_REF] provided by Saleeb, Chang and Arnold in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF].

Triple eigenvalues

Let us consider now the case of triple eigenvalues (λ 1 = λ 2 = λ 3 ). In this case, C is obviously proportional to I:

C = λ 1 I. (30)
It results in

n 1 = 1 0 0 , n 2 = 0 1 0 , n 3 = 0 0 1 . ( 31
)
We then deduce from Eqs. ( 17), (31) that:

M 1 = 100 000 000
, M 2 = 000 010 000

, M 3 = 000 000 001

. ( 32
)
By reporting Eq. (32) in Eq. ( 19), we obtain the final result for this section:

S 1 = S(λ 1 )I, (33)
which is again consistent with the result provided by Saleeb, Chang and Arnold (refer to Eq. (26c) in [START_REF] Saleeb | On the development of explicit robust schemes for implementation of a class of hyperelastic models in large-strain analysis of rubbers[END_REF]).

Incremental Ogden's law

In order to construct the tangent stiffness matrix for the analysis of non-linear structures, it is necessary to determine the stress-strain tangent operator D which appears in the incremental constitutive law:

dS = D dE. (34)
This fourth-order tensor D is obtained from the derivative of S with respect to E in Eq. ( 3) and from the Ogden's strain energy density equation ( 10):

D = ∂S ∂E = 4 ∂ 2 W ∂C 2 = 4 N i=1 μ i α i ∂ 2 W i 1 ∂C 2 + N i=1 μ i α i β i ∂ 2 W i 2 ∂C 2 , ( 35
)
where the strain energy densities W i 1 and W i 2 are defined by:

W i 1 = λ α i /2 1 + λ α i /2 2 + λ α i /2 3 ,W i 2 = I 3 -α i β i /2 . ( 36
)
For the sake of simplicity, we note α = α i /2 and β =-α i β i /2. The superscript and subscript i then vanish and Eq. (36) becomes

W 1 = λ α 1 + λ α 2 + λ α 3 ,W 2 = I 3 β . (37) 4.1. Computation of ∂ 2 W 1 ∂C 2
By deriving twice Eq. (37) 1 and using Eq. ( 18), we obtain

∂ 2 W 1 ∂C 2 = α (α -1) 3 a=1 λ α-2 a M a ⊗ M a + 3 a=1 λ α-1 a ∂M a ∂C . ( 38
)
The derivatives of the eigenvalue-bases M a with respect to the strain tensor C involved in the above equation are computed in the following subsections. Because the case of coalescent eigenvalues needs specific computation, these subsections will be divided into three parts, respectively dealing with distinct, double and triple eigenvalues. If C admits distinct eigenvalues, we start first from a formula provided by Miehe (see Eq. ( 8) in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF]). In order to obtain a formula more suitable for a numerical implementation, we next simplify this formula by using as much as possible the spectral decomposition and eigenvalue-bases. If C admits coalescent eigenvalues, this formula is not available anymore because it involves a null denominator. To overcome this difficulty, we will use a limiting procedure. The particular situations where the eigenvalues are double or triple then appear as a limiting case of the situation where the eigenvalues are distinct.

Distinct eigenvalues

We start from Eq. ( 8) provided by Miehe in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF] by adapting his notations to ours:

∂M a ∂C = λ a (λ a -λ b )(λ a -λ c ) I ⊗ I - I 3 λ a C -1 ⊗ C -1 + 3 b=1 I 3 λ a λ 2 b -1 M b ⊗ M b . (39)
To simplify Eq. (39), we rewrite the two first terms by using their spectral decompositions:

I ⊗ I - I 3 λ a C -1 ⊗ C -1 = 3 d=1 M d ⊗ 3 e=1 M e - I 3 λ a 3 d=1 1 λ d M d ⊗ 3 e=1 1 λ e M e (40)
or equivalently:

I ⊗ I - I 3 λ a C -1 ⊗ C -1 = 3 d=1 3 e=1 1 - I 3 λ a λ d λ e M d ⊗ M e . (41)
It is easy to see that the terms in the bracket are equal to different values depending on the summation indices d and e:

1 - I 3 λ a λ d λ e = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 -I 3 λ a λ 2 a if d = e = a,
1 -I 3 λ a λ 2 b if d = e = b,
1 -I 3 λ a λ 2 c if d = e = c, 1
-I 3
It then follows from Eqs. (41), (42) that:

I ⊗ I - I 3 λ a C -1 ⊗ C -1 = 1 - I 3 λ a λ 2 a M a ⊗ M a + 1 - I 3 λ a λ 2 b M b ⊗ M b + 1 - I 3 λ a λ 2 c M c ⊗ M c + 1 - I 3 λ 2 a λ b (M a ⊗ M b + M b ⊗ M a ) + 1 - I 3 λ 2 a λ c (M a ⊗ M c + M c ⊗ M a ). ( 43
)
It is noticed that:

M a ⊗ M a = M a ⊗ M a (44)
because:

(M a ⊗ M a ) ij kl = (M a ) ij (M a ) kl = (n a ) i (n a ) j (n a ) k (n a ) l = 1 2 (n a ) i (n a ) k (n a ) j (n a ) l + (n a ) i (n a ) l (n a ) j (n a ) k = 1 2 (M a ) ik (M a ) jl + (M a ) il (M a ) jk = (M a ⊗ M a ) ij kl . ( 45
)
The first three terms in the right-hand side of Eq. (43) thus appear as the opposite to the summation term into the bracket of Eq. (39). Consequently, by reporting Eq. (43) in Eq. (39), the derivative of the eigenvalue-bases M a with respect to C is reduced to:

∂M a ∂C = λ a {(1 -λ c /λ a )(M a ⊗ M b + M b ⊗ M a ) + (1 -λ b /λ a )(M a ⊗ M c + M c ⊗ M a )} (λ a -λ b )(λ a -λ c ) (46)
or equivalently:

∂M a ∂C = 3 b=1,b =a M a ⊗ M b + M b ⊗ M a λ a -λ b . (47)
It is noted that Eq. ( 47) is more compact than Eq. (39) and it does not require the calculation of C -1 . So, this formula is well adapted for a numerical implementation. Now, we are achieving the tangent operator calculation by reporting Eq. (47) in Eq. (38):

∂ 2 W 1 ∂C 2 = α (α -1) 3 a=1 λ α-2 a M a ⊗ M a + 3 a=1 λ α-1 a 3 b=1,b =a M a ⊗ M b + M b ⊗ M a λ a -λ b . (48)
By sorting and factorizing terms in Eq. ( 48), we obtain finally:

∂ 2 W 1 ∂C 2 = α (α -1) 3 a=1 λ α-2 a M a ⊗ M a + 3 a=1 λ α-1 a -λ α-1 b λ a -λ b [M a ⊗ M b + M b ⊗ M a ] . (49)
To the best of the authors' knowledge, this tangent operator expression of Ogden's law has not been published in the literature. It should be noticed that Eq. ( 49) is suitable for a computational implementation because of its very simple and compact form. Moreover, it only involves eigenvalues of C, computed by Eqs. ( 8), ( 9), and eigenvalue bases computed by Eq. ( 18). So, the implementation of Eqs. ( 8), ( 9), ( 18) will be useful for calculation of both the stress tensor and the stiffness tangent operator. Computational time will thus be optimized. In addition, Eq. (49) will turn out to be very suitable to deal with the case of coalescent eigenvalues. In fact, it involves a ratio, (λ α-1 a -λ α-1 b )/(λ a -λ b ), which is well adapted to the limiting procedure as it will be shown in the two next subsections.

Double eigenvalues

In this subsection, we consider the special case of double eigenvalues (say λ 1 = λ 2 = λ 3 ). As outlined in the above section, Eq. ( 49) is no longer valid because it involves a ratio whose denominator is becoming equal to zero. Consequently, we are considering Eq. ( 49) by assuming that all eigenvalues are distinct and studying the limiting case λ 1 → λ 2 . By this way, it is noticed that the ratio making trouble in Eq. (49) reaches the following limit:

lim λ 2 →λ 1 λ α-1 1 -λ α-1 2 λ 1 -λ 2 = (α -1)λ α-2 1 . (50)
Eq. (49) can thus be written as:

∂ 2 W 1 ∂C 2 = α(α -1)λ α-2 1 (M 1 ⊗ M 1 + M 2 ⊗ M 2 + M 1 ⊗ M 2 + M 2 ⊗ M 1 ) + α λ α-1 1 -λ α-1 3 λ 1 -λ 3 (M 1 ⊗ M 3 + M 3 ⊗ M 1 + M 2 ⊗ M 3 + M 3 ⊗ M 2 ) + α(α -1)λ α-2 1 λ α-2 3 M 3 ⊗ M 3 . (51)
By using Eqs. ( 25), (44), it is easy to show that the bracketed terms in Eq. ( 51) can be simplified as:

M 1 ⊗ M 1 + M 2 ⊗ M 2 + M 1 ⊗ M 2 + M 2 ⊗ M 1 = (I -M 3 ) ⊗ (I -M 3 ), (52)
M 1 ⊗ M 3 + M 3 ⊗ M 1 + M 2 ⊗ M 3 + M 3 ⊗ M 2 = (I -M 3 ) ⊗ M 3 + M 3 ⊗ (I -M 3 ). (53)
Eq. ( 51) is thus reduced to:

∂ 2 W 1 ∂C 2 = α(α -1) λ α-2 1 (I -M 3 ) ⊗ (I -M 3 ) + λ α-2 3 M 3 ⊗ M 3 + α λ α-1 1 -λ α-1 3 λ 1 -λ 3 (I -M 3 ) ⊗ M 3 + M 3 ⊗ (I -M 3 ) . (54)
As we can see, Eq. ( 54) is obtained in a very convenient way from Eq. (49). Its compact and simple form is well adapted to a numerical implementation because only λ 1 ,λ 3 and M 3 are involved and these quantities have been already computed for stress evaluation. That means that no additional computation is required to apply Eq. (54).

Triple eigenvalues

In this subsection, we consider the special case of triple eigenvalues (λ 1 = λ 2 = λ 3 ). For analogous reasons as in the preceding subsection, Eq. ( 49) is not valid. Eq. ( 49) is thus considered by assuming first that all eigenvalues are distinct and by studying secondly the limiting case λ 2 → λ 1 and λ 3 → λ 1 . We have to study consequently the limiting behavior of the three ratios involved in Eq. ( 49), i.e., (λ α-1

1 -λ α-1 2 )/(λ 1 -λ 2 ), (λ α-1 1 -λ α-1 3 )/(λ 1 -λ 3 ) and (λ α-1 2 -λ α-1
3 )/(λ 2 -λ 3 ). The first one has been already studied in the preceding subsection. Its limit is provided by Eq. ( 50). The second ratio can be studied in the same way as the first:

lim λ 3 →λ 1 λ α-1 1 -λ α-1 3 λ 1 -λ 3 = (α -1)λ α-2 1 . ( 55
)
The third ratio is deduced from a first order development of λ α-1 2 and λ α-1 3 :

λ α-1 2 = λ α-1 1 + (λ 2 -λ 1 )(α -1)λ α-2 1 +•••, ( 56
)
λ α-1 3 = λ α-1 1 + (λ 3 -λ 1 )(α -1)λ α-2 1 +•••. ( 57
)
By subtracting Eq. ( 56) from Eq. ( 57):

λ α-1 2 -λ α-1 3 = (λ 2 -λ 3 )(α -1)λ α-2 1 +••• (58)
it follows that:

lim λ 2 →λ 1 λ 3 →λ 1 λ α-1 2 -λ α-1 3 λ 2 -λ 3 = (α -1)λ α-2 1 . (59)
We then deduce from Eqs. ( 50), ( 55) and (59) that Eq. ( 49) can be written as:

∂ 2 W 1 ∂C 2 = α(α -1)λ α-2 1 3 a=1 M a ⊗ M a + 3 a=1 [M a ⊗ M b + M b ⊗ M a ] . (60)
To simplify Eq. (60), we use Eqs. ( 25), (44):

3 a=1 [M a ⊗ M b + M b ⊗ M a ]=I ⊗ I - 3 a=1 M a ⊗ M a ( 61
)
and we obtain finally a simple formula to calculate the second derivative in case of triple eigenvalues:

∂ 2 W 1 ∂C 2 = α(α -1)λ α-2 1 I ⊗ I. (62) 4.2. Computation of ∂ 2 W 2 ∂C 2
We now deal with the potential of the form W 2 = I β 3 in Eq. ( 37) 2 . As this potential does not have problem of regularity encountered for the potential W 1 , the calculation is much simpler. By using the following standard result (see for example Eq. (3.37) b in [START_REF] Schröder | Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions[END_REF]):

∂I 3 ∂C = I 3 C -1 (63)
we obtain then:

∂W 2 ∂C = βI β-1 3 ∂I 3 ∂C = βI β 3 C -1 (64)
and

∂ 2 W 2 ∂C 2 = β βI β-1 3 ∂I 3 ∂C ⊗ C -1 + I β 3 ∂C -1 ∂C . ( 65
)
We then use another standard result (refer for example to p. 892 in [START_REF] Miehe | Computation of isotropic tensor functions[END_REF]):

∂C -1 ∂C =-C -1 ⊗ C -1 (66)
to obtain finally:

∂ 2 W 2 ∂C 2 = βI β 3 βC -1 ⊗ C -1 -C -1 ⊗ C -1 . ( 67
)
It should be noticed that a numerical computation of C -1 can be avoided by using the spectral decomposition of C -1 :

C -1 = 3 a=1 λ -1 a M a . (68)
Actually, in case of distinct eigenvalues, the eigenvalue-bases M a are calculated by Eq. ( 18) and C -1 can be directly deduced from Eq. ( 68). In case of double eigenvalues, Eq. ( 68) can be simplified as:

C -1 = λ -1 1 (M 1 + M 2 ) + λ -1 3 M 3 (69)
or, equivalently, by using Eq. ( 25):

C -1 = λ -1 1 (I -M 3 ) + λ -1 3 M 3 . (70)
M 3 being calculated by Eq. ( 26), C -1 is directly deduced from Eq. (70). Finally, in case of triple eigenvalues, Eq. ( 68) is reduced to:

C -1 = λ -1 1 I (71)
which gives again a direct way to determine C -1 . It could be noticed that an alternative method to calculate C -1 is to use the Cayley-Hamilton theorem:

C -1 = C 2 -I 1 C + I 2 I I 3 . (72)
However, Eq. (72) requires the computation of C 2 which is more costly than the matrix operations involved in Eqs. (68), (70), (71).

Finite element implementation

For the finite element implementation, the total Lagrangian formulation is adopted in this work to describe nonlinear behavior. It is well known that the strain tensor E and the stress tensor S are both symmetric. Thus, we note hereafter E and S in vector form as

E = ( E 11 E 22 E 33 2E 12 2E 13 2E 23 ) T , S = ( S 11 S 22 S 33 S 12 S 13 S 23 ) T .
In the context of the finite element method and with Eqs. ( 1) and ( 2), the virtual Green-Lagrange strain on an integration point can be formally written with linear and non-linear contributions in terms of virtual nodal displacements δu e [START_REF] Simo | Computational inelasticity[END_REF]:

δE = (B L + B NL )δu e , (73)
where B L is the matrix which relates the linear part of the strain term to the nodal displacements, and B NL , the matrix which relates the non-linear strain term to the nodal displacements. According to the principle of virtual displacement, the virtual work

δU = V 0 δE T S dV -δu T F ext (74)
has to vanish for every virtual displacement δu. In Eq. ( 74), V 0 is the domain of the initial configuration and F ext is the vector of external loads. In view of Eqs. (34), (73), it comes

δS = DδE = D(B L + B NL )δu e , (75)
where D denotes the usual material secant tangent matrix, deduced from the stress-strain tangent operator D (Eq. ( 35)) due to its major and minor symmetry. More explicitly, we write

D = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ D 1111 D
Substituting δE from Eq. (73), Eq. ( 74) results in

δU = δu T V 0 (B L + B NL ) T S dV -δu T F ext . (77)
The vector of internal forces of a finite element is defined by

F e int = V e 0 (B L + B NL ) T S dV. (78)
The total vector of internal of n elem elements is obtained by a standard assembling procedure

F int = n elem Λ e=1 F e int . ( 79
)
Since δu is arbitrary, the following set of non-linear equations is obtained:

F int -F ext = 0. ( 80
)
This equation is strongly non-linear with respect to the nodal displacements u, because of finite strains and large displacements of solid. A typical solution procedure for this type of non-linear analysis is obtained by using the Newton-Raphson iterative procedure [START_REF] Simo | Computational inelasticity[END_REF]:

K i T u = F ext -F i int , u i+1 = u i + u (81)
where i and i + 1 are the iteration numbers at which the equations are computed.

K T = Λ n elem e=1 K e
T is the tangent stiffness matrix and u, the vector of nodal displacements correction. Taking the derivative of F e int with respect to the nodal displacements u e gives the elementary tangent stiffness matrix as

K e T = ∂F e int ∂u e = K e 0 + K e σ + K e u , (82)
where K e 0 , K e σ and K e u stand respectively for the elastic stiffness matrix, the geometric stiffness (or initial stress stiffness) matrix and the initial displacement stiffness matrix:

K e 0 = V e 0 B T L DB L dV, (83)
K e σ = V e 0 ∂B T NL ∂u S dV, (84)
K e u = V e 0 B T L DB NL + B T NL DB L + B T NL DB NL dV. (85)

Numerical results

The formulation described above has been implemented in the finite element software FER developed by the authors. In order to validate the formulation, we propose here three test examples.

Homogeneous deformation

The first example is concerned with the homogeneous deformation of a cube composed of an elastic material obeying a three-term Ogden model (N = 3) with the material constants: μ 1 = 0.66, μ 2 = 0.0012, μ 3 =-0.01, α 1 = 1.3, α 2 = 5, α 3 =-2, β 1 = 10, β 2 = 10, β 3 = 10. The size of the square cube is 0.5 × 0.5 × 0.5 mm (Fig. 1). Three cases are considered:

• Case A: only F z is applied on the cube. It corresponds to the case of double eigenvalues. In this case (simple tension test), Ogden [START_REF] Ogden | Non-linear Elastic Deformations[END_REF] (p. 499) has given an analytical solution as

F z = A 0 N i=1 μ i λ α i -1 z -λ α i /2-1 z , (86)
where A 0 is the initial section and λ z the stretch in the z direction. • Case B: F z and F x are applied on the cube such that u x = 0. It corresponds to the case of distinct eigenvalues.

• Case C: F z , F y and F x are applied on the cube such that u x = u y = u z . It corresponds to the case of triple eigenvalues.

Figs. 234show the evolution of the applied force versus the displacement u z . It is shown from Fig. 2 that the computed force is identical with the analytical solution. It is still noted that in Cases A and B, small forces (< 1.1N)are enough to achieve large deformations (u z = 3 mm). However, in Case C, much larger forces (> 110 N) are necessary to obtain small deformations (u z = 0.1mm).

Non-homogeneous deformation -pure torsion of a cylinder

In this section, the efficiency of the limiting technique implemented in FER is tested with a 3D example involving non-homogeneous deformations. This example concerns the pure torsion problem of a circular cylinder (Fig. 5). If compressible materials are used, the situation is generally complicated because radial extension can occur and the cross-section of the cylinder does not remain circular. However, in the case of the Blatz-Ko model [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF], it has been demonstrated in [START_REF] Polignone | Pure torsion of compressible nonlinearly elastic circular cylinders[END_REF] that a pure torsion is possible, i.e. a deformation with zero radial displacement. The Blatz-Ko model is recovered by an appropriate selection of the Ogden's model parameters. By setting N = 1, μ 1 =-G, α 1 =-2 and β 1 = 0.5, the strain energy density function [START_REF] Gent | A new constitutive relation for rubber[END_REF] is actually reduced to By using Eq. ( 6), the summation over the eigenvalues of C -1 can be related to the second and third invariants of C by

W = G 2 λ -1 1 + λ -1 2 + λ -1 3 -3 + G I 3 -1 . (87)
λ -1 1 + λ -1 2 + λ -1 3 = I 2 I 3 . (88)
Eq. ( 87) is thus equivalent to the Blatz-Ko energy density function

W = G 2
I 2 I 3 + 2 I 3 -5 . (89)
In the case of a pure torsion deformation, the displacement u is written as follows

u = x -X, x = x y z = cos(θ ) -sin(θ ) 0 sin(θ ) cos(θ ) 0 00 1 X, X = X Y Z . ( 90
)
Where X and x respectively denote the Lagrangian and the Eulerian positions. The rotation angle θ is related to the twist angle τ (rad/m) and to the third Lagrangian coordinate Z by

θ = τZ. (91)
By using Eq. (90), the deformation gradient F and the right Cauchy-Green deformation tensor C are written as

F = cos(τ Z) -sin(τ Z) -τy sin(τ Z) cos(τ Z) τ x 00 1 , (92)
C = F T F = 10 -τY 01 τX -τY τX τ 2 R 2 + 1 . ( 93
)
Where R represents the radius (R 2 = X 2 + Y 2 ). Eq. (92) shows that the deformation is non-homogeneous. Besides, it results from Eq. (93) that the eigenvalues of C are respectively equal to 1, 1 and τ 2 R 2 + 1. A double coalescence thus occurs, except for the cylinder axis (R = 0) where a triple coalescence arises. It is demonstrated in [START_REF] Polignone | Pure torsion of compressible nonlinearly elastic circular cylinders[END_REF] that the stress response corresponding to Eq. (92) satisfies the equilibrium equation in the absence of body forces, gives a traction-free lateral surface satisfying the boundary condition and includes an axial stress component. The only nonzero Cauchy stresses are thus σ zz and σ θz . These stress components both radially expand in a quadratic form for σ zz and in a linear form for σ θz A 0.1 m height cylinder with a radius equal to 0.05 m was considered for the numerical application. This cylinder was meshed with 11616 eight-node brick elements (Fig. 5(a)). The bottom of the cylinder is clamped and the displacements are restrained on its top face according to Eq. ( 90). The twist angle τ was selected to 900 • /m in order to provide a rotation of the top face equal to 90 • (Fig. 5(b)). The shear modulus G of the material is equal to 220711 Pa according to the measurements made by Blatz and Ko [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF]. Figs. 6 and7 show respectively the distribution of stresses σ zz and σ θz . Table 1 gives the maximum values of stresses obtained by numerical and analytical solutions. A good agreement was found as we can see from this table.

σ zz =-Gτ 2 R 2 ,σ θz = Gτ R. (94)

Non-homogeneous deformation -internally pressurized hollow sphere

This example concerns the pure radial deformation problem of an internally pressurized hollow sphere. This problem was analytically solved by Chung et al. with the compressible Blatz-Ko material model [START_REF] Chung | The finite deformation of internally pressurized hollow cylinders and spheres for a class of compressible elastic materials[END_REF]. Since it is not possible to obtain a closed form solution with the Ogden model, a finite element analysis is required. The gradient deformation matrix can however be expressed in a simple manner by [START_REF] Ogden | Non-linear Elastic Deformations[END_REF]: where r and R respectively represent the radius in deformed and reference configurations. To obtain the solution, only half of the sphere was modeled with 8-node axisymmetric elements. The material constants are the same as in the example of Section 6.1. The internal and external radii of the sphere R 1 and R 2 are respectively 1 and 2. Fig. 8 shows the initial mesh and the deformed shape with von Mises stress distribution with a maximum value of 0.761 corresponding to an internal pressure of 0.2. It is easily deduced from Eq. ( 95) that ( dr dR ) 2 and ( r R ) 2 are simple and double eigenvalues of C. The deformation is non-homogeneous as shown on Fig. 9 where the double eigenvalue is represented vs. R for different values of the internal pressure p. It is worth noting that the double eigenvalue of C depends on the spatial coordinate, contrarily to the example of Section 6.2 where the double eigenvalue was constant and equal to 1.

F = f(R)I + f ′ (R) R X ⊗ X with f(R)= r R (95)

Conclusions

In this paper, we have presented detailed calculations of stress and stress-strain tangent operator for the finite element implementation of the Ogden's hyperelastic model. To avoid eigenvectors computation, we have selected a spectral decomposition expressed in terms of eigenvalue-bases. If the eigenvalues of the right Cauchy-Green strain tensor are distinct, a simple and compact tensor formula, well adapted for numerical implementation, is proposed. A limiting technique is used to account for the special case of coalescent eigenvalues in which non-differentiability occurs.

To improve numerical results and computational efficiency, we have used tensor calculus and keep as much as possible eigenvalue-bases expression. The tensor formula obtained offers two advantages. Firstly, its compact and simple form appears to be suitable for a numerical implementation. Secondly, it provides a very convenient way to perform the limiting procedure. Implementation in the university finite element code FER is presented and three test examples are performed, including homogeneous and non-homogeneous deformations, to validate the developed algorithms.

  if d = a and e = b or d = b and e = a, 1 -I 3 λ 2 a λ c if d = a and e = c or d = c and e = a, 0i f d = b and e = c or d = c and e = b.

  1122 D 1133 D 1112 D 1113 D 1123 D 1122 D 2222 D 2233 D 2212 D 2213 D 2223 D 1133 D 2233 D 3333 D 3312 D 3313 D 3323 D 1112 D 2212 D 3312 D 1212 D 1213 D 1223 D 1113 D 2213 D 3313 D 1213 D 1313 D 1323 D 1123 D 2223 D 3323 D 1223 D 1323 D 2323

Fig. 1 .

 1 Fig. 1. Square cube.

Fig. 2 .

 2 Fig. 2. Case A: double eigenvalues.

Fig. 3 .

 3 Fig. 3. Case B: distinct eigenvalues.

Fig. 4 .

 4 Fig. 4. Case C: triple eigenvalues.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Pure torsion of a circular cylinder. (a) Mesh; (b) Deformed geometry.

Fig. 8 .

 8 Fig. 8. Deformed shape and von Mises stress in a pressurized hollow sphere.

Fig. 9 .

 9 Fig. 9. Evolution of the double eigenvalue of C.

Table 1

 1 Comparison between analytical and numerical solutions.

	Stress	Analytical	Numerical	Relative error (%)
	σ zz (Pa)	-136146	-133281	2.1
	σ θz (Pa)	173346	171208	1.2