Constitutive model for monotonic and cyclic responses of loosely cemented sand formations

Mojtaa Rahimi,Dave Chan,Alireza Nouri

aDepartment of Mechanical Engineering,Azad University,Khomeinishahr Branch,Isfahan,Iran

bDepartment of Civil and Environmental Engineering,University of Alberta,Edmonton,T6G 2W2,Canada

cChina Three Gorges University,Yichang,China

Keywords:Cyclic loading Monotonic loading Cemented sand Plasticity Constitutive model

ABSTRACT This paper presents a model to simulate the monotonic and cyclic behaviours of weakly cemented sands.An elastoplastic constitutive model within the framework of bounding surface plasticity theory is adopted to predict the mechanical behaviour of soft sandstone under monotonic and cyclic loadings.In this model,the loading surface always passes through the current stress state regardless of the type of loading.Destruction of the cementation bonds by plastic deformation in the model is considered as the primary mechanism responsible for the mechanical degradation of loosely cemented sands/weak rock.To model cyclic response,the unloading plastic and elastic moduli are formulated based on the loading/reloading plastic and elastic moduli.The proposed model was implemented in FLAC2D and evaluated against laboratory triaxial tests under monotonic and cyclic loadings,and the model results agreed well with the experimental observations.For cyclic tests,hysteresis loops are captured with reasonable accuracy.

1.Introduction

This paper focuses on the constitutive model and degradation behaviour of cemented sand/soft sandstone under monotonic and cyclic loadings.Two main approaches for the study of weak sandstone degradation behaviour caused by repeated loading include laboratory testing by conducting cyclic loading tests and the development of cyclic plasticity theories.The majority of the studies in the area of cyclic loading of sandstone or cemented sand are limited to earthquake(dynamic)type of loading.Few studies have characterized the deformation properties of sandstone under slow cyclic loading both theoretically and experimentally.In general,the behaviour of geomaterials under cyclic loading is remarkably complex.This may be due to the dependence of the constitutive relationship on pressure and void ratio as well as the nonlinear behaviour of the sand matrix(Russell and Khalili,2004;Khalili et al.,2005,2006).Even the most sophisticated models cannot provide accurate predictions under general cyclic loading(O’Reilly and Brown,1991).Therefore,development of a reliable model to capture the cyclic behaviour of geomaterials has become one of the most challenging issues in constitutive modelling(Vermeer and de Borst,1984).

Most models for slow cyclic loading have been proposed for cohesionless soils and few studies have been conducted on sandstone behaviour in response to slow cyclic loading.Thus a model that could capture soft sandstone response under slow cyclic loading would be a significant advancement.The particular dif ficulty in the integration of the critical state concept in cyclic modelling is noted in the literature(Imam and Chan,2008).In this paper,a critical state constitutive model is presented for slowcyclic loadings for soft sandstone.

The traditional plasticity theory appears to be unsuitable for modelling cyclic loading since it predicts a purely elastic response during unloading and reloading within the yield surface.That is,no plastic deformation is predicted for unloading and reloading unless the stress path reaches the yield surface again(Chen and Han,2007).This is not suitable for modelling cyclic loading because,in reality,all unload-reload cycles result in the gradual accumulation of plastic strain and energy dissipation(Khong,2004;Lenart,2008)as shown schematically in Fig.1.In other words,the response in Fig.1 suggests that the loading and unloading stress paths are not the same.This is known as hysteresis,and it shows that the material fails to recover all the energy it receives in the loadingunloading cycle(O’Reilly and Brown,1991).This is attributed to energy dissipation due to plastic deformation(Lenart,2008).Hysteresis is the result of non-uniform deformation of the material in which different parts of the material are undergoing different stages of loading and unloading.

The effect of non-uniform deformation at different stages of loading and unloading can be illustrated using a friction block model.Uniform deformation is analogous to a single block as shown in Fig.2.In this case,there is only one displacement in the system,which is represented by u.There is no slipping until the horizontal force(T)reaches the maximum frictional force(Fk)when the block starts to move.At some displacement(u1),if T decreases below Fk,movement will cease immediately and the force will vary between zero and Fk.There is no movement until T reaches Fkagain and the displacement will continue from u1.There is no hysteresis effect,and there is no study reported during the unloading and reloading cycle at u1.

