Gradient-dependent plasticity model and dynamic strain localisation analysis of saturated and partially saturated porous media: one dimensional model

Eur. J. Mech. A/Solids 19 (2000) 503–524  2000 Éditions scientifiques et médicales Elsevier SAS. All rights reserved S0997-7538(00)00177-7/FLA

Gradient-dependent plasticity model and dynamic strain localisation analysis of saturated and partially saturated porous media: one dimensional model Hong Wo Zhang a , Bernard A. Schrefler b, * a Research Institute of Engineering Mechanics, Department of Engineering Mechanics, State Key Laboratory of Structural Analysis

and Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China b Department of Structural and Transportation Engineering, University of Padua, Via F. Marzolo, 9-35131 Padua, Italy

(Received 4 March 1999; revised and accepted 20 January 2000) Abstract – Dynamic strain localisation in saturated and partially saturated porous media is investigated with a one-dimensional model in this paper. The porous medium is treated as a multiphase continuum, with the pores filled by water and air, this last one at atmospheric pressure. A gradient-dependent plasticity model is introduced to describe the plastic behaviour of the solid skeleton. Material instability due to the softening behaviour of the solid skeleton and the well-posedness of the initial value problem are studied. The advantages of the enhanced model are that the governing equations remain hyperbolic even in the softening regime and convergent solutions with mesh refinements are obtained. Moreover, the influence of permeability in the seepage process for the development of the localised zones is discussed. We find that the permeability plays an important part in the compressive wave propagation, but not in the shear wave cases. For numerical implementation of the present method, a parametric variational principle is introduced by which the original problem is reduced to a standard linear complementary problem in mathematical programming. The results of a one dimensional example are given to illustrate the efficiency of the techniques presented here.  2000 Éditions scientifiques et médicales Elsevier SAS strain localisation / gradient / dependent plasticity / porous media

1. Introduction Shear bands are strain accumulations in well defined zones which are observed in a wide class of engineering materials, including solids like metals and multiphase porous geomaterials. The problem has been analysed in the past as a material instability of a single phase rate-independent solid model (see, e.g. (Hadamard, 1903; Thomas, 1961; Hill, 1962 and Rice, 1975)). It appeared that such a model is deficient both in physical and numerical aspects. The reason has been elucidated, e.g., in (Read and Hegemier, 1984; Lasry and Belytschko, 1988) and (De Borst et al., 1993) and is due to the local change of the character of the governing equations, which results in a loss of their well-posedness when instability occurs. This means that the wave speed becomes imaginary and consequently a solution without physical significance is obtained in dynamic localisation problems. Further, the model does not contain an internal length scale and the corresponding numerical finite element solutions show a pathological mesh dependence. Various kinds of modifications and generalisations of standard continuum plasticity have been proposed to avoid the difficulties in localisation simulation of single phase materials. For instance, Needleman (1988), Asaro (1983), Loret and Prevost (1990), Perzyna (1992) and Sluys (1992) introduced rate-dependence in their constitutive models which results in a length-scale being included. For granular materials, (Mühlhaus and Vardoulakis, 1987) has developed a Cosserat model and the mean grain size diameter has been introduced into their constitutive equations. On the other hand, gradient-dependent models have been recognised by several * Correspondance and reprints.

504

H.W. Zhang, B.A. Schrefler

authors, such as Mühlhaus and Aifantis (1991), De Borst and Mühlhaus (1992) and Sluys et al. (1993) to provide a satisfactory framework for the analytical and numerical analysis of strain localisation in single phase solids. In fact, the use of a higher-order gradient model results in a well-posed set of partial differential equations. From a dispersion analysis it is clear that the governing equations remain hyperbolic under strain softening states and the continuum is capable of transforming a travelling wave into a stationary localisation wave. An internal length scale exists in the model and the numerical results do not suffer from a pathological mesh dependence. As far as multiphase materials are considered, Rice (1975), Rudnicki (1984) and Vardoulakis (1986) investigated the stability of saturated inelastic porous media with idealised initial and boundary conditions. Loret and Prevost (1991) studied strain localisation in a large mass of saturated soil in a high frequency situation, where a velocity jump was applied at the top and bottom surface of the specimen. Viscoplasticity was added to the model to regularise the numerical solution. A multiphase material model based on averaging procedures was used in (Schrefler et al., 1995) and (Schrefler et al., 1996) to study shear bands in saturated and partially saturated porous media. Ramp loading was applied to the same example of soil mass as studied in (Loret and Prevost, 1991) and the presence of negative water pressures during localisation of undrained dilatant material was evidenced. Also cavitation, which is experimentally observed (Mokni and Desrues, 1999), was modelled in (Schrefler et al., 1996). Moreover, no regularisation mechanism was introduced in the model because it was expected that the Laplacian in the mass balance equation when Darcy’s law is used would provide some sort of regularisation. Zhang et al. (1999) have shown that in multiphase materials a wave number domain exists for which the material model is well-posed when instability due to softening behaviour occurs for the solid skeleton and that an internal length scale can be derived. This scale, however, in one-dimensional dynamic problems works only for axial waves and not for ideal shear waves because, in this last case, the multiphase material model behaves like a single phase one in which no length scale exists without an additional regularisation mechanism. Regularisation of the localisation problem with a multiphase model also for ideal shear waves can be achieved only through the inclusion of one of the above mentioned regularisation methods. The purpose of the present work is to extend the gradient-dependent plasticity model to the one-dimensional analysis of shear banding in two- and three phase porous media subject to dynamic transient loading, thus obtaining a well-posed localisation problem both for axial and shear waves. This enhanced model has been chosen for its good numerical performance, while a link with experimental observation is still an open question, to the authors knowledge. An inverse analysis could be done to relate the gradient parameter to the band width. A first link between the two can be found in Vardoulakis and Sulem (1995), where the material length introduced directly in a gradient-dependent model was estimated by analysing some experimental results from a biaxial test on a fine grained sand. The material length was related to the mean grain diameter d50% . The present paper is organised as follows. First, the governing equations of a fully and partially saturated porous medium are summarised in section 2 and the associated stability problem is investigated in section 3 by means of a dispersion analysis for compressive and shear waves, respectively. Next, a description is given of the multiphase model with gradient-dependent plasticity, in which the yield function not only depends on stress and a hardening parameter, but also on the Laplacian of the latter quantity. The present model leads to a prediction of the width of the localisation zone in multiphase porous media similarly to the single phase material (De Borst and Mühlhaus, 1992). In section 5, we discuss the influence of the seepage process on the localisation results, in particular for the compressive wave propagation. It will be shown that the length scale introduced by the gradient dependent plasticity model and that naturally contained in the multiphase material model are competing ones. The interaction between the two length scales is only briefly mentioned because it is still an open question and will be analysed in a future work. Finally, a parametric variational principle and the corresponding mathematical programming algorithm are proposed in section 6 to overcome the difficulties

