G Chagnon

G Marckmann

E Verron

A COMPARISON OF THE HART-SMITH MODEL WITH ARRUDA-BOYCE AND GENT FORMULATIONS FOR RUBBER ELASTICITY

The present paper demonstrates that the Hart-Smith constitutive model and the more recent Arruda and Boyce eight chains and Gent constitutive models are closely related. The ability of these three models to predict both small and large strain responses of rubbers is highlighted and equations that relate their material parameters are established.

INTRODUCTION

Elastomers exhibit a very complex mechanical behavior, that includes large strain, timedependent response, hysteresis, and stress-softening (Mullins effect). In the last few years, there was a significant effort to propose realistic constitutive equations for these materials. 1-10 These models focus on the whole behavior of elastomers or on only one of the previously mentioned phenomena. Nevertheless, all constitutive equations are based on the general theory of hyperelasticity. 11 So, consider now the non-linear elastic response of rubber-like material. It is characterized by the existence of a strain energy potential W. Assuming that elastomers are both isotropic and incompressible, this strain energy function only depends on the two first invariants of the left Cauchy-Green stretch tensor B: (1) with:

(2)

The true stress tensor is defined by the differentiation of W with respect to B: (3) This general theory being established, the major difficulty remains the choice of the function W. Many developments were proposed in the bibliography (see for example the review of Boyce). 12 An efficient hyperelastic model can be defined by four major qualities:

• First, it should be able to accurately reproduce the whole "S" shaped response of rubbers, • Second, the change of deformation modes should not be problematic, i.e. if the model behaves satisfactorily in uniaxial tension, it should also be quite accurate in simple shear , or in equibiaxial extension, • Third, the number of relevant material parameters must be as small as possible, in order to reduce the number of experimental tests needed for their identification, • Finally, the mathematical formulation has to be quite simple to render possible the numerical implementation of the model. Taking into account these four prerequisites, three constitutive equations can be selected: the recent models of Arruda and Boyce, 13 and Gent, 14 and the previous formulation of Hart-Smith. 15 Through the rest of the paper, they will be referred to as AB, G and HS theories, respectively.

σ = -+ ∂ ∂ = -+ ∂ ∂ + ∂ ∂       - ∂ ∂ p W p W I I W I W I I B B I B B 2 2 2
The aim of the present work is to establish the equations that relate these three constitutive equations. More precisely, their small and large strain characteristics will be thoroughly examined. As Boyce 16 recently compared the AB and G approaches, the emphasize is laid on the equivalence of the HS model with the two other ones. In the next Section, the three strain energy functions are briefly recalled. Then, the theoretical relationships between material parameters of the models are established and the relevance of these results is demonstrated. The relevance of these results is demonstrated in Section 3. Finally, some concluding remarks are proposed.

STRAIN ENERGY FUNCTIONS

In this section, the strain energy functions that define the HS, AB and G models are briefly recalled.

THE HS MODEL

In the 60's, Hart-Smith 15 proposed an empirical form of W, that satisfactorily reproduces the response of rubberlike materials for the entire range of deformation endured by rubbers: (4) where C 1 , C 2 , C 3 are the three material parameters. In this equation, the first right-hand side part of W describes the global response of the material. It only depends on the first strain invariant I 1 . The second term that involves I 2 is based on the work of Gent and Thomas 17 and improves the accuracy of the model for moderate strain (less than 150%). It can be related to the Mooney-Rivlin theory that involves both I 1 and I 2 and is widely recognized for its ability to satisfactorily reproduce rubberlike behavior for small and moderate strain. 11

THE AB MODEL

The model of Arruda and Boyce 13 is a molecular based constitutive equation, also known as the 8-chains model. Following Kuhn and Grün, 18 the behavior of an individual chain is considered non-gaussian. In order to describe the response of the polymer network, this individual response should be integrated in an unit sphere defined by its chain density. Nevertheless, the corresponding constitutive equation necessitates numerical integration. 19 So, in order to overcome this difficulty, privileged directions for chains orientation are defined by the half diagonals of a unit cube contained in the unit sphere. The corresponding strain energy function is given by:

W C C I dI C I HS = - ( ) [ ] + ∫ 1 3 1 2 1 2 2 3 3 exp ln (5)
where with L being the Langevin function defined as . In the previous equation, C R and N are the two material parameters. Due to the non-gaussian nature of the chain modelling, 3N stands for the maximum value that I 1 can reach and represents the extensibility limit of chains.

This model only depends on the first strain invariant I 1 , and it is well-known that its performances at moderate strain are not very good. 19 Very recently, Meissner and Matejka 20 proposed the improvement of chains models by adding a function of I 2 to the strain energy. Moreover, Fried 21 explained the dependance of W on I 2 invoking assumptions on the network deformation.

THE G MODEL

Recently, Gent 14 proposed an empirical relation for the strain energy function: [START_REF] Beatty | [END_REF] where E and J m are the two material parameters. J m represents the maximum value of I 1 -3 that can be undergone by the material. It is the empirical counterpart of the extensibility limit of chains previously defined for the AB model.