Non-uniform deformation in a material can be conceptually represented bya series of blocks as shown in Fig.3.In this case,four blocks are connected by three springs with stiffness(K),and each block also is subjected to a normal force(N).The horizontal force(T)is slowly increased until the first block on the right,block 4,starts to move.There is no horizontal force applied on the other blocks until the spring between blocks 4 and 3 starts to compress.The force that will transmit to block 3 will be equal tothe difference between T and Fk,where Fkis the frictional resistance at the base of each block.Again,there is no force applied on block 2 until the force in the spring between blocks 4 and 3 exceeds Fkin block 3.The process continues until all of the blocks start to move when T is equal to or exceeds 4Fk.Since the movement of all the blocks is in the direction of T,Fkwill be acting in the opposite direction(Fig.3a).

Fig.2.Single frictional block on a flat surface and corresponding force-displacement response.

If T decreases below 4Fkafter some movement of the blocks,the frictional force on block 4 will start to decrease until the direction is reversed as shown in Fig.3b.There is no movement of block 4 until T decreases below 2Fk.In this case,the force in the spring between blocks 4 and 3 will also decrease until the direction of the frictional force under block 2 reverses its direction.In this case,T will become zero,representing a fully unloaded state as shown in Fig.3c.

If T increases again,there is no movement in block 4 until its frictional force changes its direction and the value of T is equal to 2Fkas shown in Fig.3d.This represents that the reloading stage after T has been fully unloaded.It is clear that the reloading path is different from the unloading path since the mobilization of the frictional force under the blocks is different,unlike the case of a single block.Movement of block 4 occurs when T increases above 2Fk.When T is equal to 4Fk,the frictional forces under all of the blocks point in the same direction and movement will continue in the direction of T.Fig.4 shows force-displacement response of block 4.

As demonstrated in this simple system of blocks,the hysteresis effect is a result of non-uniform mobilization of frictional forces under the blocks since theyare connected by deformable springs.If the blocks are connected by rigid springs,the hysteresis effect will disappear.

In the case of a real material,since the stresses and strains in the material are generally non-uniform at the mesoscopic scale,it will give rise to the hysteresis effect much like the series of blocks connected by deformable springs.Therefore,energy dissipation occurs during the unloading and reloading process below the latest yield point.

The shortcomings of the classical plasticity theory led to extensive research,beginning in the 1960s,on developing more sophisticated plasticity models to capture the cyclic behaviour of geomaterials(Yu,2006).Advanced constitutive models that have been introduced within the plasticity framework include multisurface plasticity(Iwan,1967;Mroz,1967;Mroz et al.,1978,1979),bounding surface plasticity(Dafalias and Popov,1975;Krieg,1975;Dafalias,1982,1986;Bardet,1986;Khong,2004;Khalili et al.,2005,2006;Yang et al.,2011),generalized plasticity(Zienkiewicz et al.,1985;Pastor et al.,1985,1990;Ling and Yang,2006; Chung, 2010), and subloading surface plasticity(Hashiguchi,1989;Hashiguchi and Chen,1998).

Note that for our target material(i.e.naturally or artificially cemented sand),few constitutive models have been developed to simulate cyclic behaviour.For instance,Weng and Ling(2012)and Weng(2014)proposed their constitutive models based on generalized plasticity theory and verified their models against laboratory one-way cyclic loading of sandstone.Tariq and Maki(2012)conducted one-way cyclic loading tests on artificially cemented sand.However,they did not develop any constitutive model to simulate their experimental observations.Zhang et al.(2013)conducted one-way cyclic loading tests on a specific sandstone known as red sandstone.However,their proposed elastoplastic constitutive model was not verified against cyclic shear stress tests.Fu et al.(2014)developed a constitutive model within the framework of generalized plasticity and verified their model against experimental observations of one-way cyclic loading of rufous sandstone.Liu et al.(2016)conducted several laboratory one-way cyclic loading tests on structured(i.e.cemented)soils.However,they did not suggest any constitutive model to simulate experimental results.