Gradient-dependent plasticity model and dynamic strain localisation analysis

505

caused by the incorporation of the gradient effect, thus obtaining the numerical model used for the computation. Results of numerical analyses for a one-dimensional soil bar under the loading conditions of pure shear and compression are then presented. 2. Governing equations of fully and partially saturated porous media Dynamics of wave propagation in fully and partially saturated porous media may be studied following Lewis and Schrefler (1998), where the macroscopic balance equations of mass and momentum are obtained by systematic application of averaging procedures to the relevant balance equations of the constituents at microscopic level. For our purpose the following simplifying hypotheses are introduced: isothermal conditions are assumed and gas phase (air) remains at atmospheric pressure which is taken as reference pressure (passive air phase assumption). Small strain analysis is adopted, hence the convective acceleration terms are neglected. Also the fluid (water) density variation in space is neglected and the air density is assumed to be zero because it is much smaller than the water density. The constituents are compressible. These assumptions are common in soil mechanics (Zienkiewicz et al., 1999). The general field equations are hence the linear momentum balance equation for a multiphase medium


div σ + ρ g − a s − nSw ρ w a ws = 0,

(1)

the continuity equation for water with zero mass rate of water evaporation 1 − n ∂ρ s n ∂ρ w n ∂Sw s + + div v + s w ρ ∂t ρ ∂t Sw ∂t
1 w ws + div nS ρ v =0 w Sw ρ w and the linear momentum balance equation for water

v ws =


 k 0 k wπ  − grad p w + ρ w g − a s − a ws , w µ

(2)

(3)

where σ is the total stress tensor, ρ π (π = s, w) is the phase averaged density, ρ is the average density of the multiphase medium, g is the external momentum supply related to gravitational effects and ρ π a π is the volume density of the inertia force. v s is the mass averaged velocity of the solid phase, v ws the intrinsic mass averaged velocity of the water relative to the solid phase and a ws the water acceleration relative to the solid; n is the porosity, k 0 is the intrinsic permeability tensor, µw is the dynamic viscosity and k wπ is the relative permeability, function of the degree of saturation Sw . For an isotropic material the Biot’s constant is α = 1 − Kd /Ks , where Kd and Ks are the bulk moduli of the medium and of the solid material, respectively (Zienkiewicz et al., 1999). α is taken equal to 1 in the following analysis. The capillary pressure p c is defined from the pressure equilibrium equation (Lewis and Schrefler, 1998) as p c = −p w ,

(4)

in which p w is the water pressure. Because of isothermal conditions, the saturation of liquid water is an experimentally determined function of capillary pressure only


Sw = Sw p c .

(5)

506

H.W. Zhang, B.A. Schrefler

The constitutive law for the solid phase in fully and partially saturated conditions is introduced through the concept of modified effective stress σ 00 (Bishop stress) (see, e.g., Lewis and Schrefler, 1998): σ 00 = σ + αpI ,

(6)

in which I is the unit second-order tensor and p is an average pressure in the solid phase which, by neglecting the dependence of Helmholtz free energies on void fraction (Gray and Hassanizadeh, 1991), is given as: p = Sw p w .

(7)

The classical plasticity constitutive relationship for the solid skeleton, the flow rule, the consistency and the Kuhn–Tucker complementary conditions have the incremental form of:


dσ 00 = D dε − dε p ,

(8)

dε = γ m,

(9)

p

00

f (σ , χ) 6 0,

γ > 0,

γf = 0,

(10) p

where D is the elastic constitutive tensor (or modulus for one-dimensional problem), ε and ε the total and plastic small strain tensor, respectively, γ the consistency parameter, m the plastic flow direction, f the yield function and χ the equivalent plastic strain. 3. Linear stability analysis First a fully saturated porous medium is considered. The dynamic governing equations can be expressed in the following u-U form (see also (Zienkiewicz and Shiomi, 1984)): 

00 σij,j + (α − n)Q nUj,j + (α − n)uj,j



nQ (α − n)uj,j + nUj,j



,i


 ,i




− (1 − n)ρ s u¨ i + n2 k −1 U˙ i − u˙ i = 0,


− nρ w U¨ i − n2 k −1 U˙ i − u˙ i = 0

(11)

in which ui and Ui are the solid skeleton and fluid displacements, respectively. k = kw /(ρ w g) is the isotropic total permeability of the saturated material, kw the hydraulic conductivity and Q−1 = n/Kw + (α − n)/Ks , where Ks and Kw are the bulk moduli of the fluid and solid material, respectively. To study the stability of the dynamic governing equation (11) in the plastic softening range, we consider a one-dimensional wave propagation problem for which the governing equations of motion for compressive wave (ux 6= 0, Ux 6= 0, uy = Uy = uz = Uz = 0) and shear wave (uy 6= 0, Uy 6= 0, ux = Ux = uz = Uz = 0) can, respectively, be written in the following forms: 

00 σxx,x + (α − n)Q nUx,x + (α − n)ux,x



nQ (α − n)ux,x + nUx,x



,x


 ,x




− (1 − n)ρ s u¨ x + n2 k −1 U˙ x − u˙ x = 0,


− nρ w U¨ x − n2 k −1 U˙ x − u˙ x = 0

and

(12)




00 σxy,x − (1 − n)ρ s u¨ y + n2 k −1 U˙ y − u˙ y = 0,




− nρ w U¨ y − n2 k −1 U˙ y − u˙ y = 0.

(13)

Both equations (12) and (13) can be rewritten in a standard form KY (x, t),xx = M Y¨ (x, t) + C Y˙ (x, t),

(14)

Gradient-dependent plasticity model and dynamic strain localisation analysis

507

where the coefficient matrices and vector Y are, respectively, expressed in the following forms: for the case of equation (12) 

Y (x, t) = 

