Giulio Sciarra

NICOLETTA IANIRO Francesco Dell'isola

Angela Madeo

A VARIATIONAL DEDUCTION OF SECOND GRADIENT POROELASTICITY PART I: GENERAL THEORY

Keywords: poromechanics, second gradient materials, lagrangian variational principle

Second gradient theories have to be used to capture how local micro heterogeneities macroscopically affect the behavior of a continuum. In this paper a configurational space for a solid matrix filled by an unknown amount of fluid is introduced. The Euler-Lagrange equations valid for second gradient poromechanics, generalizing those due to Biot, are deduced by means of a Lagrangian variational formulation. Starting from a generalized Clausius-Duhem inequality, valid in the framework of second gradient theories, the existence of a macroscopic solid skeleton Lagrangian deformation energy, depending on the solid strain and the Lagrangian fluid mass density as well as on their Lagrangian gradients, is proven.

Introduction

Poroelasticity stems from Biot's pioneering contributions on consolidating fluid saturated porous materials [START_REF] Biot | General theory of three-dimensional consolidation[END_REF]] and now spans a lot of different interrelated topics, from geo-to biomechanics, wave propagation, transport, unsaturated media, etc. Many of these topics are related to modeling coupled phenomena (for example, chemomechanical swelling of shales [START_REF] Dormieux | [END_REF]Coussy 2004], or biomechanical models of cartilaginous tissues), and nonstandard constitutive features (for instance, in freezing materials [Coussy 2005]). In all these cases, complexity generally remains in rendering how heterogeneities affect the macroscopic mechanical behavior of the overall material.

It is well known from the literature how microscopically heterogeneous materials can be described in the framework of statistically homogeneous media [Torquato 2002] considering suitable generalizations of the dilute approximation due to Eshelby [Nemat-Nasser and Hori 1993;Dormieux et al. 2006]; however, some lack in the general description of the homogenization procedure arises when dealing with heterogeneous materials, the characteristic length of which can be compared with the thickness of the region where high deformation gradients occur. This could be due, for example, to external periodic loading, the wavelength of which is comparable with the characteristic length of the material, or to phase transition, etc.

From the macroscopic point of view the quoted modeling difficulties, arising when high gradients occur, are discussed in the framework of so called high gradient theories [Germain 1973], where the assumption of locality in the characterization of the material response is relaxed. In these theories, the momentum balance equation reads in a more complex way than the classical one used for Cauchy continua. As a matter of fact, it is the divergence of the difference between the stress tensor and the divergence of so-called hyperstresses that balance the external bulk forces. Stress and hyperstress are introduced by a straightforward application of the principle of virtual power, as those quantities working on the gradient of velocity and the second gradient of velocity, respectively [Casal 1972;[START_REF] Casal | Equations du mouvement des fluides thermocapillaires[END_REF] 508 GIULIO SCIARRA, FRANCESCO DELL'ISOLA, NICOLETTA IANIRO AND ANGELA MADEO Gouin 1988]. Even the classical Cauchy theorem is, in this context, revised by introducing dependence of tractions not only on the outward normal unit vector but also on the local curvature of the boundary [dell' Isola and Seppecher 1997]; moreover symmetric and skew-symmetric couples (the actions called "double-forces" by Germain) must be prescribed on the boundary in terms of the hyperstress tensor together with contact edge forces along the lines where discontinuities of the normal vector occur.

Following the early papers on fluid capillarity [Casal 1972;Casal and Gouin 1988], the second gradient model can indeed be introduced by means of a variational formulation where the considered Helmholtz free energy depends both on the strain and the strain gradient tensors.

In the case of fluids, second gradient theories are typically applied for modeling phase transition phenomena [de Gennes 1985] or for modeling wetting phenomena [de Gennes 1985], when a characteristic length, say the thickness of a liquid film on a wall, becomes comparable with the thickness of the liquid/vapor interface [Seppecher 1993], annihilation (nucleation) of spherical droplets, when the radius of curvature is of the same order of the thickness of the interface [dell 'Isola et al. 1996], or topological transition [Lowengrub and Truskinovsky 1998].

In the case of solids, second gradient theories are applied, for instance, when modeling the failure process associated with strain localization [Elhers 1992;[START_REF] Vardoulakis | [END_REF]Chambon et al. 2004]. To the best of our knowledge, second gradient theories are very seldom applied in the mechanics of porous materials [dell 'Isola et al. 2003] and no second gradient poromechanical model, consistent with the classical Biot theory, is available except the one presented in [Sciarra et al. 2007]. As gradient fluid models, second gradient poromechanics will be capable of providing significant corrections to the classical Biot model when considering porous media with characteristic length comparable to the thickness of the region where high fluid density (deformation) gradients occur. We refer, for instance, to crack/pore opening phenomena triggered by strain gradients or fluid percolation, the characteristic length being in this case the average length of the space between grains (pores).

Several authors have focused their attention on the development of homogenization procedures capable of rendering the heterogeneous response of the material at the microlevel by means of a second gradient macroscopic constitutive relation [START_REF] Pideri | [END_REF]Camar-Eddine and Seppecher 2003]; however, very few contributions seem to address this problem in the framework of averaging techniques [START_REF] Drugan | [END_REF]Gologanu and Leblond 1997;Koutzetzova et al. 2002]. The present work does not investigate the microscopic interpretation of second gradient poromechanics, but directly discusses its macroscopic formulation. It is divided into two papers: in the first paper the basics of kinematics, Section 2; the physical principles, Section 3; the thermodynamical restrictions, Section 4; and in Section 5 the variational deduction of the governing equations for a second gradient fluid filled porous material are presented.

In particular, in Section 2 a purely macroscopic Lagrangian description of motion is addressed by introducing two placement maps in χ s and φ f (Equation ( 1)). We do not explicitly distinguish which part of the current configuration of the fluid filled porous material is occupied at any time t by the solid and fluid constituents, this information being partially included by the solid and fluid apparent density fields, which provide the density of solid/fluid mass with respect to the volume of the porous system (Equation ( 5)).

The deformation power, or stress working (Equation ( 12)), following Truesdell [1977] is deduced in Section 3 starting from the second gradient expression of power of external forces (Equation ( 9)) Cauchy theorem (Equation ( 10)) and balance of global momentum (see ( 11)).

In the spirit of Coussy et al. [1998] and Coussy [2004] thermodynamical restrictions on admissible constitutive relations are stated in Section 4, finding out a suitable overall potential, defined on the reference configuration of the solid skeleton. This last depends on the skeleton strain tensor and the fluid mass content, measured in the reference configuration of the solid, as well as on their Lagrangian gradients, in Equation ( 18).

Finally a deduction of the governing equations is presented in Section 5, based on the principle of virtual works, by requiring the variation of the internal energy to be equal to the virtual work of external and dissipative forces (see ( 19)). A second gradient extension of the two classical Biot equations of motion [Coussy 2004;Sciarra et al. 2007], endowed with the corresponding transversality conditions on the boundary, is therefore formulated (see Equations ( 30)-( 33)). Generalizing the treatment developed, for example, by [START_REF] Baek | [END_REF] for first gradient theories, one of the equations of motion found by means of a variational principle is interpreted as the balance law for total momentum, when suitable definitions of the global stress and hyperstress tensors are introduced (see ( 34)).

In a subsequent paper (Part II, to be published in a forthcoming issue of this journal), an application of the second gradient model to the classical consolidation problem will be discussed. Our aim is to show how the present model enriches the description of a well-known phenomenon, typical of geomechanics, curing some of the weaknesses of the classical Terzaghi equation [von Terzaghi 1943]. In particular we will figure out the behavior of the fluid pressure during the consolidation process when varying the initial pressures of the solid skeleton and/or the saturating fluid. From the mathematical point of view, the initial boundary value problem will be discussed according with the theory of linear pencils.

Kinematics of fluid filled porous media and mass balances

The behavior of a fluid filled porous material is described, in the framework of a macroscopic model, adopting a Lagrangian description of motion with respect to the reference configuration of the solid skeleton. At any current time t the configuration of the system is determined by the maps χ s and φ f , defined as

χ s : s × → , φ f : s × → f , (1)
where α (α = s, f ) is the reference configuration of the α-th constituent, while is the Euclidean place manifold, and indicates a time interval. The map χ s ( • , t) prescribes the current (time t) placement x of the skeleton material particle X s in s . The map φ f ( • , t), on the other hand, identifies the fluid material particle X f in f which, at time t, occupies the same current place x as the solid particle X s . Therefore the set of fluid material particles filling the solid skeleton is unknown, to be determined by means of evolution equations. Both these maps are assumed to be at least diffeomorphisms on . The current configuration t of the porous material is the image of s under χ s ( • , t). In accordance with the properties of χ s and φ f it is straightforward to introduce the fluid placement map as

χ f : f × → , such that χ f ( • , t) = χ s ( • , t) • φ f ( • , t) -1 ,  t x  s  s  f  t * X s x *   s   f  f X f X * f  f *  f -1  f *-1 Figure 1.
Lagrangian variations of the placement maps χ s , φ f , and χ f .

where χ f ( • , t) is still a diffeomorphism on . Figure 1 shows how the introduced maps operate on the skeleton particle X s ∈ s ; admissible variations of the two maps χ s ( • , t) and φ f ( • , t) are also depicted, in Section 5. In this way the space of configurations we will use has been introduced. Independently of t ∈ , the Lagrangian gradients of χ s and φ f are introduced as

F s ( • , t) : s → Lin (V ) , f ( • , t) : s → Lin (V ) , X s → ∇ s χ s (X s , t) , X s → ∇ s φ f (X s , t) , (2)
with V being the space of translations associated to the Euclidean place manifold. In Equation (2) ∇ s indicates the Lagrangian gradient in the reference configuration of the solid skeleton; analogously, the gradient of χ f is given by

F f X f , t = F s (X s , t) . f (X s , t) -1 , where X f = φ f (X s , t). 1
In the following the fluid Lagrangian gradient of χ f will be indicated both by F f or ∇ f χ f when confusion can arise. Moreover, the time derivatives of χ s and χ f , say the Lagrangian velocities of the solid skeleton and the fluid, can be introduced as

for all X α ∈ α , α (X α , • ) : → V , t → dχ α dt (X α ,t)
.

