Influence of particle contact models on soil response of poorly graded sand during cavity expansion in discrete element simulation

Yng Dong,Behzd Fthi,*,Hdi Khz,Henry Zhng

aSchool of Civil and Environmental Engineering,University of Technology Sydney(UTS),Sydney,Australia

bWSP in Australia&New Zealand,Sydney Of fi ce,Australia

Keywords:Discrete element method(DEM)Cavity expansion Contact models Behaviour of geomaterials

A B S T R A C T The discrete element method(DEM)has been extensively adopted to investigate many complex geotechnical related problems due to its capability to incorporate the discontinuous nature of granular materials.In particular,when simulating large deformations or distortion of soil(e.g.cavity expansion),DEM can be very effective as other numerical solutions may experience convergence problems.Cavity expansion theory has widespread applications in geotechnical engineering,particularly to the problems concerning in situ testing,pile installation and so forth.In addition,the behaviour of geomaterials in a macro-level is utterly determined by microscopic properties,highlighting the importance of contact models.Despite the fact that there are numerous contact models proposed to mimic the realistic behaviour of granular materials,there are lack of studies on the effects of these contact models on the soil response.Hence,in this study,a series of three-dimensional numerical simulations with different contact constitutive models was conducted to simulate the response of sandy soils during cylindrical cavity expansion.In this numerical investigation,three contact models,i.e.linear contact model,rolling resistance contact model,and Hertz contact model,are considered.It should be noted that the former two models are linear based models,providing linearly elastic and frictional plasticity behaviours,whereas the latter one consists of nonlinear formulation based on an approximation of the theory of Mindlin and Deresiewicz.To examine the effects of these contact models,several cylindrical cavities were created and expanded gradually from an initial radius of 0.055 m to a final radius of 0.1 m.The numerical predictions confirm that the calibrated contact models produced similar results regarding the variations of cavity pressure,radial stress,deviatoric stress,volumetric strain,as well as the soil radial displacement.However,the linear contact model may result in inaccurate predictions when highly angular soil particles are involved.In addition,considering the excessive soil displacement induced by the pile installation(i.e.cavity expansion),a minimum distance of 11a(a is the cavity radius)is recommend for practicing engineers to avoid the potential damages to the existing piles and adjacent structures.

1.Introduction

The rapid advancement in computational technology has facilitated the application of discrete numerical analysis for many complex geotechnical problems, especially investigations involving large deformations or displacements.For instance,analysing the installation mechanism of driven piles and vertical drains adopted for ground improvement requires cavity expansion theory to better interpret the smear zone characteristics.Huang and Ma(1994)and Jiang et al.(2006)have employed discrete element simulation to study the mechanism of deep penetration in granular materials.Arroyo et al.(2011)and Ciantia et al.(2016)successfully modelled the cone penetration tests in a calibration chamber while taking the particle crushing behaviour into consideration.Geng et al.(2013)simulated the cylindrical cavity expansion and pressuremeter test using discrete element analysis.All of these studies have indicated that the discrete element method(DEM)is capable and reliable for simulating the behaviour of granular materials,while comprehending the microscopic contact properties employed by the discrete numerical studies is of paramount importance.

In DEM simulations,particles interact at pairwise contacts by means of internal forces and moments.In addition,contact mechanics are embodied in particle-interaction laws that are responsible for recognising particle interfaces,updating the internal forces and moments,computing the induced displacements and rotations of all particles involved in the system(Cundall and Strack,1979;Potyondy and Cundall,2004).Hence,it can be explicitly concluded that only by adopting the most appropriate contact model,the simulation can attain the most accurate outcome.There are mainly three types of contact models widely used by researchers that are capable of simulating granular materials such as sand,i.e.linear,rolling resistance,and Hertz contact models.Linear contact model,developed by Cundall and Strack(1979),describes the constitutive behaviour in normal and tangential directions between particles at their contact interface adopting linear springs with normal and shear stiffnesses.This model has been extensively employed to investigate granular materials in many disciplines ranging from chemical to geotechnical engineering.Its capabilities in mimicking the behaviour of sandy soils under various loading conditions have been successfully verified by many researchers(Jenck et al.,2009;Zhao and Evans,2009;Chen et al.,2011;Zhou et al.,2012;Lai et al.,2014;Falagush et al.,2015;Tran et al.,2016).Despite the fact that classical theories believe that it is mainly the sliding between adjacent particles that plays the dominant role in controlling the strength and dilatancy,Oda and Kazama(1998)observed that the effect of rolling between contacting particles is very significant.Therefore,the rolling resistance contact model has been developed based on the linear contact model by introducing a rolling resistance mechanism that incorporates a torque acting on the contacting pieces to counteract the rolling motion(Iwashita and Oda,1998;Oda and Kazama,1998;Irazábal et al.,2017).It is indispensable that the rolling resistance model is suitable to simulate granular materials particularly when most grains are angular,intensifying the interlocking of particles(Tordesillas and Walsh,2002;Jiang et al.,2005;Ai et al.,2011;O’Sullivan,2011).In addition to these two linear based models,Hertz contact model also demonstrates its competency in capturing the behaviour of granular materials in discrete numerical modelling.The most distinct feature differentiating Hertz contact model from others is that it has nonlinear formulations and utilises the simplified theory of Mindlin and Deresiewicz as force-displacement laws to compute contact forces(Mindlin,1953;Tsuji et al.,1992).As the Hertz contact model requires more accessible and macroscopic parameters that can be directly correlated to the macroscopic parameters measured in the laboratory,such as Poisson’s ratio and shear modulus,it is also considered as a simple yet efficient constitutive contact model by many researchers(Gu and Yang,2013;Ng and Meyers,2015;Ning et al.,2015;Zeghal and Tsigginos,2015;Ciantia et al.,2016).Although those three contact models have already well demonstrated their own competency in simulating the behaviour of sandy soils,there is very limited research to investigate the effects of these contact models on soil response.Furthermore,the cavity expansion theory,an effective technique to interpret the mechanism of pile installation and in situ tests,has been extensively employed by many researchers.However,obtaining closed-form solutions may become difficult when more realistic constitutive models are incorporated,and therefore numerical approaches can offer a better solution(Carter et al.,1979a).Hence,there are numerous analyses available adopting numerical techniques and simulations based on continuum method for cavity expansion(Carter et al.,1979b;Carter and Yeung,1985;Mo et al.,2014;Li et al.,2017).However,the discrete element studies may provide researchers with a microscopic level of insight to study the complex cavity expansion related problems.Therefore,the purpose of this study is to compare the results obtained from the cavity expansion simulation adopting different contact models and evaluate the effectiveness of each contact model in discrete element simulations.