M=





ux (x, t) , Ux (x, t)

K=

E + Q(α − n)2 Qn(α − n)



(1 − n)ρ s 0



0 , nρ w

C=



Qn(α − n) , Qn2 

n2 k −1 −n2 k −1

−n2 k −1 ; n2 k −1

for the case of equation (13), M, and C have the same forms as above, and 

Y (x, t) =





uy (x, t) , Uy (x, t)



E K= 0

0 . 0

Because shear wave governing equation (13) can be obtained by setting Q = 0 in (12) we concentrate now our attention on equation (12), which corresponds to the compressive wave propagation problem. For the investigation of stability, a small perturbation is applied, which induces plasticity around a ground state of fully homogeneous deformation. During such a small perturbation, the material behaviour is assumed to be described by the constitutive modulus E, which is assumed to be the same for all possible deformation continuations including harmonic perturbations. The stress–strain relationship of the solid skeleton is hence given by, σ˙ 00 = E ε˙ p , where E is a negative value to study the stability due to softening behaviour. In the following, we carry out a dispersion analysis for the one-dimensional bar, and apply a single harmonic perturbation wave propagating through the soil bar with a displacement field of the form: 

du dU





=

Au AU



ei(Kx−ωt ) = A eiKx+ζ t ,

ζ = −iω.

(15)

ω is the angular frequency and K is the wave number counting the number of wave length λ in the bar over 2π . Substitution of equation (15) into the wave equation (12) gives the following equation 

EK 2 + (α − n)2 QK 2 + n2 k −1 ζ + m1 ζ 2 (α − n)nQK 2 − n2 k −1 ζ

(α − n)nQK 2 − n2 k −1 ζ n2 QK 2 + n2 k −1 ζ + βm1 ζ 2



Au AU



=0

and m1 = (1 − n)ρ s ,

β = nρ w /m1 .

The determinant equation for eigenvalue ζ is then obtained as Det(ζ ) = a0 ζ 4 + a1 ζ 3 + a2 ζ 2 + a3 ζ + a4 = 0

(16)

in which a0 = βm21 ,

a1 = m1 (1 + β)n2 k −1 ,




a2 = m1 K 2 βE + β(α − n)2 Q + n2 Q ,


a3 = n2 k −1 K 2 E + Qα 2 ,

a4 = Qn2 EK 4 .

To study the stability of the problem, we use the Routh–Hurwitz criterion. Stability requires that the roots of equation (16) have negative real parts. This holds if and only if the coefficients of the characteristic polynomial

508

H.W. Zhang, B.A. Schrefler

in (16) satisfy a0 > 0,

11 > 0,

12 > 0,

13 > 0,

14 > 0,

(17)

where 11 = m1 (1 + β)n2 k −1 , 12 = m21 n2 k −1 K 2 β 2 E + Q β(α − n) − n


2


13 = m21 n4 k −2 K 4 βE + αQ β(α − n) − n

,



2

(18) ,

14 = m21 n6 k −2 K 8 QE βE + αQ β(α − n) − n



2

.

Obviously, the determinant conditions (17) are always satisfied when the condition E > 0 holds. On the other hand, from the condition of E < 0, follows the loss of stability because 14 becomes negative. This implies that the initial homogeneous state is unstable and shear band instability can develop. Further when E < 0 we have Det(0) < 0,

Det(+∞) > 0

(19)

which means that then there is always a positive root in equation (16). To consider fully and partially saturated conditions we use the following simpler u-p form of equations, which can be obtained from the above governing equations by neglecting the fluid acceleration relative to the solid phase. This form of equation is valid for dynamic behaviour in the low frequency range (Zienkiewicz and Shiomi, 1984), and can be expressed as 00 σij,j − αSw pj − ρ u¨ i = 0,

(20)

αSw u˙ i,i − (kp,i ),i + p/Q ˙ = 0, where Q is now 



n∂Sw nSw Sw (α − n) ∂Sw Q= + Sw + + p ∂p Kw Ks ∂p

−1

.

By using the same procedure as shown in the saturated case, we obtain the determinant equation of the problem: Det(ζ ) = b0 ζ 3 + b1 ζ 2 + b2 ζ + b3 = 0,

(21)

where b0 = m = ρ,

b1 = mQkK 2 ,

b2 = K 2 E + α 2 Sw2 K 2 Q,

b3 = kK 4 EQ.

(22)

The determinant conditions in the Routh–Hurwitz criterion are 10 = b0 ,

11 = mQK 2 k,

12 = mα 2 Sw2 kK 4 Q2 ,

13 = mα 2 Sw2 k 2 K 8 Q3 E.

(23)

Again, the stability is controlled by the sign of E, as in the saturated case. In fact when E < 0 the coefficient 13 is negative and hence loss of stability occurs.

Gradient-dependent plasticity model and dynamic strain localisation analysis

509

4. Fully and partially saturated porous media with gradient-dependent plasticity model The framework of the gradient-dependent continuum is set by the following equations:


dσ 00 = D dε − dε p ,

dε p = γ m,




f σ 00 , χ, ∇ 2 χ 6 0,

f γ = 0,

(24)

γ > 0.

To demonstrate the salient features of the gradient-dependent model presented here, we consider a onedimensional model for which the constitutive equation for a stress point on the yield surface can be written in rate form as (Sluys et al., 1993): σ˙ 00 = h˙ε p − c∗

∂ 2 ε˙ p ∂x 2

(25)

in which c∗ = ∂f/∂∇ 2χ is a constant and h is the strain softening modulus. Combining equations (24), (25) and (11) or (20), and taking D = E (elastic modulus), we can obtain the dynamic wave propagation equations for a one-dimensional gradient-dependent strain softening saturated and partially saturated porous medium in the following form: for u-U saturated model: 





E + Q(α − n)2 nQ(α − n) u c∗ E+h 0 0 U ,xxxx       s ∗ (1 − n)ρ 0 u¨ −n2 k −1 n2 k −1 u˙ c c∗ − + ¨ E+h 0 0 U ,xx E + h U˙ ,xx 0 0  Eh     2 −1   + Q(α − n)2 nQ(α − n) u nk −n2 k −1 u˙  − E+h + 2 −1 2 −1 U ,xx U˙ −n k n k nQ(α − n) n2 Q    (1 − n)ρ s 0 u¨ + =0 w 0 nρ U¨

