B Nedjar
email: boumediene.nedjar@enpc.fr

Formulation of a nonlinear porosity law for fully saturated porous media at finite strains

Keywords: Nonlinear porosity law, High pore pressure, Compressible fluid, Poroelasticity, Poroplasticity, Finite strain

In this work, we develop a mathematical formulation of a nonlinear porosity law suitable for finite strain and high pore pressure conditions in porous media. The approach is built around the physical restriction that, by definition, the actual porosity is bounded in the interval [0, 1] for any admissible process. Specifically, the model is motivated by elementary considerations that have been extended to the nonlinear range and, at the limiting case of an infinitesimal approximation, it reaches the porosity law of the classical linear poromechanics. In a next step, the formulation is integrated within the unified framework of continuum thermodynamics of open media which is crucial in setting the convenient forms of the constitutive relations and evolution equations to fully characterize the behavior of porous materials.

Finite strain poroelasticity as well as poroplasticity are considered in this work where, furthermore, a generalized constitutive law for the saturating fluid has been introduced such that both the incompressible fluid and ideal gas are embedded as particular cases. Parametric studies are conducted

Introduction

Porous materials belong to a very important class of materials which are frequently employed in a wide variety of industrial applications, for instance structural applications employing polymeric foams and elastomeric gels. For geomaterials, the knowledge of the poromechanical behavior is very important for the mineral exploration and mining, and in soil mechanics, the so-called consolidation phenomena are perhaps the most studied problems in geomechanics. Moreover, and from the scientific point of view, much interest has been devoted for the modeling of porous media in a wide range of challenging domains such as biomechanics whose study is nowadays of growing interest, see among others the recent review paper [START_REF] Ambrosi | Perspectives on biological growth and remodeling[END_REF].

Most of these applications involve large deformations and/or large variations of the pore pressure level.

Historically, two approaches have been used for the modeling of porous media: mixture theories, see for example [START_REF] Bowen | Compressible porous media models by use of the theory of mixtures[END_REF], and the macroscale Biot's theory, see for example [START_REF] Biot | Theory of finite deformation of porous solids[END_REF]. The former approach is mostly used to model species migration where the mixture equations for mass balance are used in combination with classical equations for linear momentum balance in terms of rule-of-mixture relations for the stress response. As recent examples of application, see [START_REF] Duda | A theory for species migration in a finitely strained solid with application to polymer network swelling[END_REF][START_REF] Baek | Inhomogeneous deformation of elastomer gels in equilibrium under saturated and unsaturated conditions[END_REF] among many others. The present work is based on the latter approach.