We also introduce the Eulerian velocities v α as the push-forward of α into the current domain

v α ( • , t) = α ( • , t) • χ α ( • , t) -1 .
In the following we do not explicitly distinguish the map χ s from its section χ s ( • , t) if no ambiguity can arise. Moreover we will distinguish between the Lagrangian gradient (∇ s ) in the reference configuration of the solid skeleton and the Eulerian gradient (∇) with respect to the current position x. Analogously, the solid Lagrangian and the Eulerian divergence operations will be noted by div s and div, respectively.

All the classical transport formulas can be derived both for the solid and the fluid quantities; in particular,

1 From now on we will indicate single, double and triple contraction between two tensors with ., :, and . . . respectively.

those ones for an image volume and oriented surface element turn to be

d t = J α d α , nd S t = J α F -T α .n α d S α ,
where d t and d S t represent the current elementary volume and elementary oriented surface corresponding to d α and d S α , respectively, where J α = det F α , and where n and n α are the outward unit normal vectors to d S t and d S α . As far as only the solid constituent is concerned, we can understand that deformation induces changes in both the lengths of the material vectors and the angles between them.

As it is well known, the Green-Lagrange strain tensor ε measures these changes, and is defined as

ε := 1 2 F T s .F s -I , (3)
where I clearly represents the second order identity tensor.

The balance of mass both for the solid and the fluid constituent are introduced as

α = t ρ α d t = const = α ρ 0 α d α , (α = s, f ) , ( 4
)
where α is the total mass of the α-th constituent, ρ α is the current apparent density of mass of the α-th constituent per unit volume of the porous material, while ρ 0 α is the corresponding density in the reference configuration of the α-th constituent. When localizing, Equation (4) reads

ρ α J α = ρ 0 α , (α = s, f ) ,
or, in differential form,

d α ρ α dt + ρ α div (v α ) = 0, (α = s, f ) , (5)
where d α ρ α /dt represents the material time derivative relative to the motion of the α-th constituent. In other words,

d α dt := d dt X α =const
.

The macroscopic conservation laws could also be deduced in the framework of micromechanics [Dormieux and Ulm 2005;Dormieux et al. 2006] starting from a refined model, where the solid and the fluid material particles occupy two disjoint subsets of the current configuration, and considering an average of the solid and fluid microscopic mass balances. The macroscopic laws do involve the so called apparent density of the constituents and suitable macroscopic velocity fields. For a detailed description of the procedure which leads to averaged conservation laws we refer to the literature [Coussy 2004].

2.1. Pull back of continuity equations. It is clear that Equation (5) consists of Eulerian equations, meaning that they are defined on the current configuration of the porous medium. Following Wilmanski [1996] and Coussy [2004] we want to write both these equations in the reference configuration of the solid skeleton. With this purpose in mind let us define the relative fluid mass flow w as w := ρ f v fv s . The use of this definition allows us to rearrange the fluid continuity (5) in the form

d s ρ f dt + ρ f div v s + div w = 0. ( 6
)
We want now to rewrite the continuity equation for the fluid constituent in the reference configuration of the solid skeleton. The Lagrangian approach to the fluid mass balance can be carried out by introducing the current Lagrangian fluid mass content m f , defined as

m f := J s ρ f • χ s . ( 7
)
Furthermore, let M be the Lagrangian vector referred to the reference configuration of the solid and related to the flow w through the relations

M := J s F -1 s . (w • χ s ) , J s (div w • χ s ) = (div s M) . ( 8
)
By using the definitions from Equations ( 7) and ( 8) in ( 6) the fluid Lagrangian mass balance takes the form dm f dt + div s M = 0.

Power of external forces