(26)

and for u-p saturated and partially saturated model: −

∗



E c E+h 0 

0 0



u p





+ ,xxxx





∗



0 αSw c E+h 0 0 

ρ 0 u¨ 0 −αSw c∗ + E+h 0 0 p¨ ,xx 0 0       0 0 u˙ ρ 0 u¨ + + = 0. 0 1/Q p˙ 0 0 p¨

+

u p



u p





+ ,x



 Eh − E +h

,xxx

0 αSw

0 0





0

0 −k −1  u˙ p˙ ,x





u p

 ,xx

(27)

4.1. Linear stability analysis Again we consider a linear solid (equations (26) and (27)), with the stress–strain relationship given by equation (25). As in section 3 we then investigate the general solution for a single harmonic perturbation wave propagating through a one-dimensional bar. The eigenvalue determinant equations for the two different models are then, respectively, changed into:

510

H.W. Zhang, B.A. Schrefler

for saturated case (u-U model): Det(ζ ) = a0 ζ 4 + a1 ζ 3 + a2 ζ 2 + a3 ζ + a4 = 0,

(28)

where a0 = βm21 , a1 = m1 (1 + β)n2 k −1 ,





a2 = m1 K 2 β h + c∗ K 2 + Q (α − n)2 + n2 ,

(29)




a3 = n2 k −1 K 2 h + c∗ K 2 + Qα 2 , a4 = n2 K 4 Q h + c∗ K 2




m1 = (1 − n)ρ s and for fully and partially saturated case (u-p model): Det(ζ ) = b0 ζ 3 + b1 ζ 2 + b2 ζ + b3 = 0,

(30)

in which b0 = m, b1 = mkK 2 Q,




b2 = α 2 Sw2 K 2 Q + K 2 h + c∗ K 2 , b3 = kK 4 Q h + c∗ K

(31)


2

m = ρ. Consequently, the corresponding conditions for the two models in the Routh–Hurwitz criterion can be expressed as 10 = βm21 , 11 = m1 (1 + β)n2 k −1 ,




12 = m21 n2 k −1 K 2 β 2 h + c∗ K 2 + Q β(α − n) − n 13 = m21 n4 k −2 K 4 β h + c∗ K

2



2


+ αQ β(α − n) − n

14 = a4 13 = m21 n6 k −2 K 8 Q h + c∗ K


2

β h + c∗ K

,



2


2

(32)

,

+ αQ β(α − n) − n



2

and 10 = b0 , 11 = mkK 2 Q, 12 = mα 2 Sw2 kK 4 Q2 ,

(33)


13 = b3 12 = mα 2 Sw2 k 2 K 8 Q3 h + c∗ K 2 . If h + c∗ K 2 6 0 the Routh–Hurwitz criterion is violated and hence loss of stability due to softening can develop. It is interesting to note that this instability is dependent on the wave number K. In fact the unstable

Gradient-dependent plasticity model and dynamic strain localisation analysis zone appears because of waves with wave number K 6 hence a cut-off limit given by K>

p

−h/c∗

and

511

√ −h/c∗ . The wave numbers and the wave lengths have

λ < 2π lg ,

where lg =

p

−c∗ / h

(34)

from which it can be observed that only wave lengths smaller than 2π/K can propagate in the shear band. The strain softening regions remain small however and no wave lengths larger than 2π lg can occur because they do not fit in this region (Sluys, 1992). Similarly to single phase material model (De Borst et al., 1997), the parameter lg can be defined as the internal length scale of the model associated with the gradient effect. In the examples section, we will show that the width of the localisation zone is the value of wave length λ = 2π lg . For this wave length the damping term vanishes from the equation of the problem. This was shown by (Sluys, 1992; Sluys et al., 1993) for a single phase material and holds also for our case, see equation (52), when the length scale parameter depends only on E1 defined in equation (36). For compressive waves, however, condition (34) may not be the unique condition to predict the length scale in multiphase materials. In fact, in this case, the parameter lg of equation (34c) can not be taken directly as the length scale because it may be necessary to take into account the influence of the seepage process through permeability and other parameters, i.e. of another competing length scale, as we will explain in the following. 5. Loss of hyperbolicity — the influence of permeability k Seepage process is known to introduce, in compressive wave propagation, an internal length through Darcy’s equation (Zhang et al., 1999). Dispersion analysis is now carried out to take into account its influence in the localisation phenomenon. In the following, we limit our analysis to the simpler u-p model whose governing equations are (20). The determinant equation (30) can be written in the standard form: ζ 3 + aζ 2 + bζ + c = 0

(35)

in which a = a 0 Kx,

a 0 = Q,

b = b0 K 2 ,

b0 = E1 + α 2 Sw2 Q m−1 ,




c0 = QE1 m−1 , E(c∗ K 2 + h) E1 = . E + c∗ K 2 + h

c = c0 K 3 x, x = kK,

(36)

Remarks The difference between multiphase problems such as the saturated and unsaturated problems studied here and a classical (strain softening) single phase solid is that here there is not only a positive real root when strainsoftening occurs, but also two other roots whose properties are controlled by parameters such as permeability k and wave number K besides the constitutive modulus E1 . When k → 0, equation (35) will reduce to:



ζ ζ 2 + K 2 m−1 E1 + α 2 Sw2 Q

= 0.

(37)

512

H.W. Zhang, B.A. Schrefler

In this case, the loss of hyperbolicity will be controlled by the following parameter: φ = E1 + α 2 Sw2 Q

(38)

which could be understood as the constitutive modulus of the gradient dependent porous medium in undrained condition. For shear wave propagation, Q = 0, hence E1 will be the control parameter to ensure that the equation maintains hyperbolicity when strain localisation occurs in the continuum. We consider now the case when all the eigenvalues in equation (35) are non-complex (i.e. ω purely imaginary): the problem will lose hyperbolicity, the positive real solution will further induce the loss of stability of the dynamic equations and mesh dependence of the finite element analysis of the localisation problem will appear. To ensure that all the roots in (35) will be non-complex, we have the determinant condition  3

p W= 3 where p=−

a2 + b, 3

 2

+

q 2

6 0,

 3

q=2

a 3

−

(39)

ab + c. 3

Substitution of (36) into (39) yields W (x) = wx 4 + rx 2 + z 6 0,

(40)