2.Material contact models

2.1.General characteristics of adopted contact constitutive models

In this study,three contact models have been employed to mimic the interaction between particles at their contacting points,i.e.linear,rolling resistance and Hertz contact models.It is noteworthy to state that all adopted contact models are elasto-perfectly plastic,in which the former two contact models are linear elasto-perfectly plastic whereas the latter is elas to perfectly plastic model with nonlinear formulations.Therefore,all three contact models can capture both elastic and plastic behaviours.In the linear contact model,the behaviour is linearly elastic in normal direction with zero tension cut-off,while in the tangential direction,a frictional behaviour using linear springs and plastic sliders is utilised,as illustrated in Fig.1.In this contact model,slip is accommodated by imposing a Coulomb limit on the shear force using the friction coefficientμ(i.e.slider).Both normal and shear forces are computed linearly using the linear force-displacement law based on the linear springs,and the maximum shear force at the contact is the frictional force(μFn),whereFnis the normal force acting at the contact interface.

Fig.1.Schematic diagrams of(a)linear contact model and(b)rheological model.gsis the surface gap(non-tension joint);knand ksare the normal and shear stiffnesses of the linear springs,respectively;and are the normal and shear dashpots,respectively;R1and R2are the radii of two contacting particles;O1and O2are the centres of the contact interface between balls A and B and balls A and B′,respectively;and Δδsis the shear displacement increment in a time step.

Rolling resistance contact model,as depicted in Fig.2,developed from linear contact model,has a similar behaviour,except that the internal moment is incremented linearly with the accumulated relative rotation of the contacting pieces at the contact points.The limiting torque in the contact is proportional to the normal contact force and acts in the direction opposite to the rolling direction.

Hertz contact model,utilising a nonlinear formulation based on an approximation of the theory of Mindlin(1953),incorporates both normal and shear forces based on the theoretical analysis of the deformation of elastic spheres.Fig.3 describes the rheological model for Hertz contact mode.The main difference between Hertz and linear models is that the Hertz contact model employs nonlinear formulations to compute the normal and shear forces at the contact.The maximum shear force at the contact is also determined based on the frictional force,whereis the Hertz normal force.

2.2.Formulations of adopted contact constitutive models

Fig.2.Schematic diagrams of(a)rolling resistance contact model and(b)rheological model.kris the rolling stiffness of the linear spring,L is the distance between the centres of balls C and D,Δφris the rotation increment in a time step,and Mrris the rolling resistance moment.

Fig.3.Rheological model for Hertz contact model.Hnand Hsare the Hertz normal coefficient and tangent shear stiffness of the springs,respectively;and and are the normal dashpot and shear dashpot,respectively.

In discrete element simulations,the behaviour of granular materials is fundamentally determined by the adopted contact properties at the contact interfaces following the Newton’s second law of motion and force-displacement law.Despite the fact that the particles in discrete element simulation may be rigid bodies with finite mass that move independently of one another,the behaviour of the contacts is characterised using a soft contact approach allowing particles to overlap in the vicinity of the contact point.The contact between particles is detected prior to the force-displacement computation based on the contact overlapping(gs)detected(Cundall and Strack,1979;Cundall and Hart,1992;Itasca,2016).The linear contact model adopted in this study corresponds to the model developed by Cundall and Strack(1979),which utilises two linear springs in normal and tangential directions to mimic the linearly elastic and frictional behaviours between contacting particles.The rheological model and the operating mechanism are illustrated in Fig.1.The magnitudes of the normal force(Fn)and the shear force(Fs)are determined by the normal stiffnessknand shear stiffnessks,respectively,based on the following equations:

wheregsis the surface gap,measuring the overlapping between particles(i.e.gs=0 when particles contact yet not overlap,andgs＜0 when particles overlap);dis the distance between two contacting particles;andis the initial shear force.

In DEM simulations,the magnitudes of the normal stiffnessknand shear stiffnesskscan be specified by the users directly.In this circumstance,referring to Fig.4a and b,applying the same compressive force(F),the normal contact forces between two particles(e.g.A1 and B1 or A and B)will be doubled if the particles are enlarged twice.Therefore,it is obvious that the magnitude of the normal contact force increases with the increase inparticle size.Therefore,it is evident that the normal stiffnessknshould be upscaled,accordingly to response to the increase in normal contact forces.However,besides the direct approach,the linear model deformability method can be alternatively employed to control the normal and shear stiffnesses in real time by introducing the effective Young’s modulusE*and the stiffness ratiok*of the normal and shear stiffnesses of the linear springs based on the following equations(Itasca,2016):

