AKNOWLEDGEMENT

I would like to extend my sincere gratitude to my supervisors, Franc ¸ois PEYRAUT, for his instructive advice, useful suggestions and his patience to help with the revision and improvement on my thesis. I am also grateful to my co-director, Zhi-Qiang FENG, for providing me with valuable advice on my thesis and helping me in the finite element implementation. At the same time, I would like to express my gratitude to Fr éd éric HOLWECK helped me to overcome mathematical problem and valuable comments on the writing of this thesis. I am deeply grateful of their help in the completion of this thesis. Without their invaluable help and generous encouragement, the present thesis would not have been accomplished.

I thank Jena JEONG and Erwan VERRON for accepting to be the reviewers of my thesis, as well as Sabine CANTOURNET and Julie DIANI for accepting to be the examiners of this thesis. It is also an honor for me that Julie DIANI has accepted to be the President of my thesis Jury.

Besides, I wish to thank my colleagues in University Bourgogne Franche-Comt é, UTBM.

Here I cannot list their names one by one to express my gratitude to thank them for giving me help and care. I am particularly grateful B éatrice ROSSEZ for her help in life, her smile and good humor. I thank my friends for their accompany and encourage during my PhD periods. I greatly appreciate my parents and boyfriend for their continuous support and endless love. My heart swells with gratitude to all the people who helped me. The People's Republic of China is acknowledged for its financial support through a grant on the China Scholarship Council (CSC) and the National Natural Science Foundation of China (No. 11372260).

We also warmly thank the Assistant Professor Kamenskiy (University of Nebraska Medical Center) to have kindly provided us the numerical data corresponding to the measurements included in [1]. v

CONTENTS INTRODUCTION

The goal of this thesis is to propose advanced finite element tools applied to the modeling of anisotropic hyperelastic materials. This kind of material, which are typically fiberreinforced rubbers, composites and soft biological tissues such as ligament, tendons or arterial wall, has very extensive applications in engineering science and health. For example, fiber-reinforced rubbers are used in manufacturing [2] and textile applications. Material deformation, instability, destruction and limited service life, and structure stability problem of finite deformation have become very important attributes of hyperelastic material. The behavior of anisotropic hyperelastic materials is also the keystone of scientists' research because the modeling of soft biological tissues has a wide range of applications in pharmaceuticals, therapeutic, medical prosthesis, ergonomics and so on. For instance, the basic problem in virtual surgery simulation is to obtain a realistic rendering of biological soft tissue behavior under real-time constraints including collision detection, interactive operation, visual rendering and tactile feedback [3, 4, 5]. An example of application of anisotropic hyperelastic materials in medicine is artificial heart. To design an artificial heart, scientists need to have precise simulations about a person's blood circle system, especially hydrodynamics of blood and solid mechanics of vessels and cardiac muscles in heart. Because anisotropic hyperelastic materials don't have universal properties and regular structures in all directions, we need to use finite element method to simulate their behaviors. To describe behaviors of those materials, we need to use constitutive equations. One of the most important constitutive equations is the stress-strain relationship which describes the mechanical properties of materials. The conventional anisotropic linear elasticity may be used to describe anisotropic hyperelastic materials under small deformations. However, the behavior of anisotropic hyperelastic materials exhibits nonlinear elasticity when they undergo large deformations [6, 7]. In this situation their stress responses can be derived from a given strain energy function (SEF) leading to highly nonlinear problems in structural mechanics.

The properties of anisotropic hyperelastic materials are directionally dependent. Unlike isotropic materials that have material properties identical in all directions, anisotropic material's properties change with directions. In the literature, it is widely accepted that anisotropy is due to the collagen fibers and the ground substance, or matrix performance behaves isotropically [8]. The main idea to study the properties of this type of material is to build specific strain energy functions which are invariant under a group of transformations in accordance with the symmetric properties of matter. Usually a structural tensor is introduced to account for the effect of the fiber directions of the anisotropic materials and to establish a link between anisotropy and isotropy [9]. In this way, anisotropic constitu-CONTENTS tive functions can be transformed into an isotropic one. Additionally, the structural tensors lead to invariant basis by using representation theorems for anisotropic tensor functions [10, 11]. The invariants of the deformation tensor and the additional structural tensors are necessary for constructing the strain energy function of anisotropic hyperelastic materials.

The theory of invariant has been used as early as 1950 to build the constitutive equations for isotropic and anisotropic materials [10, 12, 13, 14, 15]. By using invariant theory, we can deduce the form of the response function which is invariant under the considered material symmetry group. In other words, the material symmetry group is used to characterize the symmetries of isotropic and anisotropic materials. At the beginning, the researchers thought that material symmetry would impose certain restrictions on the forms of the response functions which appear in constitutive equations. Most form-invariance problems arising from specific constitutive equations and material symmetry groups have been solved by making use of the assumption that the response functions are polynomials [12, 16, 17, 18, 19, 20, 21, 22]. Pipkins et al. [23] expounded the restrictions of material symmetry on non-polynomial constitutive equations and proved that the polynomial assumption is not essential. They demonstrated that the form of the response function cannot be subjected to any a priori restrictions with a finite group of symmetries. Then, Wineman et al. [13] extended this theorem that there is no a priori restriction of any kind on the response function. That leads to a large variety of proposals in the literature to construct strain energy function, such as logarithmic, exponential, polynomial or power forms [24, 25, 26, 27].

Up to now, several strain energy functions have been presented for anisotropic hyperelastic materials with one or several fiber families. One of the first model for representing collagenous soft tissues behavior is based on a structural tensor and on a multiplicative decomposition of the deformation gradient tensor and was proposed by Weiss et al. [28]. Zulliger et al. [29] proposed a SEF for arteries that account for the wall composition and structure. A transversely isotropic viscohyperelastic constitutive law was suggested in [30] to describe the mechanical characteristics of the human anterior ligament. This method depends on the right Cauchy-Green deformation tensor. Pioletti et al. [31] introduced a constitutive law for the human cruciate ligament and patellar tendons. It can precisely fit the non-linear stress-strain curves at different strain rate. For fiber-reinforced materials, Qiu et al. [32] proposed a standard reinforcing model and Merodio et al. [33] also have presented a simple reinforcing model which takes into account the influence of reinforced fiber on the shear response. Although Qiu and Merodio's models are simple, they exhibit monotonic behavior during extension in the fiber direction and non-monotonic behavior during compression in the fiber direction. Finally, Guo et al. [34] proposed another reinforcing model based on the multiplicative decomposition of the deformation gradient tensor introduced earlier in [27]. They found that using neo-Hookean material to model the fiber can express the monotonicity in the stress-strain response during compression in the fiber direction.

On the other hand, research on experimental aspects of mechanical behaviors of anisotropic hyperelastic materials combined with theories have also been conducted by several scientists. For experiments carried out to study the mechanical behavior of arterial walls, living cells, and organs, see the review of Ar é et al. [35]. Ohashi et al. [36] proposed a pipette aspiration technique for characterizing nonlinear and anisotropic mechanical properties of blood vessel wall. This technique can eliminate the difficulty in the specimen mounting on the experimental apparatus and measure the homogeneity and heterogeneity of the specimen. Davarani et al. [37] have offered an interesting insight to discuss the optimal number of fiber families for describing the collagen behavior. Groves et al. [38] used tensile tests on circular skin (human skin and murine skin) and choose the transversely isotropic hyperelastic constitutive model of Weiss et al. [28] previously described. Tests of mechanical behaviors of two different fiber-reinforced rubbers with a one fiber family have been done by Ciarletta et al. [39]. They used a non classical measure of the strain to build two different types of strain energy functions for both tensile and shear deformation. In the same vein, Fereidoonnezhad et al. have built later a model using this kind of strain for two types of rubbers [40]. They also investigated the torsion of a circular fiber-reinforced rubber.

At this stage, we must mention that the notion of polyconvexity, originally introduced by Ball [41], constitutes a key issue for discussing the existence of solutions of hyperelastic problems. Holzapfel et al. [26] used a polyconvex exponential function for the description of the strain energy stored in the collagen fibers. The biomechanical behavior of the arterial wall and its numerical characterization have been analyzed and discussed by using this model. Later, by using again an exponential form function, Holzapfel, Ogden and Gasser introduced the so-called HGO model [42]. This model includes a Neo-Hookean density for determining the isotropic response. This Neo-Hookean density is actually one of the most commonly used to describe the mechanical behavior of the isotropic ground substance and also constitutes a simple polyconvex function. Based on the HGO model, the mechanical behaviors of different anisotropic materials have been studied in the literature [8, 43, 44, 45]. Additionally, a very large range of polyconvex functions was studied by Schr öder et al. [46]. However, the development of accurate computational models of anisotropic hyperelastic materials is challenging because of the difficulty for choosing an appropriate SEF among the numerous proposed in the literature and of the fact that a few of them are implemented in finite element codes. It is yet noticed that a finite element model was presented in [47] for three-dimensional (3-D) nonlinear analysis of soft hydrated tissues such as articular cartilage in diarthrodial joints under physiologically relevant loading conditions. A sample problem of unconfined compression is used to further validate the finite element implementation. Weiss et al have described a three-dimensional constitutive model for biological soft tissues and its finite element implementation for fully incompressible material behaviors [28]. The well known HGO model is available through the commercial code ABAQUS and was also implemented in the university code FER by Peyraut et al. [48]. Some techniques were described in [49] which can facilitate the construction, analysis CONTENTS and validation of FE models of ligaments. Finally, we can mention some problems often met with the modeling of hyperelastic material:

• If the existence and uniqueness of solutions are ensured by assuming the convexity of the SEFs, the most popular densities are not convex. This fact requires to account for polyconvexity which is not easy to do everytime.

• Some invariants associated with the SEFs are uneasy to be physically interpreted.

• If the SEFs are not elliptic or convex, some numerical problems can arise [50].

• In order to account for the shear interaction between the matrix and the fibers and for the normal interaction between fibers, sophisticated models [2, 51, 52] combine up to 4 SEFs, 9 material parameters and 11 invariants. That situation leads to handle complex and difficult-to-use models in order to be able of predicting the full behavior of soft biological tissues or reinforced-fibers rubbers.

Considering all the factors above, we choose a new set of invariants which was introduced by Ta et al. [53, 54]. This is an original approach mixing the isotropic and the anisotropic parts in a single SEF (most of the papers published in the literature propose to separate the energy density into an isotropic and an anisotropic parts) which was inspired by the pioneer work of Thionnet et al. [55]. This approach is mathematically justified by the theory of invariant polynomials, particularly by the Noether's theorem and the Reynolds operator [56]. This constructive approach allows to build an integrity basis generating all the polynomial invariants related to a specific anisotropic material. In this way, it would be possible to reduce drastically the number of invariants, material parameters and SEFs required to simulate the behavior of the material. The approach based on the polynomial invariant theory therefore constitutes a significant improvement for decreasing the complexity of models. However, up to now, and to the best of our knowledge, the mathematical foundations introduced in [53, 54] have not met a practical extension. So the first target of our work is to propose some new strain energy functions for anisotropic hyperelastic materials with different fiber structures by using the integrity basis made of the new invariants exhibited in [53, 54]. The guideline of our proposal is to combine appropriately these new invariants in order to provide a polyconvex property and a physical meaning. Following this strategy lead us to build two different SEFs, one representing the behavior of a one-fiber family material while the other is dedicated to a four-fibers family material. To confirm the accuracy and practicability of our models, the predicted results are compared with experimental ones extracted from the papers published by Ciarletta et al. [39] and Kamenskiy et al. [1] for a one-fiber family and a four-fibers family materials, respectively.

To use our proposed theoretical models into practical situations, the second target of this thesis is to implement our proposed SEFs in a finite element code. This work was conducted in close partnership with the Laboratory of Mechanics of the University of Evry (France) by using the FER university software [57]. Some sample problems are used to further validate the finite element implementation. The proposed SEFs and the finite element codes can therefore be applied for understanding the nature of behavior laws for materials with different fiber family structures and various loading cases with homogeneous or inhomogeneous deformations. This ensures that our model is effective and efficient when we need to use it for pragmatic applications.

This Phd dissertation is divided into three chapters:

The first one mainly introduces the foundations of continuum mechanics, the strain energy functions used for modeling the behavior of isotropic or anisotropic hyperelastic materials, the notion of polyconvex strain energy functions and the total Lagrangian formulation used for the finite element implementation in the context of nonlinear structural studies.

As polyconvexity is often considered as a prerequisite for ensuring the existence of solutions compatible with physical requirements [41], we provide a summary of common polyconvex functions. We finally remind some standard results as the Rivlin-Ericksen representation theorem used for isotropic hyperelastic materials [58] and we also recall the fact that an anisotropic strain energy function can be replaced by an isotropic one by including in the model an additional structural tensor [9]. Many authors used this method to build anisotropic energy functions [10, 11, 59, 60] and we briefly recall some of the most popular ones [2, 25, 44, 42, 61, 62, 6, 63, 64, 51].

The second chapter mainly consists in building a new polyconvex family of transverse anisotropic invariants and in designing a strain energy function (SEF) for incompressible fiber-reinforced materials. Only materials made with a one-fiber family are considered in this chapter. Unlike most papers published in the literature, that proposed to separate the energy density into an isotropic part and an anisotropic part, we introduced a new energy density mixing these two parts in a single function. As the invariants defined in [53] are well appropriate for this purpose, we use them and their polyconvexity and physical meaning are discussed extensively in this chapter. Several polynomial expressions were tested for the SEF but none is satisfying for properly describing the material behavior, particularly in the case of a shear testing. This is the reason why we have finally expressed the SEF by a combination of a polynomial with a power form function. In order to validate the usability and creativeness of the proposed model, two different fiber-reinforced rubber materials studied in [39] under uniaxial and shear testing are considered. In these two testing cases, the predicted results of our model show a fair agreement with experimental and predicted results extracted from [39] and from [40] which is a sequel of [39]. Finally, we provide all the details of the tensor calculation for the determination of the explicit expression of the tangent stiffness matrix. This matrix is needed for the finite element implementation of the model in the context of a total Lagrangian formulation. The practical implementation of our one-fiber family model was performed with C++ language in the university code FER [57]. In order to assess and check the validity of the FE implementation, several numerical simulations were successfully compared to experimental data extracted from the paper published by Ciarletta et al. [39]. We also use this model in 3D configuration to simulate more complex situations, namely an inhomogeneous extension loading. Finally, we have to mention that the research work presented in this second chapter has been published in an international journal for the SEF construction [65].

We introduce in the third chapter a new model for predicting the nonlinear mechanical properties of anisotropic hyperelastic materials with four-fibers families. The proposed strain energy function (SEF) can be applied for understanding the behavior of materials which have potential applications in biomechanics, surgical and interventional therapies for peripheral artery disease (PAD). For example ischemia, angina pectoris, myocardial infarction, stroke, or heart failure and other fatal diseases are consequences of atherosclerosis [66]. Like in our first model, this SEF is built with a recent and new invariant system based on the mathematical theory of invariant polynomials [53]. By recombining these invariants in an appropriate manner, we demonstrate that it is possible to build a polyconvex integrity basis. The next step is to associate properly this polyconvex invariants in order to build a SEF consistent with experimental data extracted from the literature. To reach this goal, a very simple expression, namely a linear polynomial form, appears to be efficient. Actually, based on this simple polynomial form, the corresponding model was validated by a comparison with experimental and numerical results extracted from [1]. These results concerned diseased superficial femoral (SFA), popliteal (PA) and tibial arteries (TA) from one patient under planar biaxial extension. For each kind of arteries tested with five combinations of different biaxial stretch, the predicted results of the proposed model and the experimental data are consistent. Our model includes seven material parameters and one additional advantage of the model is related to the parameters identification. Actually, the identification procedure provides a single solution because of the linear polynomial form of the SEF. Based on this energy function, a finite element program has been implemented inside the FER software in the same conditions as the ones described in the second chapter for the implementation of the one-fiber family model. The aspects related to the building of our SEF is included in a paper accepted for publication in an international journal [67] .

NOMENCLATURE AND ABBREVIATION

•

ON LINEAR ALGEBRA

We introduce in tgohis section some standard notations for tensor, matrix and vector calculus as well as some classical results on linear algebra. These notations and results will be extensively used further in the manuscript. 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, while we use lowercase letters for scalars. The standard Euclidean inner product is symbolized by a double bracket:

Aa, a = 3 i=1 3 j=1 A i j a j a i (1)
and the related Euclidean norm is noted . :

u = u, u 1 2
(2)

The tensor product between two vectors a and b is defined by:

(a ⊗ b) i j = a i b j (3)
If (a, b, c) forms an orthonormal vector basis:

a ⊗ a + b ⊗ b + c ⊗ c = I (4)
where I stands for the unity tensor.

The tensor product between two matrix A and B keeps the same notation as the one used between two vectors but is defined by:

( A ⊗ B) i jkl = A i j B kl (5)
The tensor notation stands for:

(A B) i jkl = 1 2 A ik B jl + A il B jk . (6)
The superscripts T and -1 respectively stand for the transpose and the inverse of a matrix:

(A T ) i j = A ji (7) CONTENTS A A -1 = A -1 A = I (8)
The cofactor matrix of an invertible matrix A is defined by:

Co f (A) = det(A)A -T (9)
where det symbolizes the determinant of a matrix.

We recall below a very standard result which is often used to describe the square of the stretch of the fibers for fiber-reinforced materials:

T r( A T Aa ⊗ a) = Aa, Aa = Aa 2 (10)
where Tr represents the trace of a matrix.

In equation [START_REF] Boehler | Applications of tensor functions in solid mechanics[END_REF] , A is frequently set to the gradient deformation matrix F and a often represents the direction of fibers.

The following equation is also a common property related to the trace operator:

T r(AB) = T r(B A) (11)
It is finally reminded the standard inner product operating on matrix:

A, B Fr = T r( AB T ) (12)
which induces the so-called Frobenius norm operating on matrix:

A Fr = (T r(A A T )) 1/2 (13)
The link between the Euclidian and Frobenius norms is classical and given by:

Aa ≤ A Fr a ∀A , ∀a (14)
CHAPTER 1

STATE OF THE ART 1 STATE OF THE ART 1.1/ INTRODUCTION
The first chapter of this Phd dissertation concerns the state of the art in the field of the modeling of anisotropic hyperelastic materials. This chapter is divided into six separate sections.

Section 1.1 introduces the fundamentals of continuum mechanics with emphasis on the field of deformations, strain and stress tensors.

Section 1.2 focuses on the principle of material frame-indifference and deduces from this principle the fundamental result that the second Piola-Kirchoff stress tensor can be expressed with respect to the right Cauchy-Green strain tensor instead of only the deformation gradient matrix.

Section 1.3 presents some classical and general results on the notions of isotropy and anisotropy while section 1.4 is related to more practical aspects with the description of some standard strain energy functions (SEFs). We also focus on the fact that some of these SEFs [2, 51, 52] include numerous material parameters and invariants in order to simulate coupled effects as the shear interaction between the matrix and the fibers or normal interaction between fibers. We finally present recent invariants proposed in [53, 54] because they form an integrity basis generating all the other invariants. Using them could therefore allow to drastically reduce the complexity of the above mentioned SEFs. We aim to use these recent invariants in the forthcoming chapters 2 and 3.

Section 1.5 gives an overview of the polyconvexity conditions ensuring the existence of solutions in compatibility with physical requirements [41]. The definitions of the convexity and polyconvexity are introduced in details and the relations between them are also discussed. The polyconvexity of the strain energy function is emphasized in the hyperelastic case. In addition, this section enumerates some commonly used polyconvex terms.

Section 1.6 presents the standard total Lagrangian formulation [68, 69] because we plan to use it for the finite element implementation of the two SEFs introduced in chapters 2 and 3.

1.2/ CONTINUUM MECHANICS

Continuum mechanics is the branch of mechanics dealing with the analysis of the kinematics and the mechanical behaviors of materials in terms of strain and stress. The purpose of continuum mechanics is to provide a macroscopic model for fluids, solids and organized structures. A fundamental assumption is the "continuous medium hypothesis": namely, that the real space occupied by a fluid or a solid can be approximately regarded as continuous without voids between the particles of matter.

1.2.1/ DEFORMATION AND STRAIN

The reference frame we introduce is made of rectangular coordinate axes at a fixed origin O with orthonormal basis vectors e a , a = 1, 2, 3 as shown in figure 1.1. Consider a continuum body Ω that moves in space from one instant of time to another. It occupies a continuum space denoted by Ω, . . . , Ω t . The region occupied by the continuum body B at the reference time is known as reference (or undeformed or Lagrangian) configuration. An initial region Ω at time t = 0 is referred as the initial configuration. A region Ω t at time t (t > 0) is referred as the current (or deformed or Eulerian) configuration. We set a typical point X (X ∈ Ω) occupied by a particle P ∈ B at the time t = 0. The particle P moves to the corresponding point x (x ∈ Ω t ) at the time t > 0. The map X = κ(P, t) is a one to one correspondence between a particle P ∈ B (see Figure 1.1) and its coordinates in the reference configuration. The map x = κ(P, t) acts on B to produce the region Ω t at time t.

The relation between the coordinates x in the current configuration and X in the reference configuration is described by:

x = κ κ -1 (X, t), t = φ(X, t) (1.1)
The motion φ is suitably regular and carries points X located in Ω to x in the current configuration Ω t . In terms of components, the vectors X and x can be described as:

X = 3 i=1 X i e i x = 3 i=1 x i e i (1.2)
The deformation (or displacement) vector field in the Lagrangian description is denoted by

U: U(X, t) = x(X, t) -X (1.3)
Replacing equation (1.1) into equation (1.3) , the deformation can be represented as: The deformation vector field is called u in the current configuration:

U(X, t) = φ(X, t) -X (1.4)
u(x, t) = x -X(x, t) = x -φ -1 (x, t) (1.5)
The two descriptions can be related by φ, namely:

U(X, t) = U(φ -1 (x, t), t) = u(x, t) (1.6)
From the map φ, we introduce the classical deformation gradient tensor F:

F = ∂φ(X, t) ∂X = Gradx(X, t) (1.7)
where Gradx(X, t) represents the gradient operator applied to the map φ(X, t).

It is well known that the determinant of F, commonly noted J, represents the volume change between the reference and the current configurations:

J = det(F) = dV dV 0 > 0 (1.8)
This determinant can also be interpreted as the determinant of the Jacobian matrix between the two configurations:

dx = F(X, t)dX (1.9)
Replacing equation (1.4) into the equation (1.7) leads to

F = ∂(U + X) ∂X = I + ∂U ∂X (1.10)
The classical right Cauchy-Green deformation tensor C is an important measure of strain in continuum mechanics. It can be defined from F by:

C = F T F (1.11)
Sometimes, the left Cauchy-Green strain tensor B, which plays a similar role as C, is also used:

B = FF T (1.12)
Finally, the Green-Lagrange strain tensor E is defined by:

E = 1 2 (C -I) (1.13)
The use of E is popular because it corresponds to a zero strain for a material at rest

(U = 0 ⇒ F = I ⇒ C = I ⇒ C -I = 0).
Additionally, by a linearization of E in the framework of a small displacement assumption, it is found that E is equivalent to the symmetric part of the gradient of the displacement:

E ≈ 1 2 (gradU + grad T U) (1.14)
Remark 1.1 . C, B and E are symmetric matrices.

We finally recall the relationship between two outward unit vectors N and n related to two infinitesimally small areas dS and ds in the undeformed and deformed configurations respectively (figure 1.2):

nds = JF -T NdS (1.15)
Equation (1.15) , also known as the Nanson's formula, can be reformulated by:

n = F -T N F -T N (1.16)

1.2.2/ STRESS TENSORS

There are three different kind of stress measures that are widely used in the framework of nonlinear continuum mechanics, namely the Cauchy, the first Piola-Kirchhoff and the second Piola-Kirchhoff stress tensors. If the assumptions of small displacement and small strain are considered (which is not the case in our framework), all of these three tensors are equal. The use of these three tensors is very common [70] and the two first allow to calculate the infinitesimal resultant force d f acting on a surface element as described on figure 1.2: The Cauchy stress tensor σ is often simply called the true stress while the first Piola-Kirchhoff stress tensor P is known as the engineering stress. Note that the transpose of P is frequently called the nominal stress tensor.

d f = σnds = PNdS (1.17)
Reporting the Nanson's formula of equation (1.15) in the equation (1.17) yields to:

P = JσF -T (1.18)
Or equivalently:

σ = J -1 PF T (1.19)
Note from equation (1.18) that, even if the Cauchy stress tensor is symmetric (σ T =σ), the first Piola-Kirchhoff stress tensor P is not. A symmetrization of P from equation (1.18) leads to the definition of the second Piola-Kirchhoff stress tensor S:

S = JF -1 σF -T (1.20)
In the framework of continuous thermo-elasticity, S is assumed to be derived from a po-tential W frequently called the strain energy density (SEF):

S = 2 ∂W ∂C (1.21)
It is noticed that this derivative can be calculated indifferently with respect to C and E by using the equation (1.13) :

S = 2 ∂W ∂C = 2 ∂W ∂E ∂E ∂C (1.22)
It should also be noticed that W contains all the informations related to the organization of the matter (number of fibers families, direction of the fibers etc.). These informations are expressed through the dependence of W with respect to physical invariants and their related material parameters. The choice of these invariants and the best way to combine them for building W, and then computing S, σ and P, constitutes a key point for the modeling of fibers-reinforced materials. Some standard choices for W will be presented in the forthcoming section 1.5. However, before introducing these standard models, we have to recall the classical properties that W must satisfy, for example the material frame indifference principle.

1.3/ MATERIAL FRAME INDIFFERENCE

Material frame indifference requires the invariance of the constitutive equation under different configurations. That means that, whatever is the selected basis for the evaluation of physical quantities, these quantities must remain invariant. In terms of stress, this concept of material frame indifference leads to [70]:

σ(QF) = Qσ(F)Q T ∀Q ∈ S O(3) (1.23)
where σ represents the Cauchy stress tensor in a given basis, says B, while σ is the Cauchy stress tensor in the original basis B. B is deduced from B by applying a positive orthogonal transformation Q, that is to say a rotation. It should be noticed that we have adopted here the classical and intuitive assumption stating that the Cauchy stress σ depends on the deformation gradient matrix F. SO(3) represents the positive orthogonal group, that is say the set of any matrix satisfying:

QQ T = Q T Q = I ; det(Q) = 1 (1.24)
We deduce easily from equations (1.20) and (1.23) that the second Piola-Kirchhoff stress tensor S satisfies:

S(F) = Ŝ(QF) ∀Q ∈ S O(3) (1.25)
By using next a standard mathematic theorem regarding the polar decomposition of any matrix in a single product made of an orthogonal matrix with a positive definite symmetric one, it is possible to demonstrate [70] that the second Piola-Kirchhoff stress tensor S (which depends on F according to equation (1.25) ) can also be expressed only with C:

S(F) = S(C) (1.26)
1.4/ ISOTROPIC AND ANISOTROPIC MATERIALS 1.4.1/ ISOTROPIC MATERIAL Although we are interested in this PhD work by the behavior of anisotropic materials made of a single or several fibers-family, we focus in this section on isotropic materials for two reasons:

• Liu et al. [9] have established a deep link between isotropic and anisotropic properties through the introduction of a structural tensor representing the material symmetry group.

• Most of the hyperelastic anisotropic models proposed in the literature associate through a summation an isotropic density with an anisotropic one.

A material is said to be isotropic if its behavior is the same in all directions. In terms of stress, this property leads to the following invariant equalities for the Cauchy stress and for the second Piola-Kirchhoff stress, respectively [70]:

σ(FQ) = σ(F) ∀Q ∈ S O(3) (1.27) S(FQ) = Q T S(F)Q ∀Q ∈ S O(3) (1.28)
Combining the equation (1.28) , corresponding to the isotropic property, with the materialframe indifference principle, represented by equation (1.26) , leads to the following lemma:

Lemma 1.1 . Consider an isotropic material and the material frame-indifference principle, the Piola-Kirchhoff stress tensor S(C) satisfies.

S(QCQ

T ) = Q S(C)Q T ∀Q ∈ S O(3) (1.29)
The proof can be found in [70].

We now introduce a very classical and famous theorem, known as the Rivlin-Ericksen representation theorem. It gives an explicit quadratic relation between S and C: Theorem 1.1 . Let us consider an isotropic material satisfying the material frameindifference principle. Then the second Piola-Kirchhoff stress tensor S adopts a quadratic representation with respect to C:

S(F) = S(C) = B 0 (I 1 , I 2 , I 3 )I + B 1 (I 1 , I 2 , I 3 )C + B 2 (I 1 , I 2 , I 3 )C 2 (1.30)
where B 0 , B 1 and B 2 are scalar-valued functions of the three principal invariants of the tensor C:

I 1 = tr(C) I 2 = 1 2 [tr(C) 2 -tr(C 2 )] I 3 = det(C) (1.31)
The proof of the theorem can be found in [70]. The main steps are a spectral analysis of C in terms of eigenvalues and eigenvectors, an appropriate choice of a change of basis matrix P for diagonalizing the symmetric matrix C, a spectral decomposition of S, similar to the one used for C and, finally, the resolution of a Vandermonde system in order to determine the matrix terms of the spectral decomposition. Note that the occurrence of a double or a triple coalescence of the eigenvalues of C leads to two particular results: Double coalescence (for example:λ 1 = λ 2 λ 3 ):

S(C) = B 0 I + B 1 C (1.32) Triple coalescence (λ 1 = λ 2 = λ 3 ): S(C) = B 0 I (1.33)
1.4.2/ ANISOTROPIC MATERIALS Anisotropic materials are materials whose properties are directionally dependent. In the most complicated situation, the matter is randomly distributed and there are no specific directions to characterize the organization of matter. This case is for example encountered with abradable materials used for aeronautics application [71]. However, many anisotropic materials can be characterized by specific directions or planes through unit vectors m 1 , ..., m a or tensors M 1 , ..., M b . This leads to introduce the subgroup g of O(3) (the full group of orthogonal transformations), which is called the symmetry group of the material, and which is defined by:

g = {Q ∈ O(3), Qm = m, QMQ T = M} (1.34)
The main property of this subgroup is to preserve the geometric characteristics of the material. In general, transversely isotropic materials, orthotropic materials and some classes of crystalline solids can be specified by a symmetry group of the type described by equation (1.34) [9]. That corresponds to our topic because we will study transversely isotropic materials in the chapter 2 and orthotropic materials in the chapter 3 of this manuscript. These two types of materials are considered as anisotropic due to the collagen fiber embedded in the matrix and are endowed with a natural symmetry structure:

• a material with one fiber-family of direction a is said to be a transversely isotropic material, see figure 1.3. The corresponding material symmetry group g is the set of all rotations around the fiber direction a:

g = {Q θ , ∀θ ∈ [0, 2π]} ; Q θ =              1 0 0 0 cos(θ) -sin(θ) 0 sin(θ) cos(θ)              (1.35)
where Q θ is the rotation of angle θ around the axis a. • a material with two fibers-family of directions a and b is said to be an orthotropic material. As illustrated on Figure 1.4, it is first observed that the planes P 1 and P 2 , respectively perpendicular to the directions e 1 and e 2 , which are the bisector of a and b and the co-bisector of a and b, form two planes of symmetry for the material. Moreover, since the fibers lie in the plane P 3 generated by e 1 and e 2 , it is obvious that P 3 is also a plane of symmetry for the material. The material properties remain therefore invariant under the action of the three reflections related to the three planes P 1 , P 2 and P 3 and the three rotations by an angle π around e 1 , e 2 and e 3 . The six orthogonal tensors related to these reflections and rotations are:

R(e 1 ) = -Q π (e 1 ) = I -2e 1 ⊗ e 1 R(e 2 ) = -Q π (e 2 ) = I -2e 2 ⊗ e 2 R(e 3 ) = -Q π (e 3 ) = I -2e 3 ⊗ e 3 (1.36)
The material also remains invariant under the action of I and -I. The material symmetry group denoted by S 8 therefore contains the 8 invariant matrix operators: The two material symmetry groups defined by equations (1.35) and (1.37) are included in the widest framework described by the following theorem:

S 8 = {T m , m = 1, ..., 8} (
Theorem 1.2 . : Let g be a subgroup of O(3) and S a matrix function depending on the symmetric matrix C. The anisotropic property given by equation (1.39) is therefore equivalent to the isotropic property given by equation (1.40) :

S(QCQ T ) = QS(C)Q T ∀Q ∈ g ⊂ O(3) (1.39) ⇐⇒ S(QCQ T , QPQ T ) = QS(C, P)Q T ∀Q ∈ O(3) ∀P ∈ N (1.40)
where the link between the anisotropic matrix S and the isotropic matrix S is:

S(C) = S(C, M) (1.41)
N is the set of matrix defined by:

N = {QMQ T , ∀Q ∈ O(3)} (1.42)
and M is the structural tensor satisfying:

M = QMQ T ∀Q ∈ g. (1.43)
The proof of this theorem is given in [9]. It means that the anisotropic property described by equation (1.39) can be reformulated in an isotropic form corresponding to equation (1.40) provided that a structural tensor M defined by equation (1.43) is introduced in the formulation. The link between the anisotropic one-argument function S and the corresponding isotropic two-arguments function S is given by equation (1.41) . The choice of M depends on which kind of anisotropy we are interested in. For the two kinds understudy here, the common choices are:

• For transversely isotropic materials with a one fiber-family of direction a (figure 1.3):

M = a ⊗ a (1.44)
• For orthotropic materials with a two fibers-family:

M 1 = e 1 ⊗ e 1 , M 2 = e 2 ⊗ e 2 , M 3 = e 3 ⊗ e 3 (1.45)
where e 1 , e 2 and e 3 are the three unit perpendicular vectors shown on the figure 1.4.

At this stage, it should be underlined that the choice corresponding to equations (1.44) and (1.45) is not unique. It depends in fact on the scalar valued function used to build the strain energy function W. If the standard mixed invariants J 4 = T r(Ca ⊗ a) is for example selected to build W, that naturally leads to equation (1.44) by deriving W with respect to C. It is therefore mandatory to discuss first the choices of the scalar invariants and the subsequent strain energy function (SEF).

1.5/ COMMON STRAIN ENERGY FUNCTIONS

Usually, an hyperelastic material behavior law is a type of constitutive model for ideally elastic material for which the stress-strain relationship derives from a strain energy density function W [70]. This strain energy density function depends on F, or equivalently on C, but, practically, this dependence is expressed through invariants scalar valued functions related to C. Some of these invariants are considered as classical and we provide in the following a brief introduction to the most classical strain energy functions using them.

1.5.1/ SEFS FOR ISOTROPIC HYPERELASTIC MATERIAL

In this section, we present one of the most renowned isotropic SEFs because a standard way to model anisotropic hyperelastic materials is to combine an isotropic SEF with a fully anisotropic one. To do that, we first need to introduce the classical isotropic invariants through the following theorem:

Theorem 1.3 . : If the material frame-indifference is assumed, the strain energy density W related to an hyperelastic isotropic material only depends on the three principal invariants I 1 , I 2 and I 3 of the matrix C introduced by equation (1.31) :

W(C) = W(I 1 , I 2 , I 3 ) (1.46)
CHAPTER 1. STATE OF THE ART Additionally, the second Piola-Kirchhoff stress tensor is given by:

S = 2 ( ∂W ∂I 1 + I 1 ∂W ∂I 2 + I 2 ∂W ∂I 3 )I -( ∂W ∂I 2 + I 1 ∂W ∂I 3 )C + ∂W ∂I 3 C 2 (1.47)
The proof of this theorem is very easy to establish by using first an appropriate choice of a change of basis diagonalizing C, by employing next the derivative chain rule with equation (1.21) and by applying finally the well known Cayley-Hamilton theorem.

Because the Mooney-Rivlin model is one of the most popular isotropic density (used in conjunction with an anisotropic one) to describe the behavior of the ground substance surrounding the fibers, we now focus on this model. It consists in a linear combination of the two first main invariants of the right Cauchy-Green deformation tensor C. It was proposed by Mooney in 1940 [72] and expressed in terms of invariants by Rivlin in 1948 [62]. The authors gave to the strain energy function the following polynomial form:

W(I 1 , I 2 , I 3 ) = ∞ p,q=0 c pq ( I 1 -3) p ( I 2 -3) q + 1 d (J -1) 2 (1.48)
where the c pq are material constants related to the distortional response and d is a material constant related to the volumetric response. The reduced invariants I 1 , I 2 and I 3 are introduced instead of the principal ones defined by equation (1.31) by using the volume preserving part F of the deformation gradient matrix:

F = F vol F F vol = J 1/3 I F = J -1/3 F (1. 49
)

I 1 = I 1 J -2/3 I 2 = I 2 J -4/3 I 3 = 1 (1.50)
This separation of the dilational and volume-preserving parts of the deformation gradient matrix can overcome numerical difficulties, such like numerical ill-conditioning of the stiffness matrix [28]. It is noticed that:

det(F vol ) = det(J 1/3 I) = J det( F) = I (1.51)
So only F vol contributes to the volume change of the material and is called the dilational part of the deformation.

Remark 1.2 .

1.

The three constant coefficients -3, -3 and -1 contained in each bracket of equation (1.48) come from the fact that a material at rest must give a strain energy function equal to zero. Let us assume indeed that the material is at rest, that is to say a displacement field equal to zero:

U = O (1.52)
We therefore deduce from equations (1.10) and (1.11) that:

F = I + ∇U = I ⇒ C = F T F = I (1.53)
The three main invariants and the three reduced invariants are then evaluated from equations (1.31) and (1.50)

I 1 = 3 I 2 = 3 I 3 = 1 ; I 1 = 3 I 2 = 3 I 3 = 1 (1.54)
It is finally obvious from the definition (1.48) of W that the SEF is equal to zero:

W(3, 3, 1) = 0 (1.55)

2.

To enforce the incompressibility condition J = det(F) = 1, a penalty function is introduced at the end of equation (1.48) by:

W vol = 1 d (J -1) 2 (1.56)
where d = 2 k and k, the initial bulk modulus, is related to the Poisson ratio v by:

k = 2(C 10 + C 01 ) (1 -2v) (1.57)

3.

The first order Mooney-Rivlin model is reduced to the two first terms of the summation included in equation (1.48) :

W( I 1 , I 2 , I 3 ) = c 10 ( I 1 -3) + c 01 ( I 2 -3) (1.58)
If c 01 = 0, we recover the commonly used neo-Hookean [6] strain energy function:

W( I 1 ) = c 10 ( I 1 -3) (1.59)
It should be noticed that we have not examined in this section more isotropic densities than the Mooney-Rivlin one because we have chosen in our thesis work to mix in a single SEF the isotropic and anisotropic effects. This choice will be explained later with more details.

1.5.2/ SEFS FOR ANISOTROPIC HYPERELASTIC MATERIAL

As mentioned previously, there exists an extraordinary variety of anisotropic SEFs in the literature with a large scope of applications using biological soft tissues, textile tissues or reinforced rubbers. We aim first in this section to present some of these SEFs and to introduce next the very recent invariants proposed in [53, 54]. We actually wish to combine these invariants in order to produce two new original SEFs (see the forthcoming Some well-organized soft tissues, such as a single layer of the annulus fibrosus of the human intervertebral disc [27], are reinforced with collagen fiber embedded in the matrix ground substance. Their mechanical behaviors are considered as transversely isotropic [73]. To take into account this kind of anisotropy due to the fiber direction a, the structural tensor M defined by equation (1.44) is introduced. This leads to build the so-called mixed invariants I 4 and I 5 :

I 4 = T r(CM) = T r(Ca ⊗ a) I 5 = T r(C 2 M) = T r(C 2 a ⊗ a) (1.60)
It is well known that I 4 represents the square of the stretch in the fiber direction [8] and that I 5 is related to the shear fiber-matrix interaction through the coefficient I 5 -I 4 2 [27, 74].

The energy function W for transversely isotropic materials can be therefore represented as a function of five invariants [75]:

W(C, M) = W(I 1 , I 2 , I 3 , I 4 , I 5 ) (1.61)
Based on the invariants dependence described by equation (1.61) , it is quite easy to demonstrate the following representation theorem [76]: If we consider the case of orthotropic materials made of two-fibers family of directions a and b respectively, additional mixed invariants can be introduced:

Theorem 1.
I 4 a =T r(Ca ⊗ a); I 4 b = T r(Cb ⊗ b); I 6 = T r(Ca ⊗ b); I 5 a =T r(C 2 a ⊗ a); I 5 b = T r(C 2 b ⊗ b); I 7 = T r(C 2 a ⊗ b) (1.63)
The following does not aim to give an exhaustive review of the existing anisotropic SEFs but we will try to present some of the most popular ones. The Fung-type model has for example inspired a lot of work in the literature [77, 78, 79, 80] because Fung and his co-authors have done a lot of pioneer research works to understand and model the mechanical behavior of soft tissues [63, 81, 82, 83, 84]. One of their more attractive result [25] is to have demonstrated that an exponential form was preferable than a polynomial mathematical expression:

W = c(e Q -1) (1.64)
where c is a material parameter and the argument Q is defined by:

Q = c 1 E 2 11 + c 2 E 2 22 + 2c 3 E 11 E 22 + c 4 E 2 12 (1.65)
where the E i j are the components of the Green-Lagrange strain tensor and c 1 , c 2 , c 3 and c 4 are additional material parameters.

To model the multiaxial mechanical behavior of human saccular aneurysms, Seshaiyer et al. [78] ignored the in-plane shear effect and proposed a similar model but without the last part c 4 E 2 12 . Bass et al. [85] performed uniaxial and biaxial tests on specimens of annulus and expanded Fung-type model to describe the healthy human annulus fibrosus. The convexity of the Fung model was finally discussed in [86, 87] and we will examine the convexity issue from a more general point of view in the next section.

The exponential form appearing in equation (1.64) is very popular for modeling hyperelastic materials. For example, in the case of an isotropic behavior, Hart-Smith [88] has proposed the following SEF:

W = C 1 exp{C 2 (I 1 -3) 2 }dI 1 + C 3 ln( I 2 3 ) (1.66)
where C 1 , C 2 and C 3 are the material parameters of the model.

One other very popular model is the HGO model proposed by Holzapfel et al. in [42]. As the Fung-type model, the HGO model uses an exponential type function to describe the behavior of a one-fiber family soft biological tissue but by introducing the physically motivated invariant I 4 (see equation (1.60) ) and by adding the neo-Hookean isotropic density (for modeling the non-collagenous matrix of the media) to the anisotropic density:

W = W ani + W iso (1.67) W ani =            k 1 2k 2 {exp[k 2 (I 4 -1) 2 ] -1} f or I 4 ≥ 1 0 f or I 4 < 1
(1.68)

W iso = c 1 (I 1 -3) (1.69)
where c 1 , k 1 and k 2 are material parameters.

It must be underlined that the anisotropic part W ani of the strain energy function is case sensitive with the value of I 4 . If I 4 is lower than 1, meaning that the fibers are shortened in a compressible state, W ani is actually assumed to be null. It is indeed considered that a compressed fiber generates no stress [89].

The proof of the convexity of the HGO model with respect to the deformation gradient matrix F is given in Schr öder et al. [76] and its implementation in a university finite element code is detailed in [90]. Peyraut et al. also exhibited in [8] a closed form solution in the special case where uniaxial unconstrained tension loading is applied to a biological soft tissue modeled by the HGO density.

Many researchers use the HGO model to study the mechanical behavior of different soft tissues and composite materials. Cabrera et al. [91] have for example performed equibiaxial tensile tests on four adult ovine pulmonary artery walls and compared the experiment data to the predicted results obtained by the HGO model. Merodio and Goicolea [92] used the HGO model to study potential viscoelastic effect on soft biological material. Maher et al. [93] accounted for inelastic phenomena, such as softening and unrecoverable inelastic strains, into the construction of constitutive models to describe stress softening and permanent deformation in arteries tissue. They used the HGO model to represent the anisotropic components of the strain energy function. Bass et al. [85] indicated that there are few hyperelatic models that can accurately predict the mechanical behavior of biological soft tissue with uniaxial and biaxial tests at the same time. But all the abovementioned models consider that the mechanical contribution comes from the fiber and the matrix, but ignore the interaction between them. In the following, we will introduce some strain energy functions which account for this interaction effect.

Wu et al. [94] have for example proposed an hyperelastic model including the fibermatrix interaction for modeling annulus fibrosus in the case of a tension test along the circumferential direction. However, this model cannot be applied to the case of a uniaxial test along the axial direction. Gasser et al. [44] have proposed an extension of the HGO model by accounting for the distribution of the collagen fiber orientation in the arterial layer. This distribution and the fiber-matrix interaction were obtained by replacing in the HGO model formulation the standard fourth invariant I 4 by:

I 4 = kI 1 + (1 -3k)I 4 (1.70)
So the strain energy function W is written as:

W = 1 2 µ(I 1 -3) + k 1 2k 2 exp k 2 [kI 1 + (1 -3k)I 4 -1] 2 -1 (1.71)
where µ is the shear modulus of the ground substance or matrix, k 1 and k 2 are the material parameters representing the mechanical behavior of the collagen fibers and k is the material parameter which expresses the collagen fibers distribution. It is noted that if k = 0, the collagen fibers are ideally aligned and the HGO model is recovered.

Guo et al. [27] presented a composite-based hyperelastic constitutive density including shear effect to model the human annulus fibrosus. The deformation gradient is decomposed into a uniaxial contribution along the fiber direction and another one including shear effect. Guo et al. [73] verified the significance of the interaction between the fiber and the matrix for the human annulus fibrosus by analysing the experimental data obtained by Bass et al. [85] in the case of uniaxial and biaxial tests. Peng et al. [64] also proposed a strain energy function to describe the behavior of the annulus fibrosus by decomposing this function into three parts representing the contributions coming from the matrix W M , the fiber W F and the fiber-matrix shear interaction W F M , respectively:

W = W M + W F + W F M (1.72)
The energy W M related to the matrix combines the neo-Hookean model with a penalty contribution for accounting the incompressibility condition:

W M = C 10 ( I 1 -3) + 1 D 1 (J -1) 2 (1.73)
where C 10 and D 1 are the material parameters.

The strain energy function W F related to the fibers is a case-sensitive polynomial depending on the fiber stretch:

W F =        C 2 (I 4 -1) 2 + C 3 (I 4 -1) 4 I 4 > 1 0 I 4 ≤ 1 (1.74)
where C 2 and C 3 are material parameters.

Finally, the fiber-matrix shear interaction energy W F M can be expressed as:

W F M = f (I 4 , ϕ) = γ 1 + exp[-β(λ F -λ F * )] [ I 4 I 3 (I 5 -I 1 I 4 + I 2 ) -1] 2 (1.75)
where β, γ, λ F and λ F * are material parameters. The quantity λ F * may be related to the transition point between the toe region and the linear region in the uniaxial tensile stress-strain curves and the angle ϕ represents the shearing between the matrix and the fiber.

Hollingsworth et al. [52] considered that the combination of different effects coming from the proteoglycan matrix, the collagen fiber and the interaction between the constituents allows to model properly the annulus fibrosus. Contrary to Peng et al. [64], they assumed that the intra-lamellar fiber-fiber crosslinking dominates the interaction terms due to the relatively weak fiber-matrix interaction. The strain energy function is therefore assumed to be made of an isotropic proteoglycan matrix density (W m ), a primary collagen fiber families density (W f ), a shear interaction density W and a fiber-fiber interaction density

W ⊥ : W = W m + W f + W + W ⊥ (1.76)
With:

W m = a 1 (I 3 - 1 I 3 ) 2 + a 2 (I 1 I -1/3 3 -3) 2 (1.77) W f = a 3 b 3 (exp[b 3 (I 4 a + I 4 b -2)] -b 3 (I 4 a + I 4 b ) + 2b 3 -1) (1.78) W = a 4 (γ a 2 + γ b 2 ) (1.79) W ⊥ = a 5 b 5 (exp[b 3 (I 4 c + I 4 d -2)] -b 3 (I 4 c + I 4 d ) + 2b 5 -1) (1.80)
where a i (i = 1, ..

I 4 a = T r(Ca ⊗ a); I 4 b = T r(Cb ⊗ b); I 4 c = T r(Cc ⊗ c); I 4 d = T r(Cd ⊗ d) (1.81)
where a and b represent the fiber directions and c and d the directions perpendicular to a and b respectively in the plane contained by a and b. The introduction of c and d allows to account for the normal interaction, that is to say the ability to resist to the forces normal to the fiber directions. Finally γ a and γ b represent the shear strain along the fiber directions a and b, respectively:

γ a 2 = I 5 a -(I 4 a ) 2 ; γ b 2 = I 5 b -(I 4 b ) 2 (1.82)
In a similar spirit as Hollingsworth et al. [52], but with the objective of modeling a ground rubber reinforced by steel cords for manufacturing tire belt layers, Peng et al. [2] have proposed a strain energy function made of four contributions:

W = W M + W F + W shear + W normal (1.83)
The strain energy contribution W M of the isotropic rubber is defined by the Yeoh model [95]:

W M = 3 i=1 C i0 ( Î1 -3) i + 3 i=1 1 D i (J -1) 2i (1.84)
where Î1 is the first reduced invariant (see equation (1.50) ) and D i reflects the material compressibility.

The strain energy of the cords structure W F adopts the same form as the case sensitive polynomial introduced by equation (1.74) which was used for modeling the behavior of the annulus fibrosus.

The shear interaction strain energy W shear , caused by the angle variation between cord orientation and rubber normal, is defined as:

W shear = exp[c(I 4 -1)](aχ 2 + bχ) (1.85)
where a, b and c are material parameters and χ = tan 2 ϕ, ϕ being the angle variation between cord orientation and rubber normal.

Finally, in order to include the normal interaction in the model, a last strain energy W normal is introduced:

W normal = g(χ) k 1 k 2 [e -k 2 (I 6 -1) + k 2 (I 6 -1) -1] (1.86) g(χ) = exp(-lχ) + mχ + nχ 2 (1.87)
where k 1 , k 2 , l, m and n are material parameters and I 6 stands for the squared stretch in the direction perpendicular to the original cord orientation b:

I 6 = T r(Cb ⊗ b) (1.88)
The SEFs mentioned above use classical invariants directly related to the structural ten-sor through the trace. In the following of this section, we will introduce an example of strain energy function using non-classical invariants [51]. In this reference, the annulus fibrosus is actually considered as an isotropic ground substance reinforced by two families of collagen fibers and modeled by the following strain energy function with a single exponential function but a rather complicated argument:

W = 1 I 3 n βexp[β 1 (I 1 -3) + β 2 (I 2 -3) + β 3 (I 1 -3) 2 + β 4 (I 9 -2) + β 5 (I 9 -3) 2 + β 6 (I 1 -3)(I 9 -2) + β 7 (I 11 -2) + β 8 (I 10 -1) + β 9 (I 8 -cos 2 2φ) + β 10 (I 3 -3)(I 8 -cos 2 2φ) + β 11 (I 8 -cos 2 2φ)(I 9 -2)] (1.89)
where n and β i (i = 0, ..., 11) represent the 13 material parameters of the model, φ is the half angle between the two families of fibers, I 1 , I 2 and I 3 are the three classical invariants described by equation (1.31) , I 4 (resp. I 6 ) and I 5 (resp. I 7 ) are the classical mixed invariants related to the first fiber direction a (resp. the second fiber direction b) introduced by equation (1.60) and I 8 to I 11 are new mixed invariants combining the previous ones:

I 8 = (cos 2 2φ)T r(Ca ⊗ b), I 9 = I 4 + I 6 , I 10 = I 4 I 6 , I 11 = I 5 + I 7 (1.90)
By asking that there exists a stress free configuration corresponding to the material at rest, the authors reduce the number of independent material parameters from 13 to 11, which still remains a high number. Moreover, the authors having experienced some difficulties to identify these parameters, they placed the dependence of the invariants in separate exponentials instead of a single one as in the equation (1.89) :

W =α 0 {exp[α 1 (I 1 -3)] + exp[α 2 (I 2 -3)] + exp[α 3 (I 1 -3)] + exp[α 4 (I 9 -2)] + exp[α 5 (I 8 -cos 2 2φ)] + exp[α 6 (I 10 -1)] + exp[α 7 (I 11 -2)] + exp[α 8 (I 3 -3)(I 8 -cos 2 2φ)] + exp[α 9 (I 1 -3)(I 9 -2)] + exp[α 10 (I 8 -cos 2 2φ)(I 9 -2)]} (1.91)
where α i (i = 0, ..., 10) are the material parameters of this new energy density. By using again the fact that the configuration must be free of stress if the material is at rest, the authors reduced the number of material parameters from 11 to 9. In this way, they were able to fairly predict the mean response of the annulus fibrosus when measured in four different experimental deformations. However, it should be underlined that this model is rather complicated because it involves 9 material parameters and 11 invariants. That suggests that a rigorous mathematics approach, based on the polynomial invariant theory [56] rather than empirical considerations, could be an alternative in order to produce models containing all the relevant informations but no more. This is the main strategy followed by Ta and his co-authors in [53, 54]. In this way, they succeeded to prove that all the polynomial invariants are generated by only 6 of them (resp. 7) for a one (resp. two) fiber family material. Such sets of invariants having the property of generating all the others are called integrity basis. The objective of this thesis work is to build SEFs by using these basis. That will be done in the two forthcoming chapters by accounting for polyconvexity which is a prerequisite property for finding solutions in compatibility with physical requirements [41]. This property is reminded in the next section.

1.6/ POLYCONVEXITY

The mathematical treatment of structural boundary value problem in the framework of computer simulation is based on the calculus of variation, like finding a minimal deformation of the elastic free energy. This requires that the constitutive behavior law not only reflects the material properties, but also meets some convexity conditions as illustrated by figure 1.5 and corresponding to definitions 1 and 2:

Definition 1: convex set (Figure 1.6)
A set K is said to be convex if and only if the following relation holds: Let K be a convex set and f a scalar-valued function.

λx 1 + (1 -λ)x 2 ∈ K ∀x 1 , x 2 ∈ K ∀λ ∈ [0, 1] (1.92
Then f : K → R is said to be convex if and only if: Unfortunately, as highlighted by Ciarlet in [70], the convexity of SEFs must be ruled out because this property is incompatible with the physical insight that the energy must tend towards infinity if the volume of the material is reduced to zero. Additionally, the gradient deformation matrix F belongs to the set U of all the 3 × 3 real matrix defined by:

f (λx 1 + (1 -λ)x 2 ) ≤ λ f (x 1 ) + (1 -λ) f (x 2 ) ∀x 1 , x 2 ∈ K, ∀λ ∈ [0, 1] (1.93)
U = {M, det(M) > 0} (1.94)
because the determinant of F represents the positive volume change of the material.

Unfortunately, it is easy to demonstrate, by choosing an appropriate counter example [96], that U is not a convex set in the sense of definition 1 and figure 1.6. As a consequence, it does not make sense to include the invariant det(F) in the SEF expression if convexity is a mandatory requirement despite the fact that this invariant is widely used in the literature to describe the volume change. Observing these evidences, Ball [41] suggested to replace the convexity property by a weaker requirement called polyconvexity:

Definition 3: Polyconvexity

Let W: F → W(F) be a scalar-valued energy function. Then W is polyconvex if and only if there exists a function T :

M 3×3 × M 3×3 × R → R so that: W(F) = T (F, Co f (F), det(F)) (1.95)
and the function T : R 19 → R defined by:

(F, Co f (F), det(F)) → T (F, Co f (F), det(F))
is convex with respect to the set of arguments (F, Co f (F), det(F)).

By introducing this interesting property, it immediately results that a function of the form (det(F)) 2 is a convex real-valued function of det(F), that is to say a polyconvex function, even if it can not be a convex function of F as previously explained. It must be also noticed that, in the definition of polyconvexity, some authors like in [46] uses the adjugate of F, that is to say the transpose of the cofactor matrix of F, instead of Co f (F) but it does not constitute an essential difference. Note finally that, if convexity implies polyconvexity, the reverse is not true.

In the context of hyperelastic problems, the polyconvexity of the strain energy density is often considered as a prerequisite for ensuring the existence of solutions in compatibility with physical requirements as explained by Ball in [41]. Subsequently to the pioneer work of Ball, Marsden et al. [97] and Ciarlet [70] have provided more results related to this concept. For isotropic materials, some well-known strain energy functions, such as Ogden model [6], Mooney Rivilin model [62, 72] and Neo-Hookean model [6] are polyconvex. Schr öder et al. extended the concept of polyconvexity for anisotropic materials [98] and a wide survey with many proofs on polyconvexity of isotropic and transversely isotropic functions is given in [46].

In fact, the polyconvex and convex properties are deeply linked together and it is simply a matter of arguments to distinguish them (F on the one hand and (F, Co f (F), det(F)) on the other hand). From this point of view, it is therefore useful to recall a very practical condition to determine if a twice-differentiable function ψ of x is convex or not:

ψ (x)(dx)(dx) ≥ 0 ∀ dx (1.96)
To end this section, we will discuss the polyconvexity of the standard invariants I 1 , I 2 , I 3 , I 4 and I 5 introduced by equations (1.31) and (1.60) . It follows from equations (1.11) and (1.31) that:

I 1 = T r(F T F) (1.97)
We then use the property ( [START_REF] Spencer | Isotropic polynomial invariants and tensor functions[END_REF] ) on the trace operator and the definition [START_REF] As Wineman | Material symmetry restrictions on constitutive equations[END_REF] of the Frobenius norm to conclude that the invariant I 1 can be considered as a quadratic function on F:

I 1 = F, F Fr (1.98)
The second derivative of I 1 with respect to F is directly deduced from (1.98) :

I 1 (F)(H)(H) = 2 H, H Fr = 2 H 2 Fr (1.99)
That proves that I 1 is convex with respect to F and thus polyconvex.

To establish the polyconvexity of I 2 , we start from the Cayley-Hamilton theorem applied to the matrix C:

C 3 -I 1 C 2 + I 2 C -I 3 I = 0 (1.100)
Multiplying equation (1.100) by C -1 yields to :

C 2 -I 1 C + I 2 I = I 3 C -1 (1.101)
The trace operator applied to equation (1.101) thus gives:

I 2 = I 3 T r(C -1 ) (1.102)
We next remark that:

I 3 T r(C -1 ) = (det(F)) 2 T r(F -1 F -T ) = T r((Co f (F)) T Co f (F)) = T r(Co f (F)(Co f (F)) T )
(1.103) By using both equations (1.102) and (1.103) , as well as like the Frobenius scalar product [START_REF] Adkins | Symmetry relations for orthotropic and transversely isotropic materials[END_REF] , it is obtained that I 2 is a quadratic function of Co f (F):

I 2 (Co f (F)) = Co f (F), Co f (F) Fr (1.104)
The second derivative of I 2 with respect to Co f (F) is directly deduced from (1.104) :

I 2 (Co f (F))(H)(H) = 2 H, H Fr = 2 H Fr (1.105)
That proves that I 2 is convex with respect to Co f (F) and thus polyconvex.

The third invariant is obviously polyconvex because it is clear that I 3 is a quadratic function with respect to det(F):

I 3 (det(F)) = det(C) = (det(F)) 2 (1.106) =⇒ I 3 (det(F)) = 2 (1.107)
I 3 is consequently convex with respect to de f (F) and thus polyconvex.

Let us turn now our attention to the first mixed invariants I 4 defined by equation (1.60) . Thanks to equation [START_REF] Boehler | Applications of tensor functions in solid mechanics[END_REF] , I 4 can be written as a quadratic function of F:

I 4 (F) = Fa, Fa (1.108)
The second derivative of I 4 with respect to F is thus:

I 4 (F)(H)(H) = 2 Ha, Ha = 2 Ha 2 (1.109)
That proves that I 4 is convex with respect to F and thus polyconvex.

Finally, it has been demonstrated in [46] that I 5 is not polyconvex. Anyway, there are at least two standard ways to build polyconvex invariants containing I 5 . The first consists in using again the Cayley-Hamilton theorem described by equation (1.101) and by multiplying it by the structural tensor M defined by equation (1.44) :

C 2 M -I 1 CM + I 2 M = I 3 C -1 M (1.110)
Taking the trace of each member of equation (1.110) yields to:

I 5 -I 1 I 4 + I 2 = T r(I 3 C -1 M) (1.111)
Developping the right hand-side of equation (1.111) leads to:

T r(I 3 C -1 M) = T r((det(F)) 2 F -1 F -T a ⊗ a) = T r((Co f (F)) T Co f (F)a ⊗ a) (1.112)
Equation (1.112) can be simplified to a quadratic form with respect to Co f (F) by appliying the property [START_REF] Boehler | Applications of tensor functions in solid mechanics[END_REF] :

T r(I 3 C -1 M) = Co f (F)a, Co f (F)a (1.113)
The second derivative of T r(

I 3 C -1 M) with respect to Co f (F) is thus: T r(I 3 C -1 M) (Co f (F))(H)(H) = 2 Ha, Ha = 2 Ha 2 (1.114)
That proves that I 5 -I 1 I 4 + I 2 is convex with respect to Co f (F) and thus polyconvex.

The second way to obtain a polyconvex invariant depending on I 5 is to first multiply the Cayley-Hamilton theorem represented by equation (1.110) by I -M:

C 2 (I -M) -I 1 C(I -M) + I 2 (I -M) = I 3 C -1 (I -M) (1.115)
Taking the trace of equation (1.115) yields to:

I 1 I 4 -I 5 = T r(I 3 C -1 (I -M)) (1.116)
The right hand-side of equation (1.116) can be developed into:

T r(I 3 C -1 (I -M)) = T r(Co f (F) T Co f (F)) -T r(Co f (F) T Co f (F)a ⊗ a) (1.117)
By using [START_REF] Boehler | Applications of tensor functions in solid mechanics[END_REF] , [START_REF] Spencer | Isotropic polynomial invariants and tensor functions[END_REF] and [START_REF] Adkins | Symmetry relations for orthotropic and transversely isotropic materials[END_REF] , equation (1.117) can be reformulated in terms of two scalar products:

T r(I 3 C -1 (I -M)) = Co f (F) T , Co f (F) Fr -Co f (F)a, Co f (F)a (1.118)
That proves that I 1 I 4 -I 5 depends on Co f (F). To demonstrate that this dependence is convex, we derive twice I 1 I 4 -I 5 with respect to Co f (F):

(I 1 I 4 -I 5 ) (Co f (F))(H)(H) = 2{ H, H Fr -Ha, Ha } (1.119)
To prove that the right hand-side of equation (1.119) is positive, we use the inequality [START_REF] Rs Rivlin | Stress-deformation relations for isotropic materials[END_REF] and the fact that a is a unit vector:

Ha, Ha = Ha 2 ≤ H 2 Fr a 2 = H, H Fr (1.120)
The proof of the convexity of I 1 I 4 -I 5 with respect to Co f (F) is then complete and the polyconvexity of this invariant is demonstrated.

It is of course possible to combine the standard polyconvex invariants studied above in many different ways and an excellent overview on this issue is provided in [46].

1.7/ FINITE ELEMENT METHOD FOR STRUCTURAL NONLINEAR ANALYSIS

In the past few years, the role of computational modeling is becoming increasingly important in all fields of physics and particularly for biomechanics research. Computational analysis of the biomechanics of soft biological tissues provides a framework for quantitative description of biomedical material, which has a large potential of applications in medical science, biology simulation and robotics for real-time surgery simulation. As the focus of this thesis is not only to propose new strain energy functions for anisotropic materials, but also to implement them in a finite element code, we will describe below the standard total Lagrangian formulation used for this implementation. The geometrically nonlinear analysis may actually be described by using the total or the updated Lagrangian formulations [68, 69]. The total Lagrangian formulation is derived with respect to the initial configuration while the updated Lagrangian formulation is derived with respect to the current configuration. In other words, the total Lagrangian formulation constructs the tangent stiffness matrix with respect to the initial configuration. This simplifies the computation. Therefore, the total Lagrangian formulation was selected in this work for the finite element discretization. Using the symmetry of the strain tensor E and the stress tensor S (equations (1.11) and (1.21) ), we start by denoting 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 (1.121)
In the context of the finite element method and with equations (1.10) , (1.11) and (1.13) , the Green-Lagrange strain can be formally written with linear and nonlinear contributions in terms of nodal displacements u:

E = B L + 1 2 B NL (u) u (1.122)
where B L is the matrix which relates the linear part of the strain terms to the nodal displacements, and B NL (u), the matrix which relates the nonlinear strain terms to the nodal displacements. From equation (1.122) , the incremental form of the strain-displacement relationship is

δE = B L + B NL (u) δu . (1.123)
In static analysis of solids, the virtual work δU is

δU = V 0 δE T S dV 0 -δu T F ext = 0 (1.124)
where V 0 is the domain of the initial configuration and F ext the vector of external loads. In view of equations (1.21) and (1.123) , it comes

δS = DδE = D B L + B NL (u) δu (1.125)
where D denotes the matrix deduced from the fourth-order stress-strain tangent operator

D: D = ∂S ∂E = 2 ∂S ∂C = 2 ∂ 2 W ∂C 2 (1.126) D =                                 D1111 D1122 D1133 D1112 D1113 D1123 D1122 D2222 D2233 D2212 D2213 D2223 D1133 D2233 D3333 D3312 D3313 D3323 D1112 D2212 D3312 D1212 D1213 D1223 D1113 D2213 D3313 D1312 D1313 D1323 D1123 D2223 D3323 D2312 D2313 D2323                                 . (1.127)
Substituting δE from equation (1.123) into equation (1.124) results in:

δU = δu T V 0 B L + B NL (u) T S dV 0 -δu T F ext = 0. (1.128)
The vector of internal forces, appearing in the first term of equation (1.128) , is defined by:

F int = V 0 B L + B NL (u) T S dV 0 . (1.129)
Since δu is arbitrary, the following set of nonlinear equations is obtained:

F int -F ext = 0. (1.130)
The nonlinear equation (1.130) is solved numerically by using a classical Newton-Raphson scheme:

K i ∆u = F ext -F i int u i+1 = u i + ∆u (1.131)
where i and i + 1 refer to the current and to the next iterations. The displacement u i+1 is updated by the incremental nodal displacement ∆u and the tangent stiffness matrix K i is evaluated at each iteration by taking into account the internal forces vector F int . Deriving F int with respect to the nodal displacements u and using equation (1.125) yield to the tangent stiffness matrix:

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

K e = V 0 B T L DB L dV 0 (1.133) K σ = V 0 ∂B T NL ∂u S dV 0 (1.134) K u = V 0 B T L DB NL + B T NL DB L + B T NL DB NL dV 0 . (1.135)
The practical implementation of the two SEFs proposed in the two next chapters was performed in the university code FER [57] following the total Lagrangian approach described above.

1.8/ CONCLUSIONS

This first chapter was mainly devoted to introduce essential results needed for the two forthcoming chapters and also to put in perspective our work with the literature in order to better understand the following.

We have particularly explained why the selection of appropriate invariants is not an easy task in view of building a strain energy function. Actually, this selection often conducts to elaborate models involving many invariants, many material parameters and many densities to account for complex phenomena such as the fibers-fibers interaction or the shear interaction between the fibers and the matrix [2, 51, 52].

In order to construct in the next chapters the simplest models as possible, but capable of embedding all of the complex mechanical effects, we have explained in the end of section 1.5.2 why the invariants introduced by Ta et al. [53, 54] could be an interesting alternative to more standard invariants.

In order to select the best way for combining these invariants, we have finally explained in section 1.6 why the concept of polyconvexity, originally introduced by Ball [41] and developed later by Ciarlet [70], can serve as a guideline.

The forthcoming chapter 2 deals with the construction of a SEF for modeling a one-fiber family material. This second chapter will include a discussion on polyconvexity as well as the description of the finite implementation of the model by using the approach presented in section 1.7 of this chapter. The last chapter 3 of this manuscript follows the same logic of presentation as chapter 2 but this time for a four-fibers family material.

CHAPTER 2

A NEW SEF FOR ONE-FIBER FAMILY MATERIALS 2 A NEW SEF FOR ONE-FIBER FAMILY MATERIALS 2.1/ INTRODUCTION
The main goal of this chapter is to design a strain energy function (SEF) for incompressible fiber-reinforced materials by using the family of transverse anisotropic invariants proposed by Ta et al. [54]. These invariants derived from the application of Noether's theorem and we will organize them in order to form a polyconvex integrity basis. Based on this new constitutive model, we developed a finite element program in the FER software [57].

These past twenty years, many strain energy functions have been proposed for transversely isotropic materials to investigate the mechanical behavior of biological soft tissues. As mentioned in the first chapter, these materials are considered as anisotropic due to the collagen fiber behavior [44]. The number of fiber families is set to 1 to model tissues such as ligament, tendons or fiber-reinforced rubber materials, while it is set to 2 to represent the arterial wall [8, 65]. Several constitutive finite element models were built for biological soft tissues, such as ligament, tendons and the fiber-reinforced rubber materials [28, 47]. Shearer [99] built a new strain energy function for the hyperelastic modelling of ligaments and tendons based on the geometrical arrangement of their fibrils. Limbert et al. [100] proposed a phenomenological constitutive law to describe the anisotropic viscohyperelastic behaviour of the human posterior cruciate ligament (PCL) at high strain rates.

In general, it is assumed that the mechanical behavior of the material is not affected if the fibers are in a compressive state [101, 102]. Taking advantage of this situation, most of the papers published in the literature propose to separate the energy density into an isotropic part and an anisotropic part. The first part accounts for the low strain behavior of the ground matrix and the second part captures the behavior of the fibers at higher strain [28, 103]. More recently, an original approach mixing the isotropic and the anisotropic parts in a single SEF was introduced by Ta et al. [53, 54]. This approach was inspired by the pioneer work of Thionnet et al. [55] and is mathematically justified by the theory of invariant polynomials. It provides an alternative to the classical method found in the literature for building invariants and allows to exhibit an integrity basis made of six invariants (K 1 ,. . .,K 6 ), some of them being original, in the case of a one-fiber family material.

We have adopted the same approach as the one developed in [54] but with the following complementary results:

• one of the six invariants exhibited in [54] can be excluded from the integrity basis by adding the appropriate transformation in the material symmetry group;

• three of the six invariants are well known polyconvex functions;

• the last two invariants are original, physically motivated and directly connected to shear effects. Additionally, those two invariants shed a new light on the classical mixed invariant I 5 = T r(C 2 a ⊗ a) (where a represents the fiber direction) and allows to link it with shear strain while it is often reported in the literature the difficulty to provide a physically-based motivation for I 5 .

However, up to now and to the best of our knowledge, the mathematical foundations introduced in [54] have not met a practical extension. The new strain energy function proposed in this chapter by using the integrity basis made of five of the six invariants exhibited in [54] constitutes a first attempt in this direction. This choice is motivated by the fact that these invariants do not require a separation of the SEF into an isotropic and an anisotropic part. Another motivation is the rigorous mathematical foundations used by Ta et al. to define those invariants. In the same spirit as the Mooney-Rivlin models in the framework of isotropic hyperelasticity [72, 62], we have introduced some original SEFs as polynomial functions of these new invariants. The conclusions about those original behaviors laws are as follows:

• A linear or a quadratic expansion of the invariants is not sufficient to well describe the material behavior with the four experimental set-up considered, particularly with the shear test. In fact we prove that any polynomial SEF in the invariants will not be suitable to fit the experimental data.

• A quadratic expansion of the invariants combined with an appropriate power-law form provides accurate predictions of all the experimental results.

In most situations, the finite element method is used as a foundation for modeling the mechanical response of anisotropic materials. As an illustrative example, we can cite the work of Weiss et al. [28] who have implemented a finite element formulation for incompressible hyperelastic materials in the general purpose finite element code NIKE3D developed by Maker [104]. Noted that the famous HGO model [44] was also implemented inside the ABAQUS commercial code. Thus, the second primary focus of this chapter is to perform the finite element implementation of our model by following the total Lagrangian formulation describe in section 1.7 of chapter 1. To assess the appropriateness of this new density, numerical simulations are compared with experimental results extracted from the paper published by Ciarletta et al. [39]. For our purpose, the interest of the work of Ciarletta et al. is threefold:

• It provides a large variety of experimental results by testing two different materials, each in four different situations (tensile and pure shear loadings parallel and transverse to the rubber-reinforcement direction), covering a large scope of the material behavior. It therefore constitutes a good trial for the assessment of models because a single set of material parameters should have to match all the four experimental tests.

• If tensile tests prevail in the literature, shear tests are uncommon although they can be considered as a severe benchmark case for rubber material models. As outlined by Horgan et al. in [105]: "The classical problem of simple shear in nonlinear elasticity has played an important role as a basic pilot problem involving a homogeneous deformation that is rich enough to illustrate several key features of the nonlinear theory, most notably the presence of normal stress effects. (• • • ) Since shearing is one of the dominant modes of behavior of rubbers in applications, this raises major concerns. Put another way, simple shear is not so simple after all".

• A new hyperelastic model using a non classical measure of strain is also proposed in [39]. In the same vein, Fereidoonnezhad et al. [40] have built later a model using this kind of strain, reporting the non-linearity aspect from the form of the SEF to the strain invariant, and have used the experimental data provided by Ciarletta et al. to assess their model. Our new model can therefore be compared not only with experimental results but also with numerical simulations.

Finally, it must be noted that the part of this research work related to the construction of the new SEF has been published in the International Journal of Solids and Structures [65].

The chapter is organized as follows:

• In section 2.2, some preliminaries on the material model and on the material symmetry group are introduced.

• The polyconvexity and physical interpretation of the new invariants proposed in [54] are investigated in section 2.3.

• Four different tests are briefly introduced in section 2.4. Those tests have been performed by Ciarletta et al. [39]. They include uniaxial tension and simple shear loadings.

• A new hyperelastic model based on the new invariants is presented in section 2.5.

The predicted results of our proposed model is compared to the experimental data extracted from [39].

• The finite element implementation of our new model is given in section 2.6. Firstly, a penalty function procedure is introduced to extent the constitutive model from the compressible to the incompressible range. Secondly, the calculations of the first and second derivatives of the strain energy density is performed. These calculations have been implemented in the finite element software FER [57].

• Finally, numerical results obtained thanks to the finite element software FER are presented in section 2.7. These results concern homogeneous deformations (with simple tension and shear tests) as well as inhomogeneous deformations (with a 3D tension test). In the case of homogeneous deformations, several numerical simulations were successfully compared to experimental and theoretical results extracted from [39]. This allows us to validate the finite element implementation.

2.2/ PRELIMINARIES

In this chapter, we focus on a fiber-reinforced material with a one-fiber family of direction a as depicted on Figure 2.1. We assume that a lies in the plane (E 1 , E 2 ) and forms an angle θ with E 1 :

a =              c s 0              , b =              -s c 0              with c = cos(θ), s = sin(θ) (2.1)
Practically, we will only consider the following two cases where the fibers are parallel

(θ = 0) or transverse (θ = π 2 ) to E 1 : Parallel: a =              1 0 0              , b =              0 1 0              , c =              0 0 1              (2.2) Transverse: a =              0 1 0              , b =              -1 0 0              , c =              0 0 1              (2.3)
To model this material, Ta et al. considered in [54] the group S O(3) of all the proper orthogonal transformation (that is to say the set of the 3 × 3 real matrix satisfying equation (1.24) ) and the material symmetry group G containing all the orthogonal transformations of S O(3) leaving invariant the material structure. This group G can be described as the group of all rotations around the fiber direction a because it is consistent with the geometric symmetries described in figure 1.3. Using a mathematical argument based on an

K 1 = ρ 1 ; K 2 = ρ 2 + ρ 3 ; K 3 = ρ 2 5 + ρ 2 4 ; K 4 = ρ 2 6 -ρ 2 ρ 3 K 5 = (ρ 2 5 -ρ 2 4 )ρ 6 + ρ 4 ρ 5 (ρ 2 -ρ 3 ) ; K 6 = (ρ 2 4 -ρ 2 5 )(ρ 2 -ρ 3 ) + 4ρ 4 ρ 5 ρ 6 (2.4)
where the coefficients ρ i stand for:

ρ 1 = Ca, a ; ρ 2 = Cb, b ; ρ 3 = Cc, c ρ 4 = Ca, b ; ρ 5 = Ca, c ; ρ 6 = Cb, c (2.5)
However, at this stage, it is important to notice that the invariant K 5 can be excluded from the integrity basis provided that three orthogonal symmetries through three orthogonal planes (one perpendicular to a and the two others containing a) are added to the material symmetry group. Indeed, if we consider for example the orthogonal symmetry S 2 related to the plane P 2 (see figure 2.2), we obtain:

S 2 (a) = a ; S 2 (b) = b ; S 2 (c) = -c (2.6)
Figure 2.2: Three orthogonal planes of symmetry It follow the next from equation (2.5) that:

S T 2 CS 2 a, a = CS 2 a, S 2 a = Ca, a = ρ 1 S T 2 CS 2 b, b = CS 2 b, S 2 b = Cb, b = ρ 2 S T 2 CS 2 c, c = CS 2 c, S 2 c = C(-c), -c = ρ 3 S T 2 CS 2 a, b = CS 2 a, S 2 b = Ca, b = ρ 4 S T 2 CS 2 a, c = CS 2 a, S 2 c = Ca, -c = -ρ 5 S T 2 CS 2 a, a = CS 2 b, S 2 c = Cb, -c = -ρ 6 (2.7)
Equation (2.7) proves that, under the action of S 2 , the coefficients ρ 1 , ρ 2 , ρ 3 and ρ 4 are unchanged while ρ 5 and ρ 6 are transformed to their opposite. Using this result in equation (2.4) gives:

((-ρ 5 ) 2 -ρ 2 4 )(-ρ 6 ) + ρ 4 (-ρ 5 )(ρ 2 -ρ 3 ) = -K 5 (2.8)
Therefore K 5 is not invariant under the new material symmetry group extended with the three orthogonal symmetries. A general discussion on the different ways to define the material symmetry group can be consulted in [9]. We therefore define a strain energy density W only depending on five of the six invariant polynomials given by equation (2.4) :

W = W(K 1 , K 2 , K 3 , K 4 , K 6 ) (2.9)
Additionally, for accounting for the incompressibility condition J =det(F)= 1, we introduce the extra pressure p (which plays the role of a Lagrange multiplier) into the formulation of the second Piola-Kirchhoff stress tensor S in equation (1.21) :

S = 2 ∂W ∂C -pC -1 (2.10)
The corresponding Cauchy stress tensor σ is obtained by combining equation (1.20) with equation (2.10) :

σ = 2F ∂W ∂C F T -pI (2.11)
By reminding that the nominal stress is the transpose of the engineering stress P introduced by equation (1.18) , we deduce from equation (2.11) that:

P T = JF -1 σ = 2 ∂W ∂C F T -pF -1 (2.12)
or equivalently by using equation (2.9) and the chain derivative rule:

P T = 2 6 i=1,i 5 ω i ∂K i ∂C F T -pF -1 (2.13)
where ω i represents the derivative of the SEF with respect to the invariants:

ω i = ∂W ∂K i (2.14)
The derivatives ∂K i ∂C of the invariants with respect to C are calculated straightforwardly from equations (2.4) and (2.5) :

∂K 1 ∂C = a ⊗ a ; ∂K 2 ∂C = b ⊗ b + c ⊗ c ∂K 3 ∂C = ρ 4 [a ⊗ b + b ⊗ a] + ρ 5 [a ⊗ c + c ⊗ a] ∂K 4 ∂C = ρ 6 [b ⊗ c + c ⊗ b] -ρ 2 c ⊗ c -ρ 3 b ⊗ b ∂K 6 ∂C = (ρ 4 [a ⊗ b + b ⊗ a] -ρ 5 [a ⊗ c + c ⊗ a])(ρ 2 -ρ 3 ) + (ρ 2 4 -ρ 2 5 )[b ⊗ b -c ⊗ c] + 2(ρ 4 ρ 5 [b ⊗ c + c ⊗ b] + ρ 4 ρ 6 [a ⊗ c + c ⊗ a] + ρ 5 ρ 6 [a ⊗ b + b ⊗ a]) (2.15)
Note that the Lagrange multiplier p involved in equation (2.13) will only be used later for analytical calculations (section 2.4). For the finite element computation (section 2.6), we will prefer to introduce a penalty function because it allows to reduce the number of unknowns.

To conclude this section, we summarize in the following table the main mechanical properties of the invariants obtained by Ta et al. in [54]. This table will be helpful to understand the role played by those invariants when we will write explicit SEF models. In this table, "purely non-linear" means the corresponding invariants tend to zero when the strains are small. Purely non-linear Tensile behavior Shear behavior

K 1 K 2 , K 4 K 3 , K 6

2.3/ POLYCONVEXITY AND PHYSICAL INTERPRETATION OF THE IN-VARIANTS

As explained in Chapter 1, the polyconvexity of the strain energy density is a prerequisite for finding solutions in compatibility with physical requirements [41]. We thus investigate in this section the polyconvexity of the new five invariants K i (i = 1, 2, 3, 4, 6) introduced by equation (2.4) . To this end, let us recall (see definition 3, section 1.6, chapter 1) that a function is said to be polyconvex, if it can be expressed as a convex function of the three arguments F, Co f (F) and det(F). Because some combinations of classical invariants are known to be polyconvex [46], we first remind that the new invariants given by equation (2.4) are related to the classical ones by [54]:

K 1 = I 4 ; K 2 = I 1 -I 4 ; K 3 = I 5 -I 4 2 ; K 4 = I 1 I 4 -I 5 -I 2 (2.16
)

K 6 = (I 1 -I 4 )(I 5 -I 4 2 ) + 2[I 3 + I 4 (I 1 I 4 -I 5 -I 2 )] (2.17)
where the classical invariants I i (i = 1, ..., 5) have been introduced in chapter 1 by equations (1.31) and (1.60) .

Conversely, the classical invariants can be expressed with respect to the new ones by:

I 1 = K 1 + K 2 ; I 2 = K 1 K 2 -K 3 -K 4 ; I 3 = -K 1 K 4 + 1 2 (K 6 -K 2 K 3 )
(2.18)

I 4 = K 1 ; I 5 = K 2 1 + K 3 (2.19)
We first note that K 1 , K 2 and -K 4 are well known polyconvex functions (see [46] for details). These three invariants respectively represent (see [54], [98] and [103]):

• the elongation squared in the fiber direction,

• the elongation squared in the isotropic plane perpendicular to the fiber direction,

• the deformation of an area element with a unit normal parallel to the fiber direction.

It is also noted that the expression (2.16) of K 3 is close to I 1 I 4 -I 5 which is known to be polyconvex [46]. That implies that K 1 K 2 -K 3 is a polyconvex combination of K 3 with the other K i . This is straightforward to prove from equations (2.18) and (2.19) :

I 1 I 4 -I 5 = (K 1 + K 2 )K 1 -(K 2 1 + K 3 ) = K 1 K 2 -K 3 (2.20)
Another way to build a polyconvex combination including K 3 is to start with I 2 1 -K 3 and use equation (2.16) :

I 2 1 -K 3 = I 2 1 + I 2 4 -I 5 = 1 2 [I 2 4 + I 2 1 + (I 1 -I 4 ) 2 ] + I 1 I 4 -I 5 (2.21)
The two terms I 1 -I 4 and I 1 I 4 -I 5 are known to be polyconvex [46] and the square and the sum of polyconvex functions is still polyconvex. We then conclude from equations (2.21) that I 2 1 -K 3 is polyconvex. It can be expressed only with the new invariants K i by using equations (2.18) :

I 2 1 -K 3 = (K 1 + K 2 ) 2 -K 3 (2.22)
We have thus obtained two polyconvex functions including K 3 (namely is observed that the second one is less conventional and has been guessed from the particular expression (2.16) of the new invariant K 3 . Additionally, this new invariant K 3 can be physically interpreted by noting first its positivity. The proof directly derives from the application of the Cauchy-Schwarz inequality:

K 1 K 2 -K 3 = I 1 I 4 -I 5 and (K 1 + K 2 ) 2 -K 3 = I 2 1 + I 2 4 -I 5 ). If
I 4 = Ca, a ≤ Ca a = Ca ⇒ I 2 4 ≤ Ca 2 = C 2 a, a = I 5 ⇒ K 3 ≥ 0 (2.23)
It is then noted that the inequality (2.23) is strict except if the vectors Ca and a are parallel, that is to say if a is an eigenvector of C:

K 3 = 0 ⇔ ∃λ ∈ R/Ca = λa (2.24)
From a physical point of view, Eq. (2.24) indicates that K 3 is equal to zero if, and only if, the material is submitted to a pure axial loading in the fiber direction. K 3 can therefore serve as an indicator of the amount of shear in the fiber direction which is null if K 3 is equal to zero. It is possible to specify more precisely this shear indicator by using a convenient orthonormal basis, says B (see Figure 2.3), for the strain calculation in the fiber direction a. In this view, we select a as a first vector basis. The second vector basis in the isotropic plane is defined as the part of Ca orthogonal to a (remind that Ca represents the strain in the fiber direction):

b = Ca -Ca, a a Ca -Ca, a a = Ca -I 4 a Ca -I 4 a (2.25)
That ensures naturally the orthogonality between a and b. It is also noticed that the norm of Ca -I 4 a adopts a remarkable form directly related to the invariant K 3 :

Ca -I 4 a = Ca -I 4 a, Ca -I 4 a = I 5 -I 2 4 = K 3 (2.26)
The cross product between a and b is used to calculate the third vector c completing the orthonormal basis B (Figure 2.3):

c = a ∧ b = a ∧ Ca -I 4 a I 5 -I 2 4 = a ∧ Ca I 5 -I 2 4
(2.27)

The strain tensor in the fiber direction is represented by Ca and its components in the basis B take a very simple form by applying Eqs. (2.25) and (2.27) :

Ca =              Ca, a Ca, b Ca, c              = 1 I 5 -I 2 4                I 4 I 5 -I 2 4 Ca, Ca -I 4 a Ca, a ∧ Ca                =                I 4 I 5 -I 2 4 0                =              I 4 √ K 3 0              (2.28)
As expected, the first component of Ca expressed in the basis B represents the elongation in the fiber direction and is equal to I 4 . The two other components, which are related to the shear effect between the fiber direction and the isotropic plane, prove that the total amount of shear can be estimated by

I 5 -I 2 4 = √ K 3 .
One other important result is the fact that this total amount of shear strain, which is of course invariant by any rotation in the isotropic plane, is concentrated in the direction b and equal to zero in the direction c (2.25) and (2.27) . In other words, the maximum shear effect between the fiber direction a and the isotropic plane arises with the orthogonal projection of Ca in the isotropic plane (equation (2.25) ).

These results demonstrate that the classical invariant I 5 plays a key role to estimate the shear effect between the fiber direction and the isotropic plane, provided that I 5 is combined with I 4 through the definition (2.16) of the new invariant K 3 . This new invariant sheds a new light on I 5 while it is often reported in the literature that a physical interpretation of A complementary interpretation of K 3 can be given with the coefficient β 4 defined below and introduced in [106] as the magnitude of the along-fiber shear strain (corresponding to a fiber sheared along an adjacent fiber):

β 4 = I 5 I 2 4 -1 = √ K 3 K 1 (2.29)
In the light of equations (2.28) , the coefficient β 4 introduced by equation (2.29) can also be interpreted as a function of the ratio between the shear strain and the axial strain in the fiber direction. That means that the total amount of shear strain between the fiber direction and the isotropic plane can be estimated by K 3 (equation (2.16) ), as the absolute difference between I 5 and I 2 4 , or by β 4 (equation (2.29) ), as the relative ratio between I 5 and I 2 4 . It is less evident to consider the invariant K 6 because equation (2.17) does not reveal any polyconvex known functions, except of course I 3 = det(F) 2 which is convex with respect to det(F). Even if it is an obvious result, we can state that the following combination including

POLYCONVEXITY AND PHYSICAL INTERPRETATION OF THE INVARIANTS 57

K 6 is polyconvex: K 6 -2K 1 K 4 -K 2 K 3 = I 3 (2.30)
Besides, it is possible to link K 6 with along-fiber shear effect by reporting equation (2.29) in the first term of equation (2.17) :

(I 1 -I 4 )(I 5 -I 2 4 ) = K 2 (K 1 β 4 ) 2 (2.31)
For its part, it is remarked that the second term of equation (2.17) is related to the shear angle ϕ between the matrix and the fiber:

2[I 3 + I 4 (I 1 I 4 -I 5 -I 2 )] = -2I 3 tan 2 ϕ (2.32)
Where the shear angle ϕ is introduced in [64] by:

tan 2 ϕ = (I 5 -I 1 I 4 + I 2 )I 4 I 3 -1 (2.33)
Combining equations (2.17) , (2.32) and (2.33) allows to connect the invariant K 6 to two different types of shear effects:

K 6 = K 2 (K 1 β 4 ) 2 -2I 3 tan 2 ϕ (2.34)
The first term in equation (2.34) takes into account the along-fiber shear effect between two adjacent fibers while the second term controls the shear interaction between the matrix and the fiber.

To conclude this section, we can say that three invariants used in this chapter (among the five included in the integrity basis introduced by equations (2.16) - (2.17) ), are well known mixed polyconvex invariants of the literature. This is a notable result because the methodology used to calculate these invariants [54]:

• is based on a generalized Reynolds operator built in the framework of the theory of polynomial invariants;

• deeply differs from classical approaches used in the literature;

• does not account, a priori, for any aspects related to polyconvexity.

We have also established that the two other invariants of the integrity basis can be combined with the others to build polyconvex functions. However, two of them (equations (2.20) and (2.30) ) are classical ones and further investigations will be needed to examine deeply whether additional original combinations are possible or not.

In the following sections, the five invariants K i defined by equations (2.16) - (2.17) will be kept for three reasons:

• they arise naturally by applying a rigorous mathematics approach based on the theory of invariant polynomials [54];

• they are all physically motivated as it has been explained above, particularly the a priori non polyconvex invariants K 3 and K 6 which are related to two different types of shear effects;

• the invariants K 3 and K 6 play a key role for nonlinear analysis as reported in [54] because they tend to zero in the case of small infinitesimal strain, meaning in that that they are purely nonlinear.

2.4/ UNIAXIAL TENSION AND SIMPLE SHEAR TESTS

The experimental data obtained by Ciarletta et al. [39], concern: i) a simple tension test parallel to the fiber direction, ii) a simple tension test transverse to the fiber direction, iii) a simple shear test parallel to the fiber direction, iv) a simple shear test transverse to the fiber direction.