where 1 03 0 1 2 1 1 02 02 1 3 r = c0 − a 0 b0 c0 − z = b0 , a c , a b , 27 4 6 108 27 are the coefficients which do not depend on permeability and wave number when c∗ = 0. The general solution of the equality of (40) will be √ −r ± r 2 − 4wz 2 x = . (41) 2w w=

It should be noticed that (41) is a non-linear equation due to the inclusion of wave number K in the variable E1 , and just holds explicitly when c∗ = 0. Due to the fact that w < 0 always holds when instability occurs, and E1 > 0 means that the problem will be stable, without loss of generality, we only consider the following two cases (1) z > 0, E1 < 0, and (2) z 6 0, E1 < 0, respectively, with c∗ = 0 for simplicity. Also, it is equivalent to substitute z > 0 or z 6 0 by φ > 0 or φ 6 0 (see (38)). Case 1. φ > 0, E1 < 0 From equation (40), we find that the condition (39) will hold with the following solution √ −r − r 2 − 4wz 2 x > (42) = ξ 2 > 0. 2w Hence, when x > ξ,

K>

ξ k

(43)

Gradient-dependent plasticity model and dynamic strain localisation analysis

513

the wave speed will become purely imaginary, and mesh sensitivity will appear. In this case, the wave number domain where the hyperbolicity still holds can be defined as: 




R1S = K | K ∈ K1S0 , K1S1 ,

K1S0 = 0,

K1S1 = ξ /k.

(44)

Case 2. φ 6 0, E1 < 0 It is easy to prove that r > 0 and r 2 − 4wz > 0 always hold. In this case, the problem will have not only an upper bound solution, but presents also lower bound limitations which are: √ −r − r 2 − 4wz ξ1 2 x > K> , = ξ12 , 2w k (45) √ −r + r 2 − 4wz ξ2 2 2 x 6 K6 . = ξ2 , 2w k Hence the wave number domain where hyperbolicity of the problems still holds can be expressed as follows: 

R2S = K | K ∈ K2S0 , K2S1




K2S0 = ξ2 /k,

,

K2S1 = ξ1 /k.

(46)

It can be observed that, for case 1, when k is very large, the critical wave number K1S1 will be very small and hence the finite element results may be mesh dependent because all the wave numbers K may exceed K1S1 . For case 2, the upper bound K2S1 will be very small for very large permeability, while the lower bound K2S0 will be very large for very small k. So there must be two critical values K0 and K1 which limit the admissible values of wave number, i.e.: 




R0S = K | K ∈ K0S0 , K0S1 ,

K0S0 = K0 ,

K0S1 = K1

(47)

and mesh sensitivity should occur both for very large and for very small permeability. For partially saturated problems, sometimes Sw can become very small when the soil is at low saturation. This will lead to φ < 0 and the solution of (45) holds in case 2. For the general case of any value of φ a wave number domain where hyperbolicity of the governing equations is guaranteed when instability occurs can finally be defined as: 

Ri = RiS ∩ R0S = K | K ∈ S0


Ki0 = max K0S0 , Ki

,



Ki0 , Ki1




,


Ki1 = min K0S1 , KiS1 ,

(48)

where subscripts i = 1, 2 denote case 1 and case 2, respectively. When Ri = 0, there will be loss of hyperbolicity, and the gradient-dependent model proposed will be useful to overcome the mesh sensitivity. To predict the internal length scale in the case of Ri 6= 0, it is necessary to consider the damping effects by the negative part of the complex root. For case 1, in a recent work (Zhang et al., 1999) an approximation is derived for the length scale prediction for compressive waves (Q 6= 0) controlled by the parameter k when R1 is not zero and c∗ = 0, i.e.: lapprox =

2cm η kα 2 Sw2 Q

(49)

514

H.W. Zhang, B.A. Schrefler

and

s

cm =

φ , ρ

η=

φ , Q

φ = α 2 Sw2 Q + E1

(50)

in which cm can be understood as the wave speed of the mixture. As already noticed, R2 will be zero when k tends to a small value or large value, i.e. the domain of R2S will be occupied by very small or large wave numbers K. In this case, the governing equation without gradient dependent plasticity will lose hyperbolicity. The gradient dependent model will ensure that the dynamic equation still remains hyperbolic, and the upper and lower bounds of the length scale corresponding to the two limit permeability cases (small and large) can be, respectively, predicted by the wave lengths of the following wave numbers: 

Kup = −

h + h0 c∗

1/2



,

Klower = −

h c∗

1/2

,

h0 =

Eα 2 Sw2 Q , E + α 2 Sw2 Q

(51)

where Kup and Klower have been obtained from (35). We must admit that the problem will be more complicated when the permeability is located at intermediate values where the gradient-dependent parameter and permeability both matter. The region bounds of strain localisation can not be predicted by (51) as shown in a numerical example. We do not pursue this problem further. The interaction of two competing length scale is still an open problem. When Q → 0, equation (35) changes into the following determinant equation for shear wave analysis:


ζ ζ 2 + E1 m−1 = 0

(52)

hence E1 (equation (36)) will be the control parameter to ensure that the equation maintains hyperbolicity when strain localisation occurs in the continuum. Permeability hence does not matter and (34) holds for length scale prediction is shear wave propagation. In particular when c∗ = 0 and h < 0, the dynamic governing equations will lose hyperbolicity. The results presented in this paper are in agreement with those of Vardoulakis and Sulem (1995) and Vardoulakis (1996), where the stability and the regularisation of the particular case of undrained deformations is analysed for deviatoric softening poro-elastoplastic soil, assuming incompressible constituents. In fact, in case of zero-order flow theory (classical theory), Vardoulakis and Sulem (1995) and Vardoulakis (1996) have found the ill-posedness of the set of governing equations, while the second-grade flow theory resulted in instability with regularised governing equations. In the present paper these results are obtained removing the incompressible restraints of the constituents and are also extended to the case of drained deformations. 6. Parametric variational principle and mathematical programming solution The parametric variational principle applicable to elastic-plastic coupled field problems in consolidation analysis of saturated porous media was developed in Zhang (1995). This principle can be used to solve problems where materials do not satisfy Drucker’s postulate of stability, such as in non-associated plasticity or softening problems. In the following, we will apply the method presented to the gradient-dependent model in strain localisation analysis of saturated and partially saturated porous media. A small modification is needed to extend the formulation developed in Zhang (1995) to the problems here because the Laplacian of the equivalent plastic strain variable is introduced in the yield function. The parametric variational principle corresponding to the governing equations (20) for the u-p model can be described as