In this section, starting from the statement of the power of external forces for a second gradient solidfluid mixture, we deduce its corresponding reduced form, accounting for the extended Cauchy theorem valid for second gradient continua [Casal 1972;Germain 1973;[START_REF] Dell'isola | Edge contact forces and quasi balanced power[END_REF], and the balance of global momentum. The external power ext v s , v f for a second gradient porous medium can be defined as a continuous linear functional of the velocity fields v α ; in particular

ext v s , v f := t b s .v s + b f .v f d t + ∂ t t s .v s + t f .v f d t + ∂ t τ s . ∂v s ∂n + τ f . ∂v f ∂n d t + m k=1 k f k s .v s + f k f .v f dl, ( 9
)
where t is the current volume occupied by the porous medium, ∂ t its boundary, and m is the number of edges k (if any) of the boundary. In addition, b α , t α , τ α , and f k α represent the body force density, the generalized traction force (Cauchy stress vector), the double force vector, and the force per unit line acting on the k-th edge of the boundary, respectively.

The physical meaning of the double force τ α can be described in a way similar to that used in different contexts in [Germain 1973] and [dell' Isola and Seppecher 1997]. It can be regarded as the sum of two different contributions, the first of which works on the rate of dilatancy along the outward unit normal n (∇v α : (n ⊗ n)), and the second being a tangential couple working on the vorticity; this nomenclature is due to Germain [1973].

Let σ α and α be the apparent Cauchy stress and hyperstress tensors per unit volume of the porous material relative to the α-th constituent [Germain 1973;[START_REF] Dell'isola | Edge contact forces and quasi balanced power[END_REF]]. The Cauchy theorem can be extended for a second gradient continuum, and in particular for a second gradient porous continuum, in order to specify how the generalized external tractions appearing in (9) can be balanced by the internal forces when considering any subdomain of the current volume. In particular we have 2

t α = (σ α -div α ).n -div S ( α .n) , τ α = ( α .n).n, f k α = [] ( α .n).ν[] k , ( 10
)
where ν is the binormal unit vector which form a left-handed frame with the unit normal n and the unit vector being tangent to the k-th edge. We note that, since n is not continuous through the edge k, the vector ν is also discontinuous when passing from one side of the edge k to the other. It is for this reason that the edge force f k α is balanced by the jump of the internal force ( α .n).ν through the edge k. Extending classical poromechanics [START_REF] Biot | General theory of three-dimensional consolidation[END_REF]Coussy 2004] we now define the overall stress and hyperstress tensors as σ := σ f + σ s and := f + s , so that the momentum balance for the porous medium as a whole reads div (σdiv

) + b = 0, ( 11
)
where b = b s + b f is the overall body force. Bearing in mind both the extended Cauchy theorem, Equation (10), and the overall balance of momentum, ( 11), together with the principle of virtual powers ( ext = def. ), ( 9) leads to the expression for the deformation power

def. (v, ω) = t σ : ∇v + div σ T f .ω + . . . ∇∇v + f . . . ∇∇ω -div div f .ω d t . ( 12
)
Here and later on v := v s and ω := v fv s . Moreover, it must be remarked that ( 12) is obtained under the hypothesis of absence of volume forces (b α = 0) so that no inertia is taken into account in our model. We refer to [Coussy 2004] for the complete form of the deformation power in the case of first gradient porous continua. From now on, we also assume that the structure of the hyperstress tensors α (α = f, s) takes the particular form

α = I ⊗ c α , α = f, s, ( 13
)
where I is the second order identity tensor and c α is a kind of hyperstress vector related to the α-th constituent. The use of this assumption restricts second gradient external forces just to vector fields τ α , which only work on the stretching velocity of the α-th constituent; in other words, no contribution to the vorticity on the boundary of t comes from τ α . The aforementioned hypotheses indeed restrict the second gradient model; however, solid microdilatancies and capillarity effects can be still described by this second gradient model. According to (13) the external power due to second gradient effects, Equations ( 9) and ( 10), reduces to

τ α . ∂v α ∂n = {[(I ⊗ c α ) .n] .n} . ∂v α ∂n = (c α .n) n. ∂v α ∂n = (c α .n) [∇v α : (n ⊗ n)] .
3.1. Pull-back operations. Let us now consider the solid reference configuration pull-back of the deformation power; in order to do so, we will introduce the Piola-Kirchhoff like stress and hyperstress tensors for the overall body and for the fluid constituent. Thus, Piola-Kirchhoff stress (S) and hyperstress (γ ) are defined so that J s σ : ∇v =: S :

dε dt ⇒ S = J s F -1 s .σ .F -T s , (14)
2 Fixed a basis (e 1 , e 2 , e 3 ), where e 1 and e 2 span the plane tangent to the surface ∂ t at x, and the surface divergence of a second order tensor field A is defined as div

S A := 2 α=1 (∂ A/∂ x α ) e α . J s . . . ∇∇v =: γ . ∇ s dε dt T : C -1 -(∇ s C) T : C -1 . dε dt .C -1 ⇒ γ = J s F -1 s .c,
where C = F T s .F s is the Cauchy-Green strain tensor and c is the total hyperstress vector defined by c = c s + c f . Moreover, the fluid ones are S f =: J s F -1 s .σ f .F -T s , and γ f =: J s F -1 s .c f . The deformation power def. can be finally written in the Lagrangian form def. = s def. d s ,

where def. = S -C -1 . ((∇ s C).γ ) .C -1 : dε dt + C -1 ⊗ γ . . . ∇ s dε dt + div s 1 m f S T f .M + ∇ s div s M m f .γ f - M m f .∇ s div s γ f + div s J -1 s M m f .∇ s J s γ f .
4. Thermodynamics: deduction of a macroscopic second gradient strain energy potential

In this section, starting from the first and second principles of thermodynamics, we will prove that a suitable macroscopic strain potential can be identified depending both on the solid strain and on the fluid mass density as well as on their Lagrangian gradients. Let e α be the Eulerian density of internal energy relative to the α-th constituent, and the corresponding energy density of the porous medium is defined as e := ρ s e s + ρ f e f . The first principle of thermodynamics can be written as [Coussy 2004]

d s dt t ρ s e s d t + d f dt t ρ f e f d t = ext + Q,
where Q :=t q.nd t is the rate of heat externally supplied, and where q is the heat flow vector. In the Lagrangian form the first principle reads

d dt s E d s = def. - s div s e f M + Q d s . ( 15
)
where we recall that d/dt is the material time derivative associated with the motion of the solid, E := J s e represents the Lagrangian density of internal energy, and Q is the Lagrangian heat flux defined by Q := J s F -1 s .q. Starting from Equation ( 15), the local Lagrangian form of the first principle is naturally given by

d E dt = def.
div s e f M + Q .

Let us now consider the second principle of thermodynamics and introduce the overall Eulerian density of entropy s as s := ρ s s s + ρ f s f . The corresponding Lagrangian entropy is S := J s s, and the Lagrangian form of the second principle can be written as [Coussy 2004]

d dt s S d s ≥ - s div s s f M + Q T d s . ( 16
)
If we now introduce the Helmholtz free energy as := E -T S, Equation ( 16) can be rewritten in the local form as

d E dt -S dT dt - d dt ≥ -T div s s f M + Q T .
Merging the local form of the first and the second principles, the extended Clausius-Duhem inequality (dissipation inequality) can be deduced. In particular, following Sciarra et al. [2007], we will distinguish different contributions to the dissipation function due to the solid and fluid motion and to thermal effects, respectively ( s , f , and th ).

We now constitutively restrict [Coleman and Noll 1963] the admissible processes only to those ones which guarantee the dissipation inequality to be satisfied, because s , f and th are separately nonnegative. In particular, the solid dissipation s reads

s = S -C -1 . [(∇ s C).γ ] .C -1 - J s C -1 .∇ s J -1 s m f ⊗ γ f m f : dε dt + C -1 ⊗ γ . . . d dt (∇ s ε) + g f - 1 + 1 tr I γ f .∇ s 1 m f - J -1 s tr I γ f m f .∇ s J s . dm f dt - γ f m f . d dt ∇ s m f -S dT dt - d dt . ( 17
)
Assuming nondissipative processes occurring in the solid skeleton ( s = 0), Equation ( 17) allows for regarding the internal energy as a state function

= ε, m f , ∇ s ε, ∇ s m f , T . (18)
From now on we will treat an isothermal problem and therefore assume the energy does not depend on the temperature field T .

Variational deduction of second gradient poroelastic equations

5.1. Basic concepts and first variation of the internal energy. In this section we deduce the governing equations for a second gradient poroelastic continuum by means of a variational procedure. Variational approaches to first gradient mixture models are available in the literature [Bedford and Drumheller 1978;[START_REF] Gavrilyuk | [END_REF]Gouin and Ruggeri 2003].

In our case, we introduce the varied placement maps χ * s and φ * f for all X s ∈ s as

χ * s (X s , t) = χ s (X s , t) + δχ s (X s , t) , φ * f (X s , t) = φ f (X s , t) + δφ f (X s , t) ,
where δχ s and δφ f represent arbitrary variations of the functions χ s and φ f , respectively. The physical meaning of the variation δχ s is well known in continuum mechanics, and stands for the virtual displacement (deformation) of the solid skeleton. The variation δφ f , instead, accounts for the virtual relative displacement of a fluid material particle with respect to a solid one (see Figure 1). Since these variations keep fixed X s ∈ s we label them Lagrangian variations and we note that the symbol δ commutes with the integral over s and with the Lagrangian gradient operator ∇ s .

Following the statements of classical mechanics [Gantmacher 1970;Arnold 1989], the principle of virtual works reads

δ = δ ext + δ diss , ( 19
)
where δ represents the Lagrangian variation of the internal energy of the porous material, defined as

:= s d s ,
while δ ext and δ diss are the virtual works due to the external and dissipative forces acting on the porous system. Because of the aforementioned properties of the Lagrangian variations we can write

δ = δ s d s = s δ d s . ( 20
)
Recalling now Equation ( 18), ( 20) implies

δ = s ∂ ∂ε : δε + ∂ ∂m f δm f + ∂ ∂ (∇ s ε) . . . δ (∇ s ε) + ∂ ∂ ∇ s m f .δ ∇ s m f d s . ( 21
)
The variations δε, δm f , δ (∇ s ε), and δ ∇ s m f must now be rewritten in terms of the variations of the primitive kinematical fields χ s and φ f , bearing in mind that the Lagrangian variation commutes with the operator ∇ s . We show here directly the results obtained in Appendix A, to which we refer for detailed calculations,

δε = 1 2 (∇ s (δχ s )) T .F s + F T s .∇ s (δχ s ) , (22)
and

δm f = m f ∇ s φ f -T : ∇ s δφ f . ( 23
)
Substituting ( 22) and ( 23) into ( 21), integration by parts (see Appendix B for details), allows us to write the variation of the second gradient potential as

δ = s Ad s + ∂ s a d s + m k=1 E k α dl, ( 24
)
where m is the number of edges E k of the body in the reference configuration of the solid and

A := -div s F s . ∂ ∂ε -div s ∂ ∂(∇ s ε) .δχ s + ∇ s φ f -T . -m f ∇ s ∂ ∂m f -div s ∂ ∂ ∇ s m f .δφ f , a := F s . ∂ ∂ε -div s ∂ ∂(∇ s ε) .n s -div S s F s . ∂ ∂(∇ s ε) .n s .δχ s + F s . ∂ ∂(∇ s ε) .n s .n s . ∂(δχ s ) ∂n s + ∇ s φ f -T . m f ∂ ∂m f -div s ∂ ∂ ∇ s m f n s -m f ∇ S s ∂ ∂ ∇ s m f .n s .δφ f + ∇ s φ f -T . m f ∂ ∂ ∇ s m f .n s n s . ∂ δφ f ∂n s , α := F s . ∂ ∂ (∇ s ε) .n s .ν .δχ s , + ∇ s φ f -T . m f ∂ ∂ ∇ s m f .n s ν .δφ f .

Dissipative governing equations.

In order to obtain the equations of motion for a second gradient poroelastic continuum, the form of the external and dissipation virtual works, δ ext and δ diss , formally introduced in Equation ( 19), must be stated. The virtual dissipation δ diss will account for the classical Darcy effects and for the so called Brinkman-like contributions [Brinkman 1947]. We define the dissipation δ diss in the Eulerian configuration as

δ diss := - t v f -v s . δχ f • φ f -δχ s • χ -1 s d t - t .∇ v f -v s : ∇ δχ f • φ f -δχ s • χ -1 s d t , ( 25
)
where is the symmetric, definite positive Darcy tensor and is a suitably defined symmetric, definite positive second gradient Darcy-like tensor. Moreover, from now on, we assume the following Eulerian expression for the external work δ ext ,

δ ext := - ∂ t t . δχ s • χ -1 s + t f . δχ f • φ f -δχ s • χ -1 s d t . (26)
We restrict our attention to t and t f , defined as

t := p ext n, t f := ρ f µ ext n, (27)
where p ext is the overall external pressure applied on ∂ t , and µ ext is the chemical potential of the fluid outside the porous system. By comparison of Equation ( 26) with ( 9), we are assuming vanishing double forces and edge forces on the external boundary, as well as vanishing bulk actions. Equation ( 26), the expression for the external work, states that the external force t works only on the displacement of the solid skeleton (δχ s ), while µ ext works on the fluid mass virtual relative displacement ρ f δχ fδχ s . We note that if µ ext is spatially constant then

∂ t ρ f µ ext δχ f • φ f -δχ s • χ -1 s d t = ∂ s µ ext • χ s δm f d s ,
that is, µ ext works on the fluid mass which leaves (or enters) the solid skeleton (see Appendix C for details). Equation ( 26) for δ ext can be rewritten (see Appendix C) on the reference configuration of the solid as

δ ext = ∂ s - p ext J s F -T s . n s .δχ s + µ ext m f ∇ s φ f -T . n s .δφ f d s . (28)
Finally, (25) for δ diss (see Appendix C) assumes the Lagrangian form

δ diss = s ∇ s φ f -T . J s D F T s . f • φ f -s .δφ f d s - s ∇ s φ f -T . F T s . div s J s .∇ s f • φ f -s .(F T s .F s ) -1 .δφ f d s + ∂ s ∇ s φ f -T . J s F T s ..∇ s f • φ f -s .(F T s .F s ) -1 .n s .δφ f d s . (29)
Starting from the principle of virtual works, (19), and using ( 24), (29), and (28) for δ, δ diss , and δ ext , respectively, we can write the local equations of motion on s as

-div s F s . ∂ ∂ε -div s ∂ ∂ (∇ s ε) = 0, (30)
and

∇ s φ f -T . -m f ∇ s ∂ ∂m f -div s ∂ ∂ ∇ s m f -J s D F T s . f • φ f -s + ∇ s φ f -T . F T s . div s J s .∇ s f • φ f -s .(F T s .F s ) -1 = 0. ( 31
)
Analogously the boundary conditions on ∂ s read

F s . ∂ ∂ε -div s ∂ ∂ (∇ s ε) .n s -div S s F s . ∂ ∂ (∇ s ε) .n s = -J s p ext F -T s .n s ∇ s φ f -T . m f ∂ ∂m f -div s ∂ ∂ ∇ s m f n s -m f ∇ S s ∂ ∂ ∇ s m f .n s + - ∇ s φ f -T . J s F T s ..∇ s f • φ f -s .(F T s .F s ) -1 .n s = ∇ s φ f -T . m f µ ext n s , (32)
F s . ∂ ∂ (∇ s ε) .n s .n s = 0, ∇ s φ f -T . m f ∂ ∂ ∇ s m f .n s n s = 0.
Finally, on the edges E k of the boundary (if any) the following conditions hold true:

F s . ∂ ∂ (∇ s ε) . n s .ν = 0, ∇ s φ f -T . m f ∂ ∂ ∇ s m f . n s ν = 0. (33)
The Darcy and Brinkman dissipations appearing in Equations ( 31) and ( 32) can be rewritten in terms of the Lagrangian vector M, previously defined as

M = m f F -1 s . v f -v s .
In fact, after some straightforward calculations, it can be proven that

∇ v f -v s = 1 m f (∇ s F s ) T . M T + F s .∇ s M + F s . M ⊗ ∇ s 1 m f .
We now show that (30) is in agreement with the classical second gradient balance law for the total momentum [Germain 1973;[START_REF] Dell'isola | Edge contact forces and quasi balanced power[END_REF]. In order to do so, considering assumption (13), it can be proven that the constitutive relations for S and γ (see Equations ( 14))

∂ ∂ε = S -C -1 . ((∇ s C).γ ).C -1 , ∂ ∂ (∇ s ε) = C -1 ⊗ γ , (34)
imply that (30) can be regarded as the solid-Lagrangian pull-back of (11). In other words,

J s div (σ -div ) = div s F s . S -C -1 . ((∇ s C).γ ).C -1 -F s . div s C -1 ⊗ γ = 0.

Concluding remarks

In this paper a purely macroscopic second gradient theory of poromechanics is presented, extending classical Biot poromechanics [START_REF] Biot | General theory of three-dimensional consolidation[END_REF]Coussy 2004]. Following a standard procedure, sketched in [Coussy et al. 1998], we determine a suitable representation formula of the deformation power, ( 12), for a second gradient porous medium, assuming the forces acting on solid skeleton to be balanced (using the generalized second gradient balance of momentum in the current domain and the generalized second gradient Cauchy theorem on its boundary) and the power of external forces to be that of two superposed second gradient continua [Germain 1973]. The principles of thermodynamics, together with the aforementioned representation of the deformation power, allow for deducing the existence of a suitable overall strain energy potential depending on the solid strain tensor ε and the solid Lagrangian fluid mass density m f , as well as on their Lagrangian gradients.

The Euler-Lagrange equations associated with the energy density are the governing equations of the problem. In particular, Lagrangian variations of the placement maps χ s and φ f are considered. It is worth noting that the governing equations associated with the solid Lagrangian displacement δχ s (when δφ f = 0) represents the balance of total momentum and therefore allows for the constitutive characterization of the overall stress and hyperstress tensors. This is a characteristic feature of the classical Biot model [START_REF] Baek | [END_REF], which is completely recovered in this more general framework. On the other hand, the governing equation associated with the fluid placement map δφ f represents the balance of momentum relative to the pure fluid, which, in this case, is a generalization of the classical Darcy law.

In part II, an application to the classical consolidation problem will show how the present model improve the classical ones. It is well known that second gradient theories are capable to detect boundary layer effects in the vicinity of interfaces; this is indeed what we will observe in the case of consolidation. In particular, a kind of fluid mass density increment in the neighborhood of the impermeable wall will be observed for the first time in a one dimensional problem [Mandel 1953;Cryer 1963].

Appendix A: Basic variations

We show here how to derive the variations δε and δm f in terms of the kinematical variations δχ s and δφ f . Equation (3) for the Green-Lagrange strain tensor implies

δε = 1 2 δ s F T s .F s + F T s .δ s F s ,
where by definition F s := ∇ s χ s ; the expression ( 22) for δε is easily derived. As far as the variation δm f is concerned, recalling definition (7) for m f we can write

m f = J s J -1 f ρ 0 f , (A.1)
where ρ 0 f is the fluid density in the reference configuration of the fluid. Since by definition

J f := det(∇ f χ f ) and φ f := χ -1 f • χ s , we have J f = det ∇ f χ s • φ -1 f = det ∇ s χ s .∇ f φ -1 f = J s det ∇ s φ f -1 ,
where for the sake of simplicity we neglect the dependence of the considered fields on the reference places. Equation (A.1) thus reads

m f = ρ 0 f det ∇ s φ f . (A.2)
By derivation rule of the determinant and assuming ρ 0 f = constant, we get

δm f = ρ 0 f δ det ∇ s φ f = ρ 0 f det ∇ s φ f tr ∇ s φ f -1 .δ ∇ s φ f = m f ∇ s φ f -T : ∇ s δφ f .
Appendix B: Variation of the internal energy

The procedure to calculate the variation δ of the internal energy will be here shown in detail.

According to Equations ( 21)-( 23) and recalling that ∂/∂ε is a symmetric second order tensor, while ∂/∂(∇ s ε) is a third order tensor symmetric with respect to its first two indices, we can write

δ = s A 1 + A 2 s + A 2 f d s , (B.1)
where

A 1 := ∂ ∂ε : F T s .∇ s (δχ s ) + m f ∂ ∂m f ∇ s φ f -T : ∇ s δφ f , A 2 s := ∂ ∂ (∇ s ε) . . . ∇ s F T s .∇ s (δχ s ) , A 2 f := ∂ ∂ ∇ s m f .∇ s m f ∇ s φ f -T : ∇ s δφ f ,
account for the first gradient contribution to δ and for the solid and fluid second gradient contributions respectively.

The following identities are recalled in order to perform integrations by parts in (B.1); let λ, a, A, and be scalar, first, second, and third order tensor fields respectively. (Here ∇ S a indicates the surface gradient operator of a vector field a defined -analogously to div S s -as

∇ S a = ∂ a/∂ x α ⊗ e α , α = 1, 2
, where e α belong to the tangent plane.) Then, div

A T .a = A :∇ a + a. div A, div (λ A) = A.∇λ + λ div A, div (λa) = a.∇λ + λ div a, div T : A = A : div + (∇ A) . . . , ∇ a = ∇ S a + ∂ a ∂n ⊗ n,
where transposition for third order tensors is defined so as T := a i jk e k ⊗ e i ⊗ e j if = a i jk e i ⊗ e j ⊗ e k .

Moreover, given second order tensors A, B, C, third and first order tensors and a the following identities are satisfied:

A : (B.C) = B T . A : C = A.C T : B, T : A .a = A : (.a) .
Finally, if ϕ = χ s or ϕ = φ f , the identity holds true that

0 = div s det (∇ s ϕ) (∇ s ϕ) -T = det (∇ s ϕ) div s (∇ s ϕ) -T + (∇ s ϕ) -T .∇ s [det (∇ s ϕ)] . (B.2)
We underline that this equality holds unchanged for the surface divergence operator div S s . For the sake of simplicity, we will perform integration by parts for the first and second gradient terms appearing in (B.1) separately. Integration by parts of the first gradient term, recalling Equation (A.2) for m f and using Equation (B.2) for φ f , leads to

1 s A 1 d s = - s div s F s . ∂ ∂ε . δχ s d s + ∂ s F s . ∂ ∂ε . n s . δχ s d s - s ∇ s φ f -T . m f ∇ s ∂ ∂m f . δφ f d s + ∂ s ∇ s φ f -T . m f ∂ ∂m f n s . δφ f d s . (B.3)
Integrating by parts the solid second gradient term we get

2 s A 2 s d s = - s ∇ s (δχ s ) : F s . div s ∂ ∂ (∇ s ε) d s + ∂ s ∇ s (δχ s ) : F s . ∂ ∂ (∇ s ε) .n s d s = - ∂s F s . div s ∂ ∂ (∇ s ε) . n s . δχ s d s + s div s F s . div s ∂ ∂ (∇ s ε) . δχ s d s + ∂ s ∇ S s (δχ s ) + ∂ (δχ s ) ∂n s ⊗ n s : F s . ∂ ∂ (∇ s ε)
.n s d s .

Performing a further surface integration by parts we finally get

s A 2 s d s = - ∂s F s . div s ∂ ∂ (∇ s ε) .n s .δχ s d s + s div s F s . div s ∂ ∂ (∇ s ε) .δχ s d s - ∂ s div S s F s . ∂ ∂ (∇ s ε) .n s .δχ s d s + ∂ s F s . ∂ ∂ (∇ s ε) . n s . n s . ∂ (δχ s ) ∂n s d s + n k=1 E k F s . ∂ ∂ (∇ s ε) . n s .ν . δχ s dl. (B.4)
We finally rewrite the fluid second gradient term as

2 s A 2 f d s = s ∂ ∂ ∇ s m f .∇ s m f div s ∇ s φ f -1 .δφ f d s - s ∂ ∂ ∇ s m f .∇ s m f div s ∇ s φ f -T .δφ f d s ;
recalling Equation (A.2) for m f , using Equation (B.2) for φ f and rearranging, we have

s A 2 f d s = s ∂ ∂ ∇ s m f .∇ s m f div s ∇ s φ f -1 .δφ f d s + s ∂ ∂ ∇ s m f .∇ s ∇ s φ f -1 .δφ f .∇ s m f d s = s ∂ ∂ ∇ s m f .∇ s div s m f ∇ s φ f -1 .δφ f d s .
Integrating by parts we get

Appendix C: External and dissipation works

The dissipation and external works have been defined in ( 25) and ( 26) on the Eulerian configuration of the system in terms of δχ s and δχ f • φ f . These works must then be rewritten in terms of the independent variations δχ s and δφ f . In order to do so, the relationship between (δχ f • φ fδχ s ) and δφ f must be established. We know by definition that χ f • φ f = χ s , so that δ χ f • φ f = δχ s . Moreover, by differentiation rule for composite functions we have

δχ s = δ χ f • φ f = δχ f • φ f + ∇ f χ f • φ f .δφ f . But since χ f = χ s • φ -1 f , we get ∇ f χ f • φ f = ∇ s χ s . ∇ f φ -1 f • φ f = ∇ s χ s . ∇ s φ f -1 , so that δχ s = δχ f • φ f + ∇ s χ s . ∇ s φ f -1 . δφ f , or, δχ f • φ f -δχ s = -F s . ∇ s φ f -1 . δφ f . (C.1)
We now prove that the external work due to the force t f appearing in Equation ( 26) and prescribed by ( 27) represents the external work f ext done to change the fluid mass inside the porous system when the external chemical potential µ ext is assumed to be constant. We define this work as

f ext = s µ ext • χ s δm f d s ;
according to (23) and neglecting composition operations we can write

f ext = s µ ext m f ∇ s φ f -T : ∇ s δφ f d s ,

  . δφ f d s . (B.5) Substituting (B.3), (B.4), and (B.5) into (B.1), the variation of the internal energy given in (24) has been recovered.

  Recalling again Equation (A.2) for m f , using (B.2) for φ f and rearranging we get ⊗ n s , performing a last surface integration by parts and using (B.2) for the surface divergence operator we finally get

	Decomposing ∇ s	δφ f		as ∇ s	δφ f		= ∇ S s	δφ f	+	∂	δφ f	/∂n s
	s -∂s A 2 f d s = ∇ s φ f	∂s -T .	∇ s φ f m f ∇ S s ∂	-T . ∂ ∇ s m f	m f .n s	∂	∂ ∇ s m f .δφ f d s + . n s k=1 n s n	E k	.	∂ ∇ s φ f δφ f ∂n s	d s + -T . m f	∂	∂ ∇ s m f	.n s	ν	.δφ f dl
								-		∂s		∇ s φ f	-T .	m f div s	∂	∂ ∇ s m f	n s	. δφ f d s
	s A 2 f d s =	∂s div s		m f		∇ s φ f	-1 . δφ f +	∂ s	∂ ∇ s m f ∇ s φ f .n s d s --T . m f ∇ s s div s	m f	∇ s φ f	-1 . δφ f	div s	∂	∇ s m f ∂	d s ,
	and, integrating again,
	s	A 2 f d s =	∂s	m f div s		∇ s φ f		-1 .δφ f	∂	∂ ∇ s m f	. n s	d s
	+	∂s	∇ s m f .	∇ s φ f	-1 .δφ f		∂	∂ ∇ s m f	.n s	d s -	∂s	m f div s	∂	∂ ∇ s m f	∇ s φ f	-T .n s	.δφ f d s
																	+	s	m f	∇ s φ f	-T .∇ s	div s	∂	∂ ∇ s m f	. δφ f d s .
	s	A 2 f d s =	∂s	∇ s	δφ f	:	m f			∂	∂ ∇ s m f	. n s	∇ s φ f	-T	d s
									-	∂s	m f div s	∂	∂ ∇ s m f	∇ s φ f	-T . n s	.δφ f d s
																	+	s	m f	∇ s φ f	-T .∇ s	div s	∂	∂ ∇ s m f	.δφ f d s .


Angela Madeo

NICOLETTA IANIRO Francesco Dell'isola

Giulio Sciarra

A VARIATIONAL DEDUCTION OF SECOND GRADIENT POROELASTICITY II: AN APPLICATION TO THE CONSOLIDATION PROBLEM

Keywords: poromechanics, second gradient materials, consolidation. 607

The second gradient model of poromechanics, introduced in Part I, is here linearized in the neighborhood of a prestressed reference configuration to be applied to the one-dimensional consolidation problem originally considered by Terzaghi and Biot. Second gradient models allow for the description of boundary layer effects both in the vicinity of the external surface and the impermeable wall.

The formulated differential problem involves linear pencils of ordinary differential operators on a finite interval, with boundary conditions depending on the spectral parameter. Taking into account the dependence of the differential problem on initial stresses a linear stability analysis is carried out. Finally, numerical solutions are compared with the corresponding classical Terzaghi solutions.

]. The two models collapse one into the other when considering one-dimensional problems; however, this is not the case when modeling, for instance, the behavior of a saturated porous slab [Mandel 1953] or a saturated porous sphere [Cryer 1963]. In both circumstances Biot's three-dimensional model provides time increasing values of the water pressure (and fluid mass density) at the center of the slab (or sphere) if the Lamé constant µ of the skeleton is different from zero; on the other hand the solution for µ = 0 coincides with the one derived from Terzaghi's three dimensional model. This localized pore-fluid segregation is known in the literature as the Mandel-Cryer effect.

Introduction

This paper addresses a geotechnical application of the macroscopic second gradient poroelasticity theory presented in the first part; in particular we aim to treat the well known soil consolidation problem [Terzaghi 1943]. The consolidation of a soil layer of depth L can be schematically described as follows: when an external load p ext. is applied on the surface of the layer, the fluid starts moving from the layer towards the surface, and it finally leaves the system. While the fluid keeps flowing, the external load is gradually distributed to the solid skeleton, which starts to deform.

Different theories which model consolidation have been developed in the literature [START_REF] Biot | General theory of three-dimensional consolidation[END_REF]Terzaghi 1943;Heinrich and Desoyer 1961], however the two due to Terzaghi and Biot are surely the most widespread ones. Actually, as it was noted also by de Boer [1996], the derivation of the Terzaghi differential equation [Terzaghi 1923] is obscure and essentially driven by the comparison between the phenomenon of soil consolidation and that of heat propagation, rather than from the statement of suitable mechanical principles. On the other hand Biot's theory seems to be well grounded from the mechanical The occurrence of compaction localization phenomena has recently been discovered by Mollema and Antonellini [1996], who presented evidence of so called "compaction bands" in outcrops of the Jurassic Navajo sandstone in the Kaibab monocline, in Utah. These bands are characterized by volume loss due to microfracturing, but essentially no grain crushing or comminution. Later on, laboratory experiments have been developed [Olsson and Holcomb 2000;2003] using triaxial compression tests to reproduce the formation of these bands. These tests proved that increasing axial stress σ 11 initially determines only homogeneous axial strain 11 , however, when a suitable stress threshold has been overwhelmed, tabular zones associated with nonhomogeneous strain can be detected close to the axial borders of the specimen. Nonuniform compaction also affects fluid flow in the porous material, being detrimental if permeability of the compacted material is much reduced with respect to the uncompacted zone.