Fig.3.(a)Four blocks in series where T is increased until movement occurs in all of the blocks;(b)T is decreased until the block 4 starts to move to the right;(c)T decreases to zero during unloading;and(d)T increases during reloading.

Fig.4.Force-displacement response of block 4 in the four-block system connected by springs.

The aim of this paper is to present a continuum elastoplastic constitutive model within the framework of bounding surface plasticity tomodel cyclic loadingof frictional and cohesive material.A critical state constitutive model proposed by Imam(1999)and Imam et al.(2005)is chosen to be the base model.This model was developed for simulating the behaviours of cohesionless sands under monotonic loading.Thus,modifications are incorporated into the model to include the behaviour of cemented sand/soft sandstone under monotonic loading.To predict cyclic behaviour using the bounding surface plasticity theory,normalized plastic and elastic moduli are modified under unloading conditions without using the mapping rule.Although the projection rule is not incorporated,concepts of the bounding surface plasticity theoryare used in the sense that plastic deformation is recorded for both loading and unloading conditions.This implies that the loading surface always passes through the current stress state regardless of the position of the stress path or the type of loading.Also by assuming a fixed size ratio between the bounding surface and the loading surface throughout the shearing process,the current stress always lies inside the bounding surface.However,unlike conventional bounding surface plasticity in which the plastic modulus is expressed as the summation of the additive and bounding surface plastic moduli,the plastic modulus is stated only as a function of that of the loading surface.This is similar to Imam and Chan(2008)approach for modelling the behaviour of cohesionless sand under cyclic loading.The proposed modified model ultimately is implemented in FLAC2D.

2.Cemented sand

Cemented soils and weak rocks constitute an intermediate class of geomaterials from a mechanical standpoint and are classified between soil and rock.They are often considered as non-textbook materials(Schnaid et al.,2001).Their mechanical behaviour under various conditions is not as well understood as those of hard rocks or soils.It is known,for example,that their strength and deformation characteristics show strong nonlinearity.Therefore,nonlinear elasticity models are more appropriate than linear elasticity models for soft rocks(Yoshinaka et al.,1998).

Cementation in naturally cemented sand originates from the precipitation of cementing agents like silica,carbonate minerals(especiallycalcite),clayminerals and iron oxides ontothe surface of sand particles(Clough et al.,1981;Fjaer et al.,2008).

Cementation in artificially cemented sand comes from an externally added cementing agent such as Portland cement(Consoli et al.,2012).With an increase in the degree of cementation,peak strength,cohesion,tensile strength and stiffness of the material increase,while strain at peak strengthdecreases.Although a rise in the friction angle due to the addition of cementation has been reported,there is no agreement on the effect of cementation on the peak friction angle of the material.Some researchers believe that cementation causes no change to the friction angle.A list of different researchers who concluded different results on the effect of cementation on peak friction angle can be found in Abdulla and Kiousis(1997)and Schnaid et al.(2001).Cementation also gives rise to a more brittle response(Clough et al.,1981)and is more important than friction at low confining pressures.The failure mode is believed to be brittle for weakly cemented sand at lower confining pressures and ductile at higher pressures(see Fig.5).Therefore,cemented sand shows a brittle response due to brittle failure of cementation bonds.At higher confining pressures,the frictional component is more dominant,resulting in a more ductile response.However,Fig.6 suggests that for stronglycemented sand,cementation providesnoticeable resistance even athigher confining pressures(Clough et al.,1981).Note that artificially cemented sands are usually considered as the better choice for conducting triaxial tests for constitutive model development because of their uniformity,homogeneity and reproducibility(Mohsin,2008).

Fig.5.Failure mode for weakly cemented sand.

Fig.6.Failure mode for strongly cemented sand.