COMPARISON OF MATERIAL PARAMETERS

The present section is devoted to the comparison of the three models. As these constitutive equations are supposed to be qualitatively efficient for the entire range of strains, both small and large strain responses will be compared. As mentioned above, the I 2 terms in constitutive equations can be seen as corrections of the phantom network theory. 11 So, we propose to only compare the part of functions W that involves I 1 , i.e. the phantom part of constitutive equations. Consequently, the following HS strain energy will be considered through the rest of the paper: (7) in which only two material parameters are considered, i.e. C 1 and C 3 .

Moreover, recalling that a previous work demonstrated the equivalence of the AB and the G approaches, 16 the present work emphasizes the comparison of the HS model with them.

POLYNOMIAL EXPANSION

Before explicitly comparing the models, their polynomial expansion in terms of I 1 are established. This will highly simplify the comparison of their small strain responses.

The HS strain energy function (7) can be developed into an infinite series of powers of I 1 :

W C C I dI HS = - ( ) [ ] ∫ 1 3 1 2 1 3 exp W E I J G m = - - -       6 1 3 1 ln L x x x ( ) = ( )- coth 1 β = ( ) - L I N 1 1 3 W C N I N AB R = +       1 3 β β β ln sinh (8)
It is to note that this expression is a particular case of the general expansion proposed by Rivlin 16 :

(9)
where the dependance on I 2 is not considered, said s = 0, and only odd powers of I 1 are relevant.

The first five terms of the development of W AB , given by Equation ( 5), into a series of powers of I 1 are: (10) This expansion is also a particular case of Rivlin's series (9), but, in comparison with Equation ( 8), it is expressed with both odd and even powers of I 1 .

Similarly to previous expansions, W G (see Equation ( 6)) can be expanded into a series of powers of I 1 : (11) Here, both odd and even powers of I 1 are involved.

SMALL STRAIN BEHAVIOR

The small strain stiffness of the models are respectively defined by first terms of Equations ( 8), (10) and (11). Then, the three strain energy functions reduce to the neo-Hookean expression 11 : (12) where C is a material parameter. Moreover, the first parameters of the HS, AB and G models are proportional. They are simply related to C by: (13) LARGE STRAIN BEHAVIOR As mentioned above, the three models are characterized by their ability to describe the strain-hardening of the material that takes place under large strain. This strain-hardening phe- ... nomenon is mainly due to the extensibility limit of polymer chains. According to Equation ( 3) and taking into account that the considered strain energy functions do not depend on I 2 , i.e. ∂W/∂I 2 =0, the stress-strain relationship reduces to: (14) in which only the first derivative of W with respect to I 1 is involved.

C C C E R = = = 2 3 1 W C I NH = - ( ) 1 2 3 1 W E J I m n n n = - ( ) + = ∞ ∑ 6 1 3 1 1 1 W C I I rs r s r s = - ( ) - ( ) ( )=( ) ∞ ∑ 1 2 0 0 3 3 , , W C C i i I HS i i i = + ( ) - ( ) = ∞ + ∑ 1 3 0 1 2 1 2 1 3 ! W C I N I N I N I N I AB R = - ( ) + - ( ) + - ( ) + - ( ) + - ( ) +       1 2 3 1
In the three cases, this derivative is given by:

(15) (16) (17)
Examining these equations, it is obvious that the three functions exhibit an upturn for large values of I 1 . Moreover, for both AB and G models, ∂W/∂I 1 tends asymptotically to infinity as the first invariant approaches a limiting value. Vertical asymptotes correspond with values of I 1 equal to 3N and J m +3 for AB and G constitutive equations, respectively. Thus, the second material parameters of AB and G models are obviously related by: (18) Consider now the HS model. It does not exhibit an asymptotic behavior, because it involves an exponential-like function. Nevertheless, in order to ensure that models reproduce the same behavior at large strain, it is important that they admit similar slopes in this range of deformation. Here, the HS model is only compared with the G model, the equivalence of AB and G models being already established in Equation (18). Slopes of stress-strain curves can be expressed in terms of the second derivatives of W with respect to I 1 that are:

(19) (20)
Thus, to ensure that the responses of HS and G models are similar, these two derivatives should be equal. Eliminating their first material parameters, i.e. C 1 and E, respectively, by the use of Equation ( 13), it yields to:

∂ ∂ = -- ( ) [ ] 2 1 2 1 2 6 1 1 3 W I E J I J G m m ∂ ∂ = - ( ) - ( ) [ ] 2 1 2 1 3 1 3 1 2 2 3 3 W I C C I C I HS exp J N m = - ( ) 3 1
∂ ∂ = -- ( ) W I E I J G m 1 1 6 1 1 3 ∂ ∂ =       - W I C N I L I N AB R 1 1 1 1 2 3 3 ∂ ∂ = - ( ) [ ] W I C C I HS 1 1 3 1 2 3 exp σ = -+ ∂ ∂ p W I I B 2 1 C C C E R = = = 2 3 1 (21)
which lies the second material parameters C 3 and J m , and depends on I 1 . This equation means that the equality of slopes depends on the considered deformation level. In order to solve this equation, it can be cast under the following form:

The solution of this equation can be expressed in terms of the Lambert W function, denoted here W L : (23) where is defined to be the function satisfying (see Corless et al. 23 and the references herein for details): (24) So, defining now the reduced first strain invariant α by: (25) where α∈[0,1] as J m represents the asymptotic value of I 1 -3, and considering some basic algebraic manipulations, Equation ( 23) can be transformed into: [START_REF] Ogden | Proc. R. Soc. Lond. A[END_REF] for α∈ [0,1] This equation relates the second material parameters C 3 and J m that satisfy the equality of slopes at large strain for HS and G models. As mentioned above, it depends on α, i.e. on the strain level.

Nevertheless, as the G model increases asymptotically to infinity in the neighborhood of I 1 = J m +3, i.e. α = 1, and the HS model is an exponential-like function in this range of I 1 , Equation ( 26) is not defined for α = 1. In practice, the equivalence between the two constitutive equations can only be demonstrated on a given range of α (or of I 1 ). In this way, a mean constant value of the function in the right-hand side of Equation ( 26) can be adopted to relate C 3 and J m . This permits to omit the dependence on α of C 3 J m 2 . Let us denote [α min ,1] with α min >0, the range in which

C J W m L 3 2 2 2 1 2 1 = - ( )       α α α α = - I J m 1 3 W x e x L W x L ( ) = ( ) C I W I J I J L m m 3 1 2 1 1 2 3 3 2 1 1 3 - ( ) = - -- ( ) [ ]         C I C I I J I J m m 3 1 2 3 1 2 1 1 2 3 3 3 2 1 1 3 - ( ) - ( ) [ ] = - -- ( ) [ ] exp 2 3 3 1 1 1 3 3 1 3 1 2 1 2 C I C I J I J m m - ( ) - ( ) [ ] = -- ( )
[ ] exp the HS and G models are close to. Thus, in this range, the two parameters can be linked by a constant coefficient:

(27) with So, the value of k depends on the range in which the equivalence of the two models is recognized. Figure 1 presents the evolution of k as a function of α min . As shown by the graph, k tends to infinity as α min tends to 0 or 1. The practical choice of the bound α min and the corresponding value k will be examined in the next section for two examples. 13), ( 18) and ( 27) that lie their material parameters.

FIT OF THE G THEORY

The first example consists in using the G model as a reference. In order to build fictitious materials and their corresponding data sets to be fitted, the two parameters E and J m are considered to vary in [0,10] and [7,100], respectively. Values of J m stand for limits of extensibility between 200% and 900% in uniaxial tension.

First, J m is set to 60 and E adopts the three values 0.1, 1 and 10. In these three cases, the function ∂W G /∂I 1 given by Equation ( 17) is discretized using a large number of values, and the

k W d L = - - ( )       ∫ 1 1 1 2 1 2 2 1 α α α α α α min min C J k m 3 2 =
corresponding stress functions of HS and AB models, respectively given by Equations ( 15) and ( 16), are fitted with these data. Fitting results are obtained using two different algorithms. First, an optimization program based on a genetic algorithm estimates approximately the material parameters. Second, these parameters are precisely determined with the help of a classical steepest descent algorithm using the genetic algorithm results as initial guess solutions. The linear relation Equation ( 13) between first parameters of the three models is exactly recovered for the three tests.

Second, the initial stiffness E is set to 1, and J m varies between 7 and 100. First parameters of HS and AB models are calculated according to Equation ( 13), said C 1 = 6 and C R = 3. For the considered values of J m , the method employed above for first parameters is used to determine both C 3 and N. Let us first examine the results for the AB model. Figure 2 shows the comparison between fitted results and the theoretical relation of Equation (18). The fitted results agree well with the theory and the linear relation between the asymptotic limit of the two models is satisfactorily recovered. As mentioned above, the case of the HS model is more difficult. Our goal is to verify Equation (27). Thus, the interval of α should be defined to calculate the integral in this equation. In order to choose the lower bound α min , we consider that slopes of HS and G models should not be compared for small strain, because only the large strain responses of these models are examined. As proposed by Kucherskii, 24 the neo-Hookean formulation is valid in the first third of the deformation range. Thus, only the two thirds of the deformation range has to be retained for the large strain response, i.e. α min ≥1/9. Moreover, as models responses are highly different in the neighborhood of α = 1, because the G model tends asymptotically to infinity and the HS model grows exponentially, we consider that the equivalence should be demonstrated for more than a third of the deformation range, said α min ≤4/9. Consequently, according to Equation ( 27) and noting that k does not evolve monotonously in α min ∈ [1/9, 4/9], the range in which k varies can be easily established: Let us now examine the relevance of the previously demonstrated relationships as applied to the parameters given in Table I. In order to simplify the discussion, G material parameters are considered as the reference and material parameters of AB and HS theories are compared with them.