The second gradient poromechanical model presented in Part I [Sciarra et al. 2008] is capable of describing fluid mass density boundary segregation even in the one dimensional model. Both in the vicinity of the external consolidating surface and the impermeable wall, suitable boundary layer effects can be predicted by the second gradient model. Formation of segregation bands enhances high gradients of density of the fluid entrapped in the pores of the solid skeleton. This can be explained by means of nonvanishing hyperstresses at the boundary (see [Sciarra et al. 2008, Equations (32) 2 -(32) 3 ]); these last cause the pore pressure to differ from its reference initial value in a transient period, when dissipation does not yet dominate the evolution process.

Having in mind the classical Terzaghi's consolidation problem, whose space-time evolution is governed by the same equation as that of heat conduction, replacing temperature with pore-pressure, we claim that the second gradient model, presented in [Sciarra et al. 2008], is capable of regularizing the behavior of the Darcy flow inside the porous medium. Because of second gradient effects, the fluid mass density diffusion is smoothed. It is self evident that Terzaghi's theory does not model those phenomena occurring at the boundaries which oppose the fluid flow, for example, pore closure, solid-fluid capillarity, etc. The present model tries to macroscopically account for some of them and aims to establish the preliminary theoretical framework necessary for conceiving and designing any kind of experimental activity. In this paper it is shown that the overpressure occurring at the impermeable boundary actually depends on second gradient coefficients; therefore a more detailed analysis of these effects is recognized to be necessary.