Since the pioneering work of Biot, see for example [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Biot | Theory of deformation of a porous viscoelastic anisotropic solid[END_REF], considerable progress has been made in the last decades to develop a concise framework in the domain of poromechanics. Briefly, it describes the evolution of a saturated porous material in terms of the deformation of its solid skeleton part in one hand, and in terms of the distribution of the mass of its fluid part, on the other hand. The resulting boundary value problem consists of a coupling between the balance equation and the mass conservation of the fluid. The reader is referred for example to [START_REF] Coussy | Mechanics of porous media[END_REF][START_REF] Coussy | Poromechanics[END_REF] for a detailed synthesis and extensive discussions concerning these topics.

Several extensions have been proposed in a relevant long literature, among others, partially saturated poromechanics, see for example [START_REF] Borja | Cam-Clay plasticity. Part V: a mathematical framework for three-phase deformation and strain localization analyses of partially saturated porous media[END_REF][START_REF] Coussy | Poromechanics of freezing materials[END_REF][START_REF] Coussy | Unsaturated poroelasticity for crystallization in pores[END_REF], higher order poromechanics, see for example [START_REF] Ito | Effect of pore pressure gradient on fracture initiation in fluid saturated porous media: Rock[END_REF][START_REF] Mroginski | A thermodynamical gradient theory for deformation and strain localization of porous media[END_REF], and finite strain poromechanics in different domains ranging from soil mechanics applications, see for example [START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF][START_REF] Karrech | Poromechanics of saturated media based on the logarithmic finite strain[END_REF], to deformation of soft media, see for example [START_REF] Hong | A theory of coupled diffusion and large deformation in polymeric gels[END_REF][START_REF] Serpieri | Formulation of a finite deformation model for the dynamic response of open cell biphasic media[END_REF]. In parallel, much interest has been devoted to the field of computational poromechanics in order to provide tools for structural finite element simulations, see for example [START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF][START_REF] Borja | Finite element formulation for transient pore pressure dissipation: A variational approach[END_REF][START_REF] Lewis | The finite element method in the static and dynamic deformation and consolidation of porous media[END_REF][START_REF] Ferronato | A fully coupled 3-D mixed finite element model of biot consolidation[END_REF][START_REF] Korsawe | Finite element analysis of poro-elastic consolidation in porous media: Standard and mixed approaches[END_REF][START_REF] White | Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients[END_REF] among many others.

The present work essentially deals with poromechanics at finite strains.

The modeling framework is focused on the porosity law that must be restricted by physical considerations. In fact, by its definition, the actual (Eulerian) porosity is the ratio of the volume of the connected porous space to the total volume of the porous medium. It must then be bounded in the interval [0, 1] for any admissible process, i.e. even for very large deformations and high pore pressure. By means of elementary considerations, a nonlinear porosity law is built and simplifed versions are proposed. Continuum thermodynamics developments are then adapted to take into account the newly proposed porosity law in an unified framework. In particular, the formulation is linked to the Lagrangian porosity as it is better for capturing the change in the porous space. The theoretical developments are conducted for both poroelastic and poroplastic materials consistent with the thermodynamics requirements. For each case, parametric studies highlight the characteristics of the proposed porosity law and its influence on the global response of the porous material. Furthermore, and irrespective to the above developments, a constitutive law for the saturating fluid is also introduced that encompasses the perfect gas and incompressible fluid as particular cases.

An outline of the remainder of this paper is as follows. Section 2 motivates the proposed nonlinear porosity law starting from the nowadays well known porosity law within the linear theory. The elementary considerations are conducted with the current Eulerian porosity that will be linked later on to the material Lagrangian porosity. Alternatively, two simplified versions of the porosity law are presented and commented as well. Next, these developments are embedded into the unified thermodynamics framework of open media, for finite strain poroelasticity in Section 3, and for finite strain poroplasticity in Section 4. In each case, parametric studies on a model example are conducted by means of simulated hydrostatic compression tests.

Finally, conclusions and perspectives are drawn in Section 5. For the reader's convenience, algorithmic details have been appended in two appendices.

Porosity law at finite strains and high pore pressure

Throughout this paper, we consider isotropic porous materials under isothermal conditions where the non occluded pore space is always supposed to be fully saturated by a single fluid phase. In this section, we study the change in porosity at the large strain range with the eventual presence of high pore pressure by taking into account physical restrictions imposed by the definition of the so-called porosity itself. This is accomplished in a straightforward manner by means of elementary considerations. Moreover, in the linear approximation, the resulting nonlinear law must be consistent with the nowadays classical porosity law. As a starting point, we first recall this latter and highlight some of its remarkable properties.

Motivation: the linear porosity law

Denoting by n 0 the initial porosity of the material, its actual value n is related to the actual infinitesimal strain tensor ε and the actual pore pressure p through the following law

n = n 0 + b tr[ε] + p -p 0 Q (1)
where p 0 is the initial pore pressure. The material parameters b and Q are the so-called Biot's coefficient and Biot's modulus, respectively. Here and in all what follows, tr [.] designates the trace operator of second order tensors.

Among others, two remarkable properties can be pointed out from the above expression:

• When tr[ε] = 0, that is when the volumetric strain is held fixed at its initial value, the porosity and the pore pressure are linearly linked by the Biot's modulus as

n -n 0 = p -p 0 Q when tr[ε] = 0.
(2)

• And when the porous medium is fully drained, the porosity and the volumetric strain are linearly linked by the Biot's coefficient as

n -n 0 = b tr[ε] when p = p 0 . (3)
Our first goal is to formulate the counterpart of this constitutive law together with its properties in the finite strain regime with pore pressures ranging from low to high levels.

A nonlinear porosity law at finite strains

Before going on, we consider standard kinematical notations. Let ϕ(., t) ≡ ϕ t (.) denotes the motion in a time interval [0, T ] of a porous continuum body with reference configuration Ω 0 . The material position of the solid skeleton particles are labeled by X ∈ Ω 0 . They are mapped into the current configuration Ω t ≡ ϕ t (Ω 0 ) to the position x = ϕ t (X). The deformation gradient is denoted by F = Dϕ t (X) with J = detF > 0 being its Jacobian.

For a better understanding, the following elementary considerations are conducted using the porosity in the current configuration, i.e. the Eulerian porosity. This latter is still denoted by n where no confusion is made with the infinitesimal theory. By definition, n is the ratio of the volume of the connected porous space to the total volume. Hence, the key observation is that this definition restricts the possible evolutions of n to strictly belong to the interval [0, 1] for all strain ranges and for all pore pressure levels.

Let us analyze the correspondence to the two properties mentioned above in Section 2.1.

Porosity evolution at fixed volumetric strain

Let us imagine an experience where no volumetric strain change is allowed, that is, by imposing J = 1. Then, what is required from the physical point of view is that:

• n = n 0 when p = p 0 . Here p 0 denotes again the initial pore pressure.

• The porosity increases monotonically with increasing pore pressure.

• And, at high pore pressure, the porosity value is limited by its upper bound, i.e. n → 1 when p → ∞.

Such requirements are satisfied by the 'n -p' curve shown in Figure 1 whose expression is given by

n(J = 1, p) = 1 -(1 -n 0 ) exp - p -p 0 Q(1 -n 0 ) (4)
where Q is a material parameter the inverse of which is the tangent to the curve at p = p 0 (see Figure 1). In fact, one can check from the above expression that

∂n ∂p p=p 0 = 1 Q (5)
Moreover, a first order development of the expression (4) near p = p 0 , and use of the result (5), gives the linear relation of Equation ( 2). The parameter Q can then be interpreted as the initial Biot's modulus. Notice that the relation ( 4) is still valid for pore pressures such that p < p 0 . If necessary, it could be adjusted by a prolongation when p ≪ p 0 to ensure that the requirement n ≥ 0 is always satisfied. However, as the applications we have in mind do not reach such situations, this prolongation will not be adopted in this paper.

Porosity evolution under fully drained conditions

Now let us imagine an easier experiment where the porous material is completely drained. Such a situation occurs at very slow loading velocities where the pore pressure does not exceed the ambiant initial pore pressure p 0 .

Then, only the volumetric strain influences the porosity, and what is required from the physical point of view is that:

• n = n 0 when J = 1.

• The porosity decreases with decreasing volumetic strain until the lower bound as n → 0 when J → 0 + .

• And, the porosity increases with increasing volumetric strain until the upper bound as n → 1 when J → ∞.

Such requirements are satisfied by the 'n -J' curves shown in Figure 2 whose expression is given by n(J, p = p 0 ) = h(J) such that

h(J) =      n 0 J m for J ≤ 1, 1 -(1 -n 0 ) exp - n 0 m 1 -n 0 (J -1) for J > 1, (6)
where m > 0 is a material parameter that controls the shape of the porosity evolution. For instance, two curves have been plotted in Figure 2, one with m > 1 and one with m < 1. The (6) 1 and (6) 2 parts of the above function ensure C 1 -continuity at J = 1 where a prolongation is needed since finite dilatations (J > 1) and finite contractions (J < 1) are common in finite strain poromechanics.

At the limiting case of an infinitesimal perturbation with J ≈ 1 + tr[ε],

|tr[ε]| ≪ 1, the expression (6) reduces for both cases where tr[ε] ≥ 0 and

tr[ε] ≤ 0 to n = n 0 + n 0 m tr[ε]. ( 7
)
This latter is to be compared with its counterpart given by Equation (3) within the linear theory. It clearly identifies the product 'n 0 m' as the initial Biot's coefficient, i.e. b ≡ n 0 m, see Figure 2.

Nonlinear porosity Law

We should now be able to build a constitutive law that encompasses the preceding elementary considerations. In this work, this is established as follows: the formula (4), which is valid only for J = 1, is generalized to the more common cases where volumetric changes occur. For this, we simply replace the initial porosity n 0 in Equation ( 4) by the pore pressurefree function h(J) defined in Equation ( 6). The following porosity law is then proposed

0 n 0 1 1 Eulerian porosity n volumetric strain J n 0 m m > 1 m < 1
n(J, p) = 1 -(1 -h(J)) exp - p -p 0 Q(1 -h(J)) (8)
One can check that n(J, p = p 0 ) = h(J) and n(J = 1, p) gives again the relation (4). Hence, all the above requirements are embedded in the Eulerian porosity law (8). As an illustration, a plot of this latter is given in Figure 3 for J ≤ 1 and p ≥ p 0 , i.e. where the porous material is under hydrostatic compression and at low to high pore pressure.

Again, at the limiting case of an infinitesimal approximation with J ≈ the result

n = n 0 + n 0 m tr[ε] + p -p 0 Q (9)
which is exactly the classical linear porosity law given in (1) with b = n 0 m.

Remark 1. The porosity constitutive law (8) can be simplified as follows

n(J, p) = 1 -(1 -h(J)) exp - p -p 0 Q . ( 10
)
In fact, all the physical requirements mentioned above are again fulfilled by this latter law. ✷ Remark 2. For porous materials with high deformability of the solid skeleton and with low-to-moderate pore pressure, another simplification can be proposed. A linearization of the expression (8) with respect to the pore pressure gives

n(J, p) = h(J) + p -p 0 Q (11)
Obviously, this latter leads again to Equation (9) in the linear approximation.

✷ 2.3. Lagrangian porosity and partial stress related to the pore space

In contrast to the Eulerian porosity n which refers to the current volume dΩ t , the change in the porous space is better captured by the Lagrangian porosity denoted by φ and which is defined such that φ dΩ 0 = n dΩ t , where dΩ 0 is the corresponding reference volume with dΩ t = JdΩ 0 . Then, the two porosities are connected by the following relation, see [START_REF] Coussy | Poromechanics[END_REF],

φ = Jn ≡ J n(J, p) (12)
The Lagrangian porosity is found to be the natural state variable that is work-conjugate to the pore pressure p in the sense that, see the developments below in Section 3.1,

φ = - ∂χ por ∂p (J, p) (13)
where χ por (J, p) is the partial free energy that characterizes the response of the saturated pore space excluding any dissipative phenomenon such like plasticity or viscoelasticity, and with the requirements that

χ por (J = 1, p = p 0 ) = 0 and ∂χ por ∂J (J, p = p 0 ) = 0. ( 14
)
The following form of the partial free energy is then adopted to end up with the porosity law (8)

χ por (J, p) = -J (p -p 0 ) -J 1 -h(J) 2 Q exp - p -p 0 Q(1 -h(J)) -1 . (15)
Notice that in the function h(J) given by Equation ( 6), the initial porosity n 0 is evaluated at the undeformed stress-free configuration of the solid skeleton with J = 1. So that by definition (12) one has n 0 = φ 0 where φ 0 is the initial Lagrangian porosity.

Standard arguments also relate the partial second Piola-Kirchhoff stress tensor S por relative to the pore space to the above partial free energy through the following state law

S por = 2 ∂χ por ∂C (J, p) = J ∂χ por ∂J C -1 (16)
where C = F T F is the right Cauchy-Green tensor, the notation (.) T is used for the transpose operator, and where use has been made of the well known kinematic formula ∂J/∂C = 1 2 JC -1 . This law is equivalently written as

τ por = J ∂χ por ∂J 1 and σ por = ∂χ por ∂J 1 (17)
for the corresponding partial Kirchhoff and partial Cauchy stress tensors, respectively, where the stress relations τ por = F S por F T and τ por = Jσ por have been used, and where, and in all what follows, 1 denotes the second order identity tensor.

Remark 3. If use is made of one of the simplified porosity laws (10) or ( 11), the partial free energy relative to the pore space would then be

χ por (J, p) = -J (p -p 0 ) -J 1 -h(J) Q exp - p -p 0 Q -1 or χ por (J, p) = -J h(J) (p -p 0 ) -J (p -p 0 ) 2 2Q (18)
respectively. ✷

Continuum thermodynamics and parametric study

The form of porosity law constructed in the last section must now be embedded within the continuum thermodynamics of open media. In this section, only finite strain poroelasticity is considered. The topic of finite strain poroplasticity will be addressed later on. Moreover, and irrespective of the solid skeleton, the saturating fluid is also characterized here with a generalized constitutive law encompassing both ideal gas and incompressible fluids as particular cases. However, and without being exhaustive, let us first recall nowadays well know notions in poromechanics that are useful for the forthcoming developments, see for instance [START_REF] Coussy | Mechanics of porous media[END_REF][START_REF] Coussy | Poromechanics[END_REF][START_REF] Lewis | The finite element method in the static and dynamic deformation and consolidation of porous media[END_REF] among others.

The fluid part of the porous solid is characterized by the Lagrangian fluid mass content denoted by m f and defined per unit of reference volume dΩ 0 .

It is related to the current fluid mass content per unit of current volume dΩ t as: ρ f n dΩ t = m f dΩ 0 , where ρ f is the current fluid density. One has then

m f = ρ f φ (19)
where use has been made of the relation ( 12). Excluding any volumetric source of mass production, the conservation of fluid mass in the material description is given by ṁf

= -Div Q f (20)
where Q f is the material flow vector of fluid mass, Div(.) is the divergence operator with respect to the material coordinates X and the dot operator (˙)

is the material time derivative. On the one hand, the material vector Q f is related to its spatial counterpart q f via the Piola transform

Q f = JF -1 q f ,
and on the other hand, if Darcy's law is used, the spatial flow vector of fluid mass can be given by

q f = -ρ f k grad p (21)
where k > 0 is the spatial permeability coefficient of the porous medium assumed to be isotropic for simplicity. It is given by k = Jk 0 in terms of the positive permeability k 0 specified per unit of reference volume, see [START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF]. In ( 21), grad(.) refers to the gradient operator with respect to the spatial coordinates x.

Neglecting the gravity effects, the purely mechanical dissipation is given by the following Clausius-Duhem inequality, see for example [START_REF] Coussy | Mechanics of porous media[END_REF][START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF] for more details,

D = S : 1 2 Ċ + µ f ṁf -ψ ≥ 0 ( 22
)
where S is the total second Piola-Kirchhoff stress tensor, ψ is the free energy of the whole porous solid, and µ f is the free enthalpy of the fluid. Observe the conjugate character between this latter and the fluid mass content m f .

Here we consider a barotropic fluid with Moreover, and as in [START_REF] Coussy | Poromechanics[END_REF], let ψ sk be the free energy of the solid skeleton alone. The additive character of the free energy gives

ρ f ≡ ρ f (p),
ψ sk = ψ -m f ψ f ( 23
)
where ψ f = µ f -p/ρ f is the specific free energy of the fluid. That is, ψ sk is obtained by extracting the volumetric free energy of the fluid from the total volumetric free energy ψ. Then, use of this result together with the definition (19) allows to rewrite the dissipation inequality ( 22) as

D = S : 1 2 Ċ + p φ -ψsk ≥ 0 (24)
Observe this time the conjugate character between the pore pressure p and the Lagrangian porosity φ.

Finite strain poroelasticity

In general, the reversible response of the solid skeleton is characterized by a free energy of the form ψ sk = ψ sk (C, φ) where the dependence on the deformation gradient F through C follows from the fundamental principle of material frame indifference with respect to superposed rigid body motions to the current configuration, see [START_REF] Truesdell | The nonlinear field theories of mechanics[END_REF]. However, in order to obtain a state law giving the Lagrangian porosity φ as a function of the pore pressure p and not the reverse, the free energy ψ sk is partially inverted with respect to the couple of conjugate variables (φ, p) by introducing the following free energy potential χ sk as

χ sk (C, p) = ψ sk (C, φ) -pφ (25)
Moreover, and for convenience in the developments below, this latter is in turn additively decomposed in the form

χ sk (C, p) = χ ′ sk (C) + χ por (J, p) (26)
where the partial free energy χ ′ sk stands for the effective response of the drained solid skeleton, the partial free energy χ por relative to the saturated pore space being the one already defined in Section 2.3. Notice that χ ′ sk and χ por are coupled through the deformation gradient of the solid skeleton since

J = detF ≡ (det C) 1/2 .
The mechanical dissipation related to a poroelastic solid skeleton is zero.

Then, replacing ( 25)-( 26) into ( 24) and imposing D = 0 for any admissible state, we end up with the state law of Equation ( 13) for the Lagrangian porosity φ together with the following one for the stress tensor

S = S ′ + S por with S ′ = 2 ∂χ ′ sk ∂C (C) (27)
where S ′ is the effective partial second Piola-Kirchhoff stress tensor relative to the drained solid skeleton while the partial second Piola-Kirchhoff stress tensor S por relative to the pore space is given by Equation ( 16).

Constitutive law of the saturating fluid

If the saturating fluid is an ideal gas, its constitutive equation is given by the well known relation

ρ f = M f R T p (28)
where M f is the molar mass, T the absolute temperature considered constant in this work, and where R the ideal gas constant. In particular, one has

ρ f 0 = M f R T p 0 ( 29
)
for the initial density in the reference stress-free configuration where p = p 0 .

Hence, the following linear relation is deduced from ( 28) and ( 29)

ρ f ρ f 0 = p p 0 (30)
In this work, the constitutive law adoped for the saturating fluid is a generalization of the precedent result. We choose for the fluid density to be given by the following constitutive law

ρ f (p) = ρ f 0 p p 0 g (31)
where the power g is a new fluid parameter such that g ∈ [0, 1] whose influence on the porous material response is investigated in the next section.

Meanwhile, one can immediatly notice that

• For g = 0, the fluid is incompressible with ρ f (p) = ρ f 0 .

• For g = 1, the constitutive law reduces to that of ideal gas.

• And in an intermediate case where 0 < g < 1, the fluid is, say, not ideal, but still compressible.

Fluid compressibility and porosity law study

With the preceding developments at hand, a deeper modeling study can be conducted. In this section, we investigate the influence of the fluid phase compressibility on the response of the whole porous material, on the one hand, and we compare the mechanical response obtained with the original porosity law given by Equation ( 8) with the ones obtained with the simplified versions as introduced in Remarks 1 and 2, on the other hand.

For illustrative purposes, we choose for the drained solid skeleton a compressible hyperelastic model of the neo-Hookean type whose partial free energy as introduced in Equation ( 26) is additively split into a volumetric and deviatoric parts as

χ ′ sk (C) = κ 2 1 2 (J 2 -1) -log J + µ 2 C : 1 -3 (32)
where κ > 0 and µ > 0 are the bulk and shear moduli, respectively, and where C = J -2/3 C is the volume-preserving right Cauchy-Green tensor, see for example [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Holzapfel | Nonlinear solid mechanics. A continuum appraoch for engineering[END_REF].

We consider undrained hydrostatic compression tests. In this case, from (17) 2 and ( 27), the total hydrostatic Cauchy stress, denoted here by ̟, is

given by (see Equation A.4),

̟(J, p) = κ 2 J - 1 J + ∂χ por ∂J (J, p) (33)
which is independent of the shear modulus µ, and where the expression of the partial free energy χ por (J, p) has been given in Section 2.3.

For the fluid part, since no diffusion is allowed, the conservation of fluid mass (20) reduces to ṁf = 0 (34)

where, with (19), ( 12) and ( 31), one has m f ≡ m f (J, p) = ρ f (p)Jn(J, p).

In summary, this simple problem is given by Equations ( 33)-( 34). The loading is applied at prescribed volumetric strain. For each given J, Equation ( 34) is solved first for the pore pressure p. Then, the ordered pair (J, p) is replaced in Equation ( 33) to get the corresponding hydrostatic stress. The theoretical and algorithmic details are fully developed in Appendix A.

For the solid skeleton, we fix the bulk modulus to κ = 0.164 M P a, which is of the same order as the elastic properties of polyurethane foams at room temperature, see for example [START_REF] Zhang | Constitutive modelling and material characterization of polymeric foams[END_REF][START_REF] Zhang | Constitutive modeling of polymeric foam material subjected to dynamic crash loading[END_REF]. We also choose to fix the initial porosity to φ 0 ≡ n 0 = 0.6 and the initial pore pressure to p 0 = 0.1 M P a. One can observe that the more the parameter g is low the more the extrastress due to the fluid pore pressure is high. In fact, the fluid approaches incompressibility as g → 0. As a comparison, the curves are superposed with the one of the fully drained response with p = p 0 .

Figure 5 shows the corresponding evolutions of the Eulerian porosity n predicted by the porosity law (8). One can observe that, depending on the compressibility of the saturating fluid, the porosity decreases at low pore pressure then increases for high pore pressure never exceeding the upper bound limit n = 1. The cases of almost fluid incompressibility, in these simulations for g = 0.2 and to a lesser extent for g = 0.4, represent in fact extreme situations. The porosity evolution for the fully drained test depends solely on the volumetric strain as n = n(J, p = p 0 ) ≡ h(J), see Equation ( 6),

and is given by a straight line as m = 1 in these simulations. In Section 2.2 we have introduced simplified versions of the original porosity law (8); a first simplified law given by Equation (10), see Remark 1, for which the partial free energy χ por is given by ( 18) 1 , and a second simple law given by Equation ( 11), see Remark 2, for which χ por is given by ( 18) 2 . The corresponding evolutions of the Eulerian porosity n are plotted in Figure 7. All the porosities are increasing even for low pore pressure. This is because the initial Biot's modulus Q is here ten times lower than the one used for the computations of Figures 4 and5. The curve of the simplified porosity law (10) in dashed line is here again close to the response of the original model. However, for the simple law of Equation ( 11), the Eulerian porosity is not bounded as it exceeds the value 1 for high pore pressure. This confirms the comment of Remark 2 that this latter law is adapted only for low-to-moderate pore pressure levels.

We can conclude that the simplified porosity law given by Equation ( 10) is also well adapted for finite strain poroelasticity and high pore pressure. It also fulfills the physical requirements as mentioned earlier in Remark 1.

An example with transient fluid flow

As an illustration, this section is devoted to hydrostatic compression tests allowing fluid flow so that the fluid mass conservation, Equation (20), is no longer reduced as for the undrained case and it greatly influences the overall response through the Darcy law, Equation ( 21). The material we consider is poroelastic as given above in Section 3.3.

A cubic sample is chosen with dimensions (100 × 100 × 100) mm 3 . The initial pore pressure is set to p 0 = 0.1 M P a and the same value is prescribed on the whole faces of the specimen during the loading process as a boundary condition for the fluid part of the coupled problem. The loading consists on imposing the same displacements on the faces so as to control the prescribed volumetric strain rate J.

For the material parameters, we choose (see the parameters that led to the results of Figure 4): κ = 0.164 M P a, µ = 0.246 M P a, Q = 1 M P a, φ 0 = 0.6, m = 1, g = 0.6 and the material permeability is fixed in our computations to k 0 = 100 mm 2 /M P a s. and 10 s -1 . As expected, the more the loading velocity is high the more the extra-stress due to the fluid pore pressure is high. At the limiting case, the response is always bounded by the one obtained under undrained conditions. For the sake of comparison, the curves of the fully drained and the undrained bench test of Section 3.3 have been superposed in the figure. Last but not least, there is qualitatively and, to a lesser extent, quantitatively a good accordance with the experimental results of [START_REF] Zhang | Constitutive modelling and material characterization of polymeric foams[END_REF][START_REF] Zhang | Constitutive modeling of polymeric foam material subjected to dynamic crash loading[END_REF] for polyurethane foams.

Finite strain poroplasticity modeling

The approach to poroplasticity starts with the nowadays well-accepted multiplicative decomposition of the solid skeleton's deformation gradient F into an elastic part F e and a plastic part F p : F = F e F p . For the pore space, we adopt the usual additive decomposition of the Lagrangian porosity φ into a reversible part denoted by φ e and an irreversible part denoted by

φ p φ = φ e + φ p (35)
This latter is equivalent to the additive decomposition of the Lagrangian fluid mass content m f as adopted in [START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF].

The free energy of the solid skeleton is this time of the following general form:

ψ sk ≡ ψ sk (C e , ξ, φ e ) ( 36
)
where C e = F e T F e is the elastic right Cauchy-Green tensor and where, for simplicity, we have considered a single scalar strain-like plastic variable ξ.

Notice that with this choice, ψ sk depends on the Lagrangian porosity solely through its reversible part φ e .

For later use, with the help of the following familiar stress-power relationship

S : 1 2 Ċ = S : 1 2 Ċe + C e S : ℓ p (37)
the Clausius-Duhem dissipation inequality given in Equation ( 24) is here equivalently written as

D = S : 1 2 Ċe + C e S : ℓ p + p φe + p φp -ψsk ≥ 0 ( 38
)
where S = F p SF p T is the total second Piola-Kirchhoff stress tensor with respect to the intermediate configuration defined by F p , and where ℓ p = Ḟ p F p -1 is the left-rate plastic distorsion tensor which is the work-conjugate of the Mandel stress tensor C e S, see for example [START_REF] Klinkel | An anisotropic fibre-matrix material model at finite elastic-plastic strains[END_REF][START_REF] Nedjar | An anisotropic viscoelastic fibre-matrix model at finite strains: Continuum formulation and computational aspects[END_REF][START_REF] Nedjar | On a continuum thermodynamics formulation and computational aspects of finite growth in soft tissues[END_REF] for similar developments.

Constitutive equations and plastic flow

Partial inversion of the solid skeleton's free energy with respect to the couple of conjugate fields (φ e , p) introduces the following dual free energy potential χ sk given by χ sk (C e , ξ, p) = ψ sk (C e , ξ, φ e ) -pφ e (39)

As for finite poroelasticity, see Section 3.1, this latter is in turn additively decomposed for convenience into an effective part χ ′ sk relative to the fully drained solid skeleton and a part χ por relative to the pore space as

χ sk (C e , ξ, p) = χ ′ sk (C e , ξ) + χ por (J e , p) (40)
where J e = detF e ≡ (detC e ) 1/2 > 0 is the elastic Jacobian.

Inserting ( 39)-( 40) into (38) and using standard arguments of continuum thermodynamics, see [START_REF] Coleman | Thermodynamics with internal variables[END_REF][START_REF] Germain | Continuum thermodynamics[END_REF], we end up with the following state laws

S = 2 ∂χ ′ sk ∂C e S ′ + J e ∂χ por ∂J e C e -1

S por

and φ e = -∂χ por ∂p (41)

where use has been made of the kinematic relation ∂J e /∂C e = 1 2 J e C e -1 . The dissipation inequality (38) takes then the following reduced form

D = C e S : ℓ p + q ξ + p φp ≥ 0 ( 42
)
where q is the stress-like plastic variable in the sense that q = -∂ ξ χ ′ sk .

The restriction to isotropy is assumed on the intermediate configuration.

Consequently, the free energy potential χ sk depends on the invariants of its + J e ∂χ por ∂J e 1 τ por (43) see also [START_REF] Truesdell | The nonlinear field theories of mechanics[END_REF][START_REF] Chadwick | Continuum mechanics. Concise theory and problems[END_REF] for the well-known result in isotropic elasticity concerning the above expression of the effective Kirchhoff stress tensor τ ′ . Moreover, the reduced dissipation inequality ( 42) is expressible in the following equivalent form

D = τ : - 1 2 (£ v b e ) b e -1 + q ξ + p φp ≥ 0 (44)
where £ v b e is the Lie derivative of b e = F C p -1 F T defined as

£ v b e = F d dt C p -1 F T (45)
and where C p = F p T F p is the plastic right Cauchy-Green tensor of the solid skeleton. The first term of Equation ( 44) is nowadays classical in finite strain isotropic anelasticity, see for example [START_REF] Simo | Numerical analysis and simulation of plasticity[END_REF][START_REF] Holzapfel | Nonlinear solid mechanics. A continuum appraoch for engineering[END_REF]Nedjar, 2002a[START_REF] Nedjar | An anisotropic viscoelastic fibre-matrix model at finite strains: Continuum formulation and computational aspects[END_REF] for similar developments. Now let G denotes the yield criterion when a single-surface is considered for simplicity. In view of the expression ( 44), its natural arguments are τ , q and p in the general case, that is, G ≡ G(τ , q, p). And for the plastic flow, we choose the following evolution equations

         £ v b e = -2γ ∂ τ F b e ξ = γ ∂ q F φp = γ ∂ p F ( 46
)
where γ is the plastic consistency parameter satisfying the Kuhn-Tucker loading/unloading conditions: γ ≥ 0, G ≤ 0 and γG = 0, and where

F ≡ F(τ , q, p
) is a plastic flow potential. In particular one has F = G for poroplasticity. The evolution equation ( 46) 3 for the irreversible porosity is equivalent to the one used in [START_REF] Armero | Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions[END_REF] for the irreversible fluid mass content.

Remark 5. The partial free energy relative to the pore space χ por (J e , p) is such that χ por (J e = 1, p = p 0 ) = 0 and ∂χ por ∂J e (J e , p = p 0 ) = 0 (47) Furthermore, χ por can also depend on the plastic Jacobian J p = detF p .

In this case, this latter acts as a fixed parameter. ✷ Remark 6. From the constitutive law ( 43), the partial Cauchy stress tensor relative to the pore space is given by

σ por = 1 J p ∂χ por ∂J e 1 ( 48
)
where use has been made of the important kinematic relation J = J e J p for the Jacobian of the solid skeleton. ✷ Remark 7. From the additive decomposition ( 35) and the relation ( 12), the Lagrangian partial porosities φ e and φ p are linked to their Eulerian counterparts n e and n p , respectively, by the relations φ e = Jn e and φ p = Jn p (49)

Note that one has also the additive decomposition n = n e + n p . Furthermore, whether the reversible Eulerian porosity is of the form n e (J e , p) or n e (J, p) is a modeling matter irrespective of the continuum thermodynamic developments that led to the state law (41) 2 . This point is, among others, discussed in the next sections. ✷

A modeling example: Effective stress poroplasticity

The particularity in poromechanics is that there is a class of behaviors that can be characterized solely via the fully drained part of the stresses.

Besides on soil mechanics where this class has been intensively used, such a characterization can also be applied to model other porous materials. For instance, experimental tests on some polymeric foams reveal that the plastic flow is independent of the pore pressure level, see for example [START_REF] Zhang | Constitutive modelling and material characterization of polymeric foams[END_REF][START_REF] Zhao | Testing of polymeric foams at high and medium strain rates[END_REF] among others. Such models are considered in this section as a modeling example for a deeper discussion on the extension possibilities to the plastic range of the porosity law developed in the earlier sections.

Accordingly, the yield criterion is expressed in the form G(τ ′ , q), and by extension, the plastic flow potential is analogously expressed as F(τ ′ , q). The evolution equations ( 46) reduce then to

         £ v b e = -2γ ∂ τ ′ F b e ξ = γ ∂ q F φp = 0 (50)
Hence, one of the consequences of effective stress poroplasticity is that φ p = 0 and the total porosity is reversible during the whole loading history

φ = φ e ⇒ n = n e (51)
The extension of the nonlinear porosity law developed in Section 2 to the finite strain poroplastic range is not unique but is a modeling matter. In this work we propose two alternatives among others:

• The reversible Eulerian porosity depends on the elastic Jacobian. That is, by using the porosity law given by Equation ( 8) where the function h of Equation ( 6) is evaluated with the argument J e instead of J,

n e = 1 -(1 -h(J e )) exp - p -p 0 Q(1 -h(J e )) (52)
This latter will be denoted by n e (J e , p) in the following.

• Or the reversible Eulerian porosity depends on the total Jacobian. That is by using the porosity law of Equation ( 8) as it is,

n e = 1 -(1 -h(J)) exp - p -p 0 Q(1 -h(J)) (53)
This latter will be denoted by n e (J, p) to avoid any confusion with the precedent one.

The subtle difference between the two porosity laws is better highlighted under fully drained conditions. In this case n e (J e , p = p 0 ) = h(J e ) which means that n e takes the value of the initial porosity n 0 after elastic unloading irrespective of the volumetric plastic strain. This constrasts with n e (J, p = p 0 ) = h(J) which means that n e depends on the volumetric plastic strain though J. See the sketch of Figure 9 for an illustration.

To end up with the porosity law of Equation ( 52), the following partial free energy relative to the pore space is adopted, see the analogy with the expression of Equation ( 15),

χ por (J e , p) = -J e J p (p -p 0 ) -J e J p 1 -h(J e ) 2 Q exp - p -p 0 Q(1 -h(J e )) -1 (54)
And for the porosity law of Equation ( 53), again the same expression as in Equation ( 15) is used, always keeping in mind the relation J = J e J p .

J = J e = J p = 1 J < 1 J e = 1 J p = J n 0 n e = n 0 n e < n 0 n e = h(J e ) > n e = h(J)
Figure 9: Fully drained porous material after plastification. Comparison of the resulting Eulerian porosity after elastic unloading given by the two variants of porosity law.

For both porosity laws, the derivative of the function χ por (J e , p) with respect to J e is computed at fixed value of the plastic Jacobian J p when evaluating the corresponding partial stress tensors τ por or σ por relative to the pore space with the constitutive relations ( 43) or (48), respectively.

Remark 8. For non effective stress poroplasticity with evolving irreversible porosity, i.e. with φp = 0, additional restrictions must be taken into account when elaborating a porosity law. The (Lagrangian) state law of Equation ( 41) 2 and the evolution equation of Equation ( 46) 3 have to be such that

n e ∈ [0, 1], n p ∈ [0, 1],
and

n = n e + n p ∈ [0, 1] (55)
for any admissible process. This needs further theoretical developments that will not be addressed in this work. ✷

Porosity law study

As in Section 3.3, we consider undrained hydrostatic compression tests where no (time-dependent) transient effects take place for the sake of easy n = n e (J, p). The fluid density is still given by the constitutive law (31). The undrained hydrostatic compression problem at hand is solved for prescribed total volumetric strain history. For each given J, the plastic flow is solved first for the internal variables J e and ξ. Next, the fluid mass conservation is solved for the pore pressure p, and finally the total hydrostatic Cauchy stress ̟ is deduced by mere function evaluation. The full algorithmic details are developed in Appendix B.

To make matters concrete, we fix the bulk modulus to κ = 2 M P a and the plastic material parameters to: σ y = 0.07 M P a, σ ∞ y = 0.145 M P a, δ = 40 and H = 0.05 M P a for the isotropic hardening, and α = 1 for the effective hydrostatic stress factor as no deviatoric contributions are present in this test. These parameters are such that the drained mechanical response of the solid skeleton is of the same order as for drained polypropylene foams at room temperature, see for example [START_REF] Zhang | Constitutive modelling and material characterization of polymeric foams[END_REF].

Figure 10 gives a set of simulated hydrostatic compression stresses versus volumetric strain computed with the two porosity laws n e (J e , p) and n e (J, p).

The porosity parameters are fixed and set to, see Equations ( 52) and ( 53):

Q = 1 M P a, m = 1 and φ 0 = 0.8, and the initial pore pressure is set to p 0 = 0.1 M P a. Two values have been used for the fluid density law parameter: g = 1 for an ideal gas, and g = 0.6 for a general compressible fluid.

One can observe the similar but quantitatively different results given with the use of the two porosity laws. Irrespective of the fluid compressibility through the parameter g, the use of the porosity law n e (J e , p) gives lower extra-stress than with the porosity law n e (J, p) at low fluid pore pressure. This tendency in reversed at high pore pressure. The curves are superposed with the one of the fully drained response with p = p 0 .

The difference between the two porosity laws is more noticeable from the corresponding evolutions of the Eulerian porosity n plotted in Figure 11. As the plastic flow starts at relatively low effective stress level, the elastic Jacobian J e is almost stabilized even when the total Jacobian J still decreases. This is due to the plastic part of the solid skeleton behavior as, with the selected plastic parameters, there is a low level of isotropic hardening. Hence, the law n e (J, p) gives lower Eulerian porosity than the law n e (J e , p). This difference is more evident for the drained tests also plotted in Figure 11. this time by extending one of the simplified porosity laws ( 10) or ( 11) instead of the one given by Equation ( 8). However, for the sake of clarity, this will not be done in this paper. ✷

An example with transient fluid flow

In this last part, a similar example as in Section 3.4 is presented with, this time, a poroplastic material in two dimensions under plane strain conditions.

A square sample is chosen with dimensions (100 × 100) mm 2 with initial pore pressure set to p 0 = 0.1 M P a. This latter is prescribed on the whole sides as a boundary condition for the fluid part. The mechanical loading consists on imposing the same displacement on the four sides so as to control the prescribed volumetric strain rate J.

For the drained solid skeleton, we fix the material parameters to κ = 1.417 M P a, µ = 0.654 M P a for the hyperelastic behavior, and σ y = 0.12 M P a, σ ∞ y = 0.18 M P a, δ = 40, H = 0.01 M P a, and α = 1 for the plastic flow. For the porous space and fluid behavior we fix the parameters to:

Q = 1 M P a, φ 0 = 0.7, m = 1, g = 0.4 and the material permeability k 0 = 100 mm 2 /M P a s. Figure 12 gives a set of simulated resultant hydrostatic Cauchy stress versus volumetric strain computed for different loading velocities, here with J = 0.1, 1 and 10 s -1 . Again, and as expected, the more the loading velocity is high the more the extra-stress due to the fluid pore pressure is high, and the response is always bounded by the one obtained under undrained conditions, this latter being superposed in the figure. Moreover, there is qualitatively and, to a lesser extent, quantitatively good accordance with the experimental results of [START_REF] Zhang | Constitutive modelling and material characterization of polymeric foams[END_REF][START_REF] Zhao | Testing of polymeric foams at high and medium strain rates[END_REF] for polypropylene foams.

For illustrative purposes, Figure 13 shows a typical deformed mesh and pore pressure field, for instance here at J = 0.4 for the computation with the volumetric strain rate J = 0.1 s -1 . One can observe the locally non uniform character of the response due to the transient fluid flow.

where the constant k 0 is still the initial permeability coefficient in the undeformed stress-free configuration, the power coefficient q ≥ 0 being a new material parameter.

As an illustration, Figure 14 shows the results of the same computations as the ones that led to the results of Figure 12, this time with the permeability law (59) where the value q = 2 has been used. However, further theoretical investigations and parametric studies are needed which are out of the scope of this paper. ✷

Conclusion and perspectives

In this paper, a nonlinear porosity law formulation has been presented that allows the combined influence of finite strain and high pore pressure in poromechanics.

Using elementary considerations, we have built a sound theory that takes into account the fact that the actual (Eulerian) porosity has to be bounded by the interval [0, 1] for any admissible state. This physical restriction is the consequence of the definition of the porosity itself as being the ratio between the volume of the porous space over the total volume of the porous solid. Furthermore, the presented porosity law reaches the classical one when linearized for the limiting case of an infinitesimal approximation.

By means of the continuum thermodynamics of porous media, we have embedded the above law for the general case of finite strain poroelasticity and poroplasticity leading to a concise way to characterize the whole set of constitutive relations and evolution equations to be appended to the coupled balance equation and the mass conservation of the fluid for a global structural resolution modeling. In particular, we have also proposed a generalized fluid constitutive law that encompasses both the incomplessible fluid and ideal gas as particular cases.

Parametric studies have been conducted and commented by means of hydrostatic compression tests. These tests are of course simple, but they have the merit to highlight the intrinsic characteristics of the present modeling.

In particular for poroplasticity, the examples have been restricted to effective stress models with no irreversible porosity. The extension to take into account evolving irreversible porosity is the object of a future work.

Appendix A. Poroelastic undrained hydrostatic compression

For the solid skeleton, a compressible model of the neo-Hookean type is chosen. The partial free energy function χ ′ sk in ( 26) is additively split into a volumetric part and an isochoric part as

     χ ′ sk (C) = U (J) + µ 2 C : 1 -3 U (J) = κ 2 1 2 (J 2 -1) -log J (A.1)
where κ > 0 and µ > 0 are the bulk and shear moduli, respectively, and C = J -2/3 C is the volume-preserving right Cauchy-Green tensor, see for Given (J n , p n ) at t = t n , and given J n+1 at t = t n+1 , then

1. Initialize: k = 0, p (0) n+1 = p n 2. Evaluate residual and check convergence r (k) n+1 = -g n (k) n+1 log p (k) n+1 p n -n (k) n+1 log J n+1 J n - ∂n ∂J (k) n+1 (J n+1 -J n ) - ∂n ∂p (k) n+1 (p (k) n+1 -p n ) IF |r (k) n+1 | > T OL THEN go to Step 3 ELSE set p n+1 = p (k)
n+1 and EXIT.

3. Compute the tangent modulus and the increment

D (k) n+1 = 1 + log J n+1 J n + g log p (k) n+1 p n ∂n ∂p (k) n+1 + g n (k) n+1 p (k) n+1 +(J n+1 -J n ) ∂ 2 n ∂J∂p (k) n+1 + (p (k) n+1 -p n ) ∂ 2 n ∂p∂p (k) n+1 ∆p (k) n+1 = r (k) n+1 D (k) (n+1)
4. Update the pore pressure: p

(k+1) n+1 = p (k) n+1 + ∆p (k) n+1
Set k ←-k + 1 and return to Step 2.

where b e = J e -2/3 b e , κ and µ are the bulk and shear moduli, respectively, and the function H(ξ) characterizes the isotropic plastic hardening.

The effective stress yield criterion we use is given by Equation (57):

G(τ ′ , q) = devτ ′ + α|p ′ | - 2 3 (σ y -q(ξ)) (B.2)
The effective hydrostatic Kirchhoff stress p ′ = 1 3 [τ ′ : 1] is given in our case by

p ′ = J e ∂U ∂J e (J e ) = κ 2 J e 2 -1 (B.3)
For hydrostatic loadings, the tensorial stress and strain responses can be rewritten solely in terms of volumetric quantities. In particular, the evolution equation ( 50) 1 is equivalent to the following one, see [START_REF] Simo | Numerical analysis and simulation of plasticity[END_REF] for details,

d dt [J p ] = γ ∂ τ ′ G : 1 (B.4)
where associated flow with F(τ ′ , q) = G(τ ′ , q) is considered.

Within a typical time interval [t n , t n+1 ], we assume known the variables {J n , J e n , ξ n } at time t n , J p n being also known since J p n = J n /J e n . The objective is to approximate the constrained problem given by Equations (B.4), (50) 2 appending the loading/unloading conditions γ ≥ 0, G ≤ 0 and γG = 0 to update the variables {J e n , ξ n } → {J e n+1 , ξ n+1 } for a prescribed total volumetric strain J n+1 at time t n+1 . This is accomplished via the well known elastic predictor/plastic corrector algorithmic concept.

For the elastic prediction, the trial state is evaluated with the following elastic volumetric strain

J e tr n+1 = J n+1 J p n (B.5)
which is used to compute the trial yield criterion

G tr n+1 = α |p ′ tr n+1 | - 2 3 (σ y -q(ξ n )) (B.6)
where p ′ tr n+1 = p ′ (J e tr n+1 ) and where use has been made of the fact that the deviatoric stress vanishes for hydrostatic loading. Then, (i) if G tr n+1 ≤ 0, the trial state is admissible and we set J e n+1 = J e tr n+1 and ξ n+1 = ξ n , (ii) if G tr n+1 > 0, the trial state is not admissible and a correction has to be performed.

Introducing for convenience the notations ϑ = log J, ϑ e = log J e , ϑ p = log J p and so on for similar quantities with indices, a backward-Euler scheme on Equations (B.4) and ( 50 where ∆γ = (t n+1 -t n )γ is a plastic multiplier to be determined by enforcing the consistency condition G n+1 = G(p ′ (ϑ e n+1 ), q(ξ n+1 )) = 0. Furthermore, noticing that ϑ n+1 = ϑ e n+1 + ϑ p n+1 and that ϑ n+1 -ϑ p n = log J e tr n+1 ≡ ϑ e tr n+1 , the discrete equation (B.7) 1 is conveniently rewritten as The nonlinear algebraic sub-problem in the primary variables {ϑ e n+1 , ξ n+1 } is solved by means of the scheme summarized in Table B.2. After resolution, one has J e n+1 = exp ϑ e n+1 for the updated elastic Jacobian.

For the fluid part, the conservation of fluid mass (20) reduces to ṁf = 0 which still expands in the form given by Equation (A.5). If the porosity law ( 52) is used, then n = n e (J e , p) and the equivalent form analogous to Equation (A.7) is given by g n e log p + n e log J + ∂n e ∂J e Je + ∂n e ∂p ṗ = 0 (B.9)

If the porosity law ( 53) is used, then n = n e (J, p) and the equivalent form remains as given by Equation (A.7). For both cases, the discrete version of the fluid mass conservation is solved for the pore pressure p n+1 with the help of the procedure summarized in Table A.1. This latter has to be slightly adapted when the porosity law n e (J e , p) is to be used, J e n+1 being known from the above plastic resolution.

Finally, the total hydrostatic Cauchy stress ̟ n+1 at time t n+1 is computed with the help of the constitutive relation If plastic loading occurs (G tr n+1 > 0), then perform the following correction 1. Initialize: k = 0, ϑ e (0) n+1 = ϑ e tr n+1 , ξ

̟ = 1 J κ 2 J e 2 -
n+1 = ξ n , ∆γ (0) = 0 2. Check yield condition and evaluate residuals Set k ←-k + 1 and return to Step 2.

p ′ (k) n+1 = p ′ (ϑ e (k) n+1 ) (k) n+1 = q(ξ (k) n+1 ) G (k) n+1 = G(p ′ (k) n+1 , q (k) n+1 ) R (k) n+1 =        -ϑ e (k) n+1 + ϑ e tr n+1 -∆γ (k) α p ′ (k) n+1 |p ′ (k) n+1 | -ξ (k) n+1 + ξ n + 2 3 ∆γ (k)        IF R (

Figure 1 :

 1 Figure1: Porosity evolution at fixed volumetric strain (with J = 1).

Figure 2 :

 2 Figure 2: Porosity evolution in a fully drained porous medium (with p = p 0 ).

1Figure 3 :

 3 Figure 3: Porosity law n(J, p). An illustration for J ≤ 1 and p ≥ p 0 .

  see Section 3.2 below, and such that 1/ρ f = ∂µ f /∂p. This latter fact has been taken into account in establishing the Darcy's constitutive relation (21). In (22), and in all what follows, the double dot symbol ':' is used for double tensor contraction, in particular, one has tr[(.)] = (.) : 1.

Figure 4 Figure 4 :

 44 Figure 4 gives a set of simulated hydrostatic compression versus volumetric strain curves for different values of the parameter g of the fluid density law, see Equation (31). The other porosity parameters are fixed and set to, see Equation (8): Q = 1 M P a and m = 1.

Figure 5 :

 5 Figure 5: Eulerian porosity under undrained hydrostatic compression.

Figure 6

 6 Figure 6 compares the responses of the three laws computed with the following same set of material parameters that have been chosen to mark the difference between the three model responses: κ = 0.164 M P a, Q = 0.1 M P a, m = 1 and g = 0.4 with the initial values φ 0 = 0.6 and p 0 = 0.1 M P a. One can observe that the response with the simplified law of Equation (10) (the curve with dashed line in the figure) is close to that of the original model while the simple law of Equation (11) gives a similar but

Figure 6 :Figure 7 :

 67 Figure 6: Undrained hydrostatic compression. A comparison with the simplified versions of the porosity law.

Figure 8 :

 8 Figure 8: Hydrostatic compression for different loading velocities. Finite element mesh used for the computations.

Figure 8

 8 Figure 8 gives a set of simulated hydrostatic Cauchy stress versus volu-

  argument C e which are identical to those of the elastic left Cauchy-Green tensor denoted as usual by b e = F e F e T , i.e. χ sk ≡ χ sk (b e , ξ, p). In particular one has that J e = (detb e ) 1/2 for the elastic Jacobian. The total Kirchhoff stress tensor which is related to the second Piola-Kirchhoff stress tensor in the intermediate configuration by the relation τ = F e SF e T is then

Figure 10 :

 10 Figure 10: Undrained hydrostatic compression with the two porosity laws n e (J e , p) and n e (J, p). An illustration with different values of the fluid parameter g.

Remark 9 .Figure 11 :

 911 Figure 11: Corresponding Eulerian porosity evolutions obtained with n e (J e , p) and n e (J, p).

Figure 12 :

 12 Figure 12: Hydrostatic compression for different loading velocities. Finite element mesh used for the computations.

Figure 13 :

 13 Figure13: Typical deformed mesh and pore pressure field.

Figure 14 :

 14 Figure14: Hydrostatic compression for different loading velocities using the permeability law (59) with the parameter q = 2.

  k) n+1 > T OL or |G (k) n+1 | > T OL go to Step 3 ELSE set ϑ e n+1 = ϑ e (k) n+1 , ξ n+1 = ξ (k)n+1 and EXIT.3. Compute increments of consistency parameter and plastic variablesJ e (k) n+1 = exp ϑ e (k) Update: ∆γ (k+1) = ∆γ (k) + δ(∆γ (k) ) ϑ e (k+1) n+1 = ϑ e (k) n+1 + ∆ϑ e (k)

Table A .

 A 1: Undrained pore pressure update.

  1 + J e ∂χ por ∂J e (J e , p) (B.10) evaluated at {J n+1 , J e n+1 , p n+1 }.

  Table B.2: Return mapping algorithm for {ϑ e n+1 , ξ n+1 }.

comparison of the different responses in general, and between the two variants of the nonlinear porosity law introduced above in particular. The evolution equations of the form (50) 1-2 characterize the plastic behavior of the solid skeleton irrespective of the fluid pore pressure. Hence, the plastic material parameters can be identified solely from fully drained tests.

For the drained solid skeleton, the following partial free energy is used

where b e = J e 2/3 b e is the volume preserving left Cauchy-Green tensor, κ and µ are the bulk and shear moduli, respectively, and the function H(ξ)

characterizes isotropic plastic hardening. The effective stress yield criterion is chosen to be pressure-dependent as

where p ′ = 1 3 [τ ′ : 1] is the effective hydrostatic Kirchhoff stress, α is a material parameter, and the stress-like plastic hardening variable in the sense that q = -∂ ξ H is of the saturation type given by, see for example [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Simo | Numerical analysis and simulation of plasticity[END_REF]Nedjar, 2002b),

where σ ∞ y ≥ σ y > 0, H ≥ 0 and δ ≥ 0 are material constants; σ y being the initial flow stress. Furthermore, we consider an associated flow with

The fluid mass conservation equation ( ṁf = 0 in our undrained conditions) is solved for both the two alternative porosity laws n = n e (J e , p) and example [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Holzapfel | Nonlinear solid mechanics. A continuum appraoch for engineering[END_REF]. The effective second Piola-Kirchhoff stress tensor defined in (27) 2 is then written as

where Dev(.) = (.)-(1/3)[(.) : C]C -1 is the deviator operator in the material description. Equivalently, the effective Cauchy stress tensor such that Jσ ′ = F S ′ F T is given by

where b = J -2/3 b is the volume-preserving left Cauchy-Green tensor with b = F F T , and where dev(.) = ( .) -(1/3)[(.) : 1]1 is the deviatoric operator in the spatial description. Hence, from (17) 2 and (A.3), the total hydrostatic Cauchy stress, denoted here by ̟, is given by

For the fluid part, since no diffusion is allowed, the conservation of fluid mass (20) reduces to ṁf = 0 which, with definitions ( 19) and ( 12), expands

On the one hand, with the fluid constitutive law (31), one has

and, on the other hand, dividing Equation (A.5) by the product ρ f J, and using (A.6), leads to the following equivalent equation for the fluid mass

Even written in rate form, this latter is time-independent. Anyhow, let still consider a typical time interval [t n , t n+1 ], and assume that J n and p n are known initial data at time t n . An implicit backward-Euler scheme gives the following finite difference approximation

where n n+1 = n(J n+1 , p n+1 ), and where the time increment ∆t = t n+1 -t n has been eliminated. This nonlinear equation is to be solved for the updated pore pressure p n+1 when the actual volumetric strain J n+1 is prescribed.

This task can be accomplished by means of a Newton iterative procedure as summarized in Table A.1. The hydrostatic Cauchy stress ̟ n+1 at time t n+1 is computed by a simple evaluation of the relation (A.4) with the henceforth known ordered pair (J n+1 , p n+1 ).

This simple modeling example can be used as a benchmark test for more general numerical tools based on the finite element method.

Appendix B. Poroplastic undrained hydrostatic compression

As for the poroelastic problem, a compressible model of the neo-Hookean type is also chosen for the elastic behavior of the solid skeleton. Its partial free energy is this time given by