• Concerning the first material parameter of models, we first compute 6C 1 = 1.08 for the HS model and 3C R = 1.02 for the AB theory as proposed in Equation ( 13). Thus, the maximum difference between these data and E = 1 is about 8%. • Second, the large strain parameters J m and N of G and AB models are analyzed.

According to Equation ( 18), the theoretical value of N should be (J m + 1)/3 = 3.10. The difference between this predicted value and the fitted result is about 10%. • Finally, we consider the relationship between J m and C 3 . The corresponding theoretical relation that lies these data is given in Equation ( 27). In order to compute k, the lower bound α min of the definition range should be determined. Examining experimental results for the uniaxial test shown in Figure 4(a), the small strain range, in which the neo-Hookean theory is sufficient, can be defined as λ∈ [1,4]. This interval corresponds to I 1 ∈[3,16.5] and the use of Equation ( 25) yields to α min ≈0.14. In regard to Equation (27), k is defined by an integral from α min to 1. Nevertheless, considering that the last experimental point of the uniaxial tensile test corresponds with λ max = 7.6, i.e. I 1 = 58.0, the upper bound of the integration range is not equal to 1 but should be replaced by α max = (I 1max -3)/Jm≈0.6. In fact, it is impossible to reach experimentally the extensibility limit of chains, because samples break for smaller stretches. Finally, the theoretical value of k to be considered is given by:

(29)

It is to note that this value is not contained in the range of Equation (28). Using this result in Equation ( 3), the second parameter of the HS model C 3 should have been equal to 3.18 x10 -4 . Comparing this value with the fitted result given in Table I leads to a discrepancy of 17.7%. In regards with differences obtained between our theoretically predicted values and fitted results, the present developments are validated. The equivalence Equations ( 13), ( 18) and ( 27) can now be used to switch between HS, AB and G hyperelastic constitutive equations.

CONCLUSIONS

The present work establishes the equations that relate three hyperelastic constitutive equations: the HS, AB and G models. Here, only the I 1 part of the models were compared. As proposed elsewhere, 20 the accuracy of the three approaches for moderate strain can be improved by adding a function of I 2 to W. This can be performed with only one additional material parameter. Note that the original Hart-Smith model included this I 2 term. The HS model presents an advantage that will highly simplify the numerical implementation in finite element softwares: its large strain behavior is governed by an exponential-like function so that its tangent stiffness always exists. Moreover, the relationships established here will also simplify the identification of material parameters, because HS parameters can be now related to physical concepts through their equivalence with AB data.

FIG. 1 .

 1 FIG. 1. -Curve of Equation (27): k vs α min . EXAMPLES Finally, two examples are studied to demonstrate the equivalence of the three constitutive equations, and the relevance of Equations (13), (18) and (27) that lie their material parameters.

FIG. 2 .

 2 FIG. 2. -Comparison of second parameters for the G and AB models: (ο) fitted results, () theoretical results of Equation (18).

  FIG. 3. -Comparison of second parameters for G and HS models: (ο) fitted results, (----) theoretical results of Equation (27) for k = 2.83, () theoretical results of Equation ( 27 ) for k = 2.98.

FirstFIG. 4 .

 4 FIG. 4. -Fit of Treloar's data, (a) uniaxial tension; (b) pure shear; (c) equibiaxial tension; (d) biaxial tension; (ο) experiments; (-) HS model; (----) AB model; (…) G model.

TABLE I VALUES

 I OF MATERIAL PARAMETERS


Nicolas Chevaugeon

Erwan Verron

Bernard Peseux

Finite element analysis of non linear transversaly isotropic hyperelastic membranes for thermoforming applications

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The finite element simulation of thermoforming process for isotropic materials has been extensively examined in the past. Both axisymmetrical and three-dimensional models have been developed. In 1989, Zamani et al. [START_REF] Zamani | Status of the finite element method in the thermoforming process[END_REF] reviewed publications dealing with this subject. In the present introduction, we only focus on the more recent works on the subject. deLorenzi and Nied [START_REF] Delorenzi | Blow molding and thermoforming of plastics: finite element modeling[END_REF][START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF] examine both blow-moulding and thermoforming processes simulation. Using experimental observations, they assume that the thermoplastic parison can be modeled by a rubberlike hyperelastic membrane. In the same period, Charrier and coworkers [START_REF] Charrier | Free and constrained inflation of elastic membranes in relation to thermoforming -non-axisymmetric problems[END_REF][START_REF] Shrivastava | Large deformation finite element analysis of non-linear viscoelastic membranes with reference to thermoforming[END_REF] used the same approach for viscoelastic membranes. All these papers report on a quasi-static treatment of the problem. More recently, Bourgin et al. [START_REF] Bourgin | A first step towards the modelling of the thermoforming of plastic sheets[END_REF] and Verron et al. [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF][START_REF] Marckmann | Finite element analysis of blow-moulding and thermoforming process using a dynamic explicit procedure[END_REF] prefer to consider inertial effects and use a dynamic explicit numerical procedure to solve the problem.