Pore fluid segregation is probably the triggering effect for vertical drain or sand volcano formation, observed after liquefaction [Kolymbas 1998]; these last can indeed be interpreted as bifurcation modes of consolidation, which, in the case of the second gradient model should correspond to the boundary layer detected close to the impermeable wall becoming larger and larger. As it is not possible to identify this bifurcation mode in the case of linearized (small strain) theory, only the limit condition describing the stability/instability limit is detected. Further developments will be devoted to study of the onedimensional nonlinear problem.

From the mathematical point of view the present model, in which both second gradient conservative (relative to the behavior of the porous solid skeleton) and dissipative (relative to the flow of the saturating fluid) contributions are taken into account, implies that the newly formulated initial boundary value problem of consolidation fits in the framework of the theory of linear pencils of ordinary differential operators on a finite interval with boundary conditions depending on the spectral parameter. We refer to the general results presented, for example, in [Shkalikov 1986;Shkalikov and Tretter 1996;Marletta et al. 2003] for more details on this topic.

Linearization of the one-dimensional differential problem

We study here the aforementioned consolidation problem referring to the equations of motion for a second gradient porous medium as obtained in [Sciarra et al. 2008, Equations (30) and ( 31)], restricting

our attention to the one-dimensional case. This will allow us to compare our results with the classical ones due to Terzaghi [1923]. Of course, Equations ( 30)-( 31) can also be applied to treat three-dimensional problems, so extending classical Biot's equations.

Clearly, because of the one-dimensional hypothesis, all the gradient and divergence operations appearing in Equations ( 30) and ( 31) become simple derivatives with respect to the space variable x, and the deformation tensor ε simply reduces to its only nonzero component ε x x along the x axis. In the following we will indicate the component ε x x simply by ε.

From now on we assume the hypothesis of small deformations in the neighborhood of a suitable solid skeleton reference configuration. For the sake of simplicity we will therefore use the same notation as used in [Sciarra et al. 2008] for ε and m f to indicate the corresponding incremental quantities with respect to the considered small deformation parameter.

In accordance with the aforementioned assumptions and the introduced nomenclature, the quadratic expression for the Hemholtz free energy density , in terms of the state parameters (ε, m f , ε I , m I f ) is adopted as,

= -p ext. 0 ε + µ ext. 0 m f + 1 2 λ + 2µ + b 2 M ε 2 + 1 2 M m f m 0 f 2 -bMε m f m 0 f + 1 2 K ss + K 2 s f (ε I ) 2 + K s f m I f m 0 f ε I + 1 2 m I f m 0 f 2 , (1)
where ε I and m I f indicate the first spatial derivatives of ε and m f , respectively. The constant coefficients p ext. 0 , µ ext. 0 , and m 0 f account for the state of stress of the solid skeleton, the chemical potential of the fluid, and the initial apparent density of the fluid before any external perturbation is applied to the porous system. Moreover, λ and µ are the classical Lamé coefficients, b and M the Biot parameters, and , K ss , and K s f are the second gradient constitutive parameters.

The nonstandard energetic contributions associated with (ε I ) 2 , (m I f ) 2 , and ε I m I f are those responsible for the presence of hyperstresses in the balance equations of the overall material and the pure fluid. They allow for describing the compaction/dilatancy localization effects arising in the fluid-filled porous material when the fluid remains entrapped in the solid skeleton. In particular K ss and provide nonvanishing hyperstress on the solid skeleton and the pure fluid if second gradient coupling is negligible. Following the interpretation of double forces given in [Sciarra et al. 2008, Section 2], these two constitutive parameters allow for describing internal actions working on the rate of dilatancy along the outward unit normal. The coupling coefficient K s f is labeled as the cocapillarity coefficient in analogy to standard second gradient theories for capillarity models [Seppecher 1987]. It describes second gradient solid-fluid interactions and can be assumed as vanishing in contrast with K ss and , which will be proved to be positive when positiveness of strain energy density is required (see Equation ( 2)).

Values of the Lamé and Biot moduli can easily be recovered from the literature [Coussy 2004]; on the other hand no identification for the second gradient moduli is available up to now. It is not the purpose of this paper to set up a constitutive identification based on experiments or mathematical homogenization; conversely, our aim is that of exhibiting the capability of the model presented in [Sciarra et al. 2008] to catch compaction/dilatancy effects. Second gradient parameters therefore will be tuned so as to permit the one-dimensional model to show boundary layer effects for the solid strain and the fluid mass density in the vicinity of the external surfaces.