Since these four experiments are used in this work as a reference to assess our model, we perform in the two following sections (2.4.1 and 2.4.2) the analytical calculation of the nominal stress in those four cases. In this way, we will be able to make further comparisons between three kind of results coming respectively from an experimental set-up, from a finite element computation and from a theoretical calculation.

2.4.1/ UNIAXIAL TENSION CASE

Consider a block of material subjected to a simple tension loading as illustrated on Figure 2.4. The back side (opposite to the applied tension t) and the two lateral faces (down and back) are simply supported. These boundary conditions lead to the following homogenous deformation:

F =              λ 1 0 0 0 λ 2 0 0 0 λ -1 1 λ -1 2              ⇒ C =              λ 2 1 0 0 0 λ 2 2 0 0 0 λ -2 1 λ -2 2              (2.35)
where λ 1 and λ 2 represent the principal stretches and account for incompressibility condition J =det(F)= λ 1 λ 2 λ 3 = 1. If the tension loading is applied parallel to the fiber direction, as shown in the left part of Figure 2.4, the three vectors a, b, and c must be selected according to equation (2.2) . The six coefficients ρ i defined by equation (2.5) and the five invariants K i (except K 5 ) defined by equation (2.4) can therefore be simplified:

ρ 1 = λ 2 1 ; ρ 2 = λ 2 2 ; ρ 3 = λ -2 1 λ -2 2 ; ρ 4 = ρ 5 = ρ 6 = 0 (2.36)
K 1 = λ 2 1 ; K 2 = λ 2 2 + λ -2 1 λ -2 2 ; K 3 = K 6 = 0 ; K 4 = -λ -2 1 (2.37)
Besides, the combination of equations (2.2) , (2.15) and (2.36) yields to:

∂K 1 ∂C =              1 0 0 0 0 0 0 0 0              ; ∂K 2 ∂C =              0 0 0 0 1 0 0 0 1              ; ∂K 4 ∂C = -              0 0 0 0 λ -2 1 λ -2 2 0 0 0 λ 2 2              ; ∂K 3 ∂C = ∂K 6 ∂C = 0 (2.38)
The nominal stress P p in the case where the tension loading is applied parallel to the fiber direction is therefore obtained by using equations (2.13) , (2.35) and (2.38) :

P p =              P p 11 0 0 0 P p 22 0 0 0 P p 33              ;              P p 11 = 2ω 1 λ 1 -pλ -1 1 P p 22 = 2(ω 2 λ 2 -ω 4 λ -2 1 λ -1 2 ) -pλ -1 2 P p 33 = 2(ω 2 -ω 4 λ 2 2 )λ -1 1 λ -1 2 -pλ 1 λ 2 (2.39)
The plane stress state P p 33 = 0 can be exploited to extract the hydrostatic pressure p and to express the tensile stress by:

P p 11 = 2(ω 1 λ 1 -ω 2 λ -3 1 λ -2 2 + ω 4 λ -3 1 ) (2.40)
Using the free boundary condition P p 22 = 0 gives:

λ 2 = λ -1/2 1 (2.41)
Reporting equation (2.41) in equation (2.40) yields to:

P p 11 = 2[ω 1 λ 1 -ω 2 λ -2 1 + ω 4 λ -3 1 ] (2.42)
Equation (2.41) makes sense because we deduce from it and from the incompressibility condition that:

λ 3 = λ -1/2 1 (2.43)
That means that the free faces perpendicular to E 2 and E 3 (see figure 2.4) are subjected to the same stretch. This is consistent with the physical experience related to a family of fibers aligned with the tension loading. At this stage, it is noticed that, if the material is at rest, the configuration must be free of stress (see equation (1.53) , section 1.5.1 of chapter 1):

λ 1 = λ 2 = λ 3 = 1 ⇒ P p 11 = 0 (2.44)
And by reporting equation (2.44) in the equation (2.42) :

λ 1 = λ 2 = λ 3 = 1 ⇒ ω 1 = ω 2 -ω 4 (2.45)
In the case where the tension loading is applied transverse to the fiber direction (right part of Figure 2.4), the three vectors a, b and c must be selected according to equation (2.3) . The six coefficients ρ i defined by equation (2.5) and the five invariants K i (except K 5 ) defined by equation (2.4) can therefore be simplified:

ρ 1 = λ 2 2 ; ρ 2 = λ 2 1 ; ρ 3 = λ -2 1 λ -2 2 ; ρ 4 = ρ 5 = ρ 6 = 0 (2.46) K 1 = λ 2 2 ; K 2 = λ 2 1 + λ -2 1 λ -2 2 ; K 3 = K 6 = 0 ; K 4 = -λ -2 2 (2.47)
Besides, the combination of equations (2.3) , (2.15) and (2.46) yields to:

∂K 1 ∂C =              0 0 0 0 1 0 0 0 0              ; ∂K 2 ∂C =              1 0 0 0 0 0 0 0 1              ; ∂K 4 ∂C = -              λ -2 1 λ -2 2 0 0 0 0 0 0 0 λ 2 1              ; ∂K 3 ∂C = ∂K 6 ∂C = 0 (2.48)
By using equations (2.13) , (2.35) and (2.48) , the nominal stress P t transverse to the fiber direction adopts a diagonal form:

P t =              P t 11 0 0 0 P t 22 0 0 0 P t 33              ;              P t 11 = 2(ω 2 λ 1 -ω 4 λ -1 1 λ -2 2 ) -pλ -1 1 P t 22 = 2ω 1 λ 2 -pλ -1 2 P t 33 = 2(ω 2 -ω 4 λ 2 1 )λ -1 1 λ -1 2 -pλ 1 λ 2 (2.49)
Reporting the hydrostatic pressure p in P t 11 from the plane stress condition P t 33 = 0 yields to:

P t 11 = 2ω 2 (λ 1 -λ -3 1 λ -2 2 ) (2.50)
Note that, this time, the free stress state is automatically satisfied if the material is at rest

(λ 1 = λ 2 = λ 3 = 1)
. Also note that, contrarily to the previous parallel loading case, the free boundary condition P t 22 = 0 leads this time to the following equation where λ 2 is unknown and λ 1 is the prescribed stretch applied to the sample:

ω 1 λ 4 2 -ω 2 λ -2 1 + ω 4 = 0 (2.51)
The coefficients ω 1 and ω 4 , defined by equation (2.14) , depends on W, λ 1 and λ 2 . At this stage, it is thus not possible to solve equation (2.51) if W, and consequently ω 1 and ω 4 , are not defined. In fact, depending on the choice of W, we obtain different kind of equations. For example, in the case of a density W quadratic with respect to the invariants (this case will be studied in details later), the equation to solve is a 7 degree polynomial equation.

We have solved it thanks to the fzero function of the MATLAB software and, by reporting the numerical solution λ 2 in the equation (2.50) , we will be able to express P t 11 with respect to λ 1 .

2.4.2/ SIMPLE SHEAR CASE

The field displacement related to a block of material subjected to a simple shear deformation (Figure 2.5) is expressed in a linear form of the amount of shear deformation k:

u = kX 2 E 1 (2.52)
It follows that the corresponding strain tensors is:

F =              1 k 0 0 1 0 0 0 1              ⇒ C =              1 k 0 k k 2 + 1 0 0 0 1             
(2.53)

.5: Simple shear test

If the shear loading is applied parallel to the fiber direction, the three vectors a, b, and c must be selected according to equation (2.2) . The six coefficients ρ i defined by equation (2.5) and the five invariants K i (except K 5 ) defined by equation (2.4) can therefore be simplified:

ρ 1 = ρ 3 = 1 ; ρ 2 = k 2 + 1 ; ρ 4 = k ; ρ 5 = ρ 6 = 0 (2.54) K 1 = 1 ; K 2 = k 2 + 2 ; K 3 = k 2 ; K 4 = -(k 2 + 1) ; K 6 = k 4 (2.55)
The combination of equations (2.2) , (2.15) and (2.54) yields to:

∂K 1 ∂C =              1 0 0 0 0 0 0 0 0              ; ∂K 2 ∂C =              0 0 0 0 1 0 0 0 1              ; ∂K 3 ∂C =              0 k 0 k 0 0 0 0 0              ∂K 4 ∂C =              0 0 0 0 -1 0 0 0 -(k 2 + 1)              ; ∂K 6 ∂C =              0 k 3 0 k 3 k 2 0 0 0 -k 2              (2.56)
The nominal stress P s p in the case where the shear loading is applied parallel to the fiber direction is therefore obtained by using equations (2.13) , (2.53) and (2.56) :

P s p =           
             ;                          P sp 11 = 2(ω 1 + k 2 ω 3 + k 4 ω 6 ) -p P sp 22 = 2(ω 2 -ω 4 + k 2 ω 6 ) -p P sp 33 = 2[ω 2 -(k 2 + 1)ω 4 -k 2 ω 6 ] -p P sp 12 = 2(kω 3 + k 3 ω 6 ) + kp P sp 21 = 2k(ω 2 + ω 3 -ω 4 + 2k 2 ω 6 )
(2.57) By using the plane stress condition P sp 33 = 0 as in Fereidoonnezhad et al. [40] to obtain the hydrostatic pressure p, the shear stress corresponding to a loading parallel to the fiber direction can be expressed by:

P sp 12 = 2[(ω 2 + ω 3 -ω 4 )k -ω 4 k 3 ] (2.58)
If the shear loading is applied now transverse to the fiber direction, the three vectors a, b, and c must be selected according to equation (2.3) . The six coefficients ρ i defined by equation (2.5) and the five invariants K i (except K 5 ) defined by equation (2.4) can therefore be simplified:

ρ 1 = k 2 + 1 ; ρ 2 = ρ 3 = 1 ; ρ 4 = -k ; ρ 5 = ρ 6 = 0 (2.59
)

K 1 = k 2 + 1 ; K 2 = 2 ; K 3 = k 2 ; K 4 = -1 ; K 6 = 0 (2.60)
The combination of equations (2.3) , (2.15) and (2.59) yields to:

∂K 1 ∂C =              0 0 0 0 1 0 0 0 0              ; ∂K 2 ∂C =              1 0 0 0 0 0 0 0 1              ; ∂K 3 ∂C =              0 k 0 k 0 0 0 0 0              ∂K 4 ∂C =              -1 0 0 0 0 0 0 0 -1              ; ∂K 6 ∂C =              k 2 0 0 0 0 0 0 0 -k 2              (2.61)
The nominal stress P st in the case where the shear loading is applied transverse to the fiber direction is therefore obtained from equations (2.13) , (2.53) and (2.61) :

P st =              P st 11 P st 12 0 P st 21 P st 22 0 0 0 P st 33              ;                          P st 11 = 2(ω 2 + k 2 ω 3 -ω 4 + k 2 ω 6 ) -p P st 22 = 2ω 1 -p P st 33 = 2[ω 2 -ω 4 -k 2 ω 6 ] -p P st 12 = 2kω 3 + kp P st 21 = 2k(ω 1 + ω 3 ) (2.62)
By using the plane stress condition P st 33 = 0 as in Fereidoonnezhad et al. [40] to obtain the hydrostatic pressure p, the shear stress with the loading transverse to the fiber direction can be expressed by:

P st 12 = 2[(ω 2 + ω 3 -ω 4 )k -ω 6 k 3 ] (2.63)
It is noticed that, in the two shear cases (equations (2.58) and (2.63) ), the free stress state is automatically reached for a material at rest corresponding to the situation where k is equal to zero (see equation (2.52) ).

We have now finished the analytical study allowing to compute the nominal stress in four different loading cases. However, to obtain a full achievement of these computations, it is necessary to know the values of the quantities ω i (i = 1, 2, 3, 4, 6) involved in the analytical formulas. To reach this goal, because ω i depends on the SEF W (see equation (2.14) ), we are going now to build the SEF W introduced by equation (2.9) as a function of the invariants K i .

2.5/ A NEW HYPERELASTIC SEF

Following the strategy used by Mooney and Rivlin to build isotropic energy densities [72, 62], we adopt in this work a polynomial form for W. To identify the coefficients of each monomial, we have performed comparison between the model and experiments data extracted from the work of Ciarletta et al. [39]. These experimental data concern two different fiber-reinforced rubbers: soft silicone rubber reinforced by polyamide (referenced as material A) and soft silicone rubber reinforced by hard silicone rubber (referenced as material B). The comparison focuses on the calculated and measured nominal stresses in the cases of a tensile loading and of a shear loading. Each loading is applied first parallel and secondly transverse to the fiber direction. The equations (2.42) , (2.50) , (2.58) and (2.63) , corresponding to these four loading cases, have been presented in the previous section and are summarized below: a) Tensile stress parallel to the fiber direction:

P p 11 = 2[(ω 2 -ω 4 )λ 1 -ω 2 λ -2 1 + ω 4 λ -3 1 ]
b) Tensile stress transverse to the fiber direction:

P t 11 = 2ω 2 (λ 1 -λ -3 1 λ -2 2 )
c) Shear stress parallel to the fiber direction:

P sp 12 = 2[(ω 2 + ω 3 -ω 4 )k -ω 4 k 3 ]
d) Shear stress transverse to the fiber direction:

P st 12 = 2[(ω 2 + ω 3 -ω 4 )k -ω 6 k 3 ]
(2.64)

We notice that those four equations depend on the derivatives ω i with respect to the second, third, fourth and sixth variable of W. We also remark that the only difference between the two shear stresses comes from the cubic term k 3 , the two linear terms in k being equal. A particular attention must therefore be paid to the model in order to account for this constraint because experimental results reveal that the slope of the stress between the two shear cases is not constant (Figures 2.7 and 2.9). To evaluate if a polynomial form of the strain energy density could be appropriate, a linear and a quadratic expressions are tested in the two next sections. In the mean time, we compare in all cases our model with the one proposed by Fereidoonnezhad et al. [40].