It is believed that the shear stress-strain curve for(artificially or naturally)cemented sand will eventually approach that for uncemented sand at a critical state which is associated with zero dilatancy rate(see Fig.7).This,however,does not always occur since cementation may not be completely broken even after large amounts of shearing far beyond the initial yield point(Lee et al.,2004).That is,cemented sand after large amounts of shearing mayarrive at the ultimate void ratio instead of the critical state void ratio which is associated with constant volume and is independent of the initial state.This implies that the concept of a unique critical state line(CSL)does not apply for(artificially or naturally)cemented sand because the CSL depends on the initial condition and especially on the cement contents.This is one of the most crucial challenges in applying the critical state theory to cemented sand.Experimental determination of the CSL for cemented sand is also difficult due to its brittle behaviour and strain localization(Marri,2010).

3.Brief review of the original model

The detailed description of the original model can be found in Imam(1999)and Imam et al.(2005).A brief description is provided here.

The CSL in many proposed constitutive models for sand is used as a reference state to which various states of the soil are compared.Dif ficulty in determining the CSL position and the uncertainty of its position when the soil is loaded in different directions of shearing are the major setbacks for this approach.Errors in the determination of the CSL in these models can result in inaccuracies in the predicted behaviour(Imam,1999;Imam et al.,2005).The proposed constitutive model also relies on the CSL position to determine the soil state at large strain.However,the impact of the uncertainties in the CSL location on soil properties at smaller strains is compensated by measuring the important properties under different directions of shearing fromthe experiment and correlating these properties to the CSL(Imam,1999).The stress ratio at the peak of the undrained effective stress path(UESP)is one of these properties by which the stress ratio at the peak of the yield surface is predicted(Imam et al.,2002).Shear strains associated with the peak of the UESP are relatively small,implying that the measurement of this property is more reliable because strain in the early stages of the undrained tests is more uniform and less affected by localization(Imam,1999;Imam et al.,2005).Fig.8 illustrates typical UESP in monotonic triaxial compressive loading on loose sands.The deviatoric stress initially increases until it reaches a peak at point(pp，qpeak),where stress ratio Mpis mobilized.The deviatoric stress then decreases(strain softening)until UESP reaches the CSL(Imam et al.,2005).The yield function of the original model is expressed as(Rahimi et al.,2015)

Fig.7.Convergence of the shear stress-strain curves for uncemented and cemented sands(ideal case).

whereηis the stress ratio,Mpis the stress ratio at the peak of UESP,p is the mean effective stress,αis a scalar which represents anisotropic consolidation and its magnitude is zero for isotropic consolidation,and pcis the effective preconsolidation pressure which is considered as hardening parameter.

The parameter Mpis calculated for compressive and extensive loading,respectively,as follows(Rahimi et al.,2016):

where φp，cand φp，eare the friction angles at the peak of UESP in triaxial compression and triaxial extension,respectively.They are evaluated by(Imam et al.,2005):

whereψp=e-epis the state parameter at the peak in which e is the void ratio,and epis the critical state void ratio which is evaluated at mean effective stress corresponding to Mp(i.e.at p=pp);φμis the friction angle associated with ψp=0 in triaxial compression and is close to the interparticle friction angle;and kpand apare the model parameters.Experimental observations which support the dependencyof Mpto void ratio and the state parameter at the peak can be found in Imam et al.(2002).Note that the deviatoric stress(q)and the deviatoric strain(εq)are both positive under triaxial compression while they are both negative under triaxial extension.

When the stress path is at the peak of UESP,it is found from Eq.(1)that pp=0.64 pcfor isotropically consolidated sands regardless of void ratio and confining pressure at consolidation.This relationship is similar to that obtained by Ishihara(1993),who observed that the ratio of mean normal stressat the peak of UESP to that at consolidation is constant.By conducting undrained tests on sands consolidated to different mean normal stresses and void ratios,he obtained the ratios of 0.61 and 0.63 forclean and siltysands,respectively(Imam et al.,2005).

Isotropic nonlinear elasticity is adopted in the original model as

where K and G are the bulk and shear elastic moduli,respectively;Gaand Kaare the respective reference moduli associated with the atmospheric pressure patm;and e is the void ratio.

Following the work of Wood(1990)and Manzari and Dafalias(1997),the stress-dilatancy relationship is defined as(Rahimi et al.,2017):

where