Requiring definite positiveness of the energy density defined in Equation ( 1), the following conditions on the parameters must hold:

λ + 2µ > 0, M > 0,
and

K ss > 0, > 0. ( 2
)
The first two conditions are well known in the framework of the classical Biot poromechanics; the second ones restrict the constitutive assumptions on the second gradient parameters.

2.1. Equations of motion. We will now deduce the linearized form of the equations of motion for the second gradient consolidation problem. In order to do so, it is worthwhile to recall that in the onedimensional, linearized problem the following chains of equalities hold:

F s I + ∇ s u = (1 + ε) I, (3)
u being the infinitesimal solid displacement field,

F -1 s (1 -ε) I, (4)
and

J s det F s = 1 + tr(∇ 0 u) = 1 + ε, (5)
where we recall that all the considered fields in the right hand side of (3)-( 5) have to be regarded as incremental quantities with respect to the small deformation parameter. Taking into account (1) for the strain energy and the one-dimensional form of (3) and ( 5) , the linearized governing equations given in [Sciarra et al. 2008, Equations (30) and ( 31)], reduce to

-p ext. 0 + λ + 2µ + b 2 M ε -bM m f m 0 f - K ss + K 2 s f ε I I -K s f m I I f m 0 f I = 0, (6)
and

-m 0 f M m I f m 0 f 2 - bM m 0 f ε I - K s f m 0 f ε I I I - m 0 f 2 m I I I f -D v f -v s + α v f -v s I I = 0, (7)
for the solid skeleton and the pure fluid, respectively. We have denoted by D and α the only nonzero component of the Darcy and Darcy-like tensors and , respectively (see [Sciarra et al. 2008, Equation (25)]), and by v fv s the vertical component of the relative velocity. In the absence of inertia forces the solid momentum conservation law, Equation ( 6), can be integrated in the form

λ + 2µ + b 2 M -p ext. 0 ε -bM m f m 0 f - K ss + K 2 s f ε I I -K s f m I I f m 0 f = const. := c 0 . (8)
Moreover, considering that v fv s is related to the apparent fluid density m f by means of the linearized continuity equation ṁf + m 0 f (v fv s ) I = 0, (7) can be rewritten, performing a derivative with respect to the space variable x, as

-m 0 f    M m I I f m 0 f 2 - bM m 0 f ε I I - K s f m 0 f ε I V - m 0 f 2 m I V f    + D m 0 f ṁf - α m 0 f ṁ I I f = 0, (9)
where we have indicated by ṁf the time derivative of m f . In order to write the linearized governing equations in a dimensionless form, the following quantities are introduced:

ξ = x L , mf = m f m 0 f , t = t τ , with τ = DL 2 M ,
where L is the depth of the solid layer.

According to these definitions, ( 8) and ( 9) can be rewritten in their dimensionless form as

(λ + 2µ + b 2 M -p ext. 0 ) λ + 2µ ε - bM λ + 2µ mf - (K ss + K 2 s f ) (λ + 2µ)L 2 ε I I - K s f (λ + 2µ) L 2 m I I f = c 0 λ + 2µ , (10)
and

M L 2 m I V f + M L 2 K s f ε I V -m I I f + b ε I I - α DL 2 ṁ I I f + ṁf = 0, (11)
which represent the linearized equations of motion for the consolidation problem. For the sake of simplicity, we will no longer distinguish between m f and mf , and, if not specified, m f will indicate the dimensionless quantity. Moreover, the dimensionless variables ξ and t will be also indicated by x and t if no confusion can arise.

Boundary conditions.

The constant c 0 is deduced from the boundary condition (BC) in x = 0 given in [Sciarra et al. 2008, Equation (32) 1 ] which, in the linearized form, reads

λ + 2µ + b 2 M -p ext. 0 λ + 2µ ε - bM λ + 2µ m f - K ss + K 2 s f (λ + 2µ)L 2 ε I I - K s f (λ + 2µ)L 2 m I I f = - p ext. (λ + 2µ) ; here p ext.
represents the incremental external pressure acting on the system, deriving from the linearization process ( p ext. = p ext. 0 + p ext.

). In other words, p ext. is the perturbation in the external load applied on the surface of the soil layer. Comparing this BC with (10) it is easy to recognize that c 0 =p ext. . Equations ( 10) and ( 11) represent a differential system of the sixth order in the space variable x and of the first order in time, the integration of which requires therefore six boundary conditions and one initial condition. In classical poromechanics the Terzaghi consolidation problem does not take into account second gradient effects, and indeed it can be obtained from ( 10) and ( 11) when the second gradient parameters , K ss , K s f , and α are vanishing. Clearly, the problem reduces in this case to a second order system with respect to the space variable x.

The boundary conditions for the consolidation problem are derived from the general ones deduced in [Sciarra et al. 2008, Equation (32)]. In particular, since the given problem is one-dimensional, surface divergence and surface gradient operations do not contribute to the BCs; moreover no edge E k of the boundary exists. Extending the classical BCs stated in the Terzaghi consolidation problem we assume:

• Zero fluid traction in x = 0. This BC states that the surface of the solid layer is kept drained, meaning that the fluid reaching the surface is continuously removed from the surface itself. This BC corresponds to the one given in [Sciarra et al. 2008, Equation (32) 1 ] which, in its linearized, dimensionless form, reduces to

m f -bε - K s f M L 2 ε I I - M L 2 m I I f + α DL 2 ṁf = m 0 f µ ext. M = 0, ( 12
)
where µ ext. represents the incremental chemical potential. In other words we have assumed the linearization µ ext. = µ ext. 0 + µ ext. . Assuming that the fluid is continuously removed from the surface of the layer, this implies a restriction to the case µ ext. = 0.

• Impermeable soil in x = L. With this BC we assume that the relative velocity is equal to zero in 7), which holds everywhere in the interval [0, L], the impermeability of the layer x = L can be rewritten in its dimensionless form as