Gradient-dependent plasticity model and dynamic strain localisation analysis Y



p = 1ui,j , t ∗ Dije kl 1uk,l

− 1∗ k1p,i , t ∗ 1p,i





∗ e  − 2 1ui,j , t γ Dij kl







− 2 1p, t ∗ αSw 1ui,j (0)

+ 2 1∗ k1p, t ∗ 1p−,i ni

∗

− 2 1u,i , t 16−i



01



− t ∗ 1p, 1p/Q







∂g ∂σij00

−2 ∗



 

1ui , t ∗ 1Fi

 + 2 1P /Q, t 1p(0)



04

+ 2 1p, t ∗ 1∗ k1f−i,i



∗

− 2 1ui , t αSw 1p− ni



+ 2 1ui,j , t ∗ αSw 1p





 

515

 



+ 1ui , p1ui

 

(53)



04



02

and the corresponding control equations are υ˙ + ni σ˙ i − hχ˙ + c∗ η∇ 2 χ˙ = 0, γ , υ˙ > 0,

(54)

γ υ˙ = 0,

where η is a constant, ni is the gradient of yield surface, υ is a slack parameter complementing the control parameter γ which does not participate as an argument in the variation. The explanation of other variables and the mathematical proof of the functional (53) can be found in Zhang (1995). For the spatial discretisation of the problem, the body can be divided into a finite number of elements and in each element the continuous displacement, pressure and plastic multiplier γ fields can be interpolated by p

1ui = Nku 1uki ,

γ

1p = Nk 1p k ,

γ = Nk γ k .

(55)

The fact that second derivatives of γ enter the yield functions makes it necessary to select C 1 -continuous shape functions for the interpolation of γ . Substituting equation (55) into (53) and (54), and with the first variation of the functional (53), i.e. Q

Q

∂ p = 0, ∂1ukl

∂ p =0 ∂1p k

(56)

we obtain 

Muu 0

0 0



1u¨ 1p¨





+

0 −Qpu

0 Spp



1u˙ 1p˙





+



−Qup + Spp

Kuu 0



1u 1p





G {γ } = 0



1fu 1fp



(57)

and [C]{1u} + [M]{γ } + {υ} − {d} = 0, {γ }T {υ} = 0, where Kuu = Qup = Spp = C=

Z Z Z Z







d =

NL Q−1 Nk d, p



Muu =

u u NL,j Dije kl Nm,l d, p u NL,i αSw Nk

u Wklβ Nm,l

p

d,

(58)

{γ }, {υ} > 0,



NLu ρNmu d,

Hpp =

QTpu , G=− M=

Z

Z 

Z 

Z



p

p

NL,i kNk,i d,

u NL,j Dije kl

mβ Nmγ d +

∂gβ d, ∂σij00

Z



c∗ η∇ 2 Nmγ d,

516

H.W. Zhang, B.A. Schrefler 1fu = − 1fp = − d= Wklβ =

Z 02

Z



d0 +

p

04

Z

NLu αSw 1p− ni

NL kni 1p−j d0 + mβ =

fβ0 d,

Z 01

Z

NLu 16−ij nj

d0 +

p

04

NL k1f−i,i d0 +

Z 

Z 

NLu 1Fi d,

p

NL,i 1fi d,

∂fβ ∂gβ ∂fβ e ∂gβ − D − hβ , p ∂εij ∂σij00 ∂σij00 ij kl ∂σij00

∂fβ e D . ∂σij00 ij kl

Discretisation in time of equation (57) is carried out by means of the Newmark scheme  1  1u − 1t u(0) ˙ − 0.51t 2 u(0) ¨ , 2 β1t  γ  1u˙ = 1t u(0) ¨ + 1u − 1t u(0) ˙ − 0.51t 2 u(0) ¨ , (59) β1t  1  1p˙ = 1p − 1t p(0) ˙ , θ1t where γ , β and θ are integration parameters, and generally, have γ , θ > 0.5 and β > 0.25. Substitution (59) into (57) leads to the following equation

1u¨ =

K ∗ S = Y, where



T

S = 1uT 1p T ,  Muu + Kuu  β1t 2  ∗ K = 

Y=



β

−Qpu

β 1fp γ



−Qup Spp θ

  , + 1tHpp 

 1 Muu  1tβ  +  

γ

   1fu  

(60)



0

1tβ Spp γ From equations (58) and (60), we obtain −1tQpu

 

υ 0



+

M D

B K∗



γ S

  







=

u(0) ˙ p(0) ˙

d Y





0.5 Muu  β  +   β − 0.5 − 1t 2 Qpu γ

0 0

   





u(0) ¨ . p(0) ¨



γ T υ = 0,

,

γ , υ > 0,

(61)

where B = [C

D = [G 0].

0],

Equation (61) can be transformed into  

υ 0



U − BK ∗−1 D + K ∗−1 D

γ T υ = 0,

γ , υ > 0.

0 I



γ S





=

d − BK ∗−1 Y K ∗−1 Y



, (62)

Gradient-dependent plasticity model and dynamic strain localisation analysis

517

Then we need just to solve the linear complementary problem 





{υ} + U − BK ∗−1 D {γ } = d − BK ∗−1 Y



γ T υ = 0,

γ,υ > 0

(63)

which can easily be solved by many methods such as the Wolf’s method and Lemke’s method (see Boot (1964)) for standard linear and quadratic programming problems. 7. Numerical results To illustrate the dispersive character of the gradient-dependent model and the influence of the permeability for localisation results in porous media, a one-dimensional soil bar in pure compressive and shear deformation is selected and tested numerically in this section. The geometrical, material and loading data for the soil bar are given in figure 1. Von Mises plasticity with gradient-dependence is used for simplicity, which reduces for the one-dimensional problem to 00 f = −σxx − σs00 − hχ + c∗ ∇ 2 χ.

(64)