Fig.4.Illustration of particle assembly in DEM:(a)Particle upscaling factor=1,and(b)Particle upscaling factor=2.

The effective Young’s modulusE*and the stiffness ratiok*have direct correlations to the material macroscopic properties,such as the Young’s modulus(E)and the Poisson’s ratio(ν).Consequently,referring to Eq.(2),it can be concluded that adopting the linear deformability model can automatically take the particle size effect into account,which is especially advantageous when the particles are upscaled to reduce the computational effort so that large-scale simulation can be conducted.

Rolling resistance contact model,developed based on the linear contact model,is used to mimic the rolling effect between particles with angular shapes(Iwashita and Oda,1998;Oda and Kazama,1998;Jiang et al.,2005;O’Sullivan,2011).Hence,it is considered as a more realistic contact constitutive model for granular materials where the rolling of particles is dominant in determining the strength.The formulations of the rolling resistance contact model are similar to the linear contact model except the incorporation of the rolling resistance moment(Mrr)acting at the contacting points of particles to counteract the relative motion and therefore simulate the interlocking behaviour between contacting particles(Ai et al.,2011;Wensrich and Katterfeld,2012).The rheological model and the schematic diagram of the particle interaction are illustratedin Fig.2.The magnitude of the rolling resistance moment(Mrr)is computed based on the following equations(Iwashita and Oda,1998;Jiang et al.,2005):

In addition to these two linear based contact constitutive models,Hertz contact model,as illustrated in Fig.3,is employed to simulate the behaviour of granular materials(Gu and Yang,2013;Ng and Meyers,2015;Zeghal and Tsigginos,2015;Ciantia et al.,2016).The Hertz contact model adopts the simplified Mindlin and Deresiewicz theory(Mindlin,1953;Itasca,2016)to reproduce the nonlinear force-displacement formulation to calculate the normal and shear forces during analysis using the following equations(Holt et al.,2005;Agnolin and Roux,2007):

whereis the Hertz shear force,Gis the effective shear modulus,αHis the Hertz exponent,is the effective radius,andis the initial Hertz shear force.

It can be clearly noted that the Hertz contact model utilises the material properties such as Poisson’s ratio and effective shear modulus to define the contact behaviour.The size of the particlesis taken into account during the force-displacement computation,confirming that the Hertz contact model,similar to the linear model deformability method,can be applicable when particles are upscaled.