2.5.1/ LINEAR STRAIN ENERGY DENSITY

We introduce the following linear polynomial strain energy function:

W 1 = a 1 K 1 + a 2 K 2 + a 3 K 3 + a 4 K 4 + a 6 K 6 (2.65)
The five polynomial coefficients a 1 , a 2 , a 3 , a 4 and a 6 represent the material parameters.

We have identified them by using the classical coefficient of determination R 2 ∈ [0, 1] defined by:

R 2 = 1 - S S res S S tot (2.66)
Where S S res and S S tot are the residual sum and the total sum of squares respectively:

S S res = n i=1 (y i -f i ) 2 S S tot = n i=1 (y i -y) 2 (2.67)
where y i stands for the experimental data (extracted from the curves plotted on Figures 2.6 to 2.9 ), f i for the theoretical data (calculated with the SEF defined by equation (2.65) with the nominal stresses introduced by equation (2.64) ) and y for the mean of the experimental data:

y = 1 n n i=1 y i (2.68)
n represents the number of experimental data considered in one test. The closest to 1 R 2 is, the best the fit of the experimental data by the theoretical data will be. The derivatives

ω i = ∂W ∂K i
of this linear polynomial strain energy function with respect to K i are calculated easily from equation (2.65) :

ω i = a i i = 1, 2, 3, 4, 6 (2.69)
Additionally, it is deduced from equations (2.45) and (2.69) that the linear model is only defined by four material parameters because a 1 , a 2 and a 4 are linked by:

a 1 = a 2 -a 4 (2.70)
The data fitting was achieved through the solver fminsearch (unconstrained nonlinear minimization) of the Optimization Toolbox provided by the MATLAB commercial software by accounting for the constraint induced by equation (2.70) . The identified parameters are presented on Table 2 Figure 2.9: Comparison between numerical and experimental shear stresses -linear strain energy density (equation (2.65) )

The comparisons between the experimental and the numerical results are presented on Figures 2.6 and 2.7 for material A and on Figures 2.8 and 2.9 for material B. We can notice that the linear model fits fairly all the experimental results except for the shear test with a loading transverse to the fiber direction (Figure 2.7 right for material A and Figure 2.9 right for material B). This observation is not surprising because the experimental results show that the shear stress follows a linear law with respect to k but with a slope depending on the loading (parallel or transverse to the fiber direction). This behavior cannot be predicted by equation (2.64) because the two types of shear loading provide the same linear term. This is confirmed by the coefficient of determination R 2 presented in Tables 2.9 and 2.10. In the shear cases, for material B and with a linear SEF, R 2 is actually equal to 0.70 and 0.43. These two values are far from 1, indicating a poor match with the experimental results. To overcome this problem, note that a case sensitive shear parameter (i.e. a shear parameter taking different values depending on the loading case) has been introduced in [39] and [40]. But the results are not so satisfactory with R 2 equal to 0.69 for material B in the shear parallel case (see last line of Table 2.10). Additionally, we consider that it is preferable to not change the material parameters values depending on the considered loading case, otherwise the model will be uneasy to extend to the general situation where a complex non homogeneous loading is applied. We will therefore examine in the next section the improvement brought by a quadratic polynomial form of the strain energy density.

2.5.2/ QUADRATIC STRAIN ENERGY DENSITY

In order to improve the quality of the numerical prediction, particularly in the cases of the two shear loadings, we introduce a quadratic polynomial form of the strain energy density:

W 2 =a 1 K 1 + a 2 K 2 + a 3 K 3 + a 4 K 4 + a 6 K 6 + a 11 K 2 1 + a 12 K 1 K 2 + a 13 K 1 K 3 + a 14 K 1 K 4 + a 16 K 1 K 6 + a 22 K 2 2 + a 23 K 2 K 3 + a 24 K 2 K 4 + a 26 K 2 K 6 + a 33 K 2 3 + a 34 K 3 K 4 + a 36 K 3 K 6 + a 44 K 2 4 + a 46 K 4 K 6 (2.71)
Note that the term a 66 K 6 2 has no influence on the four loading cases studied in this work and has been therefore removed from the quadratic expression of the SEF. This term is actually only concerned by the shear case with a loading transverse to the fiber direction (see the last equation of (2.64) ), through the coefficient ω 6 = ∂W ∂K 6 . This coefficient is equal to 2a 66 K 6 in the present situation with K 6 equal to 0 according to equation (2.60) .

The nineteen polynomial coefficients of the SEF described by equation (2.71) are identified in the same way as in the previous section by fitting the classical coefficient of determination R 2 . The identified coefficients corresponding to the two materials A and B are presented on Tables 2 Note that it could have been possible to deduce from equation (2.45) a link between the material parameters of the quadratic energy density (2.71) . Actually, by combining equation (2.45) , which holds for a tension loading parallel to the fiber direction, with equation (2.37) and with the fact that the material is assumed to be at rest (λ 1 = λ 2 = λ 3 = 1), we obtain:

a 1 + 2a 11 + a 12 -a 2 -4a 22 + 3a 24 + a 4 -2a 44 = 0 (2.72)
However, we have decided to not include this equation in the identification process. We have actually considered that saving only one material parameter from this process was not of a great interest because it still leaves 18 material parameters to identify (see equation (2.71) ). Figure 2.13: Comparison between numerical and experimental shear stressesquadratic strain energy density (equation (2.71) )

The comparisons between the experimental and the numerical results are presented on Figures 2.10 and 2.11 for material A and on Figures 2.12 and 2.13 for material B. It is noted that the quadratic model improves the accuracy of the numerical results, particularly for the shear loading in the fiber direction for material B (Figure 2.13 left). This is confirmed by the improvement of the R 2 coefficient of determination from the W 1 density to the W 2 density (Table 2.10) which changes from 0.70 to 0.98 in the case of a shear loading parallel to the fiber direction. However, the problem detected with the linear model in the case of a shear loading transverse to the fiber direction still remains (Figures 2.11 and 2.13 right). The quadratic model obviously fails to find two different slopes for the two shear tests with a loading parallel and then transverse to the fibers direction. The goal of the next section is to fix this problem.

2.5.3/ LINEAR AND QUADRATIC STRAIN ENERGY DENSITIES COMBINED WITH A POWER-LAW FUNCTION

As shown on Figures 2.6 to 2.13, the predictive results of the linear and quadratic models agree with the experimental data for both materials A and B in the case of a tensile deformation. However, in the simple shear deformation case, there still exists some differences between the prediction results and the experimental data. As mentioned before, this is not surprising because equations (2.58) and (2.63) are not able to provide a different linear term in k while the experiments show two different slopes for the two shear loadings (parallel or transverse to the fiber direction). Even if we increase the degree of the polynomial SEF with quadratic terms, this will not change the fact that both equations (2.58) and (2.63) will keep the same linear terms in k. We actually obtain the following derivatives from equation (2.71)

ω 2 = ∂W 2 ∂K 2 = a 2 + a 12 K 1 + 2a 22 K 2 + a 23 K 3 + a 24 K 4 + a 26 K 6 (2.73) ω 3 = ∂W 2 ∂K 3 = a 3 + a 13 K 1 + a 23 K 2 + 2a 33 K 3 + a 34 K 4 + a 36 K 6 (2.74) ω 4 = ∂W 2 ∂K 4 = a 4 + a 14 K 1 + a 24 K 2 + a 34 K 3 + 2a 44 K 4 + a 46 K 6 (2.75)
In the shear case parallel to the fibers direction, we use equation (2.55) to simplify the three preceding formulas

ω 2 = a 2 + a 12 + 2a 22 (2 + k 2 ) + a 23 k 2 -a 24 (1 + k 2 ) + a 26 k 4 (2.76
)

ω 3 = a 3 + a 13 + a 23 (2 + k 2 ) + 2a 33 k 2 -a 34 (1 + k 2 ) + a 36 k 4 (2.77
)

ω 4 = a 4 + a 14 + a 24 (2 + k 2 ) + a 34 k 2 -2a 44 (1 + k 2 ) + a 46 k 4 (2.78)
The linear term with k corresponding to ω 2 + ω 3 -ω 4 in equation (2.58) is therefore given by We perform the same calculation in the shear case transverse to the fibers direction by using equation (2.60)

a 2 +
ω 2 = a 2 + a 12 (1 + k 2 ) + 4a 22 + a 23 k 2 -a 24 (2.80
)

ω 3 = a 3 + a 13 (1 + k 2 ) + 2a 23 + 2a 33 k 2 -a 34 (2.81
)

ω 4 = a 4 + a 14 (1 + k 2 ) + 2a 24 + a 34 k 2 -2a 44 (2.82)
By combining ω 2 , ω 3 and ω 4 from equation (2.80) , (2.81) and (2.82) to calculate the linear term with respect to k included in ω 2 + ω 3 -ω 4 , we unfortunately find the same linear term as the one given by equation (2.79) . Compared to the linear density, the quadratic one improves the accuracy of the numerical results in the shear case but not enough as it was expected because the linear term which is the only one really concerned is not affected.

In fact we can prove a more general statement: any polynomial SEF will provide the same linear expansion in k when we consider the shear tests in the parallel and transverse situations. Let us prove this fact. We denote by W m a general polynomial of degree m in the variables K 1 , K 2 , K 3 , K 4 and K 6 :

W m = γ=(γ 1 ,γ 2 ,γ 3 ,γ 4 ,γ 6 ) γ 1 +γ 2 +γ 3 +γ 4 +γ 6 ≤m a γ K γ 1 1 K γ 2 2 K γ 3 3 K γ 4 4 K γ 6 6 (2.83)
where a γ are constant material parameters.

The linear terms with k embedded in the shear stress expressions come from the following derivatives (equations (2.58) and (2.63) ):

ω 2 + ω 3 -ω 4 = ∂W m ∂K 2 + ∂W m ∂K 3 - ∂W m ∂K 4 (2.84)
Using the equation (2.83) and deriving W m with respect to K 2 , K 3 and K 4 yields to:

ω 2 + ω 3 -ω 4 = γ=(γ 1 ,γ 2 ,γ 3 ,γ 4 ,γ 6 ) γ 1 +γ 2 +γ 3 +γ 4 +γ 6 ≤m a γ A γ (2.85)
where A γ is defined by:

A γ = K γ 1 1 K γ 6 6 (γ 2 K γ 2 -1 2 K γ 3 3 K γ 4 4 + γ 3 K γ 2 2 K γ 3 -1 3 K γ 4 4 -γ 4 K γ 2 2 K γ 3 3 K γ 4 -1 4 ) (2.86)
As we are only interested in the linear terms with respect to k, we have to remove the non constant terms from equation (2.85) . To reach this goal, it is first reminded that K 1 is equal to 1 if a shear loading is applied parallel to the fiber direction (equation (2.55) ) while it is equal to k 2 + 1 in the transverse case (equation (2.60) ). That means that the constant term coming from K

γ 1
1 is always 1, whatever the value of γ 1 is. Besides, because K 6 is equal to k 4 if the shear loading is parallel to the fiber direction (equation (2.55) ) while it is equal to 0 in the transverse case (equation (2.60) ), γ 6 must be equal to 0 in the equation (2.85) in order to produce non zero constant terms. To discuss the role of the exponents γ 2 , γ 3 and γ 4 in the same manner as γ 1 and γ 6 , we need to focus on the generic term A γ included in the sum of equation (2.85) by taking γ 6 = 0:

A γ = K γ 1 1 (γ 2 K γ 2 -1 2 K γ 3 3 K γ 4 4 + γ 3 K γ 2 2 K γ 3 -1 3 K γ 4 4 -γ 4 K γ 2 2 K γ 3 3 K γ 4 -1 4 ) (2.87)
In order to go further, we going to distinguish to case of the parallel shear loading to the case of the transverse shear loading:

• Case 1: parallel shear loading

By reporting the values of the invariants K i from the equation (2.55) in equation (2.87) yields to:

A γ = γ 2 (k 2 + 2) γ 2 -1 (k 2 ) γ 3 (-1) γ 4 (k 2 + 1) γ 4 + γ 3 (k 2 + 2) γ 2 (k 2 ) γ 3 -1 (-1) γ 4 (k 2 + 1) γ 4 -γ 4 (k 2 + 2) γ 2 (k 2 ) γ 3 (-1) γ 4 -1 (k 2 + 1) γ 4 -1 (2.88)
Let us examine now one by one the constant contributions with respect to k coming from each term of equation (2.88) . The term γ 2 (k 2 + 2) γ 2 -1 (k 2 ) γ 3 (-1) γ 4 (k 2 + 1) γ 4 gives:

       0 i f γ 3 0 γ 2 2 γ 2 -1 (-1) γ 4 i f γ 3 = 0 (2.89)
The term γ 3 (k 2 + 2) γ 2 (k 2 ) γ 3 -1 (-1) γ 4 (k 2 + 1) γ 4 gives:

       0 i f γ 3 = 0 or γ 3 > 1 2 γ 2 (-1) γ 4 i f γ 3 = 1 (2.90)
The term -γ 4 (k 2 + 2) γ 2 (k 2 ) γ 3 (-1) γ 4 -1 (k 2 + 1) γ 4 -1 gives:

       0 i f γ 3 0 -γ 4 2 γ 2 (-1) γ 4 -1 i f γ 3 = 0 (2.91)
By considering the restrictions (2.89) , (2.90) and (2.91) in equation (2.88) , it follows that the constant term with respect to k in equation (2.85) is:

γ=(γ 1 ,γ 2 ,0,γ 4 ,0) γ 1 +γ 2 +γ 4 ≤m a γ (-1) γ 4 2 γ 2 -1 (γ 2 + 2γ 4 ) + γ=(γ 1 ,γ 2 ,1,γ 4 ,0) γ 1 +γ 2 +γ 4 +1≤m a γ 2 γ 2 (-1) γ 4 (2.92)
• Case 2: transverse shear loading

By reporting the values of the invariants K i from equation (2.60) in the equation (2.87) yields to:

A γ = (k 2 + 1) γ 1 γ 2 2 γ 2 -1 (k 2 ) γ 3 (-1) γ 4 + (k 2 + 1) γ 1 γ 3 2 γ 2 (k 2 ) γ 3 -1 (-1) γ 4 -(k 2 + 1) γ 1 γ 4 2 γ 2 (k 2 ) γ 3 (-1) γ 4 -1 (2.93)
Let us examine now one by one the constant contributions with respect to k coming from each term of equation (2.93) . The term (k 2 + 1) γ 1 γ 2 2 γ 2 -1 (k 2 ) γ 3 (-1) γ 4 gives:

       0 i f γ 3 0 γ 2 2 γ 2 -1 (-1) γ 4 i f γ 3 = 0 (2.94)
The term (k 2 + 1) γ 1 γ 3 2 γ 2 (k 2 ) γ 3 -1 (-1) γ 4 gives:

       0 i f γ 3 = 0 or γ 3 > 1 2 γ 2 (-1) γ 4 i f γ 3 = 1 (2.95)
The term -(k 2 + 1)

γ 1 γ 4 2 γ 2 (k 2 ) γ 3 (-1) γ 4 -1 gives:        0 i f γ 3 0 -γ 4 2 γ 2 (-1) γ 4 -1 i f γ 3 = 0 (2.96)
By considering the restrictions (2.94) , (2.95) and (2.96) in equation (2.93) , we retrieve exactly the same result as the one described by equation (2.92) .

Equation (2.92) holds therefore for both cases of parallel and transverse shear loading. Thus for any polynomial SEF W m , the linear terms of equations (2.58) and (2.63) will be the same. Increasing the degree of the polynomials W cannot improve significantly the accuracy of prediction in the shear case. In order to overcome this problem, which leads to a rather poor fitting with the shear tests, we will add an additional term to the linear and quadratic densities defined by equations (2.65) and (2.71) . It will make sense to adopt a power-law form for this additive term since it looks like an extension of a monomial with a real number exponent. But there are numerous possibilities for combining the invariants defined by equations (2.16) and (2.17) in a power form. In order to perform the appropriate combination, we remark from equation (2.64) , which gives a view of all the four load cases at a glance, that:

i) the additive term should modify the linear terms corresponding to each shear test (equations (2.58) and (2.63) ). That means that the additive density could possibly depends on K 2 , K 3 and K 4 .

ii) as the tensile tests are perfectly fitted by the linear and quadratic densities, any additional term should not affect the tensile results. That means that the additive density should not depends on K 2 and K 4 .

iii) the term ω 1 is not concerned by the four loading cases. That means that the additive density could possibly depends on K 1 . iv) the term ω 6 is not concerned by the two tensile loading cases. That means that the additive density could possibly depends on K 6 .

Based on these considerations, it is relevant to propose a new term W add adopting the following form:

W add = αK 3 K c 1 1 + βK 6 K c 2 1 (2.97)
where α, β, c 1 and c 2 are new material parameters.

To evaluate the influence of W add on the two shear cases, we first derived equation (2.97) with respect to K 3 (resp. K 6 ) and we secondly use equation (2.55) (resp. equation (2.60) ):

∂W add ∂K 3 =       
α shear stress parallel to the fiber direction α(1 + k 2 ) c 1 shear stress transverse to the fiber direction (2.98)

∂W add ∂K 6 = β(1 + k 2 )
c 2 shear stress transverse to the fiber direction (2.99)

Note that the calculation of ∂W add ∂K 6 was only performed in the case of a shear loading transverse to the fiber direction because, in the parallel case, the coefficients ω 6 is not concerned (see equation (2.58) ). It is also remarked that equation (2.97) meets the requirement of remark i) because the contribution of the constant part with respect to k of W add in the expression ω 2 + ω 3 -ω 4 is:

      
α shear stress parallel to the fiber direction 0 shear stress transverse to the fiber direction (2.100)

One could argue that the quantities α(1 + k 2 ) c 1 of equation (2.98) could provide α as a constant term with respect to k if the exponent c 1 is a positive whole number. But, fortunately, as shown by Tables 2.5, 2.6, 2.7 and 2.8, the identification of the material parameters always gives a negative value for c 1 . That explains why the quantity α(1 + k 2 ) c 1 does no contain any constant term with respect to k, providing the zero term for the shear stress transverse to the fiber direction in equation (2.100) . So we can retain equation (2.97) and combine it with the linear strain energy (2.65) :

W 3 = (a 2 -a 4 )K 1 + a 2 K 2 + a 3 K 3 + a 4 K 4 + a 6 K 6 + αK 3 K c 1 1 + βK 6 K c 2 1 (2.101)
It is noted that we have again replaced the coefficient a 1 by a 2a 4 , exactly as in the case of the purely linear density (see equation (2.65) and (2.70) ). To justify this replacement, we use equations (2.45) and (2.101) :

a 1 + αc 1 K c 1 -1 1 K 3 + βc 2 K c 2 -1 1 K 6 = a 2 -a 4 (2.102)
As the equation (2.45) only holds for a tension loading parallel to the fiber direction with (2.37) ), equation (2.102) can be simplified as expected to:

K 3 = K 6 = 0 (equation
a 1 = a 2 -a 4 (2.103)
The curve representing the prediction of the model against the experimental data and other numerical calculations extracted from [40] ). But we also remark from Figure 2.17 left that, if the predicted curve provides a correct averaged trend of the experimental points, a quadratic form of the energy density would be probably more suitable. From now, the superposition of the quadratic energy density equation (2.71) with the additive density equation (2.97) is considered:

W 4 =a 1 K 1 + a 2 K 2 + a 3 K 3 + a 4 K 4 + a 6 K 6 + a 11 K 2 1 + a 12 K 1 K 2 + a 13 K 1 K 3 + a 14 K 1 K 4 + a 61 K 1 K 6 + a 22 K 2 2 + a 23 K 2 K 3 + a 24 K 2 K 4 + a 26 K 2 K 6 + a 33 K 2 3 + a 34 K 3 K 4 + a 36 K 3 K 6 + a 44 K 2 4 + a 46 K 4 K 6 + αK 3 K c 1 1 + βK 6 K c 2 1 (2.104)
We observe on Figures 2.18 to 2.21 a good agreement between the numerical results and the experimental data. This agreement is confirmed by the coefficient of determination R 2 which is equal to 1 and 0.99 for material A and B respectively (last column of Tables 2.9 and 2.10). It is actually considered that a value greater than 0.9 typically represents a satisfactory fit to the experimental data. As awaited, in the case of a shear loading parallel to the fiber direction, the inclusion of quadratic terms in the SEF allows to fix the problem encountered with the previous linear model (Figure 2.21 left versus Figure 2.17 left). Figure 2.21: Comparison between numerical and experimental shear stressesquadratic + power form strain energy density (equation (2.104) )

The identified values of the 23 material parameters associated with the strain energy density equation (2.104) are presented on Tables 2.7 and 2.8. These material parameters where used to plot the curve of Figures 2.18 to 2.21. The only drawback of the strain energy density described by equation (2.104) is the large number of material parameters which are needed to be identified.

2.6/ FINITE ELEMENT IMPLEMENTATION

The aim of the present section is to propose a finite element implementation dealing with the strain energy density introduced by equation (2.101) . This density combines a linear and power form expression with respect to the invariants. We have selected it for the FE implementation because it provides a good balance between the accuracy of the predictions and the number of material parameters to identify. To perform the FE implementation, the total Lagrangian formulation is adopted according to the description of this formulation given in section 1.7 of chapter 1. In order to extend the constitutive model from the compressible to the incompressible range, we introduce a penalty function W, instead of the Lagrange multiplier used to for the analytical calculations in the case of homogeneous deformation (equation (2.10) ). This function, which enforces the incompressibility condition J = det(F) = 1, permits to reduce the number of unknowns by removing the Lagrange multiplier:

W = W + W(J) (2.105) W(J) = 1 d 1 2 (J 2 -1) -ln(J) + c ln(J) (2.106)
The first term of equation (2.106) is similar to the one proposed in [107, 108] while the second term is introduced for guaranteeing the reference configuration to be stress free as suggested in [109]. The numerical parameter d is set to a value of 10 -8 which is a good balance between the satisfaction of the incompressibility condition and the convergence of the Newton-Raphson scheme (equation (1.131) ). The second parameter c will be calculated in order to ensure that the material is stress free if the displacement field is zero.

To do that, we first need to calculate the nominal stress by replacing W by W in equation (2.12) and by removing from this equation the Lagrange multiplier p because we consider a penalty method for the FE implementation:

P T = 2 ∂ W ∂C F T = 2 6 i=1,i 5       ω i ∂K i ∂C + ∂ W ∂C       F T (2.107)
The derivatives ω i = ∂W ∂K i , which is part of the first term included in the bracket of equation (2.107) , are calculated straightforwardly from equation (2.101) :

ω 1 = a 2 -a 4 + αc 1 K 3 K c 1 -1 1 + βc 2 K 6 K c 2 -1 1 ; ω 2 = a 2 ω 3 = a 3 + αK c 1 1 ; ω 4 = a 4 ; ω 6 = a 6 + βK c 2 1
(2.108) By using the equation [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] , the derivatives ∂K i ∂C , which is also part of the first term included in the bracket of equation (2.107) , can be rewritten from equation (2.15) by:

∂K 1 ∂C = M a ; ∂K 2 ∂C = M b + M c = I -M a ∂K 3 ∂C = 2(ρ 4 M ab + ρ 5 M ac ) ; ∂K 4 ∂C = 2ρ 6 M bc -ρ 2 M c -ρ 3 M b ∂K 6 ∂C = 2(ρ 4 M ab -ρ 5 M ac )(ρ 2 -ρ 3 ) + (ρ 2 4 -ρ 2 5 )[M b -M c ] + 4(ρ 4 ρ 5 M bc + ρ 4 ρ 6 M ac + ρ 5 ρ 6 M ab ) (2.109)
by introducing the symmetric matrices M a , M b , M c , M ab , M ac and M bc :

M a = a ⊗ a ; M b = b ⊗ b ; M c = c ⊗ c M ab = 1 2 (a ⊗ b + b ⊗ a) ; M ac = 1 2 (a ⊗ c + c ⊗ a) ; M bc = 1 2 (b ⊗ c + c ⊗ b) (2.110)
The second term included in the bracket of equation (2.107) is obtained straightforwardly from equation (2.106) :

∂ W ∂C = ∂ W ∂J ∂J ∂C = 1 2 1 d J 2 -1 + c C -1 (2.111)
where we have used a standard derivative result [70]:

∂J ∂C = J 2 C -1 (2.112)
Because the material must be free of stress if the displacement field is null, we consider this particular case in equations (1.10) , (1.11) , (2.4) , (2.5) , (2.108) , (2.109) and (2.111) :

F = C = I ρ 1 = ρ 2 = ρ 3 = 1 ; ρ 4 = ρ 5 = ρ 6 = 0 K 1 = 1 ; K 2 = 2 ; K 3 = 0 ; K 4 = -1 ; K 6 = 0 ω 1 = a 2 -a 4 ; ω 2 = a 2 ; ω 3 = a 3 + α ; ω 4 = a 4 ; ω 6 = a 6 + β ∂K 1 ∂C = a ⊗ a ; ∂K 2 ∂C = - ∂K 4 ∂C = b ⊗ b + c ⊗ c ; ∂K 3 ∂C = ∂K 6 ∂C = 0 ; ∂ W ∂C = c 2 I (2.113)
Replacing equation (2.113) in (2.107) gives:

U = 0 =⇒ P T = 2 (a 2 -a 4 )(a ⊗ a + b ⊗ b + c ⊗ c) + c 2 I = 0 (2.114)
Or, equivalently, by using equation [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] :

U = 0 =⇒ P T = 2 a 2 -a 4 + c 2 I = 0 (2.115)
We then deduce from equation (2.115) that:

c = 2(a 4 -a 2 ) (2.116)
Equation (2.116) links some of the material parameters of the model in order to guarantee that the reference configuration specified by U = 0 is stress free.

To construct the tangent stiffness matrix for the analysis of nonlinear structures by the finite element method, one has now to determine the stress-strain tangent operator D, which is a fourth order tensor resulting from the derivation of S with respect to E (see equations (1.126) and (1.127) ). In order to calculate D, we first compute the anisotropic part of the second Piola-Kirchhoff stress tensor S related to W from equations (1.21) and (2.101) :