First, we investigate the behaviour of the model with respect to mesh refinement with a permeability k = 10−12 m3 s/kg (near to zero). For the bar, we use 20, 30, 40, 50 and 60 linear elements, respectively. We consider a block wave (t0 = 0) travelling through the bar in linearly elastic regime until reflection occurs and the stress increase causes the initiation of the localisation process. Extra boundary conditions ∂γ /∂x = 0 at both sides of the column are used. The situation Q 6= 0 and σs = 1.1 MPa is used for the compressive wave propagation, while we set Q → 0 and σs = 1.8 MPa to simulate the shear wave propagation in the bar with the same data shown in figure 1, i.e. G = E, σxx ⇒ σxy .

Figure 1. One-dimensional soil bar in pure compression and shear loading. L = 10 m, A = 1 m2 , q0 = 1.0 MPa, td = 0, E = 28.5 MPa, E1 (c∗ = 0) = −25.0 MPa or −5 MPa Q = 20 MPa, ρ = 2000 Kg/m3 , σs = 1.1 MPa or 1.8 MPa.

518

H.W. Zhang, B.A. Schrefler

Figure 2a. Plastic strain along the bar using the gradient-dependent porous media model (permeability k = 10−12 m3 s/kg) for compressive waves and different meshes.

Figure 2b. Plastic strain along the bar using the gradient-dependent porous media model (permeability k = 10−12 m3 s/kg) for shear waves (Q = 0) and different meshes.

Because a very small permeability is adopted, the width of the localisation zone is based on a gradient model and it can be obtained from the wave length whose corresponding wave numbers have been given by equation (51)1 and (34), respectively. For pure compression and shear wave simulations, we use c∗ = 0.4 MN and c∗ = 4.0 MN, respectively, and the width of the localisation zone will be equal to λ = 3.176 m and λ = 3.443 m. In fact the first value is obtained by using equation (51):

Gradient-dependent plasticity model and dynamic strain localisation analysis

Figure 3a. Development with time of the localisation zone for compressive waves using permeability of k = 10−12 m3 s/kg.

Figure 3b. Development with time of the localisation zone for shear waves (Q = 0) using permeability of k = 10−12 m3 s/kg.

EE1 28.5 × 25.0 × 106 =− = −13.318 × 106 , E − E1 28.5 + 25.0 EQ 28.5 × 20.0 × 106 h0 = =− = −11.753 × 106 , E +Q 28.5 + 20.0 h=

519

520

H.W. Zhang, B.A. Schrefler

Figure 4. Development of the water pressure for different meshes.

Figure 5. Development of the water pressure for different values of time.

s

K=

h + h0 − ∗ = c

s

13.318 − 11.753 = 1.978, 0.4

λ = 2π/1.978 = 3.176.

In figures 2a and 2b, the strain profiles are plotted for different meshes to demonstrate the convergence of the solution to a finite width of the localisation zone for compressive and shear waves. In the figures 3a and 3b, the development of the plastic zone increases after reflection but latter, the speed of extension of the zone vanishes

Gradient-dependent plasticity model and dynamic strain localisation analysis

521

Figure 6. Development of the localisation zone using permeability of k = 10−6 m3 s/kg and c∗ = 0 [N] for different values of time.

and the localisation bands of constant width w ∼ = 0.5λ (because of the symmetry) arise for the two different waves, respectively. Similar results have been obtained by Sluys et al. (1993) for one phase materials. Figure 4 shows the results of water pressure by different meshes and figure 5 depicts the development of water pressure with time for the compressive wave. Clearly, due to the decreasing stiffness of the solid skeleton, the water pressure in the softening zone increases rapidly. To illustrate the influence of the permeability, an elastic modulus E1 (c∗ = 0) = 5 MPa is used in the following computation so that the conditions of case 2 discussed in section 5 will be satisfied. At first, a permeability k = 10−6 m3 s/kg is adopted and we do not use the gradient-dependent model and just consider the influence of permeability. The plastic strain development is shown in figure 6 where the width of the fully developed softening zone is nearly half of the value lapprox ∼ = 6.5 predicted by (49). To show the influence of the two length scales, the model with both permeability and gradient dependent plasticity is used. The result obtained with k = 10−6 , k = 10−4 and different gradient dependent parameters are shown in figures 7 and 8, respectively. It can be concluded that in the case of smaller permeability (see figure 7), the plastic zone (length scale) governed by permeability will prevail in the computational process, while in the case of larger permeability (see figure 8), the length scale will be controlled by gradient dependent plasticity. The width of plastic zone in figure 8 is consistent with the prediction of (51b) (half of the value) due to the use of a larger permeability. The case corresponding to (45) and (46) is also studied here by a numerical example. In this case the gradient dependent model is still used and the parameter c∗ = 0.4 MN is adopted. The width of the localisation zone predicted by (51) will be equal to λup = 3.176 m and λlower = 1.088 m for lower and higher permeability, respectively. Figure 9 presents the plastic strain results for different permeability values. The results are consistent with the limits: the upper and lower localisation band width is close to half of the predicted value as usual in gradient dependent plasticity (Sluys et al., 1993). The results also show that for an intermediate value of the permeability the problem will be more complicated, the width of localisation zone being somewhere in the region between the lower and upper bounds.

522

H.W. Zhang, B.A. Schrefler

Figure 7. Distribution of plastic strain along the bar with different values of c∗ [N] and permeability of k = 10−6 m3 s/kg.

Figure 8. Distribution of plastic strain along the bar with different values of c∗ [N] and permeability of k = 10−4 m3 s/kg.

8. Conclusions Theoretical and computational aspects in dynamic strain localisation of fully and partially saturated porous media have been investigated with a one-dimensional model. It is shown that permeability through Darcy’s law already introduces some sort of length scale in case of compressive waves. However, to obtain true mesh independence for axial and shear waves, gradient dependent plasticity has been used. For this case, the well

Gradient-dependent plasticity model and dynamic strain localisation analysis

523

Figure 9. Distribution of plastic strain along the bar with different values of permeability [m3 s/kg].