Furthermore,it is worthy to mention that the allowed overlapping at the contacts between particles in the normal direction is directly related to the normal stiffness of the spring and the normal forces acting on the contact plane for the linear and rolling resistance contact models.This paper employs the linear model deformability method,and therefore the allowed overlapping is determined by the effective modulus,E*,specified.Referring to Eq.(1),the detected overlapping(gs)and the normal stiffness of the linear spring between two particles are used to determine the normal force(Fn).Also,according to Eq.(2),the allowed overlapping （gs）can be expressed in terms of the normal force(Fn）and effective modulus(E*）as follows:

Hence,it can be explicitly concluded that the overlapping is determined by the magnitude of the normal force acting on the contact plane,and the effective modulus of particles.Given a constant normal force,the overlapping between contacting particles decreases as the effective modulus increases.For the Hertz contact model,referring to Eq.(4),in a similar way,the allowed overlapping can be determined as below:

Hence,it can be concluded that the allowed overlapping at the contacts is determined bythe effective shear modulus(G),Poisson’s ratio(ν),as well as the Hertz normal force）.

In this study,PFC3D5.0 software(Itasca,2016)has been adopted to evaluate the performance of these three contact models for cavity expansion simulation.Moreover,it is noted that the constitutive behaviours of linear contact model,rolling resistance contact model as well as Hertz contact model discussed above contain only the essential characteristics and readers can refer to existing literature for further detailed explanations(e.g.Mindlin,1953;Cundall and Strack,1979;Iwashita and Oda,1998;Itasca,2016).

It is acknowledged that disparities may exist between the microscopic parameters of these adopted contact models and the actual behaviour that can be investigated when laboratory tests are done in microscopic scale.For instance,according to micromechanical experiments conducted by Sandeep and Senetakis(2017,2018),with the increase in the normal loads,repeated shearing causes increase in frictional force.However,the frictional forces between particles for a given normal force in many existing contact constitutive models are constant(e.g.linear,rolling resistance,and Hertz contact models).The advanced micromechanical interparticle loading apparatus can directly measure the contact properties such asknandksbetween particles,which indisputably lead to the new direction in the micromechanical studies.It also offers the possibility to study the behaviour of the contacting particles in a microscopic level,which is advantageous for the development of more comprehensive and realistic contact constitutive models adopted in the DEM analysis(O’Sullivan,2011;Senetakis et al.,2013;Senetakis and Coop,2014).

3.Numerical modelling

3.1.Calibration of contact parameters adopting triaxial test

Investigating the effects of contact models on the macroscopic soil response can never be attained if the variables are not controlled since each contact model contains multiple parameters that fundamentally determine the contact behaviour between particles.Hence,three adopted contact models were calibrated to represent realistic sandy soil by matching the results obtained from the numerical simulations of triaxial compression tests against existing laboratory results.Typically,the stress-strain relationships obtained from the experiments are used for calibration of the contact models,which is a common practice adopted by many researchers(Chareyre et al.,2002;Arroyo et al.,2011;Wang et al.,2014;Ciantia et al.,2016).In this study,the contact properties were calibrated by matching the overall(macroscopic)stressstrain and volume change-axial strain curves obtained for the soil specimen.

3.1.1.Selection of existing experimental triaxial test results for calibration exercise

Experiments adopted for calibration exercise was conducted by Cornforth(1964)to compare the drained strengths of medium to fine grained sand under the plane strain condition as well as the conventional triaxial condition.Detailed experimental procedure and specimen information can be found in Cornforth(1964).It is preferred to use comprehensive test results(e.g.triaxial test with various confining pressures)for calibration of the micromechanical properties of the soil for general simulation of soil under different confining pressures.However,it should be noted that the calibration exercise conducted in this study compared predictions and measurements for both axial stress and the volumetric strain variations with shear strains.In addition,since the test results reported by Cornforth(1964)include the required information for calibration of micromechanical properties,i.e.particle size distribution,void ratio,volume change with shear strain,stress-strain results,as well as the experimental procedure and the loading rate,these results have been used in this study.Similar calibration techniques were also adopted by many other researchers(e.g.Arroyo et al.,2011;Ciantia et al.,2015).

3.1.2.Numerical simulation of triaxial test

Considering the current computational power,it is not realistic to simulate the triaxial test using the true particle sizes.For instance,a laboratory-scale triaxial test on a sand specimen with 30 mm diameter and 60 mm height would comprise approximately 4 million particles.Directly mapping such a huge number of particles into the simulation is unrealistic and computationally infeasible.Hence,only by reducing the number of particles,the analyses can become efficient and possible(O’Sullivan,2011).Therefore,before the initial condition has been assigned to the particles,an upscaling factor should be applied to reducing the number of particles immensely so that the simulation becomes computationally feasible.An appropriate upscaling factor is of paramount importance to ensure that the soil properties will not be altered.In this case,a very small upscaling factor has been initially adopted and has been gradually increased to 16.5 so that the porosity can steadily reach the desired value.The particle size distribution adopted by Cornforth(1964)and the uniformly upscaled particle size distribution employed in the calibration exercise are presented in Fig.5.Hence,the triaxial simulations consisted of 20,000 spherical particles after upscaling,with the maximum diameter of the particle after upscaling being 12.35 mm while the minimum particle has a diameter of 1.245 mm,and the setup of the triaxial test is shown in Fig.6.As illustrated in Fig.6,the generated model for calibration exercise has a length of 406 mm,while the height and width are 101.6 mm and 50.8 mm,respectively,which are identical to specimen size adopted in the experimental test conducted by Cornforth(1964).The initial porosity of the generated DEM model for calibration exercise is 0.392,similar to the porosity measured in the laboratory,corresponding to a medium dense specimen.The common approach adopted to generate a relatively dense specimen is to apply a small frictional force at grain contacts during the specimen generation stage so that interparticle sliding can occur in a relatively effortless manner.The boundary walls are controlled using a servomechanism until the target void ratio and stress state are achieved at equilibrium.Once the initial conditions are satisfied,the frictional force can be adjusted back to the value determined from the calibration exercise(Thornton,2000;Chareyre et al.,2002;Jiang et al.,2003;O’Sullivan,2011;Tong et al.,2012).This technique is also adopted in this study to prepare a medium dense specimen.The lateral confining pressure(Ydirection)of 275 kPa was maintained while the facets in the intermediate direction(Xdirection)were fixed to simulate the plane strain condition,and then axial load was applied on the top of the specimen(Zdirection)in a constant rate.During the simulation,stress and strain variations as well as volume changes were continuously recorded.To investigate the effects of different contact models on soil response,rolling resistance,linear and Hertz contact models have been calibrated.

3.1.3.Calibration methodology

Fig.5.Particle size distribution adopted in experimental test,calibration exercise and cavity expansion simulation.

Fig.6.Triaxial compression test in plane strain condition in discrete element simulation.

The linear and the rolling resistance contact models are based on the linear deformability model,which is a method controlling the normal and shear stiffnesses of the linear springs based on the effective modulus （E*）and the normal to shear stiffness ratio （k*）,as denoted in Eq.(2).The effective modulus （E*）is related to the Young’s modulus（E）of the material,andEincreases asE*increases,confirming a direct correlation.In addition,the Poisson’s ratio(ν)is related to the normal to shear stiffness ratio （k*）,with νincreasing up to a limiting positive value ask*increases(Itasca,2016).Therefore,it can generally be concluded that the effective modulus（E*）controls the elastic part of the stress-strain curve,while the normal to shear stiffness ratio （k*）influences the volumetric behaviour.The selection of the values ofE*andk*adopted for the initial trial is based on the recommendations made by Plassiard et al.(2009),considering a comprehensive parametric study for the triaxial test simulation.Since the macroscopic Young’s modulus(E)of the sand in the calibration exercise is approximately 70 MPa,as reported in Cornforth(1964),the corresponding magnitude of the effective modulus （E*）can be estimated to be approximately 1000 MPa referring to Plassiard et al.(2009).In addition,an interparticle friction coefficient of 0.5 is selected as the initial trial value.Considering the initial value of the effective modulus,a series of sensitivity analyses is conducted to investigate the effects of the interparticle friction coefficient and the normal to shear stiffness ratio （k*）using control variable method,as illustrated in Figs.7 and 8.

According to the preliminary studies on the influence of the interparticle friction coefficient reported in Fig.7,increasing the friction coefficient can,by and large,enhance the material shear strength,especially when the friction coefficient is less than 0.5.However,when the interparticle friction coefficient is larger than 0.5,the increase of shear strength becomes rather inapparent.Huang et al.(2014)conducted a comprehensive DEM investigation exploring the influence of interparticle friction coefficient on critical state behaviour of granular materials.The study indicated that the interparticle friction coefficient can affect the critical state parameterMas well as the position of the critical state line ine-log10p′space(eis the void ratio andp′is the mean effective stress),particularly when the friction coefficient is less than 0.5,whereas such effects are insignificant when the friction coefficient increases beyond 0.5.In addition,more rolling was observed in comparison to the sliding occurring at the contact points when a large interparticle friction coefficient was adopted.In addition,Yan et al.(2009)conducted a series of numerical studies investigating the influence of micro-properties on the macroscopic stress-strain behaviour of granular material.The results confirmed that the interparticle friction coefficient affects both strength and overall volumetric dilation.The peak of the stress-strain curve increases with the increase in friction coefficient,and the volumetric dilation becomes more significant with the increase of the interparticle friction coefficient.Similar findings were also reported by Yang et al.(2012)and Hosn et al.(2017)in studying the influence of interparticle friction coefficient and dilatancy of granular materials.All of these studies carried out confirm that the interparticle friction coefficient indeed influences the soil behaviour significantly.

However,referring to Fig.8,it is observed that,in this study,the normal to shear stiffness ratio （k*）has insignificant effect on the stress-strain relationship compared with the effect induced by the interparticle friction coefficient.Hence,the normal to shear stiffness ratio could be fixed provisionally.

The interparticle friction coefficient is then adjusted to match the rest of the stress-strain curve to obtain the best match,then the interparticle friction coefficient is fixed and then the effective modulus（E*）and the normal to shear stiffness ratio （k*）are further adjusted to achieve the optimum match with the experimental results.The method discussed above is the technique adopted for the linear model calibration.The technique adopted by authors to calibrate the rolling resistance and the Hertz contact models is similar to the one reported above.In summary,the following procedure can be used for the calibration:

Fig.7.Influence of interparticle friction coefficient on the axial stress-strain relationship.

Fig.8.Influence of normal to shear stiffness ratio(k*)on the axial stress-strain relationship.

(1)Estimating the initial trial values of micromechanical properties based on the tests and parametric studies carried out by other researchers;

(2)Conducting a sensitivity analysis and then fixing the parameters influencing the predictions to less extent in the range of stresses applied;

(3)Adjusting the parameters influencing the stress-strain predictions notably until a reasonable match is achieved;and

(4)Optimising further by fine-tuning all parameters to achieve the best possible predictions.

Further information regarding the calibration procedure adopting the Hertz contact model can refer to Cheng et al.(2017).

3.1.4.Calibration results

Calibration results for stress-strain relationship as well as volumetric behaviour are plotted in Figs.9 and 10,respectively.As observed,the axial stress shows a continuous rise with the increase in axial strain until the peak stress of 1150 kPa has been reached,while strain softening has been observed when the axial strain is greater than 2.2%.Additionally,the volumetric strain demonstrates a contraction at the initial stage(i.e.axial strain＜1%)and dilation when the simulation continues,confirming the medium dense state of the specimen simulated in DEM.Referring to Figs.7 and 8,it can be concluded that the numerical predictions have a good agreement with experimental data,indicating reliability of calibrated parameters reported in Table 1.

3.2.Cavity expansion simulation

Fig.9.Comparisons of axial stress-strain relationship obtained from calibration numerical simulation and experimental results.

Fig.10.Comparisons of variations in volumetric strain with axial strain obtained from calibration numerical simulation and experimental results.

Table 1 Summary of calibrated contact parameters for poorly graded sand.

The objective of this study is to investigate the effects of different contact models on soil response during the cavity expansion process.To achieve this goal,a large-scale numerical model simulating the cylindrical cavity expansion has been proposed adopting PFC3Dsoftware and calibrated parameters.Considering the fact that many cavity expansion related problems,such as pressuremeter test and driven pile installation,are in plane strain condition(Alshibli and Sture,2000),and taking the advantage of the axisymmetric geometry,only a quarter of the geometry with plane strain condition has been simulated,as shown in Fig.11.The discrete element model established for three different contact models contains 75,010 spherical particles,generated using the radius expansion method,with its external boundary located 1 m away from the centre of the initial cavity,as shown in Fig.11.The measured soil initial porosity is approximately 0.39,similar to the porosity reported in the experiment used for calibration of medium dense sand.The initial friction coefficient assigned to particles is 0.2 at the generation stage,and the internal and external boundaries as well as the top and bottom walls are servo-controlled(as illustrated in Fig.11)to satisfy the stress state and the void ratio.Once the target void ratio and stress state are reached at equilibrium,the friction coefficient is set to that determined in the calibration exercise.The initial stress field adopted was 500 kPa(σr0=σX= σY=500 kPa),and wall servomechanism was enabled to ensure that the external pressure remains equal to the initial stress field at all time at the outer boundary during the cavity expansion process.After the initial condition has been satisfied,all three groups of analyses adopting different contact models have the same properties,including the same initial stress field and the porosity equaling 500 kPa and 0.39,respectively.In addition,the average number of active contacts per particle(i.e.coordination number)was 5.7 for all contact models,confirming that the specimen is in medium dense state according to the empirical equation to relate coordination number and porosity proposed by Mitchell and Soga(2005).

Three groups of simulations were conducted using different yet calibrated contact models(see Table 1),and the above conditions were consistently applied to all analyses.After the initial conditions were satisfied,the cylindrical cavities were expanded gradually from an initial cavity radiusa0=0.055 m to a final cavity radiusaf=0.1 m in a constant strain rate of 0.001 m/s to ensure quasistatic loading condition.Internal cavity pressure(Pa)is measured using appropriate subroutines that can obtain the contact forces acting on the internal cavity wall and then divided by the corresponding contact area.On the other hand,there were five equally spaced gauge particles selected along angular bisector,as shown in Fig.11,and the displacements and positions of these gauge particles are continuously monitored during the cavity expansion.Additionally,the radial stress(Pr)distributions at various stages of the cavity expansion are also recorded using numerous prediction spheres created in the radial direction.Furthermore,referring to Fig.11,the variations of the deviatoric stresses and volumetric strains with the shear strain during the cavity expansion are continuously recorded adopting two particular predication spheres A and B(see Fig.9),which are situated atrA/af=2.5 andrB/af=5,respectively.

4.Results and discussion

4.1.Cavity pressure variations during cavity expansion

Variations of cavity pressure during the cavity expansion process are plotted against the cavity radius ratio(a/a0)in Fig.12.Pressure variations are represented by the ratio of the measured cavity pressure(Pa)over the initial cavity pressure(Pa0).It is observed that the cavity pressure variations follow the same pattern irrespective of the adopted contact model.The peak cavity pressure acquired for all three contact models is approximately 5.5 times of the initial cavity pressure(i.e.Pmax=2.75 MPa).However,it is noticed that the cavity pressure obtained from the simulation with Hertz contact model tends to reach the peak at slightly less expansion radius.This finding complies with the soil stress-strain behaviour obtained in the triaxial test(Fig.9),in which Hertz contact model resulted in the peak axial stress at a lower axial strain.It can be readily concluded that the soil response during the cavity expansion is largely dependent on the soil behaviour during the triaxial compression test.Additionally,the simulation process takes approximately 40 h for the simulations adopting linear contact model and rolling resistance contact model on a computer with an Intel i7 processor@3.4 GHz and 16 GB RAM,while about 50 h is required for the analysis adopting Hertz contact model.This may be attributed to the fact that the Hertz contact model employs nonlinear formulations and hence complicated and time consuming computational algorithms are required.It should be noted that for all adopted contact models,strain-softening behaviour is captured,which is likely due to the dilatancy of granular material since the specimen generated in the simulation is defined as medium dense sand.

Fig.11.Cavity expansion model setup.

Fig.12.Cavity pressure variations during cavity expansion.

The dilatancy in macroscopic view is explained as the volume expansion when the specimen is subjected to the shear deformation.However,in discrete element modelling,the dilation can be fundamentally interpreted by the interactions between contacting particles adopting the equilibrating force system.As inspired by Budhu(2008),the dense assembly of sand particles in microscopic view can be simplified to the packing shown in Fig.13.With continuous shearing,particles on the top layer(e.g.particle A)tend to override the particles underneath(e.g.particle B).Hence,the potential slipping plane between two particles can,at any stage,be idealised to an inclined surface with an angle ofα to the horizontal plane,as shown in Fig.14.Referring to Fig.14,the angleα reduces to α′(α′=α+dα,and dα ＜ 0)when the particle A moves to locationA′while a new slipping plane with particle B would be formed.

Fig.13.Illustration of the dense sand assembly in DEM.

If the horizontal force(H)is applied in a slow manner and considering a constant vertical force(V)during the shearing process,quasi-static condition can be satisfied.In this case,equilibrium of forces in both horizontal and vertical directions can be expressed as

Fig.14.Free-body diagram of the particles on the slipping plane.N is the normal force acting perpendicular to the slipping plane;f is the friction force acting parallel to the slipping plane(f=μN(yùn));H and V are the horizontal and vertical forces applied due to external loading,respectively;andαis the dip of slipping plane.

Combining Eqs.(7)and(8)leads to

whereμ =tanφ′and φ′is the friction angle.

Assuming that the surface area of the specimen remains constant during the shearing process,Eq.(9)can be expressed directly in terms of the relationship between the shear and normal stresses:

Referring to Eq.(11),when the friction angle(φ′)is assumed to be constant,the predicted deviatoric stress increases(dτ＞0)as the dilation angle increases(dα ＞ 0).On the contrary,a decrease in dilation angle(dα ＜ 0)would contribute to a reduction in the deviatoric stress(dτ＜ 0),which is known as the strain-softening behaviour,observed after the peak deviatoric stress in macroscopic point of view.

4.2.Radial stress distributions during cavity expansion

A comparison for the distribution of the radial stresses measured at various distances(r/a0)during the cavity expansion is presented in Fig.15.Results have been plotted for a selection of cavity sizes.Fig.15a and b presents the results fora/a0=1.15 and 1.25,respectively.It can be generally concluded that the magnitude of the radial pressure changes significantly with the distance from the internal cavity.Asrincreases,the radial stresses(Pr)reduce more significantly,while further away from the cavity,this reduction rate becomes less pronounced.For example,the normalised radial pressure decreases immensely when the normalised cavity radius(r/a0)is less than 5,while insignificant radial pressure variations are observed whenr/a0＞11.It is projected that the measured radial stresses eventually reach the initial in situ stress(i.e.Pa/Pa0=1)at a large radial distance(i.e.r→+∞）.Furthermore,comparisons were made for radial pressure variations with the radial distance,and it is clearly noticed that there are no pronounced differences between predictions obtained from simulations with three different contact models(see Fig.15).

4.3.Soil radial displacement during cavity expansion

Fig.16 shows the soil movement obtained in radial direction due to the cavity expansion.The results are plotted based on the total displacement of five gauge particles from the beginning of the cavity expansion(i.e.a=a0=0.055 m)till the end(a=af=0.1 m).These five gauge particles are equally spaced,and the locations of these particles are represented using the normalised distancer/a(refer to Fig.11).As expected,particles located closer to the centre of the cavity experience a larger radial displacement.For instance,the gauge particle 1 positioned 0.1 m away from the internal cavity moved 0.03 m in the radial direction when the cavity expanded from 0.055 m to 0.1 m.However,the gauge particle 5 located 0.4 m from the internal cavity experienced a total radial displacement of 0.0075 m during the entire cavity expansion process.Referring to Fig.16,it can be observed that the total radial displacement shows a dramatic decrease with radius when the normalised radius(r/a)is less than 6,while the changes become almost insignificant when the normalised cavity radius is larger than 8.Furthermore,it can be concluded that the calibrated contact models express very similar results in terms of the distribution of the radial soil movement.

4.4.Deviatoric stress and volumetric variations during cavity expansion

Comparisons of the deviatoric stress and the volumetric strain variations with the shear strain during the cavity expansion are presented in Figs.17 and 18.Results obtained are reported based on two prediction spheres situated at radial distances ofrAandrB,as illustrated in Fig.11.As can be clearly observed in Figs.17a and 18a,the deviatoric stresses reveal a continuous increase with the increase of the shear strain,and the soil particles close to the internal cavity experience larger deviatoric stresses compared with those situated further away.Similar findings were also reported by Manassero(1989)and Silvestri(2001)when attempting to interpret the self-boring pressuremeter test in the calibration chamber.Additionally,as depicted in Figs.17b and 18b,the volumetric strain variations show a contraction during the initial stage of the cavity expansion and then a dilation when the cavity is expanded further,confirming the medium dense state of the sandy soil simulated.As expected,the soil close to the centre of the cavity(e.g.rA/af=1.25)experiences more volume changes while the variations of the volumetric strain predicted become less pronounced with the increase in radial distance(e.g.rB/af=2.5).It is projected that the predicted volumetric strains would eventually approach zero when the radial distance is extremely large(i.e.r→+∞).Furthermore,it can be concluded that the numerical results obtained for three different yet calibrated contact models are in a good agreement.

4.5.Dilatancy variations during the cavity expansion

Dilation angle is the measure of the change in the volumetric strain with respect to the change of the shear strain(Budhu,2008).Hence,it is possible to use the current model and results to further investigate the dilatancy variations.Fig.19 depicts the variations of dilatancy with shear strain based on the prediction spheres A and B(as shown in Fig.11).It is observed that dilatancy is negative at the initial stage of the cavity expansion and then positive when the cavity expands further,confirming a transient contraction and subsequent continuous expansion at both predictions spheres.In addition,referring to Fig.19a,the dilatancy shows a persistent reduction after the peak,and it is projected that the dilatancy will eventually approach zero when the specimen reaches the critical state.

Fig.15.Radial stress distributions at various stages:(a)a/a0=1.15 and(b)a/a0=1.25 during cavity expansion.

Fig.16.Radial displacement of five equally spaced gauge particles during the cavity expansion.

Fig.17.(a)Deviatoric stress-shear strain relationship and(b)Volumetric strain-shear strain relationship at predication sphere A.

4.6.Contact forces and displacement contour

Fig.20 shows an example of contact forces between particles and total displacement contours adopting rolling resistance contact model.As the contours obtained from three different contact models are very similar,only one case is presented.As expected,the contact forces close to the internal cavity are larger and change more dramatically with radius.Additionally,referring to Fig.20a,the maximum soil displacement of approximately 0.045 m is observed at the internal cavity,while the displacement reduces gradually with the increasing distance from the centre of the cavity.It should be highlighted that the linear contact model may not be able to simulate the behaviour of particles that are highly angular.Instead,these angular particles may be more accurately simulated using rolling resistance contact model.Consequently,it is strongly suggested that the choice of the most appropriate contact model should be deliberated based on the properties and characteristics of the materials to be studied.

5.Further discussion and limitations

Fig.18.(a)Deviatoric stress-shear strain relationship and(b)Volumetric strain-shear strain relationship at predication sphere B.

This paper has presented an investigation on the influence of particle contact models on soil response of sandy soil during cavity expansion.According to the unified soil classification system(USCS),the specimen used in this study was categorised as poorly graded sand.The calibration exercise was conducted using the same particle size distribution reported by Cornforth(1964)as shown in Fig.5.Hence,the calibrated contact parameters could mimic the shear behaviour and obviously correspond to the poorly graded sand adopted in the experiment.Obviously,the micromechanical characteristics including the particle contact properties greatly depend on the texture of the surface of particles,particle roundness and sphericity,mineralogy of the parent material,and geological processes involved(Horn and Deere,1962;Herle and Gudehus,1999;Bareither et al.,2008).For instance,Bareither et al.(2008)indicated that sand with lower shear strength is usually derived from weathering of underlying sandstones,and tends to have medium to fine and well-rounded yet poorly graded particles.However,higher shear strength is expected from the sand particles formed as a result of recent glacial activities leading to coarser,more angular and well graded particles.Obviously,if the calibration exercise is conducted for another type of soil that is well graded,then the micromechanical properties would be different,which in turn will influence the soil response during the cavity expansion process.For instance,the shear strength of the well graded sand is generally larger due to its relatively stable particle structure resisting the sliding and the collapse of the interlocking(Enomoto et al.,2015).In addition,Kirkpatrick(1965)and Leslie(1969)pointed out that the increase in the coefficient of uniformity(Cu)can enhance the peak shear strength of granular material significantly.

As mentioned earlier,this paper has employed the discrete element simulations to investigate the influence of the particle contact models on response of poorly graded sand during cavity expansion.The contact models utilised contain simplifications and assumptions that inevitably induce disparities between the numerical predictions and the realistic soil responses.Hence,it is noteworthy to highlight the limitations of this study as summarised below:

(1)Soil particle crushing and cracking have not been captured in this study.Therefore,the ploughing behaviour causing permanent plastic deformation at the contact under high normal loads has not been captured(Senetakis et al.,2013).According to Sandeep and Senetakis(2018),the ploughing behaviour is responsible for the decrease in interparticle friction coefficient during shearing,and this ploughing might also trigger the increase in particle friction coefficient during unloading process.Hence,particle crushing can contribute to strength and stiffness degradation.Similarly findings were also report by Miura et al.(1984).Hence,it is recommended for future studies to capture particle crushing during cavity expansion process to assess the extent to which the particle crushing can influence stress and strain predictions.Indeed,the particle crushing can be considered by means of the user defined contact models through a pre-defined particle crushing criterion.When the contact forces(or stresses)between particles exceed the threshold,the particles will be crushed and replaced by a certain number of smaller particles.

Fig.19.Variations of dilatancy during cavity expansion at(a)prediction sphere A and(b)prediction sphere B.

(2)The viscous behaviours of particles and interfaces have not been included.Indeed,the material behaviour can be influenced by both compression and shear creep causing stress relaxation or further deformation(Pham Van Bang et al.,2007;Wang et al.,2008).For instance,Kuwano and Jardine(2002)indicated that the shear creep deformation observed within 2 h of the end of anisotropic consolidation contributed to 30%-80%of the “primary”deformation observed during the fully drained loading stages.

(3)In this study,ideal spherical particles have been used in both the calibration exercise and the cavity expansion simulation,while angularity of particles and sharp corners may influence the results(Zhao et al.,2015;Fu et al.,2017;Sandeep and Senetakis,2018).Zhao et al.(2015)reported that the interlocking of the particles is more dramatic in a high angular assembly,resulting in high shear strength and dilatancy.In addition,the non-uniform distribution of normal contact forces near the shear plane is more significant in high angular assembly during shearing,conifrming that few particles share the majority of the shear force.

All of these limitations discussed above have certain effects on the numerical predictions.Authors are working on follow-up papers addressing the above limitations.

6.Conclusions

Fig.20.(a)Particle displacement contour(m)and(b)particle contact force contour(N)for rolling resistance contact model.

Cavity expansion theory serves as one of the most effective techniques in geotechnical engineering to interpret the soil behaviour in many practical applications such as driven piles and in situ tests.In this study,the effects of three different contact constitutive models,i.e.linear,rolling resistance and Hertz contact models,on the sandy soil response during the cylindrical cavity expansion have been investigated adopting discrete element analysis.The simulations of the cavity expansion with different contact models have been calibrated using the experimental results obtained from the triaxial test in plane strain condition on medium dense sand.The calibration process calibrating the microscopic contact parameters including the methodology adopted and parametric studies carried out is expounded.In addition,micromechanical formulations based on the contact forces have been developed to explain the soil dilation and strain softening in a microscopic point of view.The numerical predictions on cavity expansion indicate that these calibrated contact models produced similar results in terms of cavity pressure variations,radial stress distributions,deviatoric stress and volumetric strain variations as well as soil radial movement,confirming that the soil response during the cavity expansion is highly correlated to its behaviour during the triaxial test used for the model calibration.Additionally,strain-softening behaviour due to the dilatancy of medium dense specimen has been captured during the cavity expansion,and the volumetric strain predicted shows an initial slight contraction followed by dilation at larger shear strains.However,it should be highlighted that if the soil particles are highly angular(e.g.ballast particles),modelling adopting simplified linear contact model may result in inaccurate predictions.Furthermore,according to the results reported in this study,at a radial distance beyond 11aand 8a,soil displacement and radial stress variations would be insignificant,respectively.Thus,the effects of soil displacement due to installation of displacement-based inclusions(such as driven piles and vertical drains used in ground improvement)on nearby structures should be considered carefully by practicing engineers.Moreover,the limitations of this study are elaborated including the soil crushing,viscous behaviour as well as the particle angularity.All of those effects will be thoroughly considered in the future studies.

Conflicts of interest

We 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 influenced its outcome.

List of symbols

aCurrent cavity radius(m)

a0Initial cavity radius(m)

afFinal cavity radius(m)

dDistance between the centres of two contacting particles(m)

E*Effective Young’s modulus of particles(Pa)

EYoung’s modulus(Pa)

FnNormal force at the contact interface between contacting particles(N)

FsShear force at the contact interface between contacting particles(N)

Shear force at the beginning of the timestep(N)

Hertz normal force(N)

Hertz shear force(N)

Hertz shear force at the beginning of the timestep(N)

Surface gap between contacing particles(m)

GEffective shear modulus(Pa)

HHorizontal force appied to the particle(N)

knNormal stiffness of linear springs(N/m)

ksShear stiffness of linear springs(N/m)

krRolling stiffness of linear springs(N/m)

k*Stiffness ratio between normal stiffness and shear stiffness

LDistance between the centres of contacting particles

MrrRolling resistance moment(N m)

NNormal force applied to the particle(N)

paCavity pressure(kPa)

Pa0Initial cavity pressure(kPa)

pmaxPeak cavity pressure(kPa)

prRadial stress(kPa)

R1Radius of particle 1(m)

R2Radius of particle 2(m)

Particle effective radius in Hertz contact model(m)

rARadius of the prediction sphere A(m)

rBRadial of the prediction sphere B(m)

rRadial distance(m)

VVertical force appied to the particle(N)

αHHertz contact model exponent

α Dilation angle in DEM simulation(°)

δsShear displacement increment in a timestep(m)

ν Poisson’s ratio

μrRolling resistance coefficient

σr0Initial radial stress(kPa)

σXInitial stress inXdirection(kPa)

σYInitial stress inYdirection(kPa)

φrRelative rotation increment between contacting particles in a timestep(°)

φ′Friction angle(°)

τ Shear stress(Pa)