- Email: [email protected]
Numerical analysis of dynamic strain localisation in initially water saturated dense sand with a modified generalised plasticity model
Computers and Structures 79 (2001) 441±459
www.elsevier.com/locate/compstruc
Numerical analysis of dynamic strain localisation in initially water saturated dense sand with a modi®ed generalised plasticity model H.W. Zhang a, L. Sanavia b, B.A. Schre¯er b,* a
State Key Laboratory of Structure Analysis for Industrial Equipment, Research Institute of Engineering Mechanics, Dalian University of Technology, Dalian 116024, People's Republic of China b Department of Structural and Transportation Engineering, University of Padua, via F. Marzolo 9, 35131 Padua, Italy Received 6 May 1999; accepted 1 May 2000
Abstract Shear band dominated process in fully and partially saturated sand samples is simulated by means of dynamic strain localisation analysis together with a multiphase material model. The partially saturated medium is viewed as a multiphase continuum consisting of a solid skeleton and pores ®lled by water and air (vapour) which, once it appears, is presumed to remain at the constant value of cavitation pressure (isothermal monospecies approach). The governing equations are based on the general framework of averaging theories. A modi®ed generalised plasticity constitutive model for partially saturated soils, developed from the general Pastor±Zienkiewicz sand model, has been implemented in a ®nite element code and used in the computational process. This model takes into account the eects of suction in the stiness of the porous medium (solid skeleton) in partially saturated state. A case of strain localisation, which has been tested in laboratory observing cavitation of the pore water, is studied. Negative water pressures, which are of importance in localisation phenomena of initially fully saturated undrained samples of dilatant geomaterials, are obtained similarly to those observed experimentally. Ó 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction Strain localisation is strain concentration in wellde®ned narrow zones, where the material behaviour is inelastic, while the remaining zones are elastic. This phenomenon is observed in a wide class of engineering materials, including single-phase solids like metals, and multiphase ¯uid saturated and partially saturated porous geomaterials like soils and concrete. The triggering mechanisms for the formation of shear bands are
* Corresponding author. Tel.: +39-49-8275611; fax: +39-498275604.
inhomogeneities in the material and stress concentrations, which may be, for multiphase geomaterials, due to the ¯uids±solid interaction. For such kinds of materials, the role of the ¯uids can be fundamental, as evidenced by extensive laboratory tests on dense and loose sands made by Mokni [1] and Mokni and Desrues [2]. In fact, in the case of undrained monotonic biaxial compression tests of initially water saturated dense Hostun sand under constant con®ning pressure, strain localisation was only observed when pore-water pressure of the order of ±80 kPa (value relative to the atmospheric pressure) developed, irrespective of the value of the applied back pressure [1,2]. (This drop in the water pressure is due to the dilatant behaviour of the solid skeleton associated with the impervious boundary of the sample which precludes the in¯ow or out¯ow of water.) The authors of the experiments
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 4 4 - 9
442
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Nomenclature (Variables with overbar refers to the nodal values) s solid phase g gaseous phase w liquid phase p generic phase ap acceleration of p-phase aps acceleration relative to the solid B strain operator Cs speci®c capacity CT tangential constitutive tensor eective diusivity tensor Dg fp ¯ow vector fu external load vector g gravity acceleration H permeability matrix I unit tensor krp ¯uid phase relative permeability k absolute or intrinsic permeability tensor Ks solid grain bulk modulus Kw liquid phase bulk modulus KT tangential stiness matrix molar mass of constituent p Mp M mass matrix n porosity n unit normal vector pp macroscopic pressure of the p-phase P equivalent force vector Q coupled matrix R universal constant gas Swo initial water saturation Sw water saturation Sg gas saturation S compressibility matrix t time variable t surface traction tensor
observed cavitation 1 of the pore-water with a production of vapour and stated that in this case, ``the
1 Cavitation of water may occur when the absolute value of water pressure is equal or lesser than the saturation water pressure at the considered temperature (neglecting the surface tension of the interface of the arising vapour bubbles). For thermodynamic equilibrium state, which is assumed in the model used in this paper at local level, the vapour pressure at the interface between saturated and partially saturated porous material (where the capillary pressure pc is zero) is equal, at T 20°C, to the saturation value of 2.339 kPa, (i.e. ±98.986 kPa, with reference to the atmospheric pressure).
tp T u vap
partial stress tensor temperature solid displacements velocity of the a-phase with respect to the pphase vp velocity of the p-phase a Biot's coecient c1 , c2 , H Newmark's parameters convective mass transfer coecient bc lw liquid dynamic viscosity q porous medium density qp intrinsic phase averaged density of the pphase r Cauchy stress tensor r0 eective Cauchy stress tensor I1 ®rst eective stress invariant J2 second invariant of the deviatoric eective stress J3 third invariant of the deviatoric eective stress ngL direction vector of plastic ¯ow for loading direction vector of plastic ¯ow for unloadngU ing g stress ratio Mg slope of the critical state line gmax maximum stress ratio reached gu stress ratio from which unloading takes place ag ,af ,Mf ,H0 ,cDM material parameters in Pastor± Zienkiewicz model b0 ,b1 ,Hu0 ,cu ,Kevo ,Keso material parameters in Pastor± Zienkiewicz model RI , b1 , b2 material parameters in modi®ed Pastor± Zienkiewicz model
onset of localisation is delayed until cavitation takes place in the pore-water'' and that ``non-drainage can preclude localisation as long as cavitation in the pore¯uid does not relax the isochoric constraint'' [2]. A similar experimental conclusion can be found in Ref. [3], for which ``a change in the draining conditions for a fully or nearly fully saturated soil can activate unstable behaviour and this may occur due to static or dynamic disturbances'', such as that done, for the latter case, in the numerical experiments showed in this paper. The fact that cavitation happens at the onset of localisation or shortly afterwards is still a controversial point. In this paper, we are not going to discuss this
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
open question, also because Mokni and Desrues observed cavitation at a pore pressure of ±80 kPa, which is dierent from the theoretical value of ±98.986 kPa (at 20°C) (or the corresponding value in the range of ambient temperatures). In this respect, it is worth to note that McManus and Davis [4] observed strain localisation and the associated cavitation during undrained triaxial tests of dense sand at the pore pressure value of approximately ±95 kPa and that Vardoulakis [5,6] revealed cavitation in globally undrained biaxial tests of dense St. Peter sandstone at the pore pressure of ±51 kPa Ref. [5,6]. The important role that ¯uids can have in localisation modelling is con®rmed also in Ref. [7], where the bifurcation analysis of undrained medium dense St. Peter sandstone is performed. The authors found the existence of bifurcation modes only in the form of locally drained shear band, which means the necessity to include the water ¯ow in a strain localisation modelling (see also Refs. [5,6]). In the present paper, we start from the abovementioned important theoretical and experimental role of the ¯uids in strain localisation analysis to justify the use of a multiphase model, where, air, dissolved air in water and water vapour, should be taken into account together with water and the solid skeleton. Starting from the ®rst numerical analysis of fully saturated porous media made by Shuttle and Smith [8] and Loret and Prevost [9] in statics and in dynamics respectively, several authors have worked on this subject worldwide. As far as some original aspects of our research group are concerned, a numerical model based on the general frame of averaging theories and capable of simulating also the shear band dominated processes in saturated porous media was presented in Refs. [10±12]. The model was then extended in Refs. [13,14] to deal with the cavitation phenomenon of an initially water saturated and globally undrained medium, where the so-de®ned monospecies approach and the isothermal two phase ¯ow model have been respectively proposed. When simulating the behaviour of dense sand, a Mohr±Coulomb constitutive relationship with a softening branch was used in Refs. [10,13] together with associative plasticity and in Refs. [11,12] with non-associative plasticity and Pastor± Zienkiewicz model for sand. The main focus of this paper lies in the use of a modi®ed generalised plasticity model developed in Ref. [15] in the context of the numerical simulation of strain localisation problems of initially water saturated dense (dilatant) sand. This constitutive model for the solid skeleton takes into account capillary eects which arise once cavitation takes place. This paper is organised as follows: ®rst, in Section 2, a brief summary of the adopted multiphase material model for geomaterials where phase change (cavitation) can take place is presented. To this end, the monospecies approach [13,16] used for the computations is recalled, because it is a
443
simple model with a good performance. The main features of the modi®ed generalised plasticity model developed by Bolzon et al. [15] are then brie¯y described in Section 3, where capillary pressure is taken as independent variable, beyond the eective stress, because this is a necessary aspect of partially saturated soil models [17]. This soil model has been implemented in the ®nite element code SWANDYNE [18], as described in Section 4. Numerical examples of shear band localisation within dense sand samples are shown in Section 5. In particular, the strong in¯uence that ¯uid phases can have in localisation phenomena is investigated from a numerical point of view by a comparison between the results of the monospecies approach and the water saturated model performed on a classical biaxial test. Mohr±Coulomb law for the solid skeleton was used for the sake of simplicity. Then, the results of a biaxial test under symmetric and unsymmetric loading cases taken from Refs. [1,2] are shown, using the modi®ed Pastor±Zienkiewicz model. These results show the superiority of the adopted constitutive model in reproducing the experimental pore pressure values during shear banding, because they are much closer to the observations by Mokni [1], Mokni and Desrues [2], McManus and Davis [4] and Vardoulakis [5,6] than the previous attempts using a Mohr±Coulomb constitutive relationship [10,13]. It is also worth to outline that negative pore pressures were not observed in Ref. [17] when using the basic generalised plasticity model for saturated sands, from which the model used here is derived, because loose (contractive) sands were studied. Finally, we recall the well-known fact that in localisation analysis of a single phase material with softening, one obtains an ill-posed problem. This implies a strong dependence of the shear band width on the spatial discretisation used once the numerical analysis is performed. To overcome this unphysical behaviour, different regularisation strategies can be found in the literature (see e.g. Ref. [19] for a review). On the contrary, multiphase material contains a natural regularisation mechanism due to the presence of a Laplacian term in the mass balance equation of the ¯uids when DarcyÕs law is used. The numerical properties of this natural regularisation related to axial waves have been investigated in Ref. [13], while the internal length scale contained in the model for the 1-D and 2-D cases and for axial waves is presented in Refs. [20,21] respectively. Hence, this aspect will not be further pursued in this paper. 2. Simpli®ed model for simulation of cavitation: isothermal monospecies approach The mathematical model for a multiphase porous medium is now recalled. It can be derived from the
444
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
general thermo±hydro-mechanical model developed in Ref. [22] by introducing the assumptions outlined next. A short summary of the general model is presented in Appendix A for interested readers. Small strains are assumed for the development of the following equations. We consider here globally undrained situations, where cavitation has been experimentally observed [1,2,4±6]. In such conditions, neglecting air dissolved in water and air bubbles in water, which may be present in fully saturated specimens at low pressures, only two ¯uid phases are present after desaturation: liquid water and water vapour. These assumptions can be made because the authors of the experimental tests removed air bubbles as far as possible and saturated the samples using de-aired water [2,4]. This means that the air mass balance equation (A.15) in not needed. Moreover, isothermal conditions are assumed for the sake of simplicity. An isothermal approach means that physically heat is supplied or taken away with in®nite intensity to maintain constant temperature (this implies in®nite heat capacity value), thus there are no energetic restrictions for phase change (evaporation or condensation). Our isothermal model consists ®nally of the linear momentum balance equation for the multiphase medium (A.14) and the mass balance equation for water species (A.16) with the diusion term omitted (there is no diusion of water vapour in the air). Also the Clausius±Clapeyron equation is not needed in KelvinÕs law (A.7), because the water saturation pressure pgws is constant in isothermal conditions. The monospecies approach is now derived [13] as follows. Neglecting further the vapour density qgw with respect to the water density qw , from the state equation of vapour, Eq. (A.5), it follows that also the vapour pressure pgw is negligible. Hence, the average pressure (A.10) is given by p Sw pw ;
1
Further, the dynamic seepage forcing term in Eq. (A.13), connected with the solid acceleration, is very small when compared with other terms of the equation system, so it could be neglected [13]. Thus, the ¯uid mass balance equation (A.16) has now the following form: ÿ
o nSw qw Sw qw div vs ot kk rw div qw w grad pw ÿ qw g 0; l
3
where qw is the water density, vs , the solid velocity vector, k, the intrinsic permeability tensor, krw , the water relative permeability, and lw , the dynamic viscosity. Together with Eqs. (2) and (3), we need the experimentally determined constitutive equations of the medium: Sw Sw pc ;
k rw k rw Sw
4
with pc being the capillary pressure. Finally, Eq. (A.8) reads pc ÿpw ;
5
which states that negative water pressure corresponds to the development of capillary pressure. The relationship between capillary pressure and water degree of saturation is of the type depicted in Fig. 1. In fact, because of the conditions of the experiments and the simplifying assumptions, partial saturation is reached only at cavitation, which occurs at minus the atmospheric pressure using pressures relative to the atmospheric pressure, which in turn is set to zero. For the model closure, the appropriate initial and boundary conditions (A.17b,d), (A.18b,c) and (A.20) are also needed.
where Sw is the water degree of saturation and pw is the water pressure. By ignoring the relative component of the ¯uid accelerations, the primary variables may be reduced to solid displacements u and water pressure pw . After application of the eective stress equation (A.9) to the linear momentum balance equation (A.14), we obtain its simpli®ed form as div r0 ÿ a Sw pw I q g ÿ as 0;
2
where r0 is the Bishop eective stress tensor, a 1 ÿ Kd =Ks , the BiotÕs constant, with Kd and Ks represently the bulk moduli of the medium and of the solid material respectively, q, the mixture density, g and as , the gravity acceleration and the solid acceleration vector, respectively. I is the second-order identity tensor.
Fig. 1. Saturation-water pressure relationship used for cavitation modelling.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
3. The modi®ed Pastor±Zienkiewicz generalised plasticity sand model The generalised plasticity soil model was introduced by Pastor et al. [23], which was extensively validated for fully saturated sand behaviour. It has been shown by comparison with experiments that this soil model is capable of a good performance for a wide range of conditions (see Ref. [18] for a review). The interested reader can ®nd the basic features of the Pastor±Zienkiewicz model in Appendix B. The modi®ed model will be now described following Ref. [15], who extended the Pastor± Zienkiewicz model to partially saturated sands. As mentioned in Ref. [15], both BishopÕs stress and capillary pressure (or other two combination of stresses and pressures [17]) are necessary to be introduced as independent stress parameters to describe the behaviour of the partially saturated soil. In the modi®ed model, eects of capillary pressure are taken into account by the soil stiness changes induced by suction and irreversible compressive volumetric strains ev of soil upon wetting. (
Wun1 Mn1 Dun Pn1 ÿ Qn1 H Dt Dp_ n ÿ Fun1 0 Wpn1 QTn1 c1 Dt Dun Hn1 H Dt Dp_ n Sn1 Dp_ n ÿ Fpn1 0;
The initial yield limit of the fully saturated case of the Pastor±Zienkiewicz model is hence modi®ed in partially saturated conditions to introduce its dependence on suction as [15] 0 pyi0 pc py0i RI pc ;
6
0 where py0i represents the initial yield limit in saturated conditions of the mean eective stress. RI is a material parameter which expresses an increasing function of suction when water saturation is less than one and is determined by interpolation of experimental data. Because irreversible volumetric strain is assumed as the parameter which controls the volumetric hardening, the evolution of the yield surface in the (p0 , pc ) plane is given as !1apc p0 y0 0 c 0 c py p py0i RI p 7 0 py0i
445
main by means of the NewmarkÕs scheme. The unknown ®eld variables are expressed in the whole domain by global shape function matrices, Nu and Nw , and the nodal value vectors, u and p: u Nu u;
pw Nw p:
8
Once the coupling matrix Q, the mass matrix M, the strain operator B and the external load f u , are introduced (see Appendix C), the linear momentum balance can be written as Z BT r0 dX ÿ Qp Mu f u ; 9a X
while mass balance equation becomes Hp QT u_ Sp_ f p ;
9b
where H is the permeability matrix, S, the compressibility matrix and, f p , the ¯ow vector (Appendix C). After the introduction of the NewmarkÕs scheme for time integration, equations system (9) can be written at the current time tn1 as
10
where c1 , H and c2 , which appear in Eq. (11), are the NewmarkÕs parameters, and the vectors Fun1 and Fpn1 are de®ned in the Appendix C. pn1 must be evaluated by integration of the constitutive law. The non-linear coupled equation system is linearised in a standard way thus yielding the linear algebraic equations system to be solved numerically: ! M 12 KT c2 Dt2 ÿQH Dt Du ÿH ÿQT H Dt Dp_ c1 HH Dt S u ÿW ÿ H Wp : 11 c1
0 where py0 is the initial mean stress and a is a material parameter.
Since the Newton±Raphson method requires the Jacobian matrix to be evaluated and inverted at each iteration, other modi®ed schemes are also used to achieve convergence with less computational eort. In particular, the use of secant updates, like DavidonÕs and Broyden±Fletcher±Goldfarb±ShannoÕs (BFGS) methods is found advantageous in non-linear analyses [18].
4. Finite element discretisation and solution
4.1. Initial static solution for modi®ed Pastor±Zienkiewicz model
The governing Eqs. (2) and (3) are discretised in space using a GalerkinÕs procedure and in the time do-
The modi®ed Pastor±Zienkiewicz model needs the initial value of the elastic parameters depending on the
446
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
initial mean pressure. The form of the governing equations used to obtain the initial value for the modi®ed constitutive model are the following static drained equations: Z BT r00 dX ÿ Q0 p0 f u0 ; 12a X
H0 p0 f p0 ;
12b
where Q0 and H0 are matrixes depending on p0 and Sw0 . Numerical experiments have shown that inconsistency between Eqs. (9a), (9b) and (12a), (12b) will generate an unstable dynamic computational process with Eqs. (9a) and (9b), even if the dynamic load is zero. The reason is that the material parameters for static analysis in Eqs. (12a) and (12b) are dierent from those in Eqs. (9a) and (9b), the parameters in the latter ones being totally regenerated with the de®nitions in the modi®ed Pastor± Zienkiewicz model. To avoid this problem, we move the contribution of the initial state to the right-hand side of the governing equations and rewrite the dynamic equations (9a) and (9b) in the following incremental forms: Z ~ ; BT Dr0 dX ÿ Q Dp M Du Df u Qp 13a 0 X
~ 0; H Dp QT Du_ S Dp_ Df p ÿ Hp
13b
in which the incremental variables are de®ned as Dx x ÿ x0 , Q, H and S are the matrixes dependent on the current water pressure and saturation. The expres~ H ~ can be found in the Appendix C. sions of Q,
5. Numerical results Before we deal with the simulation of a laboratory test taken from Refs. [1,2], the in¯uence of the ¯uid phases in localisation phenomena is investigated from a numerical point of view. At this end, a rectangular mass of water saturated dense sand under compression is analysed by means of the isothermal monospecies approach (IMA) and the water saturated model (WSM). This example has been studied by several authors [9± 14,16]. The sample has impervious boundaries (Fig. 2) and hence the problem is globally undrained. Vertical and horizontal displacements are constrained at the bottom surface. Ramp loading is applied at the top. During computations, the self weight of the soil and water ®lling its pores are taken into account. Hydrostatic distribution of water pressure is assumed as an initial condition. The sample is subjected to plane strain condition. The relationships between capillary pressure, saturation and water relative permeability proposed by Safai and Pinder [24] for San Fernando sand are used when cavitation develops (Fig. 1). A Mohr±Coulomb yield criterion with associative ¯ow rule and isotropic linear softening is used for the solid skeleton. The material parameters used during the computations are Young's modulus E 285 MPa, Poisson's ratio m 0:4285, solid grain density qs 2000 kg mÿ3 , water density qw 998:2 kg mÿ3 , apparent cohesion c0 1:84 MPa, hardening modulus H ÿ40 MPa, angle of internal friction u 20°, initial porosity n 0:20, solid grain bulk
Fig. 2. Description of the geometrical characteristics and of the load function of the ®rst example.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
447
Fig. 3. Contours of equivalent plastic strain at dierent times with the monospecies approach model: (a) 0.245 s, (b) 0.18 s and (c) 0.22 s.
modulus Ks 6:78 GPa, water bulk modulus Kw 0:20 GPa, permeability kw 0:25 10ÿ3 m sÿ1 , BiotÕs constant a 1:0. A mesh of 396 (18 elements in width and 22 over the height) four-node isoparametric ®nite elements of equal size, the same for displacements and water pressures, has been used. Computations are performed by use of the modi®ed version [13] of the SWANDYNE code [18]. The two models present signi®cant dierences. Since the saturated model neglects the partial saturation, the generalised stiness of the multiphase medium and the coupling between solid skeleton and ¯uids, i.e. the structural ¯uids±solid interaction, are not properly described because the matrices Q, H and S of Eq. (10) are constant (see Appendix C). As a consequence, the evolution of the process and all the state variables change considerably. In fact, the numerical results reveal a big delay in the formation of the localised bands when the saturated model is used (from 0.245 to 0.40 s, see Figs. 3(a) and 4 for the IMA and the WSM, respectively), even if the shear bands start to appear at the same time
(@0.18 s) because of the same velocity of the ®rst plastic wave. In this case, the localised zones start from the corners of the (rigid) bottom surface when the stress induced by the re¯ected loading wave reach the yield
Fig. 4. Contours of equivalent plastic strain at 0.40 s with the water saturated model.
448
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
value of the material model, and propagate versus the top surface (Fig. 3(b) and (c); Fig. 3(a)±(c) are obtained with the monospecies approach model). Then also a larger value of the equivalent plastic strain (0.0448 for the WSM, 0.0441 for the IMA) is observed, while the shape of the shear bands is the same. As far as the behaviour of the ¯uid phase is concerned, an increase of the zone which should become partially saturated can be observed for the WSM (see Fig. 6 in comparison with Fig. 5), with a decreasing of the maximum water pressure (1.62 MPa instead of 2.81 MPa) and an increase of the minimum water pressure (ÿ11.1 MPa instead of ÿ26.9 MPa). Moreover, the ¯ow rate of the water is double in the saturated case, because permeability does not decrease with the capillary pressure. From this comparison between the two models, it can be concluded that ¯uid phases are important for localisation phenomena of initially saturated multiphase dense geomaterials under globally undrained conditions. This conclusion is indirectly con®rmed by the experimental results [2], for which the presence of the pore
¯uid in globally undrained conditions delays the shear band formation. The examples presented next refer to the biaxial laboratory test of Refs. [1,2] and make use of the modi®ed Pastor±Zienkiewicz model for sand to take into account the partially saturated conditions appearing during localisation at cavitation pressure. With this model, also the constitutive tangent matrix KT is dependent from the capillary pressure. The ¯uid±solid interaction is hence described also at constitutive level. The geometrical characteristics with impervious boundary conditions, which are exactly the same as in Refs. [1,2], and the load function, are shown in Fig. 7. Vertical and horizontal displacements are constrained at the bottom surface. Ramp loading is applied at the top of the sample. The permeability in saturated condition is 0:25 10ÿ3 m sÿ1 , while, once cavitation starts, the relationship between the degree of saturation and the relative permeability for S. Fernando sands [24] has been assumed. Also the capillary pressure-saturation relationship is that given in Safai and Pinder [24] for S. Fernando sands, as shown in Fig. 1 [16]. These relationships have been used here for sake of simplicity because of the lack of experimental results of Hostun sands. Both a symmetric and an unsymmetric load case are chosen here to simulate the shear band formation. The material parameters for the dense sand are water bulk modulus Kw 0:20 GPa, solid grain bulk modulus Ks 6:78 GPa, initial porosity n 0:20, solid density qs 2000 kg mÿ3 , water density qw 998:2 kg mÿ3 . The initial conditions according to Refs. [1,2] are rxx ryy ÿ100 kPa, rzz ÿ130 kPa, rxy rxz ryz 0
Fig. 5. Contours of negative water pressure at 0.245 s with the monospecies approach model.
Fig. 6. Contours of negative water pressure at 0.400 s with the water saturated model.
Fig. 7. Description of the geometrical characteristics and of the load function of the second example.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
449
Fig. 8. Contours of equivalent plastic strain at dierent times for the symmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 4.0 ´ 10ÿ3 s.
Fig. 9. Contours of equivalent plastic strain at dierent times for the symmetric load case and initial water pressure p0 200 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 2.5 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
and p0 200 or 300 kPa. The material parameters for the modi®ed Pastor±Zienkiewicz for sands are as following: Mf 1:0, Mg 1:0, af 0:45, ag 0:45, Kevo 1050 Pa, Keso 2000 Pa, b0 4:2, b1 0:2, Hu0 600 Pa, H0 350 Pa, cu 2:0, cDM 2:0, b1 0:77, b2
16:4, RI 1:0 [11]. The load values for symmetric and unsymmetric load cases are F10 153 kN, F20 5:67 kN. To study the problem in detail, the following three cases are simulated:
450
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 10. Contours of equivalent plastic strain at dierent times for the unsymmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
1. Symmetric load case with initial water pressure, p0 300 kPa; 2. Symmetric load case with initial water pressure, p0 200 kPa;
3. Unsymmetric load case with initial water pressure, p0 300 kPa. As in the laboratory tests, in the numerical simulation, the shear band starts from the top surface of the specimen, near loading area, for both symmetric and unsymmetric load cases, due to the stress concentration caused by the lower stiness of the soil in that part with respect to the remaining part of the sample. The development of shear band for the symmetric load and unsymmetric load cases are shown by the contours of equivalent plastic strain variable depicted in Figs. 8±10 for three dierent times, for cases 1±3, respectively. In the unsymmetric load case (Fig. 10), the shear band propagates in a direction controlled by the disturbing load. The pattern of this shear band is similar to what was observed in the laboratory tests [1] before the failure of the tested sample. For the symmetric load case, two shear bands appear in the upper region of the specimen, which are symmetric because of the perfect symmetry of the solved problem. The eective plastic strain along a vertical (central) cross-section and for dierent times is shown in Figs. 11±13 for the three investigated cases, respectively. These ®gures give an idea of the values involved in Figs. 8±10. The pore-water pressure distributions at dierent times are presented in Figs. 14±16 for cases 1±3, respectively. Because the uniform water pressure is applied on the lateral boundary of the sample, the results of water pressure are not regular. However, water pressure concentration with negative values of water pressure (traction) (Figs. 17±19) can be observed at the centre of the shear bands due to the plastic dilatant behaviour of the dense sand. The water pressures distribution along a vertical cross-section through the centre line for dierent time values are shown in Figs. 17±19 for cases 1±3, respectively. These help to understand their distributions in the three previous ®gures. In particular, a traction value close to ÿ100 kPa is obtained in both the symmetric and unsymmetric load cases when shear bands develop. These values of the water traction are quite similar to the vapour saturation pressure of ÿ98.986 kPa and the negative pressure observed in Refs. [4±6] at cavitation. In the numerical results, they are slightly larger because of the development of higher capillary pressure. Moreover, from Figs. 17 and 18, it can be seen that the same negative pressure develops, irrespective of the value of the applied back pressure. This is also in agreement with the experimental ®ndings of Refs. [1,2,4]. In fact, Mokni and Desrues found that strain localisation developed during biaxial tests when cavitation developed, irrespective of the value of the applied back pressure. Similarly, also McManus and Davis [4] reported the same (cavitation) ®nal value (of ±95 kPa) during their triaxial tests, starting from dierent values of the initial water pressure. Hence, with the constitutive model used, which takes into account the capillary ef-
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
451
Fig. 11. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 300 kPa.
Fig. 12. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 200 kPa.
fects, the water tractions do not exceed values which are close to the experimentally observed ones. The situation with a Mohr±Coulomb model [10,11,13±16], (see also Fig. 6) without capillary eects was quite dierent, because there, much higher capillary pressures were obtained. This points to the need to choose an appropriate constitutive model for simulating localisation and postlocalisation behaviour. Moreover, the ®gures for eective plastic strain (Figs. 11±13) with the corresponding ®gures for pressure evolution (Figs. 17±19) are useful also to analyse the shear banding interaction with the pressure evolution. In fact,
it can be clearly seen that negative pore pressures arise at the same instance when localisation appears in the strain distribution. Finally, Fig. 20 shows the deformation of the sample when the shear band occurs for the three cases studied. 6. Conclusions Numerical analysis of dynamic strain localisation in initially water saturated dense sand samples using a multiphase material model has been presented. The
452
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 13. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for unsymmetric load case and initial water pressure p0 300 kPa.
model is capable of simulating the partially saturated behaviour of globally undrained dense sand which develops because of experimentally observed cavitation. The partially saturated medium is viewed as a two-phase continuum consisting of a solid skeleton and pores ®lled with water or water and vapour at cavitation pressure. The formulation is based on the general frame of averaging theories. For partially saturated sands, a modi®ed generalised plasticity constitutive model based on Pastor±Zienkiewicz model was adopted in the computational process. This is currently one of the simplest and yet one of the most eective models for describing the full range of saturated and partially saturated sand behaviour. Using this constitutive model, the numerical simulation of multiphase geomaterials is able to describe the ¯uid±solid interaction both at the structural and constitutive levels. A dense sand sample similar to a laboratory test was studied in detail with ramp loading to demonstrate the capability of the numerical model presented here to handle such situations. An agreement between the experimental observations and numerical results was obtained as far as the form of the shear band (unsymmetric case) and the pore pressures are concerned. These pressures are the most intriguing aspect of the experiment.
Acknowledgements This work was supported by the EC International Scienti®c Cooperation Programme and the National Natural Science Foundation of China (19872016) and by the research funds MURST 020902019 from the
Italian Ministry of Scienti®c and Technological Research (Co®nanziamento Murst 1998). Appendix A A.1. Mathematical model for the description of a multiphase porous media The partially saturated porous medium is a material with an internal microstructure and is treated here as a multiphase material where the voids of the solid skeleton are ®lled with water and gas which may be either vapour alone or a mixture of dry air and vapour. The full mathematical model necessary to simulate thermo± hydro-mechanical transient behaviour of fully and partially saturated porous media was developed in Ref. [22] using averaging theories following Hassanizadeh and Gray [25±27] and Gray and Hassanizadeh [28]. It assumes non-polar and immiscible constituents (except for the vapour and the gas) and local thermal equilibrium state. The most important equations are brie¯y summarised in what follows, assuming small strains. The state of the medium will be described by water pressure pw , capillary pressure pc , temperature T and displacement vector of the solid matrix u. The material free macroscopic balance equations for the generic p-phase are the following: Mass balance equation: oqp qp div vp qp ep : ot Linear momentum balance equation:
A:1
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
453
Fig. 14. Contours of water pressure at dierent times for the symmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 4.0 ´ 10ÿ3 s.
Fig. 15. Contours of water pressure at dierent times for the symmetric load case and initial water pressure p0 200 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 2.5 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
div tp qp gp ÿ ap qp ep ^tp 0;
exchange between constituent p and the other constituents, tp , the stress tensor (called partial stress tensor in mixture theories), g, the external momentum supply related to gravitational eects, qp ap , the volume density of
A:2
where qp is the phase average density, vp , the mass average velocity of the p-phase; qp ep represents the mass
454
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
implies symmetry of the stress tensor. The energy balance equation for the p-phase may be written as qp
oEp tp : dp qp hp ÿ div qp qp Rp ; ot
A:3
where qp Rp represents the exchange of energy between the p-phase and other phases of the medium due to phase change and mechanical interaction, qp is the internal heat ¯ux, hp results from the heat sources and dp is the rate of the strain tensor. Ep accounts for the speci®c energy of the volume element and for the kinetic energy. The entropy inequality for the mixture is X okp qp q hp A:4 qp qp ep kp div p ÿ p p P 0; ot T T p where kp is the speci®c entropy of the constituent p and, ep kp , the entropy supply due to mass exchange. The moist air (g) in the pore system is assumed to be a perfect mixture of two ideal gases, dry air (ga) and water vapour (gw). The equation of state of perfect gas (ClapeyronÕs equation) is hence valid: pga qga RT =Ma ;
pgw qgw RT =Mw ;
A:5
where T is the absolute temperature, Mp , the molar mass of constituent p and, R, the universal gas constant. Further, DaltonÕs law applies, giving the gas pressure and density as pg pga pgw ;
qg qga qgw :
A:6
Due to the curvature of the meniscus separating water from gas inside the pores of the medium considered as a capillary porous body, the equilibrium vapour pressure pgw can be obtained from the Kelvin relationship: pc Mw pgw pgws exp ÿ w ; A:7 q RT where pgws is the vapour saturation pressure which can be obtained from the Clausius±Clapeyron equation [16], and pc is the capillary pressure which is de®ned from the sorption equilibrium equation as pc pg ÿ pw :
Fig. 16. Contours of water pressure at dierent times for the unsymmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
the inertia force and ep ^tp accounts for momentum exchange due to mass supply and mechanical interactions with other phases. As for the angular momentum balance equation, the assumption of microscopically non-polar constituents
A:8
pp (p g, w) is related to tp by [28] tp ÿgp pp I, with gp as the p-phase volume fraction, and I is the second-order identity tensor. The constitutive law for the solid phase in the partially saturated state is introduced through the concept of modi®ed eective stress r0 (Bishop stress): r0 r apI; s
A:9 w
g
where r t t t is the total stress tensor, a, the BiotÕs constant and, p, an average pressure in the solid phase which, by neglecting the dependence of Helmholtz free energy on void fraction [27,28] is given as
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
455
Fig. 17. Development of water pressure along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 300 kPa.
p Sw pw 1 ÿ Sw pg :
A:10
Sw is the saturation of liquid water, which is an experimentally determined function of capillary pressure and temperature. The constitutive relationship for the solid skeleton has the form dr0 CT de ÿ deT ÿ de0 ;
A:11
where CT is a tangent constitutive tensor, deT I b3s dT is the strain increment caused by thermoelastic expansion, bs is the cubic thermal expansion coecient of the solid and de0 denotes the autogenous strain increments and the irreversible part of the thermal strains tensor. A particular form of this relation will be discussed in Appendix B. Finally for a binary gas mixture, FickÕs law gives the following relative average velocities vgp of the diusing
Fig. 18. Development of water pressure along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 200 kPa.
456
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 19. Development of water pressure along the vertical direction in the central cross-section of the sample for unsymmetric load case and initial water pressure p0 300 kPa.
species, due to the existing gradients of their concentrations: ga Ma Mw p vga ÿ D grad g 2 M pg gw Ma Mw p Dg grad A:12 ÿvgw ; M2 pg where Dg is the eective diusivity tensor and M is the molar mass of the gas mixture. A.2. Governing equations Using the Macroscopic balance equations and the constitutive equations of Appendix A.1, we obtain the following governing equations, as in Lewis and Schre¯er [22]. The linear momentum balance equation for ¯uids yields the generalised Darcy equation: nSp vps
kk rp ÿ grad pp qp g ÿ as ÿ aps ; lp
A:13
where p ga, gw, w and g, vps is the (intrinsic massaveraged) velocity relative to the solid, aps , the acceleration relative to the solid, n, the porosity of the medium, k, the intrinsic permeability tensor, lp , the dynamic viscosity and krp , the relative permeability which is a function of the degree of saturation and of the temperature. The linear momentum balance equation for the multiphase medium is divr q g ÿ as ÿ nSw qw aws ÿ nSg qg ags 0;
A:14
where q is the averaged P density of the multiphase medium q 1 ÿ nqs p6s nSp qp . The dry air mass balance equation, after introduction of DarcyÕs law and FickÕs law and neglecting acceleration terms is o kk rg nSg qga Sg qga div vs ÿ div qga g grad pg l ot gw Ma Mw p div qg 0: A:15 Dg grad M2 pg The vapour mass balance equation and the water mass balance equation are summed to obtain the mass balance equation for all water species. DarcyÕs law and FickÕs law are introduced as above, yielding rg o gw gw s gw kk g nSg q Sg q div v ÿ div q grad p ot lg gw Ma Mw p ÿ div qg Dg grad M2 pg o ÿ nSw qw ÿ Sw qw div vs ot rw w kk w w s ws div q grad p ÿ q g ÿ a ÿ a : lw A:16 The macroscopic mass balance equation for the solid has already been summed with the above mass balance equations to eliminate the time derivative of the porosity n. The energy conservation equation (enthalpy balance) is not used in the simpli®ed isothermal monospecies approach model and is hence omitted.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
457
The boundary conditions can be imposed values on Cp or ¯uxes on Cqp , where the boundary is C Cp [ Cqp . The imposed values on the boundary for gas pressure, water pressure and displacements are as follows: pg p^g on Cg ; pw p^w on Cw ; u^ u on Cu for t P t0
A:18a±c
The volume average ¯ux boundary conditions for water species and dry air conservation equations to be imposed at the interface between the porous medium and the surrounding ¯uid are the following: qga vg ÿ qg vgw n qga on Cqg ; A:19a g
qgw vg qw vw qg vgw n g ÿ gw gw bc q ÿ qgw qw 1 q
on Cqc ;
A:19b
where n is the unit vector, perpendicular to the surface of the porous medium, pointing toward the surrounding gas, qgw 1 is the mass concentration of water vapour in the undisturbed gas phase distant from the interface, and bc is the convective mass transfer coecient, while qga , qgw and qw are the imposed dry air ¯ux, imposed vapour ¯ux and imposed liquid ¯ux, respectively. Eq. (A.19) are the natural boundary conditions, respectively, for the dry air conservation Eq. (A.19a) and water species conservation Eq. (A.19b), when the solution of these equations is obtained through a weak formulation of the problem, as is usually done with the ®nite element method. The traction boundary conditions for the displacement ®eld are rn t
on Cqu for t P t0 ;
A:20
where t is the imposed traction. Appendix B Fig. 20. Deformation of the sample when shear band occurs: (a) Symmetric load case, initial water pressure p0 300 kPa at time 4:0 10ÿ3 s. (b) Symmetric load case, initial water pressure p0 200 kPa at time 3:5 10ÿ3 s. (c) Unsymmetric load case, initial water pressure p0 300 kPa at time 3:5 10ÿ3 s.
A.3. Initial and boundary conditions The initial conditions specify the full ®elds of gas pressure, water (or capillary) pressure, temperature and displacements: pg p0g ; u u0
pw p0w ; T T0 ; at t t0 :
A:17a±d
B.1. Pastor±Zienkiewicz generalised plasticity model for sand The approach developed by Pastor et al. [18,23] and Bolzon et al. [15] for saturated and partially saturated soils, respectively, explicitly de®nes direction vectors in terms of the following three invariants p0 , q and h p I1 p0 ÿ ; q 3J2 ; 3 ! p 1 ÿ1 3 3 J3 h ÿ sin ; B:1 2 J23=2 3 where I1 is the ®rst eective stress invariant, J2 and J3 are the second and the third invariant of the deviatoric eective stress, respectively. The de®nitions and
458
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
parameters that need to be determined in the modi®ed Pastor±Zienkiewicz model can be summarised as follows. Loading direction vector: n np ; nq ; nh T
B:2
p
q
h
with n 1 af Mf ÿ g, n 1, n ÿ 1=2 qMf cos 3h, where af , Mf are model parameters and g is the stress ratio de®ned as q=p. Plastic ¯ow direction vector for loading: ngL npgL ; nqgL ; nhgL T
B:3
with npgL 1 ag Mg ÿ g, nqgL 1, nhgL ÿ 1=2qMg cos 3h, in which ag is a material parameter and Mg is the slope of the critical state line, related to the residual 6 sin / friction angle by Mg 3ÿ sin . / sin 3h Plastic ¯ow direction vector for unloading: ngU npgU ; nqgU ; nhgU T
B:4
with npgU ÿjnpgL j, nqgU nqgL , nhgU nhgL . This is chosen in order to ensure that densi®cation occurs in unloading. Plastic modulus for loading: B:5
with 1ÿ
g gf
4 ;
gf
1 1 af
Mf ;
g ; Mg ÿcDM g HDM ; gmax
Hv 1 ÿ
Hs b0 b1 exp ÿb0 n;
Hw 1 b1 exp ÿb2 pc ÿ 1ÿ1
B:6
where Hu0 and cu are material parameters and gu is the stress ratio from which unloading takes place. Elastic moduli: The bulk and shear moduli, respectively, are K
Kevo 0 p; p00
G
Keso 0 p; 3p00
B:7
where both Kevo and Keso are model parameters and p00 is the initial value of p0 . Appendix C
0;
0 T
Fun1 f un1 ÿ Mn1 un Qn1 pn Dt p_ n : Fpn1 f pn1 ÿ QTn1 u_ n Dt un ÿ Hn1 pn Dt Dp_ n ÿ Sn1 p_ n :
HL H0 pHf Hv Hs HDM Hw
Plastic modulus for unloading: ÿcu gu HU Hu0 ; Mg
R M R X NTu qs 1 ÿ n qw nSw Nu dX Q RX BT Sw mNw dX, where m 1; 1; 1; 0; H R X rNw T krNw dX S RX NTw Q1 Nw dX ~ BT Sw ÿ Sw0 mNw dX, Q RX ~ rNw T kw ÿ kw0 rNw dX H RX T 0 P RX B r dX R T s w C NTu t dC f u RX N ÿ u q 1T ÿ n q nSRw g dX p T T f X rNp kw qw b dX ÿ C Nw q n dC
Mass matrix Coupling matrix Permeability matrix Compressibility matrix Incremental coupling matrix Incremental permeability matrix Equivalent force vector External load vector Flow vector
Hf
where H0 , cDM , b0 , b1 , b1 and b2 are model parameters and gmax is the maximum stress ratio reached. Hw is related to partially saturated behaviour. Another possible expression for it is Hw 1 apc , where a is a material parameter.
and
1 nSw a ÿ n 2 Cs w Cs Sw 1 p ; Q Ks Kw nSw oSw where Cs n w : op
References
Z n dep ; s
[1] Mokni M. Relations entre deformations en masse et deformations localisees dans dans les materiaux granulaires. Ph.D. Thesis, Institut de Mecanique de Grenoble, France. 1992.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459 [2] Mokni M, Desrues J. Strain localisation measurements in undrained plane-strain biaxial tests on Hostun RF sand. Mech Cohes-Frict Mater 1998;4:419±41. [3] Lade PV, Bopp PA, Peters JF. Instability of dilating sand. Mech Mater 1993;16:249±64. [4] McManus KJ, Davis RO. Dilatation induced pore ¯uid cavitation in sands. Geotechnique 1997;47(1):173±7. [5] Vardoulakis I. Deformation of water-saturated sand: I. Uniform undrained deformation and shear band. Geotechnique 1996;46(3):441±56. [6] Vardoulakis I. Deformation of water-saturated sand: II. The eect of pore water ¯ow and shear banding. Geotechnique 1996;46(3):457±71. [7] Vardoulakis I, Sulem J. Bifurcation analysis in geomechanics. Blackie, Academic and Professional, 1995. [8] Shuttle DA, Smith IM. Localisation in the presence of excess pore-water pressure. Comput Geotech 1990;9: 87±99. [9] Loret B, Prevost JH. Dynamic strain localization in ¯uidsaturated porous media. J Engng Mech 1991;11:907±22. [10] Schre¯er BA, Majorana CE, Sanavia L. Shear band localization in saturated porous media. Arch Mech 1995;47:577±99. [11] Sanavia L, Zhang HW, Schre¯er BA. Strain localisation modelling in saturated sand samples. Forschungsbericht aus dem fachbereich Bauwesen 1995;66:357±66. [12] Schre¯er BA, Zhang HW, Pastor M, Zienkiewicz OC. Generalised plasticity constitutive model based strain localisation modelling in saturated sand samples. Computat Mech 1998;22:266±80. [13] Schre¯er BA, Sanavia L, Majorana CE. A multiphase medium model for localisation and postlocalisation simulation in geomaterials. Mech Cohes-Frict Mater 1996; 1:95±114. [14] Gawin D, Sanavia L, Schre¯er BA. Cavitation modelling in saturated geomaterials with application to dynamic strain localisation. Int J Numer Meth Fluids 1998;27: 109±25. [15] Bolzon G, Schre¯er BA, Zienkiewicz OC. Elastoplastic soil constitutive laws generalized to partially saturated states. Geotechnique 1996;46:279±89.
459
[16] Schre¯er BA, Sanavia L, Gawin D. Modelling of strain localisation in fully saturated soils. Proc 4th Int Worksh Localis Bifurcat Theory for Soils and Rocks, Gifu, Japan. 28 September±2 October 1997. Balkema, 1998. p. 305±12. [17] Alonso EE, Gens A, Hight DW. Special problem soils, General report. In: Hanrahan ID, Orr TLL, Widdis DF, editors. Proc 9th Eur Conf Soil Mech Fdn Engng, vol. III. Dublin. Balkema, Rotterdam, 1987. p. 1087±146. [18] Zienkiewicz OC, Chan A, Pastor M, Schre¯er BA, Shiomi T. Computational soil dynamics with special reference to earthquake engineering. Chichester: Wiley; 1999. [19] de Borst R, Sluys LJ, Muhlhaus HB, Pamin J. Fundamental issues in ®nite element analysis of localisation of deformation. Engng Comp 1993;10:99±121. [20] Zhang HW, Sanavia L, Schre¯er BA. An internal length scale in dynamic strain localisation of multiphase porous media. Mech Cohes-Frict Mater 1999;4:443±60. [21] Zhang HW, Schre¯er BA. Uniqueness and localisation analysis of elasto-plastic saturated porous media. Mech Cohes-Frict Mater, submitted for publication. [22] Lewis RW, Schre¯er BA. The ®nite element method in the static and dynamic deformation and consolidation of porous media. Chichester: Wiley; 1998. [23] Pastor M, Zienkiewicz OC, Chan AHC. Generalized plasticity and the modelling of soil behaviour. Int J Numer Anal Meth Geomech 1990;14:151±90. [24] Safai NM, Pider GF. Vertical and horizontal land deformation in a desaturating porous medium. Adv Water Resour 1979;2:19±25. [25] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 1. Averaging technique. Adv Water Resour 1979;2:131±144. [26] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 2. Mass, momenta, energy and entropy equations. Adv Water Resour 1979;2:191±201. [27] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 3. Constitutive theory for porous media ¯ow. Adv Water Resour 1980;3:25±40. [28] Gray WG, Hassanizadeh M. Unsaturated ¯ow theory including interfacial phenomena. Water Resour Res 1991; 27(8):1855±63.
www.elsevier.com/locate/compstruc
Numerical analysis of dynamic strain localisation in initially water saturated dense sand with a modi®ed generalised plasticity model H.W. Zhang a, L. Sanavia b, B.A. Schre¯er b,* a
State Key Laboratory of Structure Analysis for Industrial Equipment, Research Institute of Engineering Mechanics, Dalian University of Technology, Dalian 116024, People's Republic of China b Department of Structural and Transportation Engineering, University of Padua, via F. Marzolo 9, 35131 Padua, Italy Received 6 May 1999; accepted 1 May 2000
Abstract Shear band dominated process in fully and partially saturated sand samples is simulated by means of dynamic strain localisation analysis together with a multiphase material model. The partially saturated medium is viewed as a multiphase continuum consisting of a solid skeleton and pores ®lled by water and air (vapour) which, once it appears, is presumed to remain at the constant value of cavitation pressure (isothermal monospecies approach). The governing equations are based on the general framework of averaging theories. A modi®ed generalised plasticity constitutive model for partially saturated soils, developed from the general Pastor±Zienkiewicz sand model, has been implemented in a ®nite element code and used in the computational process. This model takes into account the eects of suction in the stiness of the porous medium (solid skeleton) in partially saturated state. A case of strain localisation, which has been tested in laboratory observing cavitation of the pore water, is studied. Negative water pressures, which are of importance in localisation phenomena of initially fully saturated undrained samples of dilatant geomaterials, are obtained similarly to those observed experimentally. Ó 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction Strain localisation is strain concentration in wellde®ned narrow zones, where the material behaviour is inelastic, while the remaining zones are elastic. This phenomenon is observed in a wide class of engineering materials, including single-phase solids like metals, and multiphase ¯uid saturated and partially saturated porous geomaterials like soils and concrete. The triggering mechanisms for the formation of shear bands are
* Corresponding author. Tel.: +39-49-8275611; fax: +39-498275604.
inhomogeneities in the material and stress concentrations, which may be, for multiphase geomaterials, due to the ¯uids±solid interaction. For such kinds of materials, the role of the ¯uids can be fundamental, as evidenced by extensive laboratory tests on dense and loose sands made by Mokni [1] and Mokni and Desrues [2]. In fact, in the case of undrained monotonic biaxial compression tests of initially water saturated dense Hostun sand under constant con®ning pressure, strain localisation was only observed when pore-water pressure of the order of ±80 kPa (value relative to the atmospheric pressure) developed, irrespective of the value of the applied back pressure [1,2]. (This drop in the water pressure is due to the dilatant behaviour of the solid skeleton associated with the impervious boundary of the sample which precludes the in¯ow or out¯ow of water.) The authors of the experiments
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 4 4 - 9
442
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Nomenclature (Variables with overbar refers to the nodal values) s solid phase g gaseous phase w liquid phase p generic phase ap acceleration of p-phase aps acceleration relative to the solid B strain operator Cs speci®c capacity CT tangential constitutive tensor eective diusivity tensor Dg fp ¯ow vector fu external load vector g gravity acceleration H permeability matrix I unit tensor krp ¯uid phase relative permeability k absolute or intrinsic permeability tensor Ks solid grain bulk modulus Kw liquid phase bulk modulus KT tangential stiness matrix molar mass of constituent p Mp M mass matrix n porosity n unit normal vector pp macroscopic pressure of the p-phase P equivalent force vector Q coupled matrix R universal constant gas Swo initial water saturation Sw water saturation Sg gas saturation S compressibility matrix t time variable t surface traction tensor
observed cavitation 1 of the pore-water with a production of vapour and stated that in this case, ``the
1 Cavitation of water may occur when the absolute value of water pressure is equal or lesser than the saturation water pressure at the considered temperature (neglecting the surface tension of the interface of the arising vapour bubbles). For thermodynamic equilibrium state, which is assumed in the model used in this paper at local level, the vapour pressure at the interface between saturated and partially saturated porous material (where the capillary pressure pc is zero) is equal, at T 20°C, to the saturation value of 2.339 kPa, (i.e. ±98.986 kPa, with reference to the atmospheric pressure).
tp T u vap
partial stress tensor temperature solid displacements velocity of the a-phase with respect to the pphase vp velocity of the p-phase a Biot's coecient c1 , c2 , H Newmark's parameters convective mass transfer coecient bc lw liquid dynamic viscosity q porous medium density qp intrinsic phase averaged density of the pphase r Cauchy stress tensor r0 eective Cauchy stress tensor I1 ®rst eective stress invariant J2 second invariant of the deviatoric eective stress J3 third invariant of the deviatoric eective stress ngL direction vector of plastic ¯ow for loading direction vector of plastic ¯ow for unloadngU ing g stress ratio Mg slope of the critical state line gmax maximum stress ratio reached gu stress ratio from which unloading takes place ag ,af ,Mf ,H0 ,cDM material parameters in Pastor± Zienkiewicz model b0 ,b1 ,Hu0 ,cu ,Kevo ,Keso material parameters in Pastor± Zienkiewicz model RI , b1 , b2 material parameters in modi®ed Pastor± Zienkiewicz model
onset of localisation is delayed until cavitation takes place in the pore-water'' and that ``non-drainage can preclude localisation as long as cavitation in the pore¯uid does not relax the isochoric constraint'' [2]. A similar experimental conclusion can be found in Ref. [3], for which ``a change in the draining conditions for a fully or nearly fully saturated soil can activate unstable behaviour and this may occur due to static or dynamic disturbances'', such as that done, for the latter case, in the numerical experiments showed in this paper. The fact that cavitation happens at the onset of localisation or shortly afterwards is still a controversial point. In this paper, we are not going to discuss this
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
open question, also because Mokni and Desrues observed cavitation at a pore pressure of ±80 kPa, which is dierent from the theoretical value of ±98.986 kPa (at 20°C) (or the corresponding value in the range of ambient temperatures). In this respect, it is worth to note that McManus and Davis [4] observed strain localisation and the associated cavitation during undrained triaxial tests of dense sand at the pore pressure value of approximately ±95 kPa and that Vardoulakis [5,6] revealed cavitation in globally undrained biaxial tests of dense St. Peter sandstone at the pore pressure of ±51 kPa Ref. [5,6]. The important role that ¯uids can have in localisation modelling is con®rmed also in Ref. [7], where the bifurcation analysis of undrained medium dense St. Peter sandstone is performed. The authors found the existence of bifurcation modes only in the form of locally drained shear band, which means the necessity to include the water ¯ow in a strain localisation modelling (see also Refs. [5,6]). In the present paper, we start from the abovementioned important theoretical and experimental role of the ¯uids in strain localisation analysis to justify the use of a multiphase model, where, air, dissolved air in water and water vapour, should be taken into account together with water and the solid skeleton. Starting from the ®rst numerical analysis of fully saturated porous media made by Shuttle and Smith [8] and Loret and Prevost [9] in statics and in dynamics respectively, several authors have worked on this subject worldwide. As far as some original aspects of our research group are concerned, a numerical model based on the general frame of averaging theories and capable of simulating also the shear band dominated processes in saturated porous media was presented in Refs. [10±12]. The model was then extended in Refs. [13,14] to deal with the cavitation phenomenon of an initially water saturated and globally undrained medium, where the so-de®ned monospecies approach and the isothermal two phase ¯ow model have been respectively proposed. When simulating the behaviour of dense sand, a Mohr±Coulomb constitutive relationship with a softening branch was used in Refs. [10,13] together with associative plasticity and in Refs. [11,12] with non-associative plasticity and Pastor± Zienkiewicz model for sand. The main focus of this paper lies in the use of a modi®ed generalised plasticity model developed in Ref. [15] in the context of the numerical simulation of strain localisation problems of initially water saturated dense (dilatant) sand. This constitutive model for the solid skeleton takes into account capillary eects which arise once cavitation takes place. This paper is organised as follows: ®rst, in Section 2, a brief summary of the adopted multiphase material model for geomaterials where phase change (cavitation) can take place is presented. To this end, the monospecies approach [13,16] used for the computations is recalled, because it is a
443
simple model with a good performance. The main features of the modi®ed generalised plasticity model developed by Bolzon et al. [15] are then brie¯y described in Section 3, where capillary pressure is taken as independent variable, beyond the eective stress, because this is a necessary aspect of partially saturated soil models [17]. This soil model has been implemented in the ®nite element code SWANDYNE [18], as described in Section 4. Numerical examples of shear band localisation within dense sand samples are shown in Section 5. In particular, the strong in¯uence that ¯uid phases can have in localisation phenomena is investigated from a numerical point of view by a comparison between the results of the monospecies approach and the water saturated model performed on a classical biaxial test. Mohr±Coulomb law for the solid skeleton was used for the sake of simplicity. Then, the results of a biaxial test under symmetric and unsymmetric loading cases taken from Refs. [1,2] are shown, using the modi®ed Pastor±Zienkiewicz model. These results show the superiority of the adopted constitutive model in reproducing the experimental pore pressure values during shear banding, because they are much closer to the observations by Mokni [1], Mokni and Desrues [2], McManus and Davis [4] and Vardoulakis [5,6] than the previous attempts using a Mohr±Coulomb constitutive relationship [10,13]. It is also worth to outline that negative pore pressures were not observed in Ref. [17] when using the basic generalised plasticity model for saturated sands, from which the model used here is derived, because loose (contractive) sands were studied. Finally, we recall the well-known fact that in localisation analysis of a single phase material with softening, one obtains an ill-posed problem. This implies a strong dependence of the shear band width on the spatial discretisation used once the numerical analysis is performed. To overcome this unphysical behaviour, different regularisation strategies can be found in the literature (see e.g. Ref. [19] for a review). On the contrary, multiphase material contains a natural regularisation mechanism due to the presence of a Laplacian term in the mass balance equation of the ¯uids when DarcyÕs law is used. The numerical properties of this natural regularisation related to axial waves have been investigated in Ref. [13], while the internal length scale contained in the model for the 1-D and 2-D cases and for axial waves is presented in Refs. [20,21] respectively. Hence, this aspect will not be further pursued in this paper. 2. Simpli®ed model for simulation of cavitation: isothermal monospecies approach The mathematical model for a multiphase porous medium is now recalled. It can be derived from the
444
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
general thermo±hydro-mechanical model developed in Ref. [22] by introducing the assumptions outlined next. A short summary of the general model is presented in Appendix A for interested readers. Small strains are assumed for the development of the following equations. We consider here globally undrained situations, where cavitation has been experimentally observed [1,2,4±6]. In such conditions, neglecting air dissolved in water and air bubbles in water, which may be present in fully saturated specimens at low pressures, only two ¯uid phases are present after desaturation: liquid water and water vapour. These assumptions can be made because the authors of the experimental tests removed air bubbles as far as possible and saturated the samples using de-aired water [2,4]. This means that the air mass balance equation (A.15) in not needed. Moreover, isothermal conditions are assumed for the sake of simplicity. An isothermal approach means that physically heat is supplied or taken away with in®nite intensity to maintain constant temperature (this implies in®nite heat capacity value), thus there are no energetic restrictions for phase change (evaporation or condensation). Our isothermal model consists ®nally of the linear momentum balance equation for the multiphase medium (A.14) and the mass balance equation for water species (A.16) with the diusion term omitted (there is no diusion of water vapour in the air). Also the Clausius±Clapeyron equation is not needed in KelvinÕs law (A.7), because the water saturation pressure pgws is constant in isothermal conditions. The monospecies approach is now derived [13] as follows. Neglecting further the vapour density qgw with respect to the water density qw , from the state equation of vapour, Eq. (A.5), it follows that also the vapour pressure pgw is negligible. Hence, the average pressure (A.10) is given by p Sw pw ;
1
Further, the dynamic seepage forcing term in Eq. (A.13), connected with the solid acceleration, is very small when compared with other terms of the equation system, so it could be neglected [13]. Thus, the ¯uid mass balance equation (A.16) has now the following form: ÿ
o nSw qw Sw qw div vs ot kk rw div qw w grad pw ÿ qw g 0; l
3
where qw is the water density, vs , the solid velocity vector, k, the intrinsic permeability tensor, krw , the water relative permeability, and lw , the dynamic viscosity. Together with Eqs. (2) and (3), we need the experimentally determined constitutive equations of the medium: Sw Sw pc ;
k rw k rw Sw
4
with pc being the capillary pressure. Finally, Eq. (A.8) reads pc ÿpw ;
5
which states that negative water pressure corresponds to the development of capillary pressure. The relationship between capillary pressure and water degree of saturation is of the type depicted in Fig. 1. In fact, because of the conditions of the experiments and the simplifying assumptions, partial saturation is reached only at cavitation, which occurs at minus the atmospheric pressure using pressures relative to the atmospheric pressure, which in turn is set to zero. For the model closure, the appropriate initial and boundary conditions (A.17b,d), (A.18b,c) and (A.20) are also needed.
where Sw is the water degree of saturation and pw is the water pressure. By ignoring the relative component of the ¯uid accelerations, the primary variables may be reduced to solid displacements u and water pressure pw . After application of the eective stress equation (A.9) to the linear momentum balance equation (A.14), we obtain its simpli®ed form as div r0 ÿ a Sw pw I q g ÿ as 0;
2
where r0 is the Bishop eective stress tensor, a 1 ÿ Kd =Ks , the BiotÕs constant, with Kd and Ks represently the bulk moduli of the medium and of the solid material respectively, q, the mixture density, g and as , the gravity acceleration and the solid acceleration vector, respectively. I is the second-order identity tensor.
Fig. 1. Saturation-water pressure relationship used for cavitation modelling.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
3. The modi®ed Pastor±Zienkiewicz generalised plasticity sand model The generalised plasticity soil model was introduced by Pastor et al. [23], which was extensively validated for fully saturated sand behaviour. It has been shown by comparison with experiments that this soil model is capable of a good performance for a wide range of conditions (see Ref. [18] for a review). The interested reader can ®nd the basic features of the Pastor±Zienkiewicz model in Appendix B. The modi®ed model will be now described following Ref. [15], who extended the Pastor± Zienkiewicz model to partially saturated sands. As mentioned in Ref. [15], both BishopÕs stress and capillary pressure (or other two combination of stresses and pressures [17]) are necessary to be introduced as independent stress parameters to describe the behaviour of the partially saturated soil. In the modi®ed model, eects of capillary pressure are taken into account by the soil stiness changes induced by suction and irreversible compressive volumetric strains ev of soil upon wetting. (
Wun1 Mn1 Dun Pn1 ÿ Qn1 H Dt Dp_ n ÿ Fun1 0 Wpn1 QTn1 c1 Dt Dun Hn1 H Dt Dp_ n Sn1 Dp_ n ÿ Fpn1 0;
The initial yield limit of the fully saturated case of the Pastor±Zienkiewicz model is hence modi®ed in partially saturated conditions to introduce its dependence on suction as [15] 0 pyi0 pc py0i RI pc ;
6
0 where py0i represents the initial yield limit in saturated conditions of the mean eective stress. RI is a material parameter which expresses an increasing function of suction when water saturation is less than one and is determined by interpolation of experimental data. Because irreversible volumetric strain is assumed as the parameter which controls the volumetric hardening, the evolution of the yield surface in the (p0 , pc ) plane is given as !1apc p0 y0 0 c 0 c py p py0i RI p 7 0 py0i
445
main by means of the NewmarkÕs scheme. The unknown ®eld variables are expressed in the whole domain by global shape function matrices, Nu and Nw , and the nodal value vectors, u and p: u Nu u;
pw Nw p:
8
Once the coupling matrix Q, the mass matrix M, the strain operator B and the external load f u , are introduced (see Appendix C), the linear momentum balance can be written as Z BT r0 dX ÿ Qp Mu f u ; 9a X
while mass balance equation becomes Hp QT u_ Sp_ f p ;
9b
where H is the permeability matrix, S, the compressibility matrix and, f p , the ¯ow vector (Appendix C). After the introduction of the NewmarkÕs scheme for time integration, equations system (9) can be written at the current time tn1 as
10
where c1 , H and c2 , which appear in Eq. (11), are the NewmarkÕs parameters, and the vectors Fun1 and Fpn1 are de®ned in the Appendix C. pn1 must be evaluated by integration of the constitutive law. The non-linear coupled equation system is linearised in a standard way thus yielding the linear algebraic equations system to be solved numerically: ! M 12 KT c2 Dt2 ÿQH Dt Du ÿH ÿQT H Dt Dp_ c1 HH Dt S u ÿW ÿ H Wp : 11 c1
0 where py0 is the initial mean stress and a is a material parameter.
Since the Newton±Raphson method requires the Jacobian matrix to be evaluated and inverted at each iteration, other modi®ed schemes are also used to achieve convergence with less computational eort. In particular, the use of secant updates, like DavidonÕs and Broyden±Fletcher±Goldfarb±ShannoÕs (BFGS) methods is found advantageous in non-linear analyses [18].
4. Finite element discretisation and solution
4.1. Initial static solution for modi®ed Pastor±Zienkiewicz model
The governing Eqs. (2) and (3) are discretised in space using a GalerkinÕs procedure and in the time do-
The modi®ed Pastor±Zienkiewicz model needs the initial value of the elastic parameters depending on the
446
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
initial mean pressure. The form of the governing equations used to obtain the initial value for the modi®ed constitutive model are the following static drained equations: Z BT r00 dX ÿ Q0 p0 f u0 ; 12a X
H0 p0 f p0 ;
12b
where Q0 and H0 are matrixes depending on p0 and Sw0 . Numerical experiments have shown that inconsistency between Eqs. (9a), (9b) and (12a), (12b) will generate an unstable dynamic computational process with Eqs. (9a) and (9b), even if the dynamic load is zero. The reason is that the material parameters for static analysis in Eqs. (12a) and (12b) are dierent from those in Eqs. (9a) and (9b), the parameters in the latter ones being totally regenerated with the de®nitions in the modi®ed Pastor± Zienkiewicz model. To avoid this problem, we move the contribution of the initial state to the right-hand side of the governing equations and rewrite the dynamic equations (9a) and (9b) in the following incremental forms: Z ~ ; BT Dr0 dX ÿ Q Dp M Du Df u Qp 13a 0 X
~ 0; H Dp QT Du_ S Dp_ Df p ÿ Hp
13b
in which the incremental variables are de®ned as Dx x ÿ x0 , Q, H and S are the matrixes dependent on the current water pressure and saturation. The expres~ H ~ can be found in the Appendix C. sions of Q,
5. Numerical results Before we deal with the simulation of a laboratory test taken from Refs. [1,2], the in¯uence of the ¯uid phases in localisation phenomena is investigated from a numerical point of view. At this end, a rectangular mass of water saturated dense sand under compression is analysed by means of the isothermal monospecies approach (IMA) and the water saturated model (WSM). This example has been studied by several authors [9± 14,16]. The sample has impervious boundaries (Fig. 2) and hence the problem is globally undrained. Vertical and horizontal displacements are constrained at the bottom surface. Ramp loading is applied at the top. During computations, the self weight of the soil and water ®lling its pores are taken into account. Hydrostatic distribution of water pressure is assumed as an initial condition. The sample is subjected to plane strain condition. The relationships between capillary pressure, saturation and water relative permeability proposed by Safai and Pinder [24] for San Fernando sand are used when cavitation develops (Fig. 1). A Mohr±Coulomb yield criterion with associative ¯ow rule and isotropic linear softening is used for the solid skeleton. The material parameters used during the computations are Young's modulus E 285 MPa, Poisson's ratio m 0:4285, solid grain density qs 2000 kg mÿ3 , water density qw 998:2 kg mÿ3 , apparent cohesion c0 1:84 MPa, hardening modulus H ÿ40 MPa, angle of internal friction u 20°, initial porosity n 0:20, solid grain bulk
Fig. 2. Description of the geometrical characteristics and of the load function of the ®rst example.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
447
Fig. 3. Contours of equivalent plastic strain at dierent times with the monospecies approach model: (a) 0.245 s, (b) 0.18 s and (c) 0.22 s.
modulus Ks 6:78 GPa, water bulk modulus Kw 0:20 GPa, permeability kw 0:25 10ÿ3 m sÿ1 , BiotÕs constant a 1:0. A mesh of 396 (18 elements in width and 22 over the height) four-node isoparametric ®nite elements of equal size, the same for displacements and water pressures, has been used. Computations are performed by use of the modi®ed version [13] of the SWANDYNE code [18]. The two models present signi®cant dierences. Since the saturated model neglects the partial saturation, the generalised stiness of the multiphase medium and the coupling between solid skeleton and ¯uids, i.e. the structural ¯uids±solid interaction, are not properly described because the matrices Q, H and S of Eq. (10) are constant (see Appendix C). As a consequence, the evolution of the process and all the state variables change considerably. In fact, the numerical results reveal a big delay in the formation of the localised bands when the saturated model is used (from 0.245 to 0.40 s, see Figs. 3(a) and 4 for the IMA and the WSM, respectively), even if the shear bands start to appear at the same time
(@0.18 s) because of the same velocity of the ®rst plastic wave. In this case, the localised zones start from the corners of the (rigid) bottom surface when the stress induced by the re¯ected loading wave reach the yield
Fig. 4. Contours of equivalent plastic strain at 0.40 s with the water saturated model.
448
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
value of the material model, and propagate versus the top surface (Fig. 3(b) and (c); Fig. 3(a)±(c) are obtained with the monospecies approach model). Then also a larger value of the equivalent plastic strain (0.0448 for the WSM, 0.0441 for the IMA) is observed, while the shape of the shear bands is the same. As far as the behaviour of the ¯uid phase is concerned, an increase of the zone which should become partially saturated can be observed for the WSM (see Fig. 6 in comparison with Fig. 5), with a decreasing of the maximum water pressure (1.62 MPa instead of 2.81 MPa) and an increase of the minimum water pressure (ÿ11.1 MPa instead of ÿ26.9 MPa). Moreover, the ¯ow rate of the water is double in the saturated case, because permeability does not decrease with the capillary pressure. From this comparison between the two models, it can be concluded that ¯uid phases are important for localisation phenomena of initially saturated multiphase dense geomaterials under globally undrained conditions. This conclusion is indirectly con®rmed by the experimental results [2], for which the presence of the pore
¯uid in globally undrained conditions delays the shear band formation. The examples presented next refer to the biaxial laboratory test of Refs. [1,2] and make use of the modi®ed Pastor±Zienkiewicz model for sand to take into account the partially saturated conditions appearing during localisation at cavitation pressure. With this model, also the constitutive tangent matrix KT is dependent from the capillary pressure. The ¯uid±solid interaction is hence described also at constitutive level. The geometrical characteristics with impervious boundary conditions, which are exactly the same as in Refs. [1,2], and the load function, are shown in Fig. 7. Vertical and horizontal displacements are constrained at the bottom surface. Ramp loading is applied at the top of the sample. The permeability in saturated condition is 0:25 10ÿ3 m sÿ1 , while, once cavitation starts, the relationship between the degree of saturation and the relative permeability for S. Fernando sands [24] has been assumed. Also the capillary pressure-saturation relationship is that given in Safai and Pinder [24] for S. Fernando sands, as shown in Fig. 1 [16]. These relationships have been used here for sake of simplicity because of the lack of experimental results of Hostun sands. Both a symmetric and an unsymmetric load case are chosen here to simulate the shear band formation. The material parameters for the dense sand are water bulk modulus Kw 0:20 GPa, solid grain bulk modulus Ks 6:78 GPa, initial porosity n 0:20, solid density qs 2000 kg mÿ3 , water density qw 998:2 kg mÿ3 . The initial conditions according to Refs. [1,2] are rxx ryy ÿ100 kPa, rzz ÿ130 kPa, rxy rxz ryz 0
Fig. 5. Contours of negative water pressure at 0.245 s with the monospecies approach model.
Fig. 6. Contours of negative water pressure at 0.400 s with the water saturated model.
Fig. 7. Description of the geometrical characteristics and of the load function of the second example.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
449
Fig. 8. Contours of equivalent plastic strain at dierent times for the symmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 4.0 ´ 10ÿ3 s.
Fig. 9. Contours of equivalent plastic strain at dierent times for the symmetric load case and initial water pressure p0 200 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 2.5 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
and p0 200 or 300 kPa. The material parameters for the modi®ed Pastor±Zienkiewicz for sands are as following: Mf 1:0, Mg 1:0, af 0:45, ag 0:45, Kevo 1050 Pa, Keso 2000 Pa, b0 4:2, b1 0:2, Hu0 600 Pa, H0 350 Pa, cu 2:0, cDM 2:0, b1 0:77, b2
16:4, RI 1:0 [11]. The load values for symmetric and unsymmetric load cases are F10 153 kN, F20 5:67 kN. To study the problem in detail, the following three cases are simulated:
450
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 10. Contours of equivalent plastic strain at dierent times for the unsymmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
1. Symmetric load case with initial water pressure, p0 300 kPa; 2. Symmetric load case with initial water pressure, p0 200 kPa;
3. Unsymmetric load case with initial water pressure, p0 300 kPa. As in the laboratory tests, in the numerical simulation, the shear band starts from the top surface of the specimen, near loading area, for both symmetric and unsymmetric load cases, due to the stress concentration caused by the lower stiness of the soil in that part with respect to the remaining part of the sample. The development of shear band for the symmetric load and unsymmetric load cases are shown by the contours of equivalent plastic strain variable depicted in Figs. 8±10 for three dierent times, for cases 1±3, respectively. In the unsymmetric load case (Fig. 10), the shear band propagates in a direction controlled by the disturbing load. The pattern of this shear band is similar to what was observed in the laboratory tests [1] before the failure of the tested sample. For the symmetric load case, two shear bands appear in the upper region of the specimen, which are symmetric because of the perfect symmetry of the solved problem. The eective plastic strain along a vertical (central) cross-section and for dierent times is shown in Figs. 11±13 for the three investigated cases, respectively. These ®gures give an idea of the values involved in Figs. 8±10. The pore-water pressure distributions at dierent times are presented in Figs. 14±16 for cases 1±3, respectively. Because the uniform water pressure is applied on the lateral boundary of the sample, the results of water pressure are not regular. However, water pressure concentration with negative values of water pressure (traction) (Figs. 17±19) can be observed at the centre of the shear bands due to the plastic dilatant behaviour of the dense sand. The water pressures distribution along a vertical cross-section through the centre line for dierent time values are shown in Figs. 17±19 for cases 1±3, respectively. These help to understand their distributions in the three previous ®gures. In particular, a traction value close to ÿ100 kPa is obtained in both the symmetric and unsymmetric load cases when shear bands develop. These values of the water traction are quite similar to the vapour saturation pressure of ÿ98.986 kPa and the negative pressure observed in Refs. [4±6] at cavitation. In the numerical results, they are slightly larger because of the development of higher capillary pressure. Moreover, from Figs. 17 and 18, it can be seen that the same negative pressure develops, irrespective of the value of the applied back pressure. This is also in agreement with the experimental ®ndings of Refs. [1,2,4]. In fact, Mokni and Desrues found that strain localisation developed during biaxial tests when cavitation developed, irrespective of the value of the applied back pressure. Similarly, also McManus and Davis [4] reported the same (cavitation) ®nal value (of ±95 kPa) during their triaxial tests, starting from dierent values of the initial water pressure. Hence, with the constitutive model used, which takes into account the capillary ef-
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
451
Fig. 11. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 300 kPa.
Fig. 12. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 200 kPa.
fects, the water tractions do not exceed values which are close to the experimentally observed ones. The situation with a Mohr±Coulomb model [10,11,13±16], (see also Fig. 6) without capillary eects was quite dierent, because there, much higher capillary pressures were obtained. This points to the need to choose an appropriate constitutive model for simulating localisation and postlocalisation behaviour. Moreover, the ®gures for eective plastic strain (Figs. 11±13) with the corresponding ®gures for pressure evolution (Figs. 17±19) are useful also to analyse the shear banding interaction with the pressure evolution. In fact,
it can be clearly seen that negative pore pressures arise at the same instance when localisation appears in the strain distribution. Finally, Fig. 20 shows the deformation of the sample when the shear band occurs for the three cases studied. 6. Conclusions Numerical analysis of dynamic strain localisation in initially water saturated dense sand samples using a multiphase material model has been presented. The
452
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 13. Development of equivalent plastic strain along the vertical direction in the central cross-section of the sample for unsymmetric load case and initial water pressure p0 300 kPa.
model is capable of simulating the partially saturated behaviour of globally undrained dense sand which develops because of experimentally observed cavitation. The partially saturated medium is viewed as a two-phase continuum consisting of a solid skeleton and pores ®lled with water or water and vapour at cavitation pressure. The formulation is based on the general frame of averaging theories. For partially saturated sands, a modi®ed generalised plasticity constitutive model based on Pastor±Zienkiewicz model was adopted in the computational process. This is currently one of the simplest and yet one of the most eective models for describing the full range of saturated and partially saturated sand behaviour. Using this constitutive model, the numerical simulation of multiphase geomaterials is able to describe the ¯uid±solid interaction both at the structural and constitutive levels. A dense sand sample similar to a laboratory test was studied in detail with ramp loading to demonstrate the capability of the numerical model presented here to handle such situations. An agreement between the experimental observations and numerical results was obtained as far as the form of the shear band (unsymmetric case) and the pore pressures are concerned. These pressures are the most intriguing aspect of the experiment.
Acknowledgements This work was supported by the EC International Scienti®c Cooperation Programme and the National Natural Science Foundation of China (19872016) and by the research funds MURST 020902019 from the
Italian Ministry of Scienti®c and Technological Research (Co®nanziamento Murst 1998). Appendix A A.1. Mathematical model for the description of a multiphase porous media The partially saturated porous medium is a material with an internal microstructure and is treated here as a multiphase material where the voids of the solid skeleton are ®lled with water and gas which may be either vapour alone or a mixture of dry air and vapour. The full mathematical model necessary to simulate thermo± hydro-mechanical transient behaviour of fully and partially saturated porous media was developed in Ref. [22] using averaging theories following Hassanizadeh and Gray [25±27] and Gray and Hassanizadeh [28]. It assumes non-polar and immiscible constituents (except for the vapour and the gas) and local thermal equilibrium state. The most important equations are brie¯y summarised in what follows, assuming small strains. The state of the medium will be described by water pressure pw , capillary pressure pc , temperature T and displacement vector of the solid matrix u. The material free macroscopic balance equations for the generic p-phase are the following: Mass balance equation: oqp qp div vp qp ep : ot Linear momentum balance equation:
A:1
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
453
Fig. 14. Contours of water pressure at dierent times for the symmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 4.0 ´ 10ÿ3 s.
Fig. 15. Contours of water pressure at dierent times for the symmetric load case and initial water pressure p0 200 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 2.5 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
div tp qp gp ÿ ap qp ep ^tp 0;
exchange between constituent p and the other constituents, tp , the stress tensor (called partial stress tensor in mixture theories), g, the external momentum supply related to gravitational eects, qp ap , the volume density of
A:2
where qp is the phase average density, vp , the mass average velocity of the p-phase; qp ep represents the mass
454
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
implies symmetry of the stress tensor. The energy balance equation for the p-phase may be written as qp
oEp tp : dp qp hp ÿ div qp qp Rp ; ot
A:3
where qp Rp represents the exchange of energy between the p-phase and other phases of the medium due to phase change and mechanical interaction, qp is the internal heat ¯ux, hp results from the heat sources and dp is the rate of the strain tensor. Ep accounts for the speci®c energy of the volume element and for the kinetic energy. The entropy inequality for the mixture is X okp qp q hp A:4 qp qp ep kp div p ÿ p p P 0; ot T T p where kp is the speci®c entropy of the constituent p and, ep kp , the entropy supply due to mass exchange. The moist air (g) in the pore system is assumed to be a perfect mixture of two ideal gases, dry air (ga) and water vapour (gw). The equation of state of perfect gas (ClapeyronÕs equation) is hence valid: pga qga RT =Ma ;
pgw qgw RT =Mw ;
A:5
where T is the absolute temperature, Mp , the molar mass of constituent p and, R, the universal gas constant. Further, DaltonÕs law applies, giving the gas pressure and density as pg pga pgw ;
qg qga qgw :
A:6
Due to the curvature of the meniscus separating water from gas inside the pores of the medium considered as a capillary porous body, the equilibrium vapour pressure pgw can be obtained from the Kelvin relationship: pc Mw pgw pgws exp ÿ w ; A:7 q RT where pgws is the vapour saturation pressure which can be obtained from the Clausius±Clapeyron equation [16], and pc is the capillary pressure which is de®ned from the sorption equilibrium equation as pc pg ÿ pw :
Fig. 16. Contours of water pressure at dierent times for the unsymmetric load case and initial water pressure p0 300 kPa: (a) 2.0 ´ 10ÿ3 s, (b) 3.0 ´ 10ÿ3 s and (c) 3.5 ´ 10ÿ3 s.
the inertia force and ep ^tp accounts for momentum exchange due to mass supply and mechanical interactions with other phases. As for the angular momentum balance equation, the assumption of microscopically non-polar constituents
A:8
pp (p g, w) is related to tp by [28] tp ÿgp pp I, with gp as the p-phase volume fraction, and I is the second-order identity tensor. The constitutive law for the solid phase in the partially saturated state is introduced through the concept of modi®ed eective stress r0 (Bishop stress): r0 r apI; s
A:9 w
g
where r t t t is the total stress tensor, a, the BiotÕs constant and, p, an average pressure in the solid phase which, by neglecting the dependence of Helmholtz free energy on void fraction [27,28] is given as
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
455
Fig. 17. Development of water pressure along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 300 kPa.
p Sw pw 1 ÿ Sw pg :
A:10
Sw is the saturation of liquid water, which is an experimentally determined function of capillary pressure and temperature. The constitutive relationship for the solid skeleton has the form dr0 CT de ÿ deT ÿ de0 ;
A:11
where CT is a tangent constitutive tensor, deT I b3s dT is the strain increment caused by thermoelastic expansion, bs is the cubic thermal expansion coecient of the solid and de0 denotes the autogenous strain increments and the irreversible part of the thermal strains tensor. A particular form of this relation will be discussed in Appendix B. Finally for a binary gas mixture, FickÕs law gives the following relative average velocities vgp of the diusing
Fig. 18. Development of water pressure along the vertical direction in the central cross-section of the sample for symmetric load case and initial water pressure p0 200 kPa.
456
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
Fig. 19. Development of water pressure along the vertical direction in the central cross-section of the sample for unsymmetric load case and initial water pressure p0 300 kPa.
species, due to the existing gradients of their concentrations: ga Ma Mw p vga ÿ D grad g 2 M pg gw Ma Mw p Dg grad A:12 ÿvgw ; M2 pg where Dg is the eective diusivity tensor and M is the molar mass of the gas mixture. A.2. Governing equations Using the Macroscopic balance equations and the constitutive equations of Appendix A.1, we obtain the following governing equations, as in Lewis and Schre¯er [22]. The linear momentum balance equation for ¯uids yields the generalised Darcy equation: nSp vps
kk rp ÿ grad pp qp g ÿ as ÿ aps ; lp
A:13
where p ga, gw, w and g, vps is the (intrinsic massaveraged) velocity relative to the solid, aps , the acceleration relative to the solid, n, the porosity of the medium, k, the intrinsic permeability tensor, lp , the dynamic viscosity and krp , the relative permeability which is a function of the degree of saturation and of the temperature. The linear momentum balance equation for the multiphase medium is divr q g ÿ as ÿ nSw qw aws ÿ nSg qg ags 0;
A:14
where q is the averaged P density of the multiphase medium q 1 ÿ nqs p6s nSp qp . The dry air mass balance equation, after introduction of DarcyÕs law and FickÕs law and neglecting acceleration terms is o kk rg nSg qga Sg qga div vs ÿ div qga g grad pg l ot gw Ma Mw p div qg 0: A:15 Dg grad M2 pg The vapour mass balance equation and the water mass balance equation are summed to obtain the mass balance equation for all water species. DarcyÕs law and FickÕs law are introduced as above, yielding rg o gw gw s gw kk g nSg q Sg q div v ÿ div q grad p ot lg gw Ma Mw p ÿ div qg Dg grad M2 pg o ÿ nSw qw ÿ Sw qw div vs ot rw w kk w w s ws div q grad p ÿ q g ÿ a ÿ a : lw A:16 The macroscopic mass balance equation for the solid has already been summed with the above mass balance equations to eliminate the time derivative of the porosity n. The energy conservation equation (enthalpy balance) is not used in the simpli®ed isothermal monospecies approach model and is hence omitted.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
457
The boundary conditions can be imposed values on Cp or ¯uxes on Cqp , where the boundary is C Cp [ Cqp . The imposed values on the boundary for gas pressure, water pressure and displacements are as follows: pg p^g on Cg ; pw p^w on Cw ; u^ u on Cu for t P t0
A:18a±c
The volume average ¯ux boundary conditions for water species and dry air conservation equations to be imposed at the interface between the porous medium and the surrounding ¯uid are the following: qga vg ÿ qg vgw n qga on Cqg ; A:19a g
qgw vg qw vw qg vgw n g ÿ gw gw bc q ÿ qgw qw 1 q
on Cqc ;
A:19b
where n is the unit vector, perpendicular to the surface of the porous medium, pointing toward the surrounding gas, qgw 1 is the mass concentration of water vapour in the undisturbed gas phase distant from the interface, and bc is the convective mass transfer coecient, while qga , qgw and qw are the imposed dry air ¯ux, imposed vapour ¯ux and imposed liquid ¯ux, respectively. Eq. (A.19) are the natural boundary conditions, respectively, for the dry air conservation Eq. (A.19a) and water species conservation Eq. (A.19b), when the solution of these equations is obtained through a weak formulation of the problem, as is usually done with the ®nite element method. The traction boundary conditions for the displacement ®eld are rn t
on Cqu for t P t0 ;
A:20
where t is the imposed traction. Appendix B Fig. 20. Deformation of the sample when shear band occurs: (a) Symmetric load case, initial water pressure p0 300 kPa at time 4:0 10ÿ3 s. (b) Symmetric load case, initial water pressure p0 200 kPa at time 3:5 10ÿ3 s. (c) Unsymmetric load case, initial water pressure p0 300 kPa at time 3:5 10ÿ3 s.
A.3. Initial and boundary conditions The initial conditions specify the full ®elds of gas pressure, water (or capillary) pressure, temperature and displacements: pg p0g ; u u0
pw p0w ; T T0 ; at t t0 :
A:17a±d
B.1. Pastor±Zienkiewicz generalised plasticity model for sand The approach developed by Pastor et al. [18,23] and Bolzon et al. [15] for saturated and partially saturated soils, respectively, explicitly de®nes direction vectors in terms of the following three invariants p0 , q and h p I1 p0 ÿ ; q 3J2 ; 3 ! p 1 ÿ1 3 3 J3 h ÿ sin ; B:1 2 J23=2 3 where I1 is the ®rst eective stress invariant, J2 and J3 are the second and the third invariant of the deviatoric eective stress, respectively. The de®nitions and
458
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459
parameters that need to be determined in the modi®ed Pastor±Zienkiewicz model can be summarised as follows. Loading direction vector: n np ; nq ; nh T
B:2
p
q
h
with n 1 af Mf ÿ g, n 1, n ÿ 1=2 qMf cos 3h, where af , Mf are model parameters and g is the stress ratio de®ned as q=p. Plastic ¯ow direction vector for loading: ngL npgL ; nqgL ; nhgL T
B:3
with npgL 1 ag Mg ÿ g, nqgL 1, nhgL ÿ 1=2qMg cos 3h, in which ag is a material parameter and Mg is the slope of the critical state line, related to the residual 6 sin / friction angle by Mg 3ÿ sin . / sin 3h Plastic ¯ow direction vector for unloading: ngU npgU ; nqgU ; nhgU T
B:4
with npgU ÿjnpgL j, nqgU nqgL , nhgU nhgL . This is chosen in order to ensure that densi®cation occurs in unloading. Plastic modulus for loading: B:5
with 1ÿ
g gf
4 ;
gf
1 1 af
Mf ;
g ; Mg ÿcDM g HDM ; gmax
Hv 1 ÿ
Hs b0 b1 exp ÿb0 n;
Hw 1 b1 exp ÿb2 pc ÿ 1ÿ1
B:6
where Hu0 and cu are material parameters and gu is the stress ratio from which unloading takes place. Elastic moduli: The bulk and shear moduli, respectively, are K
Kevo 0 p; p00
G
Keso 0 p; 3p00
B:7
where both Kevo and Keso are model parameters and p00 is the initial value of p0 . Appendix C
0;
0 T
Fun1 f un1 ÿ Mn1 un Qn1 pn Dt p_ n : Fpn1 f pn1 ÿ QTn1 u_ n Dt un ÿ Hn1 pn Dt Dp_ n ÿ Sn1 p_ n :
HL H0 pHf Hv Hs HDM Hw
Plastic modulus for unloading: ÿcu gu HU Hu0 ; Mg
R M R X NTu qs 1 ÿ n qw nSw Nu dX Q RX BT Sw mNw dX, where m 1; 1; 1; 0; H R X rNw T krNw dX S RX NTw Q1 Nw dX ~ BT Sw ÿ Sw0 mNw dX, Q RX ~ rNw T kw ÿ kw0 rNw dX H RX T 0 P RX B r dX R T s w C NTu t dC f u RX N ÿ u q 1T ÿ n q nSRw g dX p T T f X rNp kw qw b dX ÿ C Nw q n dC
Mass matrix Coupling matrix Permeability matrix Compressibility matrix Incremental coupling matrix Incremental permeability matrix Equivalent force vector External load vector Flow vector
Hf
where H0 , cDM , b0 , b1 , b1 and b2 are model parameters and gmax is the maximum stress ratio reached. Hw is related to partially saturated behaviour. Another possible expression for it is Hw 1 apc , where a is a material parameter.
and
1 nSw a ÿ n 2 Cs w Cs Sw 1 p ; Q Ks Kw nSw oSw where Cs n w : op
References
Z n dep ; s
[1] Mokni M. Relations entre deformations en masse et deformations localisees dans dans les materiaux granulaires. Ph.D. Thesis, Institut de Mecanique de Grenoble, France. 1992.
H.W. Zhang et al. / Computers and Structures 79 (2001) 441±459 [2] Mokni M, Desrues J. Strain localisation measurements in undrained plane-strain biaxial tests on Hostun RF sand. Mech Cohes-Frict Mater 1998;4:419±41. [3] Lade PV, Bopp PA, Peters JF. Instability of dilating sand. Mech Mater 1993;16:249±64. [4] McManus KJ, Davis RO. Dilatation induced pore ¯uid cavitation in sands. Geotechnique 1997;47(1):173±7. [5] Vardoulakis I. Deformation of water-saturated sand: I. Uniform undrained deformation and shear band. Geotechnique 1996;46(3):441±56. [6] Vardoulakis I. Deformation of water-saturated sand: II. The eect of pore water ¯ow and shear banding. Geotechnique 1996;46(3):457±71. [7] Vardoulakis I, Sulem J. Bifurcation analysis in geomechanics. Blackie, Academic and Professional, 1995. [8] Shuttle DA, Smith IM. Localisation in the presence of excess pore-water pressure. Comput Geotech 1990;9: 87±99. [9] Loret B, Prevost JH. Dynamic strain localization in ¯uidsaturated porous media. J Engng Mech 1991;11:907±22. [10] Schre¯er BA, Majorana CE, Sanavia L. Shear band localization in saturated porous media. Arch Mech 1995;47:577±99. [11] Sanavia L, Zhang HW, Schre¯er BA. Strain localisation modelling in saturated sand samples. Forschungsbericht aus dem fachbereich Bauwesen 1995;66:357±66. [12] Schre¯er BA, Zhang HW, Pastor M, Zienkiewicz OC. Generalised plasticity constitutive model based strain localisation modelling in saturated sand samples. Computat Mech 1998;22:266±80. [13] Schre¯er BA, Sanavia L, Majorana CE. A multiphase medium model for localisation and postlocalisation simulation in geomaterials. Mech Cohes-Frict Mater 1996; 1:95±114. [14] Gawin D, Sanavia L, Schre¯er BA. Cavitation modelling in saturated geomaterials with application to dynamic strain localisation. Int J Numer Meth Fluids 1998;27: 109±25. [15] Bolzon G, Schre¯er BA, Zienkiewicz OC. Elastoplastic soil constitutive laws generalized to partially saturated states. Geotechnique 1996;46:279±89.
459
[16] Schre¯er BA, Sanavia L, Gawin D. Modelling of strain localisation in fully saturated soils. Proc 4th Int Worksh Localis Bifurcat Theory for Soils and Rocks, Gifu, Japan. 28 September±2 October 1997. Balkema, 1998. p. 305±12. [17] Alonso EE, Gens A, Hight DW. Special problem soils, General report. In: Hanrahan ID, Orr TLL, Widdis DF, editors. Proc 9th Eur Conf Soil Mech Fdn Engng, vol. III. Dublin. Balkema, Rotterdam, 1987. p. 1087±146. [18] Zienkiewicz OC, Chan A, Pastor M, Schre¯er BA, Shiomi T. Computational soil dynamics with special reference to earthquake engineering. Chichester: Wiley; 1999. [19] de Borst R, Sluys LJ, Muhlhaus HB, Pamin J. Fundamental issues in ®nite element analysis of localisation of deformation. Engng Comp 1993;10:99±121. [20] Zhang HW, Sanavia L, Schre¯er BA. An internal length scale in dynamic strain localisation of multiphase porous media. Mech Cohes-Frict Mater 1999;4:443±60. [21] Zhang HW, Schre¯er BA. Uniqueness and localisation analysis of elasto-plastic saturated porous media. Mech Cohes-Frict Mater, submitted for publication. [22] Lewis RW, Schre¯er BA. The ®nite element method in the static and dynamic deformation and consolidation of porous media. Chichester: Wiley; 1998. [23] Pastor M, Zienkiewicz OC, Chan AHC. Generalized plasticity and the modelling of soil behaviour. Int J Numer Anal Meth Geomech 1990;14:151±90. [24] Safai NM, Pider GF. Vertical and horizontal land deformation in a desaturating porous medium. Adv Water Resour 1979;2:19±25. [25] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 1. Averaging technique. Adv Water Resour 1979;2:131±144. [26] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 2. Mass, momenta, energy and entropy equations. Adv Water Resour 1979;2:191±201. [27] Hassanizadeh M, Gray WG. General conservation equations for multi-phase system: 3. Constitutive theory for porous media ¯ow. Adv Water Resour 1980;3:25±40. [28] Gray WG, Hassanizadeh M. Unsaturated ¯ow theory including interfacial phenomena. Water Resour Res 1991; 27(8):1855±63.