posedness of the problem has been shown in the softening range and the presence of an internal length scale has been pointed out. The influence of the seepage process through the permeability is discussed and the interplay of the two length scales has been addressed. For the numerical analysis, a parametric variational principle is used, by which the analysis is transformed into a linear complementary problem. Numerical examples confirm that the gradient dependent solution converges to a finite width of the localisation zone, which agrees with the analytical value for very large and for very small values of the permeability. The numerical results also demonstrate the influence of permeability and for its intermediate values, where the two length scales complete, an interaction between the two different length scales obtained by Darcy’s law and gradient dependent model is observed. This interaction will be in different forms for the different conditions of φ > 0 and φ 6 0, where φ is understood as the modulus of the mixture in undrained condition. In the former case, for a small permeability, the length scale will be mainly controlled by Darcy’s law when the parameter c∗ in gradient dependent model is relatively small, while for a large permeability the length scale will be mainly decided by gradient dependent model. In the later case, i.e. φ 6 0, the lower and upper values of length scale have been given to predict the width of shear band corresponding to large and small permeability. An intermediate value of permeability will interact with the parameter c∗ in the gradient dependent model and influence the width of shear band, and thus the length scale will be somewhere between the lower and upper values. It is clear that for the case where both permeability and gradient dependent parameter assume intermediate values, the exact shear band width is difficult to predict, and this problem is still open. Acknowledgements This work was supported by the EC International Scientific Coorperation Programme and by research funds M.U.R.S.T. 020902019. The support of the State Educational Committee and National Natural Science Foundation of China (19872016) is also acknowledged.

524

H.W. Zhang, B.A. Schrefler

References Asaro R.J., Mechanics of crystals and polycristals, Adv. Appl. Mech. 23 (1983) 2–115. Boot J.C., Quadratic programming, Amsterdam, North-Holland, 1964. De Borst R., Mühlhaus H.B., Gradient-dependent plasticity: Formulation and algorithmic aspects, Int. J. Num. Meth. Engng. 35 (1992) 521–539. De Borst R., Wang W.M., Geers M.G.D, Sluys L.J., Material instabilities and internal length scales, in: Owen D.R.J., Onate E., Hinton E. (Eds.), Computational Plasticity – Fundamentals and Applications, CIMNE, Barcelona, 1997. De Borst R., Sluys L.J., Mühlhaus H.B., Pamin J., Fundamental issues in finite element analyses of localisation of deformation, Eng. Comput. 10 (1993) 99–121. Gray W.G., Hassanizadeh S.M., Unsaturated flow theory including interfacial phenomena, Water Resour. Research 27 (1991) 1855–1863. Hadamard J., Leons sur la propagation des ondes et les equations de l’hydrodynamique, Librairie Scientifique A. Hermann, Paris, France (in French), 1903. Hill R., Acceleration waves in solids, J. Mech. Phys. Solids 10 (1962) 1–16. Lasry D., Belytschko T.B., Localization limiters in transient problems, Int. J. Solids Structures 24 (1988) 581–597. Lewis R.W., Schrefler B.A., The finite element method in the static and dynamic deformation and consolidation of porous media, J. Wiley, Chichester, 1998. Loret B., Prevost J.H., Dynamic strain localization in elasto (visco-) plastic solids. Part 1. General formulation and one-dimensional examples, Comp. Meths. Appl. Mech. Engng. 83 (1990) 247–273. Loret B., Prevost J.H., Dynamic strain localisation in fluid-saturated porous media, J. Engng. Mech. 11 (1991) 907–922. Mokni M., Desrues J., Strain localization measurements in undrained plane-strain baxial tests on Hostun RF sand, Mech. Cohes.-Frict. Mater. 4 (1999) 419–441. Mühlhaus H.B., Aifantis E.C., A variational principle for gradient plasticity, Int. J. Solids and Structures 28 (1991) 845–858. Mühlhaus H.B., Vardoulakis I., The thickness of shear band in granular materials, Geotechnique 37 (1987) 271–283. Needleman A., Material rate dependence and mesh sensitivity on localization problems, Comp. Meth. Appl. Mech. Eng. 67 (1988) 69–86. Perzyna P., Constitutive equations of dynamic plasticity, in: Owen D.R.J., Onate E., Hinton E. (Eds.), Computational Plasticity. Fundamentals and Applications, Pineridge Press, 1992, 483–508. Read H.E., Hegemier G.A., Strain softening of rock, soil and concrete – A review article, Mech. of Mat. 3 (1984) 271–294. Rice J.R., On the stability of dilatant hardening for saturated rock masses, J. Geophys. Res. 80 (1975) 1531–1536. Rudnicki J.W., Effect of dilatant hardening on the development of concentrated shear-deformation in fissured rock mass, J. Geophys. Res. 89 (1984) 9259–9270. Sluys L.J., Wave propagation, localization and dispersion in softening solids, Ph.D. Thesis, Civil Engineering Department of Delft University of Technology, 1992. Sluys L.J., De Borst R., Mühlhaus H.B., Wave propagation, localization and dispersion in a gradient-dependent medium, Int. J. Solids and Structures 30 (1993) 1153–1171. Schrefler B.A., Majorana C.E., Sanavia L., Shear band localization in saturated porous media, Archives of Mechanics 47 (1995) 577–599. Schrefler B.A., Sanavia L., Majorana C.E., A multiphase medium model for localization and post localisation simulation in geomaterials, Mech. of Cohes.-Frict. Mat. 1 (1996) 95–114. Thomas T.Y., Plastic Flow and Fracture in Solids, Academic Press, New York, N.Y., 1961. Vardoulakis I., Dynamic stability analysis of undrained simple shear on water-saturated granular soils, Int. J. Numer. Anal. Meth. Geomech. 10 (1986) 177–190. Vardoulakis I., Sulem J., Bifurcation Analysis in Geomechanics, Blackie Academic & Professional, 1995. Vardoulakis I., Deformation of water-saturated sand: I. uniform undrained deformation and shear band. II the effect of pore water flow and shear banding, Géotechnique 46 (1996) 441–456, 457–471. Zhang H.W., Parametric variational principle for elastic-plastic consolidation analysis of saturated porous media, Int. J. Numer. Anal. Meths. Geomech. 19 (1995) 851–867. Zhang H.W., Sanavia L., Schrefler B.A., An internal length scale in dynamic strain localization of multiphase porous media, Mech. of Cohes.-Frict. Mat. 4 (1999) 443–460. Zienkiewicz O.C., Shiomi T., Dynamic behaviour of saturated porous media, the generalized Biot formulation and its numerical solution, Int. J. Numer. Anal. Meths. Geomech. 8 (1984) 71–96. Zienkiewicz O.C., Chan A., Pastor M., Schrefler B.A., Shiomi T., Computational Soil Dynamics with Special Reference to Earthquake Engineering, J. Wiley, Chichester, 1999.