All these works are focused on isotropic materials. There are just few studies concerning the thermoforming of fibers-reinforced thermoplastics. Some authors use the approach developed in injection molding process of fiber loaded material [START_REF] Bhattacharyya | Composite sheet forming[END_REF]. The material is described at a micromechanical scale (unit cell model) and is considered as a viscous fluid loaded by rigid fibers. From this microscopic point of view and using homogeneization technics, authors build macroscopic constitutive equations, formulated in terms of strain rate. The corresponding behaviour is anisotropic viscous fluid like [START_REF] Okine | Numerical modeling of sheet forming for thermoplastic composites[END_REF][START_REF] De Luca | An explicit finite element solution for the forming prediction of continuous fibre-reinforced thermoplastic sheets[END_REF][START_REF] Pickett | Numerical simulation of the forming process in long fibre reinforced thermoplastics[END_REF].

The main goal of our work consists in using the classical developments relative to thermoforming of isotropic materials and in adapting them to the general case of anisotropic materials. Following the work of Spencer [START_REF] Spencer | Continiuum Theory of the mechanics of fibre reinforced Composite[END_REF], we choose to describe the anisotropy induced by aligned fibers in the rubberlike thermoplastics at the macromechanical scale using material symmetry properties. The present paper reports on the use of transversely hyperelastic constitutive models in the thermoforming simulation. Such theory was successfully used by Kyriacou et al. [START_REF] Kyriacou | Finite element analysis of nonlinear orthotropic hyperelastic membranes[END_REF] in membrane inflation problems.

In the next section the theoretical background used to describe stress and strain for finite deformation in the Lagrangian point of view is recalled, and the general form of the strain energy function for hyperelastic transversely isotropic materials is established. Section 3 is devoted to the presentation of the finite element formulation suitable to describe the inflation of membrane under an inflating pressure. The numerical implementation of the constitutive equation is highlighted. Some numerical results obtained by our method for free and confined inflation cases are presented in section 4. Last, concluding remarks and perspectives are proposed.

HYPERELASTIC TRANSVERSELY ISOTROPIC CONSTITUTIVE EQUA-TION

In this paper, fibers-reinforced thermoplastic materials are considered. Fibers are assumed to be aligned and materials are studied far above their glass transition temperature. Therefore materials obey hyperelastic transversely isotropic constitutive equations.

Stress and strain definitions

Consider a material point P of a body B. P occupies the position X relative to a frame of reference R 0 : {O, e I } in the undeformed initial configuration C 0 . In the deformed configuration C, x is the position of P relative to an observer frame of reference R : {o, e i }. The transformation between C 0 and C is defined by:

x = ξ(X, t) (1)
where t is the present time.

In the neighborhood of P, the deformation gradient F, the right Cauchy-Green deformation tensor C and the Green-Lagrange strain tensor E are respectively defined by:

F(X, t) = ∂x(X, t) ∂X (2) C(X, t) = F T F (3) E(X, t) = 1 2 [C(X, t) -I] ( 4
)
where I is the second-order identity tensor and • T denotes the transposition. Let σ and S be the Cauchy "true" stress tensor defined in the deformed configuration and the second Piola-Kirchhoff stress tensor relative to the undeformed configuration. They are respectively defined by: df = σ . n ds (5)

dF = S . N dS (6)
in which df is the force in the deformed configuration acting on a deformed oriented area n ds and dF is the force in the undeformed configuration acting on an undeformed oriented area N dS. df and n ds are the transformed vectors of dF and N dS by the deformation gradient. With some standard manipulations the two stress tensors can be related by:

S = J F -1 σ F -T (7)
where J is the Jacobian of the transformation (for a volume-preserving deformation: J = 1).

Construction of a Strain energy function

Hyperelastic materials are defined by the existence of a scalar function W (F) called stored energy or strain energy function, from which stresses can be derived at each point. In order to satisfy the objectivity principle, the strain energy function must be invariant under changes of observer frame of reference. It is well-known that the Cauchy-Green deformation tensor is invariant under changes of observer frame of reference [START_REF] Beatty | Topics in finite elasticity: hyperelasticity of rubber, elastomers and biological tissues -with examples[END_REF]. Thus, if the strain energy function W can be written as a function of C, it automatically satisfies the objectivity principle.

We now have to take into account the material symmetries that restrict the way W depends on C. Every orthogonal transformation member of the material symmetry group must keep the strain energy unchanged. For isotropic material, the material symmetry group consists in the entire group of orthogonal transformations. In the case of unidirectional fibers-reinforced composite, the symmetry group has to be determined. The fibers orientation is defined by a local unit vector denoted a 0 (X). Therefore, the material symmetry group is the group of orthogonal transformations which keep a 0 (X) unchanged. Consequently the material is isotropic in all planes orthogonal to the fibers direction and is called transversely isotropic. The general form of the corresponding strain energy function must satisfy:

W (C) = W (Q.C.Q T ) ∀Q orthogonal with Q.a 0 (X) = a 0 (X) (8)
Another way of writing W consists in introducing explicitly the dependence on a 0 [START_REF] Spencer | Continiuum Theory of the mechanics of fibre reinforced Composite[END_REF][START_REF] Govindjee | Finite implementation of incompressible, transversaly isotropic hyperelasticity[END_REF].