where D is the dilatancy rate,MPTis the phase transformation stress ratio,φcsis the critical state friction angle,ψs=e-ecsis the state parameter(Been and Jefferies,1985;Jefferies,1993),ecsis the critical void ratio,and kPTand aPTare both model parameters.

Similar to how Mpis calculated using sinφpunder triaxial compression and extension conditions,MPTvalues under compression and extension conditions are obtained from sinφPT.Note that Eqs.(12)and(13)are similar to the relationship of Manzari and Dafalias(1997).Hardening in this model depends on the proximity to the critical state,in contrast to conventional critical state models which relate the size of the yield surface to void ratio(Jefferies,1993).Shear hardening law is expressed as(Imam,1999):

where kfis a material parameter.

4.A two-surface model for cemented sand

In line with the general tendencies to use as many simplifying assumptions as possible and to predict the stress-strain relationships using the simplest possible approach(Chen,1994),a relatively simple constitutive model is presented.For modelling the mechanical behaviour of cemented sands,the base model is modified similar to the simplified assumptions made by Gens and Nova(1993).Gens and Nova(1993)and Nova(2005)suggested that a yield surface which had been originally proposed for cohesionless soils can be used for cemented soils after some modifications are made.They suggested that bonding(cementation between soil particles)provides tensile strength(pt)and additional strength(po)to the material.In consideration of the additional strength,they suggested that the yield surface for the unbonded geomaterials must be enlarged towards the right.They also proposed that the yield surface must also be expanded towards the left to account for the tensile strength.Thus,more bonding leads to larger expansion of the yield surface towards both the right and left.Fig.9 illustrates the schematic representation of the modified yield surface/loading surface.pbcontrols the yielding of the bonded material in isotropic compression,which determines the size of the enlarged yield surface,and pocontrols the increase in the size of the initial elastic domain.The degree of bonding may be expressed as po/pc(Gens and Nova,1993).

To account for non-uniform deformation during cyclic loading,kinematic/anisotropic hardening parameters(i.e.paand qa)are also integrated into the formula of the modified loading surface.Bounding and loading surfaces,elastic properties, flow rule,and stress-strain relationships are the major components of the modified model which are discussed herein.

4.1.Bounding and loading surfaces

The bounding and loading surfaces are assumed to have the same shape for simplicity.In terms of the conventional triaxial parameters,the bounding surface(F)and loading surface(f)can be written respectively as(see Fig.10):

Fig.9.Schematic representation of the modified yield surface.

Fig.10.Schematic representation of the loading and bounding surfaces at the start of loading(solid surfaces)and immediately after elastoplastic loading starts(dashed surfaces).In addition to change of position,size of two surfaces also changes with further elastoplastic loading(isotropic hardening).Kinematic hardening has been significantly magnified for better illustration.

where q is the deviatoric stress,ptis the tensile strength,and paand qaare the components of the kinematic hardening tensor evolving with plastic strain.The kinematic hardening tensor is assumed to lie initially at the origin of the stress space,meaning the first time loading,i.e.(pa，qa)=(0，0).The superimposed bar denotes variables of the bounding surface.

To ensure that the current stress state will not cross the bounding surface,it is assumed that the initial ratio of size of the two surfaces remains constant during the shearing process.It is also assumed that the components of the kinematic hardening tensor always coincide with the two surfaces,i.e.(pa，qa)=(pa，qa).

pband ptare direct measures of the size of the loading surface and the tensile strength,respectively.They are defined using the following relationships:

where pois a measure of the increase in size of the uncemented yield surface due to strength increase by cementation,βis defined as the ratio between the tensile strength and po,and pcplays the same role as the preconsolidation pressure for uncemented sand.

whereγis a decay parameter which determines the rate of bond breakage,and dpoindicates the change in podue to the change in plastic shear strain.

Note that the origin of plastic deformation of sand under low shear stresses mainly comes from grain crushing at the particle contact.Gross slippage at the particle contact,however,is responsible for plastic deformation under high shear stresses because the mobilized shear stresses are suf ficient to overcome the resistance of the contact(Imam,1999).

The nonlinear kinematic hardening law of Armstrong and Frederick(1966)is adopted to control the evolution of the loading surface.For triaxial conditions,it is expressed as