S = 2 ∂W ∂C = 2         6 i=1,i 5 ∂W ∂K i ∂K i ∂C         = 2         6 i=1,i 5 ω i ∂K i ∂C         (2.117)
In order to obtain the fourth-order tensor D (equation (1.126) ), we derive again W with respect to C from equation (2.117) :

D = 4          6 i=1,i 5 6 j>i, j 5 ω i j ∂K i ∂C ⊗ ∂K j ∂C + ∂K j ∂C ⊗ ∂K i ∂C + i=1,i 5 ω i ∂ 2 K i ∂C 2 + ω ii ∂K i ∂C ⊗ ∂K i ∂C          (2.118)
where the coefficients ω i j stand for the second derivative of W with respect to the invariants K i and K j . They are obtained straightforwardly from equation (2.108) :

ω 11 = αc 1 (c 1 -1)K 3 K c 1 -2 1 + βc 2 (c 2 -1)K 6 K c 2 -2 1 ; ω 13 = αc 1 K c 1 -1 1 ; ω 16 = βc 2 K c 2 -1 1 ω 12 = ω 14 = ω 22 = ω 23 = ω 24 = ω 26 = ω 33 = ω 34 = ω 36 = ω 44 = ω 46 = ω 66 = 0 (2.119)
To obtain the second derivative

∂ 2 K i ∂C 2
, we derive the first derivatives contained in equation (2.109) with respect to C:

∂ 2 K 1 ∂C 2 = ∂ 2 K 2 ∂C 2 = 0 ; ∂ 2 K 3 ∂C 2 = 2 {N abab + N acac } ; ∂ 2 K 4 ∂C 2 = 2N bcbc -D bbcc ∂ 2 K 6 ∂C 2 = 2 {ρ 4 ( D abbc + 2N acbc ) + (ρ 2 -ρ 3 )(N abab -N acac ) + ρ 5 (N abbc -D acbc ) + 2ρ 6 N abac }
(2.120) where we have introduced the following fourth-order tensors:

N abab = M ab ⊗ M ab ; N acac = M ac ⊗ M ac ; N bcbc = M bc ⊗ M bc N abbc = M ab ⊗ M bc + M bc ⊗ M ab ; N acbc = M ac ⊗ M bc + M bc ⊗ M ac N abac = M ab ⊗ M ac + M ac ⊗ M ab ; D abbc = M ab ⊗ (M b -M c ) + (M b -M c ) ⊗ M ab D acbc = M ac ⊗ (M b -M c ) + (M b -M c ) ⊗ M ac ; D bbcc = M bb ⊗ M cc + M cc ⊗ M bb (2.121)
We have now finished to calculate all the quantities involved in the anisotropic fourth-order tensor D given by equation (2.118) . But, to achieve the finite element implementation, we need to compute the fourth-order tensor D vol related to the volumetric part of the strain energy density. This tensor is obtained by derivating W twice with respect to C from equation (2.106) and by using the first derivative given by equation (2.111) :

D vol = 4 ∂ 2 W ∂C 2 = 2 ∂ ∂C J 2 -1 d + c ⊗ C -1 + J 2 -1 d + c ∂C -1 ∂C (2.122)
The first term in (2.122) is easily calculated by using equation (2.112) :

∂ ∂C J 2 -1 d + c = 2J d ∂J ∂C = 2J d J 2 C -1 = J 2 d C -1 (2.123)
The second term can be expressed by using a standard result related to the derivative of the inverse of C [42, 110, 111]:

∂C -1 ∂C = -C -1 C -1 (2.124)
where the tensor notation is defined by equation [START_REF] Rw Ogden | Non-linear elastic deformations[END_REF] .

By reporting equations (2.123) and (2.124) in equation (2.122) , we finally obtain:

(D vol ) i jkl = 2J 2 d C -1 i j C -1 kl - J 2 -1 d + c C -1 ik C -1 jl + C -1 il C -1 jk (2.125)
According to equations (2.118) and (2.125) , the finite element implementation of the second derivative of the strain energy densities described by equations (2.101) , (2.105) and (2.106) was realized inside the FER code. This university code is developed by the Laboratory of Mechanics of the University of Evry (France) and many standard hyperelastic densities have already been implemented during the past 15 years [48, 111, 112]. The implementation was achieved by using C++ language and following the procedure described in

2.7/ FE SIMULATION RESULTS

In this section, in order to validate the finite element implementation, we consider the two different fiber-reinforced rubbers (material A and material B) that have been introduced in section 2.5. We recall that these two materials were originally studied in [39]. For both of them, four different loading cases are considered leading to two different values of the angle θ introduced in Figure 2.1:

1. a simple tension test parallel to the fiber direction (θ = 0), 2. a simple tension test transverse to the fiber direction (θ = π 2 ), 3. a simple shear test parallel to the fiber direction (θ = 0), 4. a simple shear test transverse to the fiber direction(θ = π 2 ).

These four loading cases provide homogeneous deformation giving closed-form solutions which can be compared with FE simulations. The next section 2.7.1 is dedicated to this comparison. In the following section 2.7.2, a 3D example including inhomogeneous deformation is presented, demonstrating the capability of FER and of our model to deal with more complex 3D problems.

2.7.1/ 2D HOMOGENEOUS DEFORMATION

In this section, we compare the finite element computations with closed form solution and also with experimental data. The closed form solutions are summarized by equation (2.64) while the experimental data are extracted from [39]. The finite element computations were performed with FER by using the values of material parameters listed in Tables 2 In order to illustrate the capability of the FE implementation of our model, a 3D example is considered in this section. The material parameters are those of material B and are given in Table 2.6. For this study, the collagen is embedded as one family of fibers that are disposed parallel to the tensile (axial) direction. The strip model has a length of L = 60 mm, a width of W = 20 mm and a thickness of T = 10 mm. Figure 2.27 represents the initial mesh which includes 96 hexahedral elements. One side of the specimen is fixed and a displacement of 5 mm is applied on the other side leading to a uniaxial tensile test. It is noted that the embedded fibers could increase the thickness at some locations of the specimen. Due to the incompressibility character, the width of the specimen decreases in order to balance the increase of thickness. This phenomenon could be clearly observed from the isovalues of U z displacement, as shown in Figure 2.31. It is also noted from Figure 2.30 that the displacements U y and U z remain almost linear with respect to the prescribed displacement U x until a value of U x = 1.5 mm is reached. After this value, the behavior becomes clearly non-linear because of the effect of the embedded fibers. This typical behavior induced by the fibers has been already mentioned in [48]. In fact, most of the classical laws are separated into an isotropic and an anisotropic parts and are case sensitive with respect to the fiber stretch [42]. This separation is based on the fact that the shortening of the fibers is assumed to generate no stress and the stiffness is only due to the ground substance in this case. The fibers are only acting if they are extended. It is remarkable that our model is able to predict this kind of behavior while our SEF combines in a single energy the isotropic and the anisotropic effects. Moreover, we have not considered any sensitive case to separate the ground substance and the fiber influences.

It is also noticed that the amount of non-linear effects depends on the location of the considered points. For example, this amount is less pronounced for C, located in the middle of the specimen, than A located near the place where the prescribed displacements is loaded. In this second chapter, a novel strain energy function (SEF) was developed for modeling hyperelastic incompressible fiber-reinforced materials with a single fiber direction. The construction of this energy is based on a family made of five new invariants combined in a polynomial form. It has been shown that three of these invariants are polyconvex while the two others are physically motivated by shear effects. We have studied several polyno-mial forms and proved in the most general case that none is able to predict properly the behavior of fibered materials. It is actually impossible to match correctly the experimental data coming from a shear loading applied first parallel and next transverse to the fiber direction if a polynomial strain energy density is considered. In order to overcome this problem, a power form function has been added to a linear and then a quadratic polynomial. These two kinds of combination provide a fair agreement with experiments and we have selected the linear option (equations (2.101) , (2.105) and (2.106) ) because it requires less material parameters than the quadratic option (8 against 23).

The finite element implementation was performed inside the FER university code [57] by using a total Lagrangian approach. All the details of the implementation, that is to say the calculations of the strain and stress incremental forms, as well as the tangent stiffness matrix, are provided in the section 2.6 of this chapter.

In order to validate the proposed finite element model, the FER results were successfully compared with four different closed-form solutions corresponding to four different loading conditions: uniaxial tensile and shear tests, each one with a loading first parallel and next transverse to the fiber direction. These closed-form solutions have been determined in section 2.4 and are summarized by equation (2.64) . A fair agreement was also found with experimental data corresponding to these four loading cases and extracted from [39].

Finally, A 3D simulation of a strip specimen under tension loading was performed successfully to show the applicability of the constitutive model in the context of a more complex finite element analysis than the one corresponding to the closed-form solutions previously mentioned.

CHAPTER 3

A NEW SEF FOR FOUR-FIBER FAMILY MATERIALS 3 A NEW SEF FOR FOUR-FIBER FAMILY MATERIALS 3.1/ INTRODUCTION
We introduce in this chapter a new hyperelastic model for predicting the nonlinear mechanical properties of anisotropic hyperelastic materials under biaxial stretching. The proposed strain energy function (SEF) can be applied for understanding the nature of behavior laws for material with four-fiber family structure, which has a large potential of applications, particularly in biomechanics, surgical and interventional therapies for peripheral artery disease (PAD).

Atherosclerosis is for example nowadays one of the most important subjects of medical and biochemical research. Actually, ischemia, angina pectoris, myocardial infarction, stroke, or heart failure and other fatal diseases are consequences of atherosclerosis [66].

Thus some treatment must be taken in order to reduce the occurrence of diseases caused by atherosclerosis [113]. Unfortunately high treatment charges of peripheral vascular operations are often unable to guarantee a good therapeutic effect and repetitive intervention are needed ( [114][115] [116][117]). However, it has been proven in [118] and [119] that the mechanical stress and strain of arterial tissue are deeply connected to atherosclerosis. These aspects are the reason why the study about mechanical properties of the arterial tissue got more and more attention in the last decades ([120] [121][122] [123]). Holzapfel et al. [42] have for example introduced structural SEFs for describing the soft biological tissues such as the arterial wall. Kamenskiy et al. [1] used a planar biaxial extension set-up to test diseased superficial femoral, popliteal and tibial arteries from 170 patients to determine their passive biaxial mechanical properties. This very complete study includes a large variety of measurements with different loading conditions (9 combinations of a biaxial test), 3 different human parts tested (superficial femoral, popliteal and tibial arteries). Using the same experiment facility, Kamenskiy et al. tested more recently [124] the fresh femoropopliteal arteries from 70 human subjects (13-79 year old) to study the effect of the age on the physiological and mechanical characteristics. The references [1] and [124] will therefore be used as a basis of our work to assess the capability and the appropriateness of our model.

It should also be noticed that, based on the work of Baek et al. ([125][126] [127]), Kamenskiy et al. also proposed in their papers a constitutive model to describe the behavior of the diseased arteries with four-fiber families of collagen. This model adopts an exponential form which extends the well-known two-fiber family model introduced by Holzapfel, Gasser and Ogden in [42]. As most of the constitutive models met in the literature ( [125][127] [42]), the strain energy function (SEF) introduced by A. V. Kamenskiy et al. is divided into an isotropic and an anisotropic parts.