x = L, implying v f -v s = 0. Using Equation (
-m I f + b ε I + K s f M L 2 ε I I I + M L 2 m I I I f - α DL 2 ṁ I f = 0. ( 13
)
• Zero double force for the overall system in x = 0 and x = L. These BCs are those ones given in [Sciarra et al. 2008, Equation (32) 2 ], and can be rewritten in their linearized dimensionless form as

m I f + (K ss + K 2 s f ) K s f ε I = 0.
We remind that the overall double forces are the contact forces introduced in the second gradient model, which work on the rate of pore opening/pore shrinkage. With the assumption of vanishing double forces on the boundary of the porous material, we claim that no external source of double force exists; internal double forces, on the contrary, allow for capturing the effects of pressure gradient concentration in the neighborhood of the external and impermeable surfaces [Holcomb and Olsson 2003].

• Zero fluid double force in x = 0 and in x = L. These BCs are those ones given in [Sciarra et al. 2008, Equation (32) 3 ]. They can be rewritten in their linearized dimensionless form as

m I f + K s f ε I = 0. ( 14
)
The assumption on fluid double forces can be interpreted similarly to that considered for the overall double forces. In this case no external double forces are exerted on the fluid boundary, but internal double forces, associated with pressure gradient concentration, account for internal capillarity and wetting phenomena.

Initial boundary value problem

The differential Equations ( 10) and ( 11) can be reduced to a unique differential equation introducing an auxiliary function V (x, t) which satisfies the relationships

ε = K s f (λ + 2µ) L 2 V I I + bM (λ + 2µ) V, (15)
and

m f = λ + 2µ + b 2 M -p ext. 0 (λ + 2µ) V - K ss + K 2 s f (λ + 2µ) L 2 V I I + p ext. bM . ( 16
)
In such a way, Equation ( 10) is identically satisfied, while (11) can be rewritten, after some straightforward calculations, as

C 1 V V I -C 2 p ext. 0 V I V -C 3 V I V + C 4 p ext. 0 V I I + C 5 p ext. 0 V I I -C 6 p ext. 0 V = 0; ( 17
)
on the other hand, the boundary conditions given in ( 12)-( 14) read

1 V I V -2 p ext. 0 V I I -3 V I I + 4 p ext. 0 V + 5 p ext. 0 V + 6 = 0 at x = 0, ( 18
) V I = 0 at x = 0, L , V I I I = 0 at x = 0, L , V V = 0 at x = L . (19)
Finally, the initial condition (corresponding to the instant in which the external load is applied) is deduced assuming that the apparent Lagrangian fluid density is vanishing for t = 0 + . For instance, m f x, t = 0 + = 0; this initial condition states, similarly to in classical consolidation, that, at the instant in which the external load is applied there is no instantaneous variation of the fluid density m f inside the soil layer. In terms of the auxiliary function V the initial datum reads as

V (x, 0 + ) := V in = - p ext. bM 1 C 4 (π 0 ) + k 6 , (20)
where ( 16) with m f = 0 has been solved using BCs given in (19).

All the coefficients appearing in the governing equation, ( 17), as well as in the initial and boundary conditions, ( 20) and ( 19), depend on the constitutive parameters, the solid initial stress p ext. 0 , and the increment of the external force p ext. ; their expressions are listed in Appendix A.

It must be remarked that the expression for the energy density assumed in (1) would not allow the linearized differential problem to explicitly depend on the initial solid stress p ext. 0 and on the initial chemical potential µ ext. 0 . In fact, considering (1) we can write, in dimensionless form,

∂ ∂ε = 1 λ + 2µ p ext. 0 + (1 + k 6 )ε -bk 5 m f , ∂ ∂m f = m 0 f λ + 2µ µ ext. 0 + k 5 m f -bk 5 ε, (21)
so that there are no linear terms (in ε and m f ) involving p ext. 0 and µ ext. 0 coming from ∂/∂ε and ∂/∂m f . On the other hand the dependence of the differential system ( 17)-( 19) on p ext. 0 is due to so called geometrical nonlinearities; as matter of fact it is the presence of F s in the balance of the total momentum (see [Sciarra et al. 2008]) which even in linearized problems implies a nontrivial dependence of the governing equations on p ext. 0 (see the term εp ext. 0 in (10)).

Considering the linearity of the differential problem and the nonhomogeneity appearing in the BC, given in Equation ( 18), we will look for a solution V (x, t) in the form

V (x, t) = V (x) + W (x, t), (22)
where V (x) is the solution of the stationary problem

C 1 V V I -C 2 V I V + C 4 V I I = 0, (23)
with nonhomogeneous BCs

1 V I V -2 V I I + 4 V = -6 at x = 0, ( 24
) V I = 0 at x = 0, L , V I I I = 0 at x = 0, L , V V = 0 at x = L , (25)
while the deviation W (x, t) is the solution of the initial boundary value problem (IBVP)

C 1 W V I -C 2 W I V -C 3 Ẇ I V + C 4 W I I + C 5 Ẇ I I -C 6 Ẇ = 0, (26)
with homogeneous BCs

1 W I V -2 W I I -3 Ẇ I I + 4 W + 5 Ẇ = 0 at x = 0, ( 27
)
W I = 0 at x = 0, L , W I I I = 0 at x = 0, L , W V = 0 at x = L , (28)
and nontrivial initial condition

W (x, 0 + ) := W in = V in -V (x) . (29)
For the sake of simplicity we will no longer specify the dependence of the coefficients C i and i on the prestress parameter p ext. 0 . It is easy to prove that the solution of the stationary problem given by ( 23)-( 25) when 4 = 0 is given by V

(x) = - 6 4 , 4 = 0, (30)
while if 4 = 0 the stationary solution V (x) exists if and only if 6 = 0; in this case a family of constant solutions of the stationary problem arises, so that we can write V (x) = K , and 4 = 6 = 0, where K is an undetermined constant. On the basis of the preliminary study of the stationary solution V (x) we can state that, according to the assumption ( 22), a solution of the given problem for the variable V exists if and only if 4 = 0 or 4 = 6 = 0. We will now restrict our attention to the case 4 = 0 and will analyze the case 4 = 0 later.

Fourier series solution

The initial boundary value problem given by ( 26)-( 29) is solved using the method of separation of variables; in other words we assume that W (x, t) = X (x)T (t). A straightforward calculation yields to the definition of the real eigenparameter λ as λ = Ṫ /T, which leads to T (t) = T 0 e λt . Consequently, the function X (x) must satisfy the eigenvalue problem

C 1 X V I -C 2 X I V + C 4 X I I = λ C 3 X I V -C 5 X I I + C 6 X , (31)
endowed with the BCs

1 X I V -2 X I I + 4 X = λ 3 X I I -5 X , at x = 0, (32)
X I = 0 at x = 0, L ,

X I I I = 0 at x = 0, L , X V = 0 at x = L . ( 33
)
This is a nonclassical spectral problem since the BCs also depend on the spectral parameter λ; in the literature this kind of spectral problem is referred to as a linear pencil (X ) = λ(X ). Many authors investigate the spectral properties of the differential operators and in suitable function spaces in order to guarantee completeness and orthonormality for the eigenfunction system and discreteness of the spectrum [Shkalikov 1986;Shkalikov and Tretter 1996;Marletta et al. 2003]. Here we rely on these general results and numerically determine a subset of the eigenfunction space so as to approach the requirements of the Parseval equality [Kolmogorov and Fomin 1975].

According to the aforementioned properties of the eigensystem, the solution of the considered IBVP can be given in Fourier series form as

W (x, t) = +∞ k=0 p k X k (x)e λ k t , ( 34
)
where p k denotes the k-th Fourier coefficient, and, in particular, p 0 the Fourier coefficient relative to the null eigenvalue λ = 0 (if any). It is easy to prove that if 4 = 0 the eigenfunction X 0 relative to the null eigenvalue is the trivial one X 0 = 0, so that in Equation ( 34) k runs now from one to infinity.

The eigenfunctions (X k ) k∈ are orthogonal with respect to the following bilinear form defined, in the Hilbert space

H 3 ([0, L]) × H 3 ([0, L]), as X k , X h := α 0 L 0 X k X h dx + α 1 L 0 X I k X I h dx + α 2 L 0 X I I k X I I h dx + α 3 L 0 X I I I k X I I I h dx, (35)
where the coefficients α i are defined as

α 0 = 4 C 6 , α 1 = 4 C 5 -C 4 5 + 2 C 6 , α 2 = 4 C 3 -3 C 4 + 2 C 5 -C 2 5 + 1 C 6 , α 3 = 2 C 3 -C 2 3 + 1 C 5 -C 1 5 .
It must be noted that expression (35) is indeed an inner product over the aforementioned function space if and only if all the coefficients α i are positive definite.

We notice that when the initial stress p ext. 0 is vanishing the α i coefficients are all positive (assuming positiveness of the energy , see (2)), so that (35) always represents an inner product. On the other hand, it is easy to verify that in the presence of prestress the positiveness of the aforementioned coefficients α i is guaranteed if and only if p ext. 0 < λ + 2µ (⇔ 4 > 0). The explicit form of the Fourier coefficients p k is determined by projecting the initial datum (29) on the k-th element of the Fourier series according to the inner product (35); in particular we can find

W in , X k = α 0 W in L 0 X k dx, (36) W in , X k = +∞ h=1 p h X h , X k = +∞ h=1 p h X h , X k = p k X k 2 , ( 37
)
where we have noted by • = • , • 1/2 the norm induced by the inner product • , • . Comparing Equation ( 36) with (37) it is easy to recognize that

p k = α 0 W in L 0 X k dx X k 2
, so that, recalling (34), the final form of the solution is

W (x, t) = α 0 W in +∞ k=1 1 X k 2 L 0 X k (ξ )dξ X k (x)e λ k t ,
and, according to ( 22) and ( 30), the solution for the variable V (x, t) is finally given by

V (x, t) = - 6 4 + α 0 W in +∞ k=1 1 X k 2 L 0 X k (ξ )dξ X k (x)e λ k t . ( 38
)
Finally, the fields ε and m f can be evaluated using ( 38) with ( 15) and ( 16) respectively.

4.1. The limit case 4 = 0. We have already mentioned that when 4 = 0 (p ext. 0 = λ + 2µ) the stationary solution V (x) exists if and only if 6 = 0 (⇔ p ext. = 0), and it is an undetermined constant K . This means that, corresponding to a critical value of the prestress p ext. 0 , no solution can be found if perturbing the porous system with an external load p ext. . The only possible solution is relative to the unloaded configuration of the porous system (p ext. = 0). In this case the solution for V (x, t) is found by solving the differential problem given by ( 17), ( 19), and ( 20) when 4 = 6 = 0. Separating the variables, V (x, t) = X (x)T (t), the solution can be found in Fourier series form as

V (x, t) = +∞ k=0 p k X k (x)e λ k t . ( 39
)
It must be noticed that when 4 = 0 the inner product (35) reduces to

X k , X h 4 =0 := α 1 L 0 X I k X I h dx + α 2 L 0 X I I k X I I h dx + α 3 L 0 X I I I k X I I I h dx,
and it is still well defined over the quotient space of the H 3 ([0, L]) functions, differing at most by a constant. It follows that the Fourier coefficients p k are now determined on the basis of the reduced form • , • 4 =0 of the inner product, according to the identities involving the initial condition

V in = constant, 0 = V in , X 0 4 =0 = p 0 X 0 + +∞ k=1 p k X k (x), X 0 4 =0 = p 0 X 0 , X 0 4 =0 , ( 40
) 0 = V in , X k 4 =0 = p 0 X 0 + +∞ k=1 p h X h (x), X k 4 =0 = p k X k 2 4 =0 , for all k ∈ , ( 41
)
where we have noted by

• 4 =0 = • , • 1/2
null eigenvalue λ 0 = 0, respectively, while (X k ) k∈ are the remaining eigenfunctions. Since X 0 = constant, Equation ( 40) reads p 0 0 = 0 ⇒ p 0 undetermined, moreover, (41) gives p k X k 2 4 =0 = 0 ⇒ p k = 0, ∀k ∈ . According to (39), the solution for V (x, t) is an undetermined constant, so

V (x, t) = p 0 X 0 := p 0 K = constant.
We want to remark that all the Fourier coefficients p k corresponding to nonvanishing eigenvalues turn to be zero only because the initial condition V in has been assumed to be constant; if it was not the case, (41) would have stated the expression for the coefficients p k , and the solution for V (x, t) would have been known except for a constant K . The fact that a family of constant solutions for V (x, t) arises can be seen as a sort of bifurcation phenomenon which is triggered when p 0 reaches the critical value λ + 2µ.

Finally, we underline that the null eigenvalue λ 0 = 0 belongs to the spectrum of the differential problem only when 4 = 0; in the following section we will show that when 4 > 0 only negative eigenvalues exist, while if 4 < 0 some positive eigenvalues appear.

Numerical results

In this section we will show the numerical solution of the differential problem, given by ( 17)-( 20), for a particular set of values of the constitutive parameters, which are listed in Table 1. The first gradient parameters are those relative to a water saturated clay, while the values of the second gradient dimensionless numbers are chosen in order to let boundary layer effects arise.

Fixing suitable values for the initial external pressure ( p ext. 0 = 4.9 GPa) and for the increment of this latter (p ext. = 1 MPa), so as to guarantee 4 > 0 and 6 = 0, we look for a numerical solution V (x, t) given by (38). In particular, we look for a numerical solution X (x) of the differential problem, given by (31)-(33), in the form

X (x) = 6 i=1 K i e β i (λ) x , (42)
where K i are the integration constants and β i (λ) are the solutions of the characteristic polynomial associated with the differential equation, (31). Consequently, BCs given by ( 32)-(33) yield

A(λ)v = 0, (43)
where A(λ) is a suitably defined 6 × 6 matrix and v := (K 1 , ..., K 6 ). We notice that the matrix A(λ) depends on the eigenparameter λ both because it appears in the differential equation, (31), and in the boundary condition, (32). It follows that the resulting eigenvalue problem cannot be classified as a standard eigenvalue problem. The system of algebraic equations, (43), has a nontrivial solution if and only if det [A(λ)] = 0, which leads to the calculation of the eigenvalues λ k (discrete spectrum). For each eigenvalue λ k an eigenfunction X k (x) is deduced in the form of Equation (42) so that a numerical solution for V (x, t) can be found according to (38).

M(G Pa) λ(G Pa) µ(G Pa) k 1 k 2 k 3 k 4 5 2 .3 1 .5 10 -2 10 -2 10 -2 10 -2
The numerical solution for V (x, t) involves a finite number of eigenfunctions N , where the choice of N is made so as to approach the condition stated by the Parseval equality [Kolmogorov and Fomin 1975].

Once the numerical solution for V (x, t) has been found, we can deduce the corresponding solutions for the fields ε and m f simply be referring to ( 15) and ( 16).

We now show the behavior of the fields ε and m f corresponding to the aforementioned values of the constitutive parameters, and initial and incremental pressures. In Figure 1 the fluid apparent density m f versus x is depicted for times in the very close neighborhood of t = 0 + . It can be noticed that a critical depth x cr 0.8 exists such that the density m f decreases for 0 ≤ x < x cr , while it increases for x cr < x ≤ 1. This means that the fluid of the upper regions actually leaves the layer, while the fluid contained in the deeper regions remains entrapped in the pores whose deformation consequently increases the apparent density m f . When increasing time (see Figure 2), the apparent density m f decreases along the whole depth of the layer, and finally approaches a constant value. This means that the fluid starts flowing out also from the deeper regions until the system reaches a new equilibrium and no fluid leaves the layer anymore. This effect is evidently related to viscosity, which dominates the evolution of the fluid density as time becomes larger and larger.

As far as the vertical deformation ε is concerned, the same qualitative behavior as that of m f is detected (see Figure 3 and Figure 4). For times close to t = 0 + (see Figure 3) the upper regions of the layer undergo to a vertical compression, which is connected to the fact that less fluid is present in the pores, while the deeper regions experience a sort of dilatancy which is connected to an over pressurization of the saturating fluid.

For increasing times (see Figure 4) a general further compression is detected along the whole depth of the layer (this is due to the fact that the fluid is uniformly flowing along the layer) until the layer does not deform anymore (equilibrium).

We remark that the chosen values of the prestress p ext. 0 are such that 4 > 0 so that the inner product Equation ( 35) is well defined. Consequently the solution for V (x, t) (and thus for ε and m f ) can be numerically evaluated.

It is interesting to notice that when the initial stress p ext. 0 is such that p ext. 0 < λ + 2µ ( 4 > 0) only negative eigenvalues λ k < 0 have been found, while in the region where 4 < 0 a finite number of positive eigenvalues arise. In Figure 5 the behavior of the first eigenvalue λ 1 is shown when varying p ext. 0 through the threshold p ext. 0 = λ + 2µ . It is worth noticing that when λ 1 passes from negative to positive values, the solution V (x, t) given in the form of Equation ( 38) blows up due to the presence of positive time exponentials; the solution thus experiences an unstable behavior related to the fact that p ext. 0 reaches a critical value. This kind of instability is known as geometrical instability since the presence of p ext. 0 in the differential problem is due to the geometry of the problem (see Equation ( 21)). The linearity of the present model does not allow us to capture solutions associated with unstable conditions; this is evident when considering that the bilinear form given in ( 35) is no longer a well defined inner product.

In order to show the influence of the solid prestress on the behavior of ε and m f we have found solutions for different values of p ext. 0 and noticed changes in the solution when approaching the threshold 4 = 0. Figure 6 shows the behavior of m f when p ext. 0 progressively approaches the critical value p ext. 0 = λ + 2µ. When increasing the value of p ext. 0 the fluid density decreases in the superficial regions of the layer, while increasing in the deeper ones. This means that the initial stress increases the capability of the fluid to flow out from the skeleton matrix close to the external surface, while pumping it in the deeper layers.

Let us now consider the second gradient constitutive parameters and the initial stresses to be vanishing. The resulting differential problem reduces to the classical Terzaghi consolidation problem. More particularly, Equation (10) reduces to

ε = bM λ + 2µ + b 2 M m f - p ext. λ + 2µ + b 2 M , which, substituted in (11), gives ṁf = am I I f , a = (λ + 2µ) λ + 2µ + b 2 M . ( 44
)
The Terzaghi consolidation problem thus reduces to the differential equation, Equation ( 44), together with the initial datum m f x, 0 + = 0 and the BCs, ( 12) and ( 13), which simplify into

m f = λ + 2µ + b 2 M λ + 2µ b p ext. λ + 2µ := c at x = 0. ( 45
)
and m I f = 0 at x = L, respectively. It is easy to notice that the BC, (45), and the initial datum, m f x, 0 + = 0, are not consistent, so the Terzaghi solution for m f exhibits the well known behavior of the classical unidimensional heat equation. As it will be shown in the following, this discontinuity between the initial datum and the BC is cured by the second gradient model. The Terzaghi solution given in terms of the fluid mass density m f reads

m f (x, t) = c + +∞ k=1 √ 2a sin π 2 + kπ x e λ k t , λ k = a π 2 + kπ 2 .
In Figure 7 we show the comparison between the Terzaghi and the second gradient solutions for m f and ε, respectively (in absence of prestresses), corresponding to the initial condition. The Terzaghi solution (blue line) tends to a step function due to the discontinuity between the initial datum and the BC; on the other hand this discontinuity is not present in the second gradient solution (red line). Moreover, we underline that the second gradient Fourier series solution converges more quickly to its limit compared with the Terzaghi one.

Figure 8 shows the comparison between the Terzaghi and second gradient solutions for increasing time.

It must be noted that, due to the continuity between the initial datum and the BC, the second gradient solution smoothly decreases with respect to the initial datum, while the Terzaghi solution is not able to describe the behavior of m f close to the external surface. The second gradient allows for describing compaction of the solid in the vicinity of the external surface, which contrasts instantaneous escape of the fluid out of the porous skeleton. This effect has been indeed recognized both in experiments and in situ measurements [Mollema and Antonellini 1996;Holcomb and Olsson 2003;2000].

Concluding remarks

In this paper an application of the second gradient theory of poromechanics to the consolidation problem is discussed. In particular, we present some results within the hypothesis of small deformations around a prestressed reference configuration of the solid skeleton. Even in the framework of the linearized theory, the considered second gradient model gives rise to several interesting questions, concerning both the mathematical formulation of the problem and the mechanical interpretation of the results.

From the mathematical point of view the problem could be studied within the framework of linear pencils of ordinary differential operators on a finite interval, with boundary conditions depending on the spectral parameter. Several applications of this theory to physics and mechanics can be found in the literature [Tretter 2000;Marletta et al. 2003]; it is our purpose to investigate in the future how the very special problem we are dealing with can fit within the general theory.

From the mechanical point of view, the results presented also look quite interesting, in particular concerning the capability of the model to describe fluid segregation. It has to be remarked that second gradient models, in general, regularize the solutions of evolutionary or equilibrium equations (see, for example, Figure 8). In the case of their application to phase transition phenomena they allow for the coexistence of different phases at equilibrium, in the case of strain concentration phenomena for the description of shear and compaction bands, and in the case of wetting for the description of drop/film stability. In the first and third instances, the second gradient is necessary to describe capillarity, and in A linear stability analysis provides the limit value of the initial stress, which forces the first eigenvalues to become at least nonnegative. The goal in the future will be that of identifying which are the buckling modes, and in particular to corroborate the idea that bifurcated modes of consolidation can interpret liquefaction phenomena and occurrence of sand boils [Kolymbas 1998].

Appendix A: Coefficients of the differential problem

The constant coefficients C i and i appearing in the differential problem given by Equations ( 17) and (19) are defined as with

C 1 = k 1 k 3 , C 2 p ext. 0 = k 1 + k 3 k 5 (k 2 + b) 2 + k 3 C 4 p ext.
k 1 = K ss (λ + 2µ) L 2 , k 2 = K s f , k 3 = M L 2 , k 4 = α DL 2 , k 5 = M λ + 2µ , k 6 = b 2 k 5 .
Moreover, the coefficients appearing in the BC, Equation (18), are defined as bM .

Figure 1 .

 1 Figure 1. Behavior of the fluid mass density m f versus depth x in the vicinity of the time t = 0 + . A segregation of the fluid arises close to the impermeable wall at x = 1.

Figure 3 .

 3 Figure 3. Profile of the vertical solid strain ε versus x for times close to t = 0 + . A dilatancy of the solid skeleton is detected in the neighborhood of x = 1.

Figure 5 .Figure 6 .

 56 Figure 5. The value of the first eigenvalue λ 1 versus the prestress p ext. 0 .

Figure 7 .Figure 8 .

 78 Figure 7. Comparison between the classical Terzaghi solution and the second gradient solution for t = 0 + .

Table 1 .

 1 The values of first gradient elasticity parameters relative to a normally consol-

	idated water saturated clay, together with trial values of second gradient dimensionless
	parameters.

  Figure 2. Profile of m f for further times. Notice that m f tends to assume a constant value for t → +∞, approaching equilibrium.

		0.0000	t 0
		0.0002
	m f	0.0004	t 0.5
		0.0006	t 1
		0.0008	t 2
				t 3
		0.0010	t 4
			0.0	0.2	0.4	0.6	0.8	1.0
						x
		0.0002	t 0
			t 0.02
		0.0003
		0.0004	t 0.1
	Ε
		0.0005	t 0.25
		0.0006
			t 0.5
			0.0		0.2	0.4	0.6	0.8	1.0
						x

  Figure 4. Behavior of ε for further times. Notice that the system tends to reach a state of equilibrium for t → +∞.

			6. 10 9	4	0
			4. 10 9 5. 10 9	4	0			p 0	ext Λ 2Μ
	p 0	ext	3. 10 9
			2. 10 9
			1. 10 9
			0
				0.6	0.4	0.2	0.0	0.2	0.4	0.6
							Λ 1
			0.0002	t 0
			0.0004
			0.0006	t 0.5
	Ε
			0.0008	t 1
			0.0010	t 2 t 2
				t 3
			0.0012	t 4
			0.0	0.2	0.4	0.6	0.8	1.0
							x

=0 the norm induced by the inner product • , • 4 =0 . Notice that p 0 and X 0 = V = K are the Fourier coefficient and the constant eigenfunction corresponding to the

It must be noticed that the constants k 1 , . . . , k 4 are introduced by the second gradient model, while k 5 and k 6 are related to the first gradient parameters M, λ, and µ, which represent the Biot bulk modulus and the Lamé coefficients of the considered material, respectively.