where c1and c2are the model constants;andare the plastic deviatoric and volumetric strain increments,respectively;and˙ωis the accumulative plastic strain increment defined as

The parameters pαand qαare the components of the kinematic hardening tensor,which are defined for triaxial conditions as

4.2.Elastic properties

The definition of the elastic moduli in the original model is modified similar to Yu et al.(2007)approach as follows:

This definition provides two important features of the bonded geomaterials:(a)extra initial elastic stiffness due to bonding;and(b)the progressive reduction of elastic moduli due to the breakage of bonding during plastic straining(Yu et al.,2007).

The elastic properties for unloading conditions with a slight modification are given as

At the time, your Prince Charles barely glanced at the tiny baby sleeping in her cot. After all, bow could a twelve-year-old boy be interested in babies?

where|η|denotes the absolute value of the stress ratio.Incorporation ofηin the expression for the elastic moduli results in the prediction of a stiffer response at the commencement of unloading and then a softer response when the stress ratio approaches the smaller deviatoric stresses.The decreasing trend is particularly more pronounced when the deviatoric stress is very small(see Fig.11).It should be noted that a decrease inηor q(usually accompanied with a decrease in deviatoric/axial strain)during compressive loading signifies the unloading conditions.Similarly,an increase inηor q(normally followed by an increase in deviatoric/axialstrain)during compressive loading indicatesthe loading/reloading conditions.

Suggestion of these expressions for the unloading elastic moduli comes from experimental observations.That is,the stiffness shows a sudden increase immediately after the inception of unloading accompanied with a sudden decrease when soil undergoes further unloading(O’Reilly and Brown,1991).The elastic moduli for reloading conditions,however,are assumed to remain the same as those of the original loading.Although assuming that the elastic properties of the reloading conditions are in the same form as those for unloading conditions leads to better predictions of the hysteresis loops(see Fig.16),it results in numerical instabilities in this case in FLAC2D model.

4.3.Flow rule

Experimental observations suggest that dilatancy is affected by the degree of cementation between sand particles.Most of studies have shown that an increase in the cement content gives rise to an increase in dilatancy and subsequent decrease in compression(Clough et al.,1981;Marri,2010).Therefore,the original Rowe’s stress-dilatancy relationship proposed for cohesionless sands cannot provide a good description for calculating the rate of dilation of cemented sands/weak sandstone.

By incorporating cohesion for compression,Rowe’s stressdilatancy relationship can be calculated by(Yu et al.,2007;Rahimi,2014;Rahimi et al.,2015):

where D is the dilatancy rate,c is the cohesion,and Mcsis the critical state stress ratio.Cohesion is assumed to degrade with the total plastic strain increment as follows:

Fig.11.Stress ratio vs.Axial strain.1-2 shows a sudden increase in stiffness after start of unloading,2-3 shows a sudden decrease when the material undergoes further unloading.3-4 shows the softer response under unloading when the material approaches lower values of the stress ratio.4-5 shows a sudden increase in stiffness after start of reloading.5-6 shows softer response when the material is subjected to further reloading.

Fig.12.Measured and predicted responses of Tehran Alluvium cemented sand with cement content of 3%.

whereξindicates the rate of cohesion degradation.

Clearly,when cohesion approaches zero,the original Rowe’s stress-dilatancy relationship for cohesionless soils is recovered.The plastic potential function(g)associated with the dilatancy relationship in Eq.(31)takes the form(Yu et al.,2007)as follows:

Fig.13.Measured and predicted responses of a slightly weathered rock under monotonic loading.

Fig.14.Measured and predicted responses of a cemented carbonate sand in a triaxial drained cyclic compression test.

where φgcan be determined for any given stress state(p，q)by solving equation Eq.(33a)(Yu,2006).

Consistent with Imam(1999)and Imam et al.(2005),the flow rule is expressed based on a variable phase transformation stress ratio instead of a constant critical state stress ratio.Thus,for triaxial compression conditions,we have

Fig.15.Measured and predicted responses of rufous sandstone under a triaxial drained cyclic compression test.

where MPTis the phase transformation stress ratio.

4.4.Stress-strain relationships

Fig.16.Prediction of rufous sandstone behaviour using 1D model(reloading elastic moduli are defined similar to those of unloading conditions).

Constitutive equations are governed by the following equations in line with incrementally linear stress-strain relationships:where mpand mqare the components of the unit vector normal to the plastic potential surface which are obtained using the dilatancy relationship.L is determined from the following relationship:

where t=1 for compressive loading and t=-1 for extensive loading.For unloading,the components of the unit vector to the plastic potential surface are obtained by those defined already:

Hnis the normalized plastic modulus which is obtained fromthe following equation:

where pfis the failure mean normal stresswhich is calculated using an iterative method from the following equation:

uQin Eq.(49)is defined as follows:

Hnfor unloading is assumed to be related to loading as follows:

where Ruis a model constant.

Due to the decrease in the stress ratio,higher plastic strains are predicted for a given total strain increment with advancement of unloading during one-way cyclic loading,which is in accordance with experimental observations.

4.5.Material parameters

The parameters kp,φμ,φcs,kPT,Ga,Ka,h,kf,and CSL are material parameters of the original model(Imam,1999).Methods of determination of these parameters have been given in Imam et al.(2005).The following material parameters were introduced in the upgraded model:

(1)Initial value of po;

(2)Initial value of tensile strength to evaluateβ;

(3)Decay parameter of bond strength(γ);

(4)Cohesion(c);

(5)Rate of degradation of cohesion(ξ);

(6)Anisotropic hardening parameters(c1，c2);and

(7)Ratio of unloading plastic modulus to that of loading(Ru).

The initial value of pois equal to the difference between the initialvalues of pband pc.The initial value of pcis obtained from the maximum mean normal stress characterizing the initial elastic domain forcohesionless soil.The initialvalue of pbcan be evaluated by conducting an isotropic compression test on weakly cemented sand to determine the initial yield value of pb(Nova,2005).

The parameterβis obtained using the initial value of the tensile strength.If there is no data regarding the tensile strength,the unconfined compressive strength can be used for an approximation of the tensile strength.The tensile strength is in order of 5%-20%of the unconfined compressive strength(Nova,2005).

The parameterγcontrols the rate at which bonds are broken.It is determined by fitting theoretical results to the experimental data in an isotropic compression test.The higher theγvalue is,the faster the compression curve will become identical to that for uncemented material(Nova,2005).

The parameter c can be evaluated using the loading function or Mohr-Coulomb yield function.The parameterξcan be evaluated by fitting theoretical results of the volumetric behaviour to those of experimental observations.

The anisotropic hardening parameters c1and c2determine the contribution of anisotropic/kinematic hardening in the combined isotropic-kinematic/mixedhardeningrule.Largerdifferencebetween c1and c2causes largercontribution of kinematichardening to overall hardening.The constants c1and c2may be determined from stressstraincurveofuniaxialtests(Araujo,2002;DunneandPetrinic,2005).

The parameter Ruis determined by fitting model predictions to unloading experimental data.A larger value of Ruresults in a stiffer unloading response and a smaller unloading-induced plastic strain increment and vice versa.

4.6.Model performance

The performance of the proposed model is examined first against two sets of triaxial monotonic loading tests.Then,it is assessed against two triaxial one-way cyclic loading tests.All simulations performed by FLAC are based on axisymmetric conditions.Fig.12 shows the calculated and observed behaviours of an arti ficially cemented gravelly sand under triaxial compression monotonic tests.Hydrated lime is used as the cementing agent for these samples(Asghari et al.,2003).Table 1 lists the material parameters used for these analyses.Due tothe lackof data,some of the material parameters such as the CSL and φcsare selected based on the best match using one set of values for these parameters under a given cement content.Isotropic hardening is assumed for monotonic loading.Thus,a zero value is allocated for kinematic hardeningrelated constants.The povalue is assigned for different confining stresses under a given cement content such that a constant initial bonding is obtained.

Table 1 Material parameters used for calibration in Figs.12 and 13(p*=p+pt).

As seen in Fig.12,the predicted and measured volumetric behaviours are in good agreement.The change of behaviour from contractive to dilative is observed in all tests.This is predicted numerically in the proposed model by a change of dilatancy rate from positive to negative.Fig.12 also indicates less dilation at higher confining pressures as expected.Unlike volumetric behaviour,however,there is small discrepancy in the predicted and observed stress-strain behaviours.Using one critical state friction angle has led to slight underestimation of the peak deviatoric strength especially for the test with the confining pressure of 110 kPa,demonstrating the difficulty in applying the critical state theory for cemented soil.This difficulty has been observed experimentally as well.For instance,Lee et al.(2004)experimentally measured a larger critical state stress ratio and thus a greater critical state friction angle at lower confining pressures for a cemented soil with a given cement percentage.

Fig.13 shows the predicted and observed behaviours of a slightly weathered rock obtained from the site of a pumped storage power station(Fu et al.,2014).The CSL and cohesion intercept are chosen based on the best match for the rock behaviour under various confining stresses due to lack of data.Isotropic hardening is assumed for the material behaviour under monotonic loading.The povalue is also allocated in a way similar to the procedure mentioned for Tehran Alluvium cemented sand in Fig.12.The material parameters used in Fig.13 are listed in Table 1.

The model was next examined against the cyclic triaxial data(Mohsin,2008)for a cemented carbonate sand with 20%gypsum content(gypsum was used as the cementing agent)and an initial dry unit weight of 13 kN/m3(see Fig.14).The model is able to predict the narrow hysteresis loops,but the loops may not be shown clearly in Fig.14.In this calculation,some of the material parameters were chosen based on the best match.Table 2 lists the material parameters used for the calibration of the triaxial cyclic and monotonic compression tests shown in Figs.14-16.

The model was then assessed against cyclic triaxial data of the rufous sandstone,which is a slightly weathered rock(Fu et al.,2014).Due to the absence of required data,some material parameters are assigned to provide the best match to the observed behaviour.As can be seen from Fig.15,the FLAC2D model is reasonably capable of capturing the observed shear and volumetric cyclic behaviours of the rufous sandstone.However,there is a discrepancy between the measured and predicted volume changes during unloading.Experimental observations show that the material contracts during unloading while the model predictionssuggest that it undergoes expansion.However,in general,the measured and predicted volumetric behaviours are in relatively good agreement.

Table 2 Material parameters used for calibration in Figs.14-16(p*=p+pt).

Fig.16 shows the model predictions using one-dimensional(1D)model(i.e.the constitutive equations are solved for one node)for the same initial conditions and material parameters as those of Fig.15,when reloading elastic moduli are assumed similar to the unloading elastic moduli(see Eqs.(29)and(30)).As observed,the hysteresis loops are captured better in this case using 1D model.However,discrepancy between the measured and predicted volumetric responses exists similar to Fig.15.Note that numerical instabilities were encountered in the FLAC2D model when the reloading moduli were defined using Eqs.(29)and(30)rather than Eqs.(27)and(28).Therefore,to avoid the instability,the reloading elastic moduli were chosen to be the same as the loading elastic moduli(Eqs.(27)and(28))in FLAC simulations.

5.Conclusions

A continuum elastoplastic constitutive model has been presented in this paper to investigate the degradation of weak arti ficially and naturally cemented sand under monotonic and one-way cyclic loadings.Destruction of bonds between sand particles due to plastic deformation has been considered as the reason for mechanical degradation under applied loads.Several triaxial monotonic and cyclic compression tests on both artificially and naturally cemented sands have been chosen to assess the model performance.A comparison between the predicted and observed behaviours shows reasonably fair agreement.Hysteresis loops have beencapturedwithreasonableaccuracytoo.However,the observed softening response cannot be reproduced by the model well.Moreover,the increase in volumetric strain(i.e.further compaction)upon unloading cannot be modelled by the proposed constitutive model.Modification of the stress-dilatancy relationship under unloading conditions appears to be a potential candidate for addressing this shortcoming.

Conflicts of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have in fluenced its outcome.

Acknowledgements

Funding for this project has been provided by BP Canada and the Nature Science and Engineering Research Council of Canada(NSERC),for which the authors are grateful.