More recently, an original approach mixing the isotropic and the anisotropic parts in a single SEF was introduced by Ta et al. [53]. This method takes advantage of the theory of polynomial invariants (namely the Noether's theorem and the Reynolds operator) to compute an integrity basis made of 7 new invariants consistent with the considered type of anisotropy. These new invariants are considered in our work instead of the classical isotropic ones and instead of the anisotropic mixed invariants found in the literature. This original approach is motivated by the fact that the new invariants do not require a separation of the SEF into an isotropic and an anisotropic part. Another motivation is the rigorous mathematical foundations used to define those invariants. Finally, we demonstrate in this chapter that the integrity basis found in [53] can be recombined in a smart way in order to form a new integrity basis made of polyconvex invariants. This is a major issue because, in the context of hyperelastic problems, the polyconvexity of the strain energy density is often considered as a prerequisite for ensuring the existence of solutions [41].

Practically, we have introduced an original SEF as a quadratic polynomial depending on these 7 new invariants. Among the 7 new invariants, 3 are linear with respect to the right Cauchy-Green deformation tensor C, 3 are quadratic with C and 1 is cubic. The cubic invariant is linked to incompressibility and will be taken into account through a Lagrange multiplier. It will therefore not be directly included in the SEF. Moreover, to enrich the quadratic part of the energy density, the squares of the 3 linear invariants are considered, leading to a total of 9 monomials: 3 are linear and 6 quadratic with respect to C. We have therefore 9 material parameters related to each of these monomials. But, as we demonstrate that 2 among the 9 material parameters are dependent, due to the fact that a zero stress corresponds to a zero strain, we have finally just 7 material parameters to identify.

To assess the appropriateness of this new density, numerical simulations were successfully compared to experimental and theoretical results extracted from [1] in the case of a biaxial testing. The main results and conclusions are:

• An excellent fit of the experimental data,

• The mean trends of the experimental curves, evaluated through the standard coef-ficient of determination (R 2 ), are better described with our model than with the ones used in [1],

• For large stretches, the model introduced in [1] provides better results than our, likely because we used polynomial functions instead of exponential functions,

• A perfect identification of the material parameters with no possible doubt on the optimal solution is offered by our model. This nice result comes from the polynomial form of the SEF which presents a linear dependence of the density with respect to the material parameters. It allows performing a linear least square identification leading to a single optimal solution.

• Compared to the model proposed in [1], which includes 8 material parameters, our model only needs 7 material parameters.

The material parameters and the strain-energy function developed in this chapter are intended to serve as a basis for a finite element implementation and to investigate the problem of atherosclerosis for possible improvements of the treatment. Therefore, the second aim of this chapter is to implement our constitutive model in a finite element code. This implementation was realized in C++ language with the FER university code [57] by adopting the total Lagrangian formulation. Several numerical examples, including homogeneous deformation (biaxial tension loading) and non-homogeneous deformation (3D uniaxial tensile loading), are presented to show the validity of the model.

Note that the research work presented in this chapter has been accepted for publication in the International Journal of Solids and Structures [67].

The chapter is organized as follows:

• Section 3.2 introduces the material understudy.

• In section 3.3, the polyconvexity and physical interpretation of the invariants are investigated.

• The material model is presented in section 3.4 with the definition of the SEF with respect to the invariants. The identification of the material parameters is performed. The closed-form solution corresponding to the biaxial stretching employed in [1, 124] is also determined and the proposed model is validated by comparison with other model and experimental data extracted from the literature [1].

• The implementation of the proposed model in the finite element code FER is presented in section 3.5.

• Numerical results obtained from the finite element software FER are presented in section 3.6. These results concern homogeneous deformation (three different arteries tested with 5 combinations of different biaxial stretch) as well as nonhomogeneous deformation (with a 3D tension test).

3.2/ MATERIAL UNDERSTUDY

In this chapter, we focus on a fiber-reinforced material with four fiber families, such as the arterial wall studied in [1] and [124]. The general opinion is that the arterial wall is incompressible, hyperelastic and anisotropic. It actually consists in a mixture of an elastindominated amorphous matrix and families of locally parallel collagen fibers. Based on the work of Kamenskiy et al.

([1][124]
), the research object (diseased superficial femoral (SFA), popliteal (PA) and tibial arteries (TA)) includes four fiber families of collagen fiber: two oriented axially and circumferential, and two symmetrically along the diagonal as depicted on Figure 3.1. We assume that the two fiber directions a and b lie in the plane (e 1 , e 2 ) and form respectively an angle θ and -θ with e 1 . The longitudinal fiber direction c and circumferential fiber direction d are parallel to e 1 and e 2 respectively.

To model this kind of materials, we adopt in this chapter the invariants introduced by Ta et al. in [53]. In the case of a fiber reinforced material with a two-fibers family of directions a and b as depicted on figure 1.4, Ta et al. have presented a systematic method to find a list of invariants associated with the material symmetry group S 8 defined by equation (1.37) . By using a mathematical argument based on the Reynolds operator and on the Noether's theorem, they have demonstrated that all the polynomial invariants can be generated by 209 invariants. Additionally, they have demonstrated that some of these 209 invariants are linked together and that the following 7 polynomial invariants form an integrity basis of the ring of invariant polynomials under the material symmetry group:

H 1 = ρ 1 H 2 = ρ 2 H 3 = ρ 3 H 4 = ρ 2 4 H 5 = ρ 2 5 H 6 = ρ 2 6 H 7 = ρ 4 ρ 5 ρ 6 (3.1)
Where the coefficients ρ i stand for:

ρ 1 = Ce 1 , e 1 ρ 2 = Ce 2 , e 2 ρ 3 = Ce 3 , e 3 ρ 4 = Ce 1 , e 2 ρ 5 = Ce 1 , e 3 ρ 6 = Ce 2 , e 3 (3.2)
However, at this stage, it is remarkable to notice that these invariants originally introduced by Ta et al. for a two-fibers family can also be employed for a four-fibers family provided that the geometric plane symmetries shown on Figure 1.4 are satisfied with the four fiber directions. Because the four fiber directions considered in this work (Figure 3.1) satisfy this requirement, we will adopt the same invariants as the ones introduced in [53]. One can argue that it is surprising that two different materials, respectively made of a two and a four fibers family, can be modeled with the same invariants. However, these two materials will not be necessarily modeled by the same SEF, even if the invariants are identical. Additionally, as the two materials behave differently, the identification process will provided different values of the material parameters.

VARIANTS

We investigate in this section the polyconvexity of the new invariants introduced by equation (3.1) . First of all, it is reminded [46] that convexity implies polyconvexity and that any function φ of the form φ(F) = Fv, Fv (where v represents any non-zero vectors) is a convex function. It is therefore immediate to conclude that H 1 , H 2 and H 3 are polyconvex functions because it is straightforward from equations (3.1) and (3.2) that: Additionally, we observe from equations (3.1) and (3.2) that:

H 1 = Fe 1 ,
C(e 1 + e 2 ), e 1 + e 2 = Ce 1 , e 1 + Ce 2 , e 2 + 2 Ce 1 , e 2 = H 1 + H 2 + 2ρ 4 (3.5) In this way, we have linked three polyconvex functions (namely C(e 1 + e 2 ), e 1 + e 2 from the left-hand side of equation (3.5) and H 1 and H 2 from the right-hand side) with ρ 4 which represents the square root of H 4 . To make H 4 appear, which is of interest for us, we square equation (3.5) :

C(e 1 + e 2 ), e 1 + e 2 2 = (H 1 + H 2 ) 2 + 4H 4 + 4ρ 4 (H 1 + H 2 ) (3.6)
In order to eliminate the double product 4ρ 4 (H 1 + H 2 ) from equation (3.6) , we proceed in the same manner but by replacing e 1 + e 2 by e 1 -e 2 in equation (3.5) :

C(e 1 -e 2 ), e 1 -e 2 2 = (H 1 + H 2 ) 2 + 4H 4 -4ρ 4 (H 1 + H 2 ) (3.7)
And we finally add equations (3.6) and (3.7) to obtain: It remains to deal with the last invariant H 7 which adopts the following form according to equations (3.1) and (3.2) :

L 4 = (H 1 + H 2 ) 2 + 4H 4 = 1
H 7 = Ce 1 , e 2 Ce 1 , e 3 Ce 2 , e 3 = Fe 1 , Fe 2 Fe 1 , Fe 3 Fe 2 , Fe 3 (3.10)
There are again no clear evidences proving that H 7 could be a polyconvex function, so we can try to exhibit a polyconvex combination of H 7 with other invariants. First we noticed that H 7 is expressed in a cubic form with respect to C (equation (3.10) ), exactly as the classical third isotropic invariant I 3 (equation (1.31) ). As we have proved in section 1.6 that I 3 is a polyconvex function, linking H 7 with I 3 could be a good strategy to obtain a polyconvex combination of invariants involving H 7 . To find this link, we first observe that the strain tensor C is expressed in the {e 1 , e 2 , e 3 } basis by:

C =              ρ 1 ρ 4 ρ 5 ρ 4 ρ 2 ρ 6 ρ 5 ρ 6 ρ 3              (3.11)
A simple algebraic calculation gives:

I 3 = det(C) = 2ρ 4 ρ 5 ρ 6 + ρ 1 ρ 2 ρ 3 -ρ 1 ρ 6 2 -ρ 3 ρ 4 2 -ρ 2 ρ 5 2 (3.12)
Finally we use equations (3.1) and (3.12) to introduce L 7 as follows:

L 7 = I 3 2 = H 7 + 1 2 {H 1 H 2 H 3 -H 1 H 6 -H 3 H 4 -H 2 H 5 } (3.13)
The new invariant L 7 introduced by equation (3.13) is thus a polyconvex function involving H 7 . We have therefore built a family {L 1 , L 2 , L 3 , L 4 , L 5 , L 6 , L 7 } of polyconvex invariants where we have set:

L 1 = H 1 L 2 = H 2 L 3 = H 3 (3.14)
This family forms an integrity basis as a polynomial combination of H 1 , H 2 , H 3 , H 4 , H 5 , H 6 and H 7 (refer to equations (3.8) , (3.9) , (3.13) and (3.14) ). Additionally, each polynomials of this integrity basis meet a physical interpretation. For example, L 1 represents the elongation squared in the direction e 1 because equations (3.1) , (3.2) and (3.14) yield to:

L 1 = Fe 1 2 (3.15)
Similarly, L 2 and L 3 represent the elongation squared in directions e 2 and e 3 , respectively.

L 7 is directly connected to the deformed volume through equation (3.13) . Finally, to give a physical meaning to the invariants L 4 , L 5 and L 6 , it is first necessary to recall that the shear angle ϕ between two directions u and v (Figure 3.2) is defined by:

cos ϕ = Fu, Fv Fu Fv = Cu, v Cu, u 1/2 Cv, v 1/2 (3.16)
The change of the shear angle between the deformed and the reference configurations can therefore be measured by:

cos ϕ -cos ϕ 0 = Cu, v Cu, u 1/2 Cv, v 1/2 - u, v u, u 1/2 v, v 1/2 (3.17)
Equation (3.17) is introduced in [128] as an invariant related to the amount of shear. However, there is no argument for proving that cos ϕcos ϕ 0 is polyconvex. Fortunately, the polyconvex invariant L 4 introduced by equation (3.8) can be linked to cos ϕcos ϕ 0 by replacing u by e 1 and v by e 2 in (3.17) : Using equations (3.1) , (3.2) , (3.8) , (3.14) and (3.18) yields to:

cos ϕ -cos ϕ 0 = Ce 1 , e 2
cos ϕ -cos ϕ 0 = 1 2 L 4 -(L 1 + L 2 ) 2 L 1 L 2 (3.19)
Therefore L 4 can be seen as a polyconvex invariant involved in the amount of shear between directions e 1 and e 2 . Similarly, L 5 and L 6 are linked to the amount of shear related to directions (e 2 , e 3 ) and (e 1 , e 3 ) respectively.

It is finally remarked that, as demonstrated in [67], the set of invariants L 1 to L 7 , described by equations (3.8) , (3.9) , (3.13) and (3.14) , is equivalent to the classical set of invariants proposed in [129] for orthotropic symmetry:

T r(C) ; T r(C 2 ) ; T r(C The Cauchy stress tensor related to a SEF W, which depends on the invariants L i , is calculated from equations (1.20) and (1.21) :

σ = 2F         7 i=1 ∂W ∂L i ∂L i ∂C         F T -pI (3.21)
where we have introduced the extra pressure p to account for incompressibility.

The derivatives of the invariants L i with respect to C are calculated straightforwardly from equations (2.112) , (3.1) , (3.2) , (3.8) , (3.9) , (3.13) and (3.14) :

∂L 1 ∂C = e 1 ⊗
∂L 7 ∂C = ∂det(C) 2∂C = 1 2 det(C)C -1 (3.22)
As seen in equation (2.10) , the second Piola-Kirchhoff stress tensor S already includes an extra pressure term pC -1 to account for incompressibility. This extra pressure term is very similar to the last line of equation (3.22) . To avoid any redundancy, it is therefore logic to exclude L 7 from the strain energy density W:

W = W(L 1 , L 2 , L 3 , L 4 , L 5 , L 6 ) (3.23)

3.4.2/ CONSTITUTIVE MODEL

Following the strategy used by Mooney and Rivlin to build isotropic energy densities ( [72], [62]), we adopt in this work a polynomial form for W. One major advantage of this choice is the extreme ease of identifying the material parameters of the model as it will be explained in the next section. The variables of the polynomial correspond to the 6 new invariants introduced by equations (3.8) , (3.9) and (3.14) . This choice is motivated by the facts that these invariants are polyconvex, linked to a physical meaning and related to the invariants exhibited with a rigorous mathematics approach by Ta et al. [53]. To obtain the best flexibility in the model, but with a moderate number of material parameters, we introduce a second order polynomial with respect to L 1 , L 2 , L 3 , L 4 , L 5 and L 6 :

W = a 1 L 1 + a 2 L 2 + a 3 L 3 + a 4 L 4 + a 5 L 5 + a 6 L 6 + a 7 L 1 2 + a 8 L 2 2 + a 9 L 3 2 (3.24)
It is noticed that the six first terms of W are linear with the invariants L i and only L 1 , L 2 and L 3 are squared (the three last terms). We have made this choice because L 4 , L 5 and L 6 are already squared functions of the invariants H i (equations (3.8) and (3.9) ) while L 1 , L 2 and L 3 are only linear functions of them (equation (3.14) ).

The nine polynomial coefficients a 1 to a 9 represent the material parameters. It is possible to reduce their number from 9 to 7 by using the fact that the stress level must be zero in the case where there is no loading. To exploit this property, we first calculate ω i (i = 1, ..., 6) the derivatives of W with respect to L 1 , L 2 , L 3 , L 4 , L 5 and L 6 from equation (3.24) :

ω 1 = ∂W ∂L 1 = a 1 + 2a 7 L 1 ; ω 2 = ∂W ∂L 2 = a 2 + 2a 8 L 2 ; ω 3 = ∂W ∂L 3 = a 3 + 2a 9 L 3 (3.25) ω 4 = ∂W ∂L 4 = a 4 ; ω 5 = ∂W ∂L 5 = a 5 ; ω 6 = ∂W ∂L 6 = a 6 (3.26)
In the case where the displacement is equal to zero, giving F = C = I, equations (3.1) , (3.2) , (3.8) , (3.9) and (3.14) are simplified to:

ρ 1 = ρ 2 = ρ 3 = 1 ; ρ 4 = ρ 5 = ρ 6 = 0 (3.27) H 1 = H 2 = H 3 = 1 ; H 4 = H 5 = H 6 = 0 (3.28) L 1 = L 2 = L 3 = 1 ; L 4 = L 5 = L 6 = 4 (3.29)
Reporting equation (3.29) in equation (3.25) yields to:

∂W ∂L 1 = a 1 + 2a 7 ; ∂W ∂L 2 = a 2 + 2a 8 ; ∂W ∂L 3 = a 3 + 2a 9 (3.30)
We finally report equations (3.22) , (3.26) , (3.27) and (3.30) in equation (3.21) :

σ = 2              a 1 + 2a 7 + 4a 4 + 4a 5 -p 2 0 0 0 a 2 + 2a 8 + 4a 4 + 4a 6 -p 2 0 0 0 a 3 + 2a 9 + 4a 5 + 4a 6 -p 2              (3.31)
Making σ equal to zero from equation (3.31) , calculating the extra pressure p from one of the three equations obtained and reporting the result in the two others lead to express a 1 and a 2 with respect to the other material parameters: a 1 = a 3 -4a 4 + 4a 6 -2a 7 + 2a 9 ; a 2 = a 3 -4a 4 + 4a 5 -2a 8 + 2a 9

(3.32)

3.4.3/ CLOSED-FORM SOLUTION FOR A BIAXIAL STRETCHING

The experimental data obtained by Kamenskiy et al. [1] are related to a large variety of samples tested quasi-statically by using a custom-made soft-tissue biaxial testing device. The arteries samples were tested under a biaxial stretching with a different ratio of loading applied to the longitudinal and circumferential directions with the following proportions: 1:1, 1:2, 1:4, 2:1 and 4:1. Since these five experiments are used in this work as a reference to assess our model, we perform below the calculation of the Cauchy stress in the case of a biaxial stretching. To reach this goal, we consider a cubic block of material subjected to a biaxial tension loading as illustrated on Figure 3.3. The different ratio of loads were applied to the right and top faces of the cube (represented by the applied loading T 1 and T 2 in the longitudinal and circumferential directions). To model the symmetric boundary conditions induced by the biaxial stretching, we consider that the bottom, left and back faces of the cube are simply supported while the front face is free. These boundary conditions, represented by arrows on figure 3.3, lead classically to the homogenous deformations described by equation (2.35) .

Because e 1 , e 2 and e 3 constitutes an orthonormal basis where C adopts the diagonal matrix expression (2.35) , the 6 coefficients defined by equation (3.2) and the invariants defined by equations (3.8) , (3.9) and (3.14) can be simplified:

ρ 1 = λ 1 2 ρ 2 = λ 2 2 ρ 3 = λ 1 -2 λ 2 -2 ρ 4 = ρ 5 = ρ 6 = 0 (3.33) L 1 = λ 1 2 L 2 = λ 2 2 L 3 = λ 1 -2 λ 2 -2
(3.34)

L 4 = (λ 1 2 + λ 2 2 ) L 5 = (λ 1 2 + λ 1 -2 λ 2 -2 ) 2 L 6 = (λ 2 2 + λ 1 -2 λ 2 -2 ) 2 (3.35)
Where λ 1 and λ 2 represent the principal stretches and where the incompressibility condition λ 3 =λ -1 1 λ -1 2 was used. In the case of a biaxial tension loading, the Cauchy stress tensor σ is finally expressed in a diagonal form by reporting equations (2.35) , (3.33) , (3.34) and (3.22) in equation (3.21) :

σ =              σ 11 0 0 0 σ 22 0 0 0 σ 33              (3.36)
with:

σ 11 = 2[ ∂W ∂L 1 + 2 ∂W ∂L 4 (λ 1 2 + λ 2 2 ) + 2 ∂W ∂L 5 (λ 1 2 + λ 1 -2 λ 2 -2 )]λ 1 2 -p σ 22 = 2[ ∂W ∂L 2 + 2 ∂W ∂L 4 (λ 1 2 + λ 2 2 ) + 2 ∂W ∂L 6 (λ 2 2 + λ 1 -2 λ 2 -2 )]λ 2 2 -p σ 33 = 2[ ∂W ∂L 3 + 2 ∂W ∂L 5 (λ 1 2 + λ 1 -2 λ 2 -2 ) + 2 ∂W ∂L 6 (λ 2 2 + λ 1 -2 λ 2 -2 )]λ 1 -2 λ 2 -2 -p (3.37)
where the derivatives of W with respect to the invariants L i are obtained from equations (3.25) , (3.26) , (3.32) and (3.34) : The free boundary condition σ 33 = 0 can be exploited from the third equation of (3.37) to extract the hydrostatic pressure p and to finally express the tensile stress, respectively in the longitudinal and circumferential directions, by:

∂W ∂L 1 = a 3 -4a 4 +
σ 11 =2a 3 (λ 1 2 -λ 1 -2 λ 2 -2 ) + 4a 4 (λ 1 4 + λ 1 2 λ 2 2 -2λ 1 2 ) + 4a 5 (λ 1 4 -λ 1 -4 λ 2 -4 ) + 4a 6 (2λ 1 2 -λ 1 -2 -λ 1 -4 λ 2 -4 ) + 4a 7 (λ 1 4 -λ 1 2 ) + 4a 9 (λ 1 2 -λ 1 -4 λ 2 -4 ) (3.40) σ 22 =2a 3 (λ 2 2 -λ 1 -2 λ 2 -2 ) + 4a 4 (λ 2 4 + λ 1 2 λ 2 2 -2λ 2 2 ) + 4a 5 (2λ 2 2 -λ 2 -2 -λ 1 -4 λ 2 -4 ) + 4a 6 (λ 2 4 -λ 1 -4 λ 2 -4 ) + 4a 8 (λ 2 4 -λ 2 2 ) + 4a 9 (λ 2 2 -λ 1 -4 λ 2 -4 ) (3.41)
The two above equations will be used in the next section to identify the seven material parameters a 3 , a 4 , a 5 , a 6 , a 7 , a 8 and a 9 by making a comparison between the theoretical and the measured stresses.

3.4.4/ MATERIAL PARAMETERS IDENTIFICATION

To identify the 7 constitutive material parameters that determine the tissue behavior, we have performed a comparison between the stress predicted by our model (equations (3.40) and (3.41) ) and experimental data extracted from the work of Kamenskiy et al. [1]. To assess the quality of the prediction, we have used the classical coefficient of determination R 2 introduced by equation (2.66) in the section 2.5.1. The closest to 1 R 2 is, the best the fit of the experimental data by the theoretical data will be. So we want the ratio S S res S S tot of equation (2.66) to be the closest possible to 0. According to the definitions of S S res and S S tot (equations (2.67) - (2.68) ), we therefore introduce the following objective function F:

F(η) = 10 k=1 n i=1 σ k exp,i -σ k th,i 2 n i=1 σ k exp,i -σ k exp 2 (3.42)
Where η = (a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 ) T represents the set of the 7 material parameters to be identified. The first summation over k (from 1 to 10) corresponds to the ten biaxial tests considered in [1] and described by Table 3 The variable σ k exp,i and σ k th,i included in equation (3.42) respectively represent the i th component of the experimental and theoretical Cauchy stress for case k.

The global minima of F satisfied the first order condition on the gradient of F with respect to η:

∇F(η) = 0 (3.43)
To calculate the gradient of F in a convenient way, we first remark from equations (3.40) and (3.41) that the theoretical Cauchy stress can be expressed in a linear form with respect to η:

σ 11 = B 1 η σ 22 = B 2 η (3.44)
where B 1 and B 2 are the 1 × 7 matrix defined by:

B 1 = 2                                       λ 1 2 -λ 1 -2 λ 2 -2 2(λ 1 4 + λ 1 2 λ 2 2 -2λ 1 2 ) 2(λ 1 4 -λ 1 -4 λ 2 -4 ) 2(2λ 1 2 -λ 1 -2 -λ 1 -4 λ 2 -4 ) 2(λ 1 4 -λ 1 2 ) 0 2(λ 1 2 -λ 1 -4 λ 2 -4 )                                       T ; B 2 = 2                                       λ 2 2 -λ 1 -2 λ 2 -2 2(λ 2 4 + λ 1 2 λ 2 2 -2λ 2 2 ) 2(2λ 2 2 -λ 2 -2 -λ 1 -4 λ 2 -4 ) 2(λ 2 4 -λ 1 -4 λ 2 -4 ) 0 2(λ 2 4 -λ 2 2 ) 2(λ 2 2 -λ 1 -4 λ 2 -4 )                                       T (3.45)
In fact, B 1 and B 2 must be indexed by two integer numbers i and k (see equation (3.42) ) but we have omitted to mention them in equation (3.45) for the sake of simplicity. The index i varies from 1 to n and refers to the number of tested values λ 1 and λ 2 while k varies from 1 to 10 and represents the load case number (see table 3.1). It therefore follows from equation (3.44) that the n stress components of σ k 11 and σ k 22 , corresponding to the k th load case, can be stored as follows:

σ k 11 =               (B 1 ) k 1 . . . (B 1 ) k n               η ; σ k 22 =               (B 2 ) k 1 . . . (B 2 ) k n               η (3.46)
σ k 11 and σ k 22 are vector of dimension n while

              (B 1 ) k 1 . . . (B 1 ) k n               and               (B 2 ) k 1 . . . (B 2 ) k n               are matrix of dimension n × 7.
In an equivalent but more compact form, we introduce the theoretical stress σ k th :

σ k th = A k η (3.47)
where the n×7 matrix A k stands indifferently for

              (B 1 ) k 1 . . . (B 1 ) k n               or               (B 2 ) k 1 . . . (B 2 ) k n              
depending on the considered stress components σ k 11 or σ k 22 . Accounting for this matrix formulation, the objective function F introduced by equation (3.42) can be reformulated by:

F(η) = 10 k=1 σ k exp -A k η 2 σ k exp -σ k mean 2 (3.48)
where . stands for the standard Euclidean norm and σ k exp and σ k mean respectively represent the vector containing all the stress measurements corresponding to the k th load case and the averaged constant vector defined by:

σ k mean = σ k exp               1 . . . 1               (3.49)
The calculation of the gradient of F is straightforward from equation (3.48) :

∇F(η) = 2 10 k=1 A k T A k η -A k T σ k exp σ k exp -σ k mean 2 (3.50)
We are thus faced to a classical linear least squares minimization with a unique solution given by:

η =           10 k=1 A k T A k σ k exp -σ k mean 2           -1           10 k=1 A k T σ k exp σ k exp -σ k mean 2           (3.51)
We thus do not need to discuss the uniqueness of the identified set of material parameters because the single solution is given by equation (3.51) . This remarkable property results from the polynomial form with respect to the invariants that we have selected for the strain energy density (equation (3.24) ), giving a linear dependence of this density with the material parameters. At this stage, it should be noticed that for the same arterial material, Kamenskiy et al. [1] also proposed a constitutive model that includes eight material parameters, but it is mentioned in [1] that there is no guarantee about the uniqueness of these parameters.

The numerical values of η deduced from equation (3.51) , by using the experimental data of Kamenskiy et al. [1] as a reference, are listed on Table 3.2. The assessment of these values is presented in the next section by comparing our model to the one proposed in [1] and by evaluating the quality of the prediction with experimental data. For all these figures, it is observed an excellent agreement between our model and the experimental data. This agreement is confirmed by the calculation of the coefficient of determination R 2 (Table 3.3). It is actually observed that R 2 is very close to 1 with our model, indicating an excellent fit with the experimental data points.

We also note that our model is very close to the numerical results obtained with the model provided by Kamenskiy et al. [1], even if it gives sometimes worse results (SFA: cases 1, 2, 4, 5, 6, 8 and 9; PA: cases 1, 2, 4, 5, 8 and 9; TA: cases 1, 2, 4, 5, 8 and 9). For all the other cases, our model offers a great improvement. That can be explained by the constitutive parts of the two models. The Kamenskiy model is actually based on the combination of four exponential functions, each one corresponding to a four-fiber family, following in that the original concept introduced by Baek et al. [125]. Logically, due to the exponential form of the SEF, the Kamenskiy model is very efficient to predict the behavior of the arterial materials in the large strain range while our model is less efficient in this kind of situation (see for example the right part of the curves plotted on Figure 3.4). But, if we consider the noisy experimental cases (Figures 3.7 left, 3.8 right, 3.11 right, 3.12 left, 3.13 right, 3.15 left, 3.16 right, 3.17 left and 3.18 right), our model appears to be more efficient with a coefficient of determination close to 0.99 while this coefficient takes in the worse situation a value of 0.932 with the Kamenskiy model. This last observation likely results from the fact that the family of invariants used in this paper, which has been proved to be an integrity basis in Ta et al. [53], was built with the material symmetry group S 8 containing all the information relative to the geometrical orientation of the collagen fibers.

We have therefore taken into account all the possible mechanical effects contained in all the possible invariants.

SFA PA TA Load case R 2 [1] R 2 (eq. (3.24) ) R 2 [1] R 2 (eq. (3.24) ) R 2 [1] R 2 (eq. ( 3

3.5/ FINITE ELEMENT IMPLEMENTATION

The aim of this section is to present the finite element implementation of the strain energy density introduced by (3.24) . As in chapter 2, the total Lagrangian formulation is adopted.

To extend the constitutive model from the compressible to the incompressible range, we use the same penalty function W as one introduced by equation (2.106) in section 2.6. By using the additive decomposition introduced by equation (2.105) to separate the anisotropic part to the volumetric part of the SEF, we deduce the Cauchy stress tensor σ from equation (3.21) :

σ = 2F         6 i=1 ∂W ∂L i ∂L i ∂C + ∂ W ∂C         F T (3.52)
Note that we have replaced the extra pressure p (which plays the role of a Lagrange multiplier) in equation (3.21) by a penalty function in order to enforce the incompressibility condition. This replacement allows to reduce the number of unknowns to be determined.

For the sake of simplicity, the derivatives ∂L i ∂C of equation (3.22) are rewritten: The second term included in the bracket of equation (3.52) has been already calculated (equation (2.111) ).

∂L 1 ∂C = M
Because the material must be free of stress if the displacement field is null (F = C = I), we use equation (3.27) in equation (3.53) and we consider this particular case in equation (2.111) :

∂L 1 ∂C = M 11 ; ∂L 2 ∂C = M 22 ; ∂L 3 ∂C = M 33 ∂L 4 ∂C = 4(M 11 + M 22 ) ; ∂L 5 ∂C = 4(M 11 + M 33 ) ; ∂L 6 ∂C = 4(M 22 + M 33 ) ∂ W ∂C = c 2 I (3.55)
Replacing equations (3.26) , (3.29) , (3.30) , (3.32) and (3.55) in (3.52) gives:

σ = 2 (a 3 + 4a 5 + 4a 6 + 2a 9 )(M 11 + M 22 + M 33 ) + c 2 I (3.56)
Or, equivalently, by using equation [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] (see the chapter related to notations and standard results at the beginning of the manuscript):

σ = 2 a 3 + 4a 5 + 4a 6 + 2a 9 + c 2 I (3.57)
From the particular case where the displacement is equal to zero, we then deduce from equation (3.57) that the material parameter c is linked to the other ones by:

c = -2(a 3 + 4a 5 + 4a 6 + 2a 9 ) (3.58)
To construct the tangent stiffness matrix for the analysis of nonlinear structures by the finite element method, one has to determine the stress-strain tangent operator D, which is a fourth order tensor resulting from the derivation of S with respect to C (see equation (1.126) ). In order to calculate D, we first compute the part of the second Piola-Kirchhoff stress tensor S related to W from equation (3.23) :

S = 2 ∂W ∂C = 2         6 i=1 ∂W ∂L i ∂L i ∂C         = 2         6 i=1 ω i ∂L i ∂C         (3.59)
Next, we derive again W with respect to C from equation (3.59) :

D = 4          6 i=1 6 j>i ω i j ∂L i ∂C ⊗ ∂L j ∂C + ∂L j ∂C ⊗ ∂L i ∂C + 6 i=1 ω i ∂ 2 L i ∂C 2 + ω ii ∂L i ∂C ⊗ ∂L i ∂C          (3.60)
The coefficients ω i j stand for the second derivative of W with respect to the invariants L i . They are obtained straightforwardly from equations (3.25) and (3.26) :

ω 11 = ∂ 2 W ∂L 1 ∂L 1 = 2a 7 ; ω 22 = ∂ 2 W ∂L 1 ∂L 2 = 2a 8 ; ω 33 = ∂ 2 W ∂L 3 ∂L 3 = 2a 9 ; ω i j = 0 otherwise (3.61)
To obtain the second derivative ∂ 2 L i ∂C 2 , we derive the first derivatives contained in equation (3.53) with respect to C:

∂ 2 L 1 ∂C 2 = ∂ 2 L 2 ∂C 2 = ∂ 2 L 3 ∂C 2 = 0 ∂ 2 L 4 ∂C 2 = 2(M 11 + M 22 ) ⊗ (M 11 + M 22 ) + 8N 1212 ∂ 2 L 5 ∂C 2 = 2(M 11 + M 33 ) ⊗ (M 11 + M 33 ) + 8N 1313 ∂ 2 L 6 ∂C 2 = 2(M 22 + M 33 ) ⊗ (M 22 + M 33 ) + 8N 2323 (3.62)
where we have introduced the following fourth-order tensors:

N 1212 = M 12 ⊗ M 12 ; N 1313 = M 13 ⊗ M 13 ; N 2323 = M 23 ⊗ M 23 (3.63)
Finally, to achieve the finite element implementation, we need to compute the fourth-order tensor D vol related to the volumetric part of the strain energy density. This computation has been already done in section 2.6 of chapter 2. The result is given by equation (2.125) .

The finite element implementation of the second derivative (equations (2.125) and (3.60) ) of the strain energy densities described by equations (2.105) , (2.106) and (3.24) was realized inside the FER code [57] with C++ language, by following the procedure described on

3.6/ FE SIMULATION RESULTS

In section 3.6.1, in order to validate the finite element implementation, we consider diseased superficial femoral (SFA), popliteal (PA) and tibial arteries (TA) from one patient under planar biaxial extension. These three materials were originally studied in [1]. For all of them, a biaxial stretching with a different ratio of loading was applied to the longitudinal and circumferential directions with the following proportions: 1 : 1, 1 : 2, 1 : 4, 2 : 1 and 4 : 1.

After the validation of the FE implementation, we present in section 3.6.2 a 3D example performed with FER. This 3D example concerns a uniaxial tensile loading involving nonhomogeneous deformation In order to demonstrate the proper implementation of the model inside FER, we make a comparison between the finite element computations and the closed form solution corresponding to equations (3.40) and (3.41) . Comparisons with experimental data extracted from [1] are also performed to prove the capability of the finite element model for predicting the actual behavior of the material. The finite element computations were performed by using the values of material parameters listed in table 3.2. The penalty factor d included in equation (2.106) has been set to a small value of 10 -9 in order to enforce the incompressibility condition. In this section, a tensile test of tibial artery (TA) materials involving non-homogeneous deformations has been processed by considering a rectangular specimen of dimension 10 × 3 × 0.5 mm. The lower part of the specimen is clamped and the upper part is submitted to a fixed displacement varying from 0 to 4 mm. The simulation was conducted in 50 loading steps with a mesh composed of 3200 cubic brick elements and 4305 nodes. Figure 3.35 shows the initial mesh and the distribution of the Von Mises stresses on the deformed mesh. One can observe a symmetrical radial stress distribution at the center of the specimen. One can also observe a decrease of the surface of the section that balances the elongation of the specimen, thus preserving the volume, which is in accordance with the incompressible nature of biological materials described by the proposed model. The proposed SEF adopts a polynomial form (equation (3.24) ) which allows to perform a least square minimization for identifying a single set of material parameters.

Based on our proposed hyperelastic model, the finite element implementation was performed inside the FER university code by using a total Lagrangian approach. All the details of the implementation, that is to say the calculations of the strain and stress incremental forms and the tangent stiffness matrix, are provided in section 3.5.

In order to validate the proposed approach, comparisons were performed between analytical results calculated with our own model and numerical and experimental data obtained by Kamenskiy et al. [1]. Each time, the numerical prediction resulting from the application of the new SEF fits nicely with the data extracted by [1].

It should be additionally underlined that:

• this nice fit was obtained with a large variety of materials (3 different kinds of arteries) and many tests (5 different tests applied to each artery),

• the same set of material parameters was considered for the same artery and applied for the 5 different tests (Table 3.2),

• only 7 material parameters are required in our model against 8 in the Kamenskiy et al. model [1],

• our model provides better results than the Kamenskiy et al. model [1] in the case of noisy experimental data. One possible explanation is the fact that our strain energy density was built by considering material symmetry group including all the properties of invariance of the fiber directions.

• the Kamenskiy et al. model [1] provides better results than ours in the large strain range. One possible explanation is the specific form of the Kamenskiy et al. model which is made of exponential functions allowing to well capture the behavior of the material at large strains.

GENERAL CONCLUSIONS AND FUTURE PROSPECTS GENERAL CONCLUSIONS AND FUTURE

PROSPECTS

This PhD thesis constitutes a first attempt to build some practical strain energy functions (SEF) for modeling different hyperelastic anisotropic materials, including one or four-fibers family by using the theoretical results obtained by Ta et al. [53, 54].

In the first chapter, the foundations of continuum mechanics and of the theory of hyperelasticity were introduced for isotropic and anisotropic materials. The most common strain energy functions based on the classical invariants have been presented, as well as the new invariant proposed by Ta et al. [53, 54]. The concept of polyconvexity, which is essential for ensuring the existence of solutions [41, 70], was next discussed and common polyconvex functions presented. The basics of the finite element method applied to nonlinear structural analysis, such as the total Lagrangian formulation, was finally introduced. It offers the essential realization means for the implementation of the new models developed in the chapters 2 and 3.

In the second chapter, an original strain energy function, based on the transverse anisotropic invariants proposed by Ta et al. [54], was built for modeling anisotropic hyperelastic materials including a one-fiber family. The model is made of new invariants forming an integrity basis derived from the application of the Noether's theorem. Three of these new invariants are well known polyconvex functions while the two others are connected to shear effects.

Two polynomial forms were used (linear first and quadratic next) to represent our strain energy function. To evaluate the relevance of these two polynomial forms, numerical examples were carried out in eight cases: uniaxial tension and shear deformations with a loading direction parallel or transverse to the fiber direction. The predicted results were compared to the experimental data extracted from the work of Ciarletta et al. [39] for two different fiber-reinforced rubbers (soft silicone rubber reinforced by polyamide and soft silicone rubber reinforced by hard silicone rubber). We found that the linear and the quadratic polynomials were not able to well describe the material behavior, particularly with the shear loading transverse to the fiber direction. We have also proved the more general result that any polynomial SEF of any degree does not allow to provide a satisfactory prediction in the case of a shear loading. To overcome this problem, we have added a power-law function to the previous polynomials. When the strain energy function consists in a linear polynomial combined with a power-law function, the predicted results are greatly improved. The quadratic polynomial combined with a power-law function also gives an excellent agreement between the numerical results and the experimental data. But, if the accuracy of the model is not a requirement, the linear option is preferable to the quadratic one because it demands less material parameters to identify (9 against 23).

In the third chapter, an original strain energy function was developed with new polyconvex invariants for modeling four-fibers biometerials. These new invariants are obtained by recombining the 7 invariants originally introduced by Ta et al. [53]. We adopt a polynomial form to express the SEF because it allows to identify a single optimal solution of material parameters through a least square minimization. Accuracy and reliability of the corresponding numerical model were validated by a comparison with experimental and numerical results extracted from [1]. These results concerned diseased superficial femoral (SFA), popliteal (PA) and tibial arteries (TA) from one patient under planar biaxial extension. For each kind of arteries tested with 5 combinations of different biaxial stretch, the predicted results of the proposed model and the experimental data are consistent.

Compared with the model proposed by Kamenskiy et al. [1], which includes 8 material parameters, our model just needs 7 material parameters. The non-linear behavior of these arterial materials can be better described by our model in most cases. But in the range of large stretches, the Kamenskiy model is more efficient than our due to the fact that they have adopted an exponential form, instead of the polynomial form we used, to express the strain energy function.

Based on the combination of the linear polynomial with a power-law function introduced in the second chapter, we have developed a finite element model. This model was implemented in C++ language in the university software FER ( [57]). Following the same strategy, the SEF representing the four-fibers family material built in the second chapter was implemented in FER as well. For both models, several numerical computations performed with FER have demonstrated their efficiency and accuracy. These computations have concerned simple loading cases leading to homogeneous deformation and providing closed-form solutions that are convenient for comparison. But more complex 3D examples involving non-homogeneous deformation were also investigated.

All the achievements completed during this thesis have been published or accepted for publication in the International Journal of Solids and Structures. [65, 67]. Moreover, we have not only built practical SEFs from the theoretical results obtained by Ta et al. [53, 54] for different hyperelastic anisotropic materials, but also proposed their finite element implementations. However, if the practical extension of the mathematical foundations introduced in [54, 53] have been achieved, some additional work could be undertaken and some open-ended questions still remind and will need further investigations:

1. Two of the five invariants used in the second chapter play a key role to predict shear effects but they are not a priori polyconvex functions. The possibility to combine them in order to obtain a polyconvex property is still questionable.

2.

The quality of the prediction of the model introduced in the second chapter, with a large number of monomials and the improvement brought by a power-law function, suggests that it could be possible to replace the proposed SEF by a transcendental function as an exponential for example, allowing in that a drastic decrease of the number of material parameters to identify. But this transcendental function, as well as the corresponding argument to use, is not so simple to guess.

3.

In chapter 3, we have proposed a strain energy function for modeling the behavior of anisotropic hyperelastic material with a four-fibers family by using the invariants introduced by Ta et al. [53]. These invariants were originally intended for modeling the behavior of materials with a two-fibers family but we have explained why it was relevant to use them in the framework of a four-fibers family material. In order to prove the versatility and feasibility of our approach, it could be interesting to apply this approach to the case of a two-fibers family material.

4.

It is remarked that some material parameters identified in chapter 3 (see Table 3.2) are negative. It could therefore be interesting to investigate the convex property of the corresponding SEF as well as the positive definite property of the related tangent stiffness matrix. Accounting for these properties in the identification process could improve the capability of the model to fit experimental data.

5.

In chapter 3, for the arterial materials with a four-fibers family, we just studied the case of a biaxial testing. It could be therefore interesting to test the efficiency of our model for other loading cases, for example in the context of a shear loading.

6.

In terms of non-homogeneous deformations, we have presented in chapters 2 and 3 two 3D FE computations but they are restricted to the case of a uniaxial tensile loading. Now that our models are validated and implemented in a FE code, we could test more complex situations involving for example contact and impact between several hyperelastic bodies.

Document generated with L AT E X and: the L AT E X style for PhD Thesis created by S. Galland -http://www.multiagent.fr/ThesisStyle the tex-upmethodology package suite -http://www.arakhne.org/tex-upmethodology/

Abstract:

This thesis has focused on the construction of strain energy densities for describing the non-linear behavior of anisotropic materials such as biological soft tissues (ligaments, tendons, arterial walls, etc.) or fiber-reinforced rubbers. The densities we have proposed have been developed with the mathematical theory of invariant polynomials, particularly the Noether theorem and the Reynolds operator. Our work involved two types of anisotropic materials, the first with a single fiber family and the second with a four-fiber family. The concept of polyconvexity has also been studied because it is well known that it plays an important role for ensuring the existence of solutions. In the case of a single fiber family, we have demonstrated that it is impossible for a polynomial density of any degree to predict shear tests with a loading parallel and then perpendicular to the direction of the fibers. A linear polynomial density combined with a power-law function allowed to overcome this problem.

In the case of a material made of a four-fiber family, a polynomial density allowed to correctly predict bi-axial tensile test data extracted from the literature. The two proposed densities were implemented in C++ language in the university finite element software FER by adopting a total Lagrangian formulation. This implementation has been validated by comparisons with reference analytical solutions exhibited in the case of simple loads leading to homogeneous deformations. More complex three-dimensional examples, involving non-homogeneous deformations, have also been studied.

Keywords:

Biomechanics, Theory of invariant polynomials, Anisotropic hyperelasticity, Finite element method, Nonlinear mechanics

R ésum é :

Cette th èse a port é sur la construction de densit és d' énergie de d éformation permettant de d écrire le comportement non lin éaire de mat ériaux anisotropes tels que les tissus biologiques souples (ligaments, tendons, parois art érielles etc.) ou les caoutchoucs renforc és par des fibres. Les densit és que nous avons propos ées ont ét é élabor ées en se basant sur la th éorie math ématique des polyn ômes invariants et notamment sur le th éor ème de Noether et l'op érateur de Reynolds. Notre travail a concern é deux types de mat ériaux anisotropes, le premier avec une seule famille de fibre et le second avec quatre familles. Le concept de polyconvexit é a également ét é étudi é car il est notoire qu'il joue un r ôle important pour s'assurer de l'existence de solutions. Dans le cas d'un mat ériau comportant une seule famille de fibre, nous avons d émontr é qu'il était impossible qu'une densit é polynomiale de degr é quelconque puisse pr édire des essais de cisaillement avec un chargement parall èle puis perpendiculaire à la direction des fibres. Une densit é polynomiale lin éaire combin ée avec une fonction puissance a permis de contourner cet obstacle. Dans le cas d'un mat ériau comportant quatre familles de fibre, une densit é polynomiale a permis de pr édire correctement des r ésultats d'essai en traction bi-axiale extraits de la litt érature. Les deux densit és propos ées ont ét é impl ément ées avec la m éthode des él éments finis et en langage C++ dans le code de calcul universitaire FER. Pour se faire, une formulation lagrangienne totale a ét é adopt ée. L'impl émentation a ét é valid ée par des comparaisons avec des solutions analytiques de r éf érence que nous avons exhib ée dans le cas de chargements simples conduisant à des d éformations homog ènes. Des exemples tridimensionnels plus complexes, impliquant des d éformations non-homog ènes, ont également ét é étudi és.

Mots-cl és :

Biom écanique, Th éorie des polyn ômes invariants, Hyper élasticit é anisotrope, M éthode des él éments finis, M écanique non lin éaire

  é de Bourgogne Franche-Comt é UTBM: Universit é de Technologie de Belfort-Montb éliard • CONTINUUM MECHANICS MOTIONS AND CONFIGURATIONS Ω 0 : continuum body in the reference configuration Ω: continuum body in the current configuration X: Lagragian reference coordinate (m) x: Eulerian current coordinate (m) U: displacement vector field (m) V 0 : reference volume (m 3 ) V: current volume (m 3 ) J: volume change between the reference and the current configurations (-) • STRESS TENSORS σ: Cauchy stress (Pa) P: first Piola-Kirchhoff stress or engineering stress (Pa) P T : nominal stress (Pa) S: second Piola-Kirchhoff stress (Pa) CONTENTS • STRAIN AND DEFORMATION TENSORS F: deformation gradient (-) F: volume preserving part of F (-) F vol : dilational part of F (-) C: right Cauchy-Green strain (-) B: left Cauchy-Green strain (-) E: Green-Lagrange strain (-) • MATERIAL SYMMETRY Q: rotation transformation R: reflections transformation M: structural tensor g: symmetry group of the material O(3): group of all orthogonal transformations S O(3): group of all positive orthogonal transformations • INVARIANTS I i (i = 1, 2, 3): classical isotropic invariants related to F I i (i = 1, 2, 3): classical isotropic reduced invariants related to F I 4 and I 5 : classical anisotropic mixed invariants K i (i = 1, ..., 6): integrity basis of invariants for a one-fiber family material H i (i = 1, ..., 7): integrity basis of invariants for a two-fibers family material L i (i = 1, ..., 7): integrity basis of polyconvex invariants for a two-fibers family material • ALL THE NOTATIONS AND OPERATORS RELATED TO LINEAR ALGEBRA ARE DESCRIBED IN THE NEXT SECTION NOTATIONS AND STANDARD RESULTS

Figure 1 . 1 :

 11 Figure 1.1: Configurations and motion of a continuum body

Figure 1 . 2 :

 12 Figure 1.2: Traction vectors acting on infinitesimal surface element with outward unit normals

Figure 1 . 3 :

 13 Figure 1.3: material with one fiber family

Figure 1 . 4 :

 14 Figure 1.4: The material plane of symmetry

CHAPTER 1 .

 1 STATE OF THE ART chapters 2 and 3).

  , 5), b 3 and b 5 are material parameters. I a 4 , I b 4 , I c 4 and I d 4 represent the square of the stretch in the directions a, b, c and d respectively:

) Definition 2 :

 2 convex function (Figure 1.5)

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: One-dimensional convex and non-convex functions

Figure 2 . 1 :

 21 Figure 2.1: A fiber-reinforced material with one fiber family

I 5

 5 is difficult to obtain, contrarily to I 4 .

Figure 2 . 3 :

 23 Figure 2.3: An appropriate orthonormal basis B=(a, b, c) for the strain calculation

Figure 2 . 4 :

 24 Figure 2.4: simple tension test -loading parallel (left) and transverse (right) to the fibers

Figure 2

 2

Figure 2 . 22 :

 222 Figure 2.22: Flow chart of the finite element implementation of the anisotropic part of the model

Figure 2 . 23 :

 223 Figure 2.23: Comparison between finite element, analytical and experimental results (tension tests with material A)

Figure 2 . 24 :Figure 2 . 25 :Figure 2 . 26 :

 224225226 Figure 2.24: Comparison between finite element, analytical and experimental results (shear tests with material A)

Figure 2 . 27 :

 227 Figure 2.27: 3D tension test: mesh

Figure 2 . 28 :

 228 Figure 2.28: Tensile force versus displacement

Figure 2 . 30 :Figure 2 . 31 :

 230231 Figure 2.30: Displacements of points A, B and C

Figure 3 . 1 :

 31 Figure 3.1: A fiber-reinforced material with four-fiber family

Ce 1 , e 1 1 / 2 2 1/ 2 ( 3

 12223 Ce 2 , e

Figure 3 . 2 :

 32 Figure 3.2: shear angle -reference (a) and current (b) configurations

Figure 3 . 3 :

 33 Figure 3.3: Boundary conditions of the biaxial tension test

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: SFA -case 1 (left) and case 2 (right)

Figure 3 . 6 :Figure 3 . 7 :Figure 3 . 8 :Figure 3 . 9 :Figure 3 . 10 :Figure 3 .Figure 3 .Figure 3 .

 36373839310333 Figure 3.6: SFA -case 5 (left) and case 6 (right)

Figure 3 .Figure 3 .Figure 3 .Figure 3 . 18 :

 333318 Figure 3.15: TA -case 3 (left) and case 4 (right)

Figure 3

 3

Figure 3 . 19 :

 319 Figure 3.19: Flow chart of the finite element implementation of the anisotropic part of the strain energy density

Figure 3 .Figure 3 .Figure 3 . 1 •Figure 3 . 25 :Figure 3 . 2 1Figure 3 .Figure 3 .Figure 3 . 1 •Figure 3 .Figure 3 .Figure 3 .Figure 3 . 33 :Figure 3 . 34 : 1 3. 6 . 2 /

 3331325323331333333334162 Figure 3.22: SFA, ratio of loading 2 : 1

Figure 3 .

 3 Figure 3.35: FER simulation -Initial mesh and Von Mises stresses on the deformed mesh

CHAPTER 3 .

 3 A NEW SEF FOR FOUR-FIBER FAMILY MATERIALS

Table 2 .

 2 1: Mechanical properties of the new invariants[54]

Table 2 .

 2 .2.

	Material parameters (MPa)	a 2	a 3	a 4	a 6
	Material A	0.06 -0.035 -0.023 -6.553
	Material B	0.567 -1.549 -1.033 -3309

2: Identified material parameters of the strain energy density W 1 (Eq.

(2.65)

)

Table 2 .

 2 .3 and 2.4.

	Linear terms	a 1	a 2	a 3	a 4	a 6
	Values (MPa) -0.5783 -0.5682 -1.985 -0.1897	48.56
	Coupled terms	a 12	a 13	a 14	a 16	a 23
	Values (MPa)	0.4229	3.58	0.842	194.3	-0.5158
	Coupled terms	a 24	a 26	a 34	a 36	a 46
	Values (MPa) -0.2068 -68.77 0.0687	202.1	114.8
	Squared terms	a 11	a 22	a 33	a 44
	Values (MPa)	0.0654	0.0058	0.4881 -0.1448
	Table 2.3: Identified quadratic material parameters of the strain energy density W 2 (Eq.
		(2.71) ) -Material A
	Linear terms	a 1	a 2	a 3	a 4	a 6
	Values (MPa) -0.5765	3.007	-11.74 -1.823 -159884
	Coupled terms	a 12	a 13	a 14	a 16	a 23
	Values (MPa)	-0.413	9.842	-7.241 163251	-1.638
	Coupled terms	a 24	a 26	a 34	a 36	a 46
	Values (MPa) -0.8456	-3586	3.019	3402	-223.1
	Squared terms	a 11	a 22	a 33	a 44
	Values (MPa)	1.23	-0.6894 0.9358 -1.568

4: Identified quadratic material parameters of the strain energy density W 2 (Eq.

(2.71)

) -Material B

  a 3a 4 + a 12 + a 13a 14 + 2(2a 22 + a 23a 24 ) -(a 24 + a 34 -2a 44 ) (2.79)

Table 2 .

 2 are presented on Figures 2.14, 2.15, 2.16 and 2.17. The values of the identified material parameters related to W 3 and used for the calculations are shown on Tables 2.5 and 2.6. They were identified by following the same procedure as one described in section 2.5.1 5: Identified material parameters of the strain energy density W 1 + W add (Eq.

	Linear terms	a 2	a 3	a 4	a 6
	Values (MPa)	0.0592	0.0478	-0.024 -0.0156
	Power form terms α (MPa)	c 1 (-)	β (MPa)	c 2 (-)
	Values	-0.0846 -984392 -1.859 -31.555
		(2.101) ) -Material A
	Linear terms	a 2	a 3	a 4	a 6
	Values (MPa)	0.5663	4.132	-1.0353	9962
	Power form terms α (MPa)	c 1 (-)	β (MPa) c 2 (-)
	Values	-5.683 -359782 -10233 -8.291

Table 2 .

 2

		0.2									0.2
			Proposed model (Eq (2.37) and Eq (2.96))			Proposed model (Eq (2.45) and Eq (2.96))
			Experiment data [39]							Experiment data [39]
			Model of Fereidoonnezhad et al. [40]				Model of Fereidoonnezhad et al. [40]
		0.15									0.15
	P 11 p (MPa)	0.1	MAT A						P 11 t (MPa)	0.1	MAT A
		0.05									0.05
				tension parallel to the fiber direction			tension transverse to the fiber direction
		1 0	1.1		1.2	λ 1	(-)	1.3	1.4	1.5	1 0	1.1	1.2	1.3 λ 1	1.4 (-)	1.5	1.6	1.7
	Figure 2.14: Comparison between numerical and experimental tensile stresses -linear +
						power form strain enerygy density (2.101)
		0.02	Proposed model (Eq (2.53) and Eq (2.96)) Experiment data [39] Model of Fereidoonnezhad et al. [40]		0.06 0.05	Proposed model (Eq (2.58) and Eq (2.96)) Experiment data [39] Model of Fereidoonnezhad et al. [40]
		0.015
											0.04
	P 12 sp (MPa)	0.01	MAT A					P 12 st (MPa)	0.03		MAT A
											0.02
		0.005
											0.01
		0 0	0.05	0.1 shear parallel to the fiber direction 0.15	0.2	0 0		0.05 shear transverse to the fiber direction 0.1 0.15
						k (-)							k (-)
	Figure 2.15: Comparison between numerical and experimental shear stresses -linear +
						power form strain enerygy density (2.101)
		1 1.2	Proposed model (Eq (2.37) and Eq (2.96)) Experiment data [39] Model of Fereidoonnezhad et al. [40]		0.6	Proposed model (Eq (2.45) and Eq (2.96)) Experiment data [39] Model of Fereidoonnezhad et al. [40]
											0.5
	P 11 p (MPa)	0.6 0.8	MAT B					P 11 t (MPa)	0.3 0.4		MAT B
		0.4									0.2
		0.2									0.1				tension transverse
		1 0	1.02	1.04 λ 1 (-) tension parallel to the fiber direction 1.06 1.08		1 0		1.05	1.1 1 (-) λ to the fiber direction 1.15	1.2
	Figure 2.16: Comparison between numerical and experimental tensile stresses -linear +
						power form strain enerygy density (2.101)

6: Identified material parameters of the strain energy density W 1 + W add (Eq.

(2.101)

) -Material B From Figures 2.14 to 2.17, it can be seen that the predicted results are greatly improved, particular for the shear tests (Figures 2.15 and 2.17) in comparison with the previous linear (Figures 2.7 and 2.9) or quadratic models (Figures 2.11 and 2.13

Table 2 .

 2

	Linear terms	a 1	a 2	a 3	a 4	a 6
	Values (MPa)	-0.1574 -0.0886 -0.3005 -0.0409 -697.7
	Coupled terms	a 12	a 13	a 14	a 16	a 23
	Values (MPa)	-0.1195 0.3649 -0.1899	90.12	-0.063
	Coupled terms	a 24	a 26	a 34	a 36	a 46
	Values (MPa)	-0.1298 -29.68	0.1416	-14.85 -33.57
	Squared terms	a 11	a 22	a 33	a 44
	Values (MPa)	0.03	0.04512	0.3139	-0.139
	Power form terms α (MPa)	c 1 (-)	β (MPa)	c 2 (-)
	Values	-0.0876 -1E11	632.1	-0.0248

If a very significant accuracy is not mandatory, the model introduced by equation (2.101) is sufficient to obtain a satisfactory correlation with measurements (see Figures 2.14 to 2.17 and Tables 2.9 and 2.10, with a R 2 coefficient equal to 0.99 and 0.97, for material A and B respectively). It requires less material parameters: 9 instead of 23. 7: Identified material parameters of the strain energy density W 2 + W add (Eq.

(2.104)

) -Material A

Table 2 .

 2 8: Identified material parameters of the strain energy density W 2 + W add (Eq.

	(2.104) ) -Material B
			tension	shear
	R 2	parallel transverse parallel transverse total
	Linear SEF Eq. ( (2.65) )	0.99	0.99	0.99	0.83	0.99
	Quadratic SEF Eq. ( (2.71) )	1	1	0.99	0.95	1
	Linear+power SEF Eq. ( (2.101) )	0.99	0.99	1	1	0.99
	Quadratic+power SEF Eq. ( (2.104) )	1	1	1	1	1
	Fereidoonnezhad et al. model [40]	0.99	0.99	0.96	0.86	0.99
	Table 2.9: Coefficient of determination R 2 for material A
			tension	shear
	R 2	parallel transverse parallel transverse total
	Linear SEF Eq. ( (2.65) )	0.96	0.98	0.70	0.43	0.84
	Quadratic SEF Eq. ( (2.71) )	0.99	0.99	0.98	0.49	0.87
	Linear+power SEF Eq. ( (2.101) )	0.96	0.98	0.70	1	0.97
	Quadratic+power SEF Eq. ( (2.104) )	0.99	0.99	0.99	1	0.99
	Fereidoonnezhad et al. model [40]	0.86	0.98	0.69	0.98	0.92

Table 2 .

 2 10: Coefficient of determination R 2 for material B

  .5 and 2.6. The comparisons are presented on Figures 2.23 and 2.24 for material A and on Figures 2.25 and 2.26 for material B. From these figures, we can observe that the finite element results match very well the closed form predictions as well as the experimental data. That proves that the finite element model has been properly implemented inside the FER code.

		0.1								0.08
			finite element model							finite element model
		0.08	closed form solution (Eq (2.37)) experimental data [39]			0.07	closed form solution (Eq (2.45)) experimental data [39]
										0.06
	P 11 p (MPa)	0.04 0.06	MAT A						P 11 t (MPa)	0.03 0.04 0.05	MAT A
										0.02
		0.02
				tension parallel to the fiber direction		0.01	tension transverse to the fiber direction
		1 0	1.05	1.1	λ 1	(-)	1.15	1.2	1.25	1 0	1.05	1.1	λ 1	(-)	1.15	1.2	1.25

  Fe 1 ; H 2 = Fe 2 , Fe 2 ; H 3 = Fe 3 , Fe 3(3.3) Unfortunately, there are no clear evidences proving that the four other invariants, namely H 4 , H 5 , H 6 and H 7 , are polyconvex. It is thus mandatory to recombine them in order to make them polyconvex. To do that, we follow the same strategy as the one used to establish the polyconvexity of H 1 , H 2 and H 3 by introducing the polyconvex quantity:C(e 1 + e 2 ), e 1 + e 2 = F(e 1 + e 2 ), F(e 1 + e 2 )(3.4)

  2 C(e 1 + e 2 ), e 1 + e 2 2 + C(e 1 -e 2 ), e 1 -e 2 It is noted that the quantity L 4 introduced by equation (3.8) is invariant as a combination of H 1 , H 2 and H 4 and is also polyconvex as a summation over squared polyconvex functions. Following the same strategies, we introduce two additional polyconvex invariants L 5 and L 6 : L 5 = (H 1 + H 3 ) 2 + 4H 5 ; L 6 = (H 2 + H 3 ) 2 + 4H 6 (3.9)

	2	(3.8)

  3 ) ; Ce 1 , e 1 ; C 2 e 1 , e 1 ; Ce 2 , e 2 ; C 2 e 2 , e 2(3.20)

  + ρ 2 )(e 1 ⊗ e 1 + e 2 ⊗ e 2 ) + 4ρ 4 [e 1 ⊗ e 2 + e 2 ⊗ e 1 ] + ρ 3 )(e 1 ⊗ e 1 + e 3 ⊗ e 3 ) + 4ρ 5 [e 1 ⊗ e 3 + e 3 ⊗ e 1 ] + ρ 3 )(e 2 ⊗ e 2 + e 3 ⊗ e 3 ) + 4ρ 6 [e 2 ⊗ e 3 + e 3 ⊗ e 2 ]

	e 1 ;	∂L 2 ∂C	= e 2 ⊗ e 2 ;	∂L 3 ∂C	= e 3 ⊗ e 3
	∂L 4 ∂C = 2(ρ 1 ∂L 5 ∂C = 2(ρ 1 ∂L 6 ∂C = 2(ρ 2

  .1.

	Ratio of loads applied in the
	longitudinal and circumferential	1 -1	1 -2	2 -1	1 -4	4 -1
	directions
	σ 11 measurement	case 1 case 3 case 5 case 7 case 9
	σ 22 measurement	case 2 case 4 case 6 case 8 case 10
	Table 3.1: Ten different biaxial loading cases [1]

  In the case of the biaxial tension loading applied to the 3 different arteries, the comparisons between the experimental data extracted from the literature[1] and the numerical results from the constitutive model are presented onFigures 3.4 to 3.8 for Superficial Femoral Artery (SFA), on Figures 3.9 to 3.13 for Popliteal Artery (PA) and on Figures 3.14 to 3.18 for Tibial Artery (TA).

		• Superficial Femoral Artery (SFA)
		200					200
		Experimental data					Experimental data
		Proposed model (Eq. (3.21))				Proposed model (Eq. (3.21))
		Model of Kamenskiy et al. [1]				Model of Kamenskiy et al. [1]
		150					150
	σ 11 (Kpa)	100				σ 22 (Kpa)	100
		50					50
		0					0
		1	1.05	1.1 λ 1 (-)	1.15	1.2		1	1.02	1.04 λ 2 (-)	1.06	1.08	1.1
	Material parameters (kPa)	a 3	a 4	a 5		a 6	a 7	a 8	a 9
		S FA		-876.97 105.4 -196.6 106.8 273.8 551.9 269.8
		PA		-2027.9 295.2 172.3 311.3 -136.5 24.4 -109.7
		T A		-3978.9 523.4 239.4	583	-149	137.2 -63.9
	Table 3.2: Identified material parameters for the Superficial Femoral Artery (SFA), the
			Popliteal Artery (PA) and the Tibial Artery (TA)
	3.4.5/ VALIDATION OF THE MODEL

Table 3 . 3

 33

	.24) )

: Coefficient of determination R 2 for 10 different load cases and three different arteries (SFA, PA and TA)

  + ρ 3 )(M 22 + M 33 ) + 8ρ 6 M 23 (3.53) where the symmetric matrix M 11 , M 22 , M 33 , M 12 , M 13 and M 23 are defined by: M 11 = e 1 ⊗ e 1 ; M 22 = e 2 ⊗ e 2 ; M 33 = e 3 ⊗ e 3 M 12 = 1 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) ; M 13 = 1 2 (e 1 ⊗ e 3 + e 3 ⊗ e 1 ) ; M 23 = 1 2 (e 2 ⊗ e 3 + e 3 ⊗ e 2 ) (3.54)

	11 ;	∂L 2 ∂C	= M 22 ;	∂L 3 ∂C	= M 33 ;	∂L 4 ∂C	= 2(ρ 1 + ρ 2 )(M 11 + M 22 ) + 8ρ 4 M 12
	∂L 5 ∂C	= 2(ρ 1 + ρ 3 )(M 11 + M 33 ) + 8ρ 5 M 13 ;	∂L 6 ∂C	= 2(ρ 2

  The comparisons are presented onFigures 3.20 to 3.24 for SFA, on Figures 3.25 to 3.29 for PA and on Figures 3.30to 3.34 for TA. On each figure, the left and right parts correspond to the σ 11 and σ 22 components of the Cauchy stress, respectively. We can observe that the finite element results provide a fair agreement with the closed form predictions as well as with the experimental data. That proves that the finite element model has been properly implemented inside the FER code and is able to well predict the behavior of real materials.

		100	experimental data [1]				200	experimental data [1]
			finite element model					finite element model
		80	closed form solution (Eq (3.37))				closed form solution (Eq (3.38))
										150
	σ 11 (Kpa)	40 60							σ 22 (Kpa)	100
		20								50
		0								0
		1		1.05		λ 1 (-)	1.1		1.15		1	1.02	1.04 λ 2 (-)	1.06	1.08	1.1
						Figure 3.21: SFA, ratio of loading 1 : 2
		250								100
			experimental data [1]					experimental data [1]
			finite element model					finite element model
		200	closed form solution (Eq (3.37))			80	closed form solution (Eq (3.38))
	σ 11 (Kpa)	100 150							σ 22 (Kpa)	60 40
		50								20
										0
		0
		1		1.05	1.1	λ 1 (-)	1.15	1.2	1.25	1	1.01	1.02	1.03 λ 2 (-)	1.04	1.05	1.06
			• Superficial Femoral Artery (SFA)
			200	experimental data [1]				200	experimental data [1]
				finite element model					finite element model
				closed form solution (Eq (3.37))				closed form solution (Eq (3.38))
			150								150
		σ 11 (Kpa)	100							σ 22 (Kpa)	100
			50								50
			0								0
			1	1.05		1.1 λ 1 (-)	1.15		1.2		1	1.02	1.04 λ 2 (-)	1.06	1.08	1.1
							Figure 3.20: SFA, ratio of loading 1 : 1

LIST OF FIGURES 2. 13 Comparison between numerical and experimental shear stressesquadratic strain energy density (equation (2.71)

2.18

Comparison between numerical and experimental tensile stressesquadratic + power form strain energy density (equation (2.104) ) . . . . . . . .

2.19

Comparison between numerical and experimental shear stressesquadratic + power form strain energy density (equation (2.104) ) . . . . . . . .

2.20

Comparison between numerical and experimental tensile stressesquadratic + power form strain energy density (equation (2.104) ) . . . . . . . .

2.21

Comparison between numerical and experimental shear stressesquadratic + power form strain energy density (equation (2.104) ) . . . . . . . .