The invariance with respect to the symmetry group is now expressed by:

W (C, a 0 ) = W (Q.C.Q T , Q.a 0 ) ∀Q orthogonal (9)
As the sign of a 0 has no physical meaning, we can introduce the orientation tensor A = a 0 ⊗ a 0 with no loss of information (⊗ stands for the tensor outer product). Then the previous equation ( 9) becomes:

W (C, A) = W (Q.C.Q T , Q.A.Q T ) ∀Q orthogonal ( 10
)
Thus W is an isotropic function of two symmetric tensors C and A. Using representation theorems [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], W can be written as a function of the invariants of C, A and their products. Following equations give a possible set of ten independent scalar invariants on which W can depend [START_REF] Sidoroff | Cours sur les grandes transformations[END_REF]:

I 1 = trC I 2 = 1 2 (trC) 2 -trC 2 I 3 = detC (11) A 1 = trA A 2 = 1 2 (trA) 2 -trA 2 A 3 = detA ( 12
)
I 4 = tr(C : A) I 5 = tr(C 2 : A) I 6 = tr(C : A 2 ) ( 13
)
I 7 = tr(C 2 : A 2 ) ( 14
)
where (I i ) i=1,3 are the invariants of C, (A i ) i=1,3 those of A and (I i ) i=4,7 are product invariants.

In the present case of transversely isotropic behaviour, A = A 2 and the invariants of A become:

A 1 = 1 A 2 = 0 A 3 = 0 (15)
Moreover, some relations can be written for the invariants of the products:

I 5 = I 7 I 4 = I 6 (16)
Therefore, the number of relevant independent invariants on which W depends is reduce to five [START_REF] Burton | A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations[END_REF]:

W = W (I 1 , I 2 , I 3 , I 4 , I 5 ) (17)
Note that in the isotropic case, the number of invariants is reduce to three [START_REF] Treloar | The mechanics of rubber elasticity[END_REF]:

I 1 , I 2 and I 3 .
Taking into account the incompressibility assumption, i.e.: I 3 = 1, the number of independent invariants is finally reduced to four.

The 2nd Piola-Kirchhoff stress tensor can be calculated by differentiation of the strain energy function with respect to the Cauchy-Green dilatation tensor C [START_REF] Beatty | Topics in finite elasticity: hyperelasticity of rubber, elastomers and biological tissues -with examples[END_REF]:

S = -pC -1 + 2 ∂W ∂C ( 18
)
where p is an arbitrary hydrostatic pressure preventing for any volume change. This pressure has no physical meaning and is determined from equilibrium equations. Using the differentiation chain rule, the stress tensor can be expressed as:

S = -pC -1 + 2 i=1,2,4,5 ∂W ∂I i ∂I i ∂C (19)
with:

∂I 1 ∂C = I (20)
∂I 2 ∂C = I 1 • I -C (21)
∂I 3 ∂C = I 3 • C -1 (22) ∂I 4 ∂C = a 0 ⊗ a 0 ( 23
)
∂I 5 ∂C = a 0 ⊗ C • a 0 + a 0 • C ⊗ a 0 (24)

FINITE ELEMENT FORMULATION

In this section, we briefly present the finite element formulation adopted to simulate the thermoforming process. For more details, reader can refer to [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF][START_REF] Marckmann | Finite element analysis of blow-moulding and thermoforming process using a dynamic explicit procedure[END_REF].

Governing equations

The polymeric sheet is supposed to be a membrane and is represented by a twodimensional continuum. Its thickness becomes a function of the position on the midsurface. Assuming that there are no body forces, the Principle of Virtual Work, is expressed on the undeformed configuration for the inertial effects and internal work, and on the deformed configuration for the pressure follower force work:

B 0 δu • ρ 0 üdV 0 = - B 0 δE : SdV 0 + δB δu.TdS ∀δu ( 25
)
where B 0 and δB stand respectively for the volume of the undeformed membrane and the surface of the deformed membrane, ρ 0 is the mass density, ü is the acceleration vector, T is the external surface force due to pressure, and δu is a virtual displacement vector compatible with displacement boundary conditions. The contact with the mold is assumed sticky: the heated parison cools down and becomes stiffer after contact [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF]. Thus the blowing pressure is not sufficient to deform the membrane. Every material points which enter in contact is fixed until the end of the simulation. This assumption simplifies the numerical procedure [START_REF] Khayat | Multiple contact and axisymmetric inflation of a hyperelastic cylindrical membrane[END_REF].

Spatial discretization

The previous equation ( 25) is discretized by the finite element method. In the membrane context, 3-nodes isoparametrical triangular elements are chosen. These elements only deform in their plane, remain flat and triangular The use of these simple elements provides analytical integration which is time saving. But it has a major flow: lot of elements are needed to obtain a good stiffness in high curvature regions of the parison. To avoid this problem, an automatic remeshing technique described in [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF] is implemented. After discretization and integration on each element surface, the problem reduces to a system of ordinary second order differential equations:

[M]{ Ü (t)} = F ext (t) -F int (t) (26)
where [M] stands for the mass matrix, U(t) is the nodal displacement vector, F ext (t),and F int (t) are respectively the external and internal nodal force vectors. By the use of the special lumping technique [START_REF] Taylor | The finite element method. Volume I: Basic formulation and linear problems[END_REF], [M] can be reduced to a diagonal matrix, splitting the system in n independent differential equations where n is the number of degrees of freedom of the model.

Time-integration scheme

In order to solve the previous system (26), we use the classical explicit second-order central difference method. The velocity of the degree of freedom i, Ui (t), and its acceleration, Üi (t), are given by:

Ui (t) = U i (t + ∆t) -U i (t -∆t) 2∆t (27) Üi (t) = U i (t + ∆t) -2U i (t) + U i (t -∆t) ∆t 2 ( 28
)
where ∆t is the time step. Consequently, the problem ( 26) is reduced to the following system:

U i (t + ∆t) = ∆t 2 M d ii F ext i (t) -F int i (t) + 2U i (t) -U i (t -∆t) (29)
for each degree of freedom i.

It is well-known that the stability of this scheme depends on the time step ∆t. ∆t must be less than the critical time step ∆t cr . ∆t cr is the smallest time step required for a wave to cross an element. As the present system is highly non-linear (large strains and material behaviour), stiffness properties change during calculations. To ensure the convergence of the scheme, ∆t cr is periodically estimated by determining the stiffest element in the whole mesh.

Constitutive equation implementation

In order to calculate the internal force term in Eq. ( 29), the relation between deformations and stresses must be computed in each element. Using the plane stress state assumption, deformation and stress matrices in a single triangular element have the following forms:

[C] =   C 11 C 12 0 C 12 C 22 0 0 0 C 33   (30)
and:

[S] =   S 11 S 12 0 S 12 S 22 0 0 0 0   (31)
where 1, 2, 3 represent the local axis in the element plane. The direction 3 is the normal to the surface.

Using nodal initial positions and displacements in the global referential axis, C 11 , C 12 and C 22 are easily computed in each element. From the incompressibility constraint detC = 1, we have:

C 33 = 1 C 11 C 22 -C 2 12 (32)
By definition, C 33 is the square of the principal stretch ratio in the membrane thickness direction:

C 33 = h h 0 ( 33
)
where h and h 0 are respectively the deformed and undeformed thicknesses of the membrane. By now, all terms of the strain tensor can be computed in each element. The components of A in element local axis are computed using the fibers direction a 0 previously defined in each undeformed element:

[A] =   a 2 01 a 01 a 02 0 a 01 a 02 a 2 02 0 0 0 0   (34)
All terms corresponding to the third direction are null because the fibers always stay in the element plane. For a given expression of the strain energy function, S can now be calculated in the element local frame. The hydrostatic pressure p is determined by using:

S 33 = 0 (35)
In order to illustrate this general method, a very simple strain energy function is chosen:

W = Φ(I 1 -3) + Ψ(I 2 -3) + ΘI 4 (36)
where Φ, Ψ and Θ are the material parameters. This constitutive equation is a tranversely isotropic generalization of the widely used Mooney-Rivlin model [START_REF] Treloar | The mechanics of rubber elasticity[END_REF]. Using Eqs ( 18), (32), ( 35) and (36), the hydrostatic pressure p is given by:

p = 2 Φ + Ψ(C 11 + C 22 ) C 11 C 22 -C 2 12 ( 37
)
And the stress components becomes:

S 11 = 2(-C 22 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(C 22 + C 33 ) + ΘA 11 ) (38) S 22 = 2(-C 11 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(C 11 + C 33 ) + ΘA 22 ) ( 39
)
S 12 = 2(C 12 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(-C 12 ) + ΘA 12 ) (40)
where C 33 is given by Eq. (32).

NUMERICAL EXAMPLES

Validation: free inflation

In order to validate the implementation of the previous constitutive equation, our results are compared with those obtained by Kyriacou et al. [START_REF] Kyriacou | Finite element analysis of nonlinear orthotropic hyperelastic membranes[END_REF]. They use a quasi-static approach with four nodes elements to solve membrane inflation problems.

The free inflation of an initially flat square membrane submitted to an internal pressure is examined. The length of the undeformed membrane is denoted l 0 and its thickness h 0 . The inflating pressure is noted P . The material obey the modified Mooney-Rivlin model previously defined. To simplify the discussion the dimensionless inflating pressure P * is defined by:

P * = P l 0 Φh 0 (41)
Φ is the first material parameter. The square membrane is first biaxially pre-stretched in its undeformed plane until the deformed square length l reaches 1.1 l 0 . Kyriacou et al. use this pre-stretching to avoid any singularity in the tangent stiffness matrix at the first loading step. It is to note that no pre-streching is needed to start the computation in our dynamical approach. Numerical values are set to:

Φ = 1. Ψ = 0. Θ = 1. H = 10 -3 L = 1. ( 42
)
In this example, the square membrane is meshed with 400 triangular elements. Kyriacou et al. use 100 four nodes finite elements. Moreover, when the fibers direction coincides with one of the global axis directions, i.e.: a 0 = X or a 0 = Y, only a quarter of the membrane is meshed using symmetry properties. First, the inflating pressure P * is increased from 1.0 to 9.5. Figure 1 shows the height of the center node of the square sheet versus the inflating pressure. At the maximum pressure, the height is equal to 0.52. Our results are identical to the previously published data. Figure 2 shows the height of the center node at P * = 6. for different fibers directions. In the figure, a relative angle between one side of the square and the fibers direction equal to zero corresponds to a 0 = X. The curve has a sinus-like shape and is symmetrical with respect to the 45 o angle because of the geometrical symmetry of the problem. The membrane is less stiff when fibers are aligned with the diagonal direction of the square. The relative difference in height between our curve and Kyriacou's one does not exceed 0.03%.

It is well-known that the pressure-displacement curve corresponding to the free inflation of a rubberlike membrane exhibits limit points [START_REF] Khayat | Inflation of an elastic cylindrical membrane: non-linear deformation and instability[END_REF]. In order to be able to reach higher deformation ratios, the load is imposed using a constant gaz flow rate. The pressure at time t is computed step by step by using the Perfect Gaz Law and the volume inside the membrane at t [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF]. Using this continuation technique, the limit point can be overcomed as shown in Figure 3. Thus the membrane reaches higher deformations.

The deformed geometry of a transversely isotropic inflated square membrane for two different fibers directions is presented in Fig. 4. The lines represent the fibers direction. The classical geometrical symmetry of the isotropic case disappears: membranes tend to exhibit a bulge perpendicularly to the fibers direction.

Finally, this example demonstrates that our results show good agreement with those of Kyriacou et al. Moreover, by the use of a gaz flow rate loading instead of a pressure loading, our method is able to go beyond limit points.

Thermoforming application

In the present example, some results concerning a simple thermoforming application are examined. The inflation of an initial circular flat sheet in a cylindrical mold is considered. Here automatic remeshing is used in order to be able to fill the corners of the mold.

As mentioned elsewhere [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF], one of the main goal of the numerical simulation of such processes is to predict the final thickness distribution of molded parts. In order to evaluate the influence of fibers-reinforcement on the final geometry, comparisons are made between the membrane thickness distributions in the isotropic and the transversely isotropic cases. Material parameters are set to:

Φ = 1. Ψ = 0. Θ = 3. ( 43
)
Figures 5 and6 present respectively four inflation stage of the moulding of a cylinder in the isotropic and the anisotropic cases. Due to geometrical and material symmetries, only a quarter of the membrane is considered. In the isotropic case, the thickness distribution is axially symmetric (see Fig. 5). Figure 6 highlights the loss of symmetry due to the presence of fibers. The fibers direction is represented by the lines on the parison. It is to note that in this case axisymmetrical calculations would be unable to predict the final thickness distribution. In the transversely isotropic case, fibers are oriented as in Figure 6 and distributions differ for the two lines. The results corresponding to lines XZ and Y Z are respectively plotted with a pointed line and a dashed line. In both cases, isotropic and transversely isotropic (for the two lines), the highest thickness is observed near the clamped edge and the smaller thickness corresponds to the corner of the mold. An important remark found in the bibliography [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF] relative to the thermoforming of isotropic material concerns the membrane thickness distribution at the end of the process. It is said that this distribution does not depend on the material model. But as pointed out in the present example, the anisotropic material behaviour due to the presence of fibers highly influence thickness distribution in thermoformed parts.

CONCLUDING REMARKS

In this paper a general form of the strain energy function for transversely isotropic hyperelastic materials has been developed. This function depends on a limited set of tensor invariants and ensure the objectivity principle. Moreover as shown by Kyriacou et al., the formulation of the constitutive equation is easy implemented in a quasi-static finite element program. Here it is shown that the implementation in a dynamic explicit finite element code is easier. Efficiency of the method is demonstrated on both free and confined inflation problems.

As a perspective, we are working on a viscoelastic counterpart of the present hyperelastic model and on models for misaligned thermoplastic composites using a more general orientation tensor. In order to determinate the material parameters, a biaxial rheometer for high strain rates is in development.

Figure 1 :Figure 2 :

 12 Figure 1: Height of the center node versus inflating pressure for a square membrane

Figure 3 :Figure 4 :

 34 Figure 3: Height of the center node versus inflating pressure for a square membrane loaded with a constant gaz flow rate

Figure 5 :

 5 Figure 5: Thermoforming stages of a circular plane sheet in a cylindrical box: isotropic case. Gray scale: relative membrane thickness

Figure 6 :

 6 Figure 6: Thermoforming stages of a circular plane sheet in a cylindrical box: transversely isotropic case. Gray scale: relative membrane thickness

Figure 7 Figure 7 :

 77 Figure 7 presents the membrane thickness evolution along two edges shown in the small figure in the corner. The results correspond to the previous simulations. In the