Numerical investigation of diffuse instability in sandy soil using discrete element method under proportional strain path loading

The Loose sandy soil can manifest static liquefaction phenomenon under undrained condition, in which the onset of instability is within the failure line, and there is no obvious shear band. This type of failure mode, very different from localized instability that occurs in dense sand, is called the diffuse instability. In this paper, a series of proportional strain tests and fully drained tests under different initial void ratio were simulated using the discrete element method. The influence of strain increment ratio and the initial void ratio affecting the instability of sandy soil were discussed in detail. The development mechanism of pore water pressure in proportional strain tests was analyzed by comparing with the volumetric curve of fully drained test. Finally, a unified mechanism of diffuse instability of sandy soil in proportional strain tests was explained. Numerical results indicate that the strain increment ratio and the initial void ratio work together affecting the instability of specimen. The occurrence of diffuse instability is the result of effective stress reduction due to the development of pore water pressure, which depends on the difference of volumetric strain between fully drained tests and proportional strain tests. The increment of pore water pressure is determined by the difference of strain increment ratio, which can be used as an index to reflect the liquefaction potential of sandy soil.


INTRODUCTION
The instability of soil can be considered as a stress state that is reached before the failure line from mechanical perspective.Any small stress increment under this state will lead to a large plastic strain.Soil instability can be roughly divided into two types: internal instability and external instability (Goddard, 2003).Internal instability refers to the loss of the inherent strength of the soil, and the external instability means that the interaction between the soil and the surrounding environment is beyond the bearing limit.The internal instability including localized instability with shear bands and diffuse instability without shear bands.The differences between these two types of instability are discussed in detail by Lade, 1990 andLancelot et al. (2004).This paper mainly discussed the diffuse instability.The traditional idea holds that soil will not unstable before the effective stress reaches the failure line.But Castro (1969) found that loose sand may manifest liquefaction phenomenon under undrained condition without obvious shear band and the onset of instability is within the failure line.Instability of liquefaction is one of the major reasons which results in the failure of earth structures such as dam."Liquefaction" is used for the first time by Hazen (1920), describing the failure of the Calaveras dam's failure.Casagrande defined the critical void ratio (CVR) concept in his early paper.Seed and Lee (1966) proposed the concept of "initial liquefaction" and Casagrande (1975) proposed the concept "steady state strength", which is used to estimate sand liquefaction failure.Later several concepts such as "steady state line", "flow structure" were proposed to describe the liquefied deformation mechanism by several researchers (Poulos, 1982;Poulos et al, 1985;Prevost,1985;Lade and Pradel, 1990;Bazier and Dobry,1995).More recently, other researchers have carried out experimental research of the instability of sand under undrained conditions (Laouafa and Darve, 2002;Daouadji et al, 2009;Daouadji et al, 2011) and developed corresponding constitutive models to explain and predict this type of instability (Mr óz and Boukpeti, 2003;Desai et al, 2005;Rahman et al, 2014;Najma and Latifi, 2017;Chandrakant and Desai, 2016;Liu et al, 2016).These studies provide an important theoretical basis for the analysis of diffuse instability.
In recent years, some studies have shown that diffuse instability can occurs under fully drained conditions (Eckersley, 1990, Olson et al., 2000).Therefore, the research on the instability mechanism of soil under different drainage conditions has attracted the attention of many scholars (Sasitharan et al., 1993;Anderson and Riemer, 1995;Chu et al., 2012;Skopek et al., 2011;Dong et al., 2016;Alipour and Lashkari, 2018).The drainage conditions have great influence on the mechanical behaviour of soils.The actual in situ drainage conditions of the soil include five situations.i.e., fully undrained, fully drained, partial drainage, excessive drainage, expansive drainage.However, most of the experimental studies were carried out under fully drained or fully undrained conditions.To investigate the diffuse instability of soils under more general conditions, several investigators have used proportional strain tests by imposing a constant strain ratio Chu et al., 1992Chu et al., , 1993;;Vaid and Eliadorani, 1998;Lancelot et al., 2004;Jrad et al., 2012;Nicot et al., 2013Nicot et al., , 2015;;Mital and Andrade, 2016;Lashkari and Yaghtin, 2018;Daouadji et al., 2013).The strain ratio α can be considered as a mathematical parameter for controlling the drainage conditions.
Limited by equipment, it is difficult to achieve accurate control of the strain path of the specimen during the experimental study.But numerical simulation can avoid these problems and several authors have used the discrete element method to study the behavior of granular soils under different loading conditions (Gu et al., 2014a;Gu et al., 2014b;Guo and Zhao, 2016;Zhou et al., 2016).Therefore, in this paper, the discrete element method was used to simulate a series of proportional strain tests of sandy soil with different initial void ratio e and different strain ratio α .The influence of drainage conditions and initial void ratio on the diffuse instability of sand was analyzed in detail.Hill's second-order work criterion was adopted to estimate the instability of soil (Hill, 1958).The stable / unstable regional map in e α − plane was obtained.Furthermore, a unified mechanism of diffuse instability of sandy soil in proportional strain tests is explored by analyzing the development of pore water pressure.

NUMERICAL MODEL
The numerical model was established by the three-dimensional discrete element program PFC 3D , as shown in Figure 1.Ni et al. (2000) found that increasing the number of particles beyond 15,000 in a 3D DEM simulation raised the computational expense with a relatively small effect on the model results.In addition, according to Wang and Li (2014), who compared the results obtained with 1000, 5000, 10,000, and 20,000 particles, there is no significant difference between the results obtained with 5000, 10,000, and 20,000 particles.In this study, we choose approximate 5000 particles to conduct simulation, which changes with the initial void ratio.The minimum and maximum particle diameters are 1.0 and 1.8mm, respectively, with a mean particle diameter of 1.4mm.
Different sizes of samples were used in previous DEM simulations or laboratory test but the ratio of height to diameter is usually taken as about 2.0.For example, the DEM specimen used by Tang (Tang et al., 2016) has a cylindrical structure with the height of 80.5 mm and the diameter of 38 mm; The triaxial specimen tested by Vaid and Eliadorani (1998) is 126 mm height × 63 mm diameter; The triaxial specimen used by Sivathayalan and Logeswaran (2008), is 125 mm height × 63 mm diameter.Therefore, we choose 2.0 as the ratio of height to diameter, to build numerical model.Then the appropriate diameter of 16 mm can be obtained ( ).We choose 20mm as the sample diameter and the height is 40mm considering the influence of initial porosity of the assembly.

CALIBRATION OF PARAMETERS AND VERIFICATION WITH TYPICAL EXPERIMENTAL RESUTLS
The linear contact model was adopted here and three micromechanical parameters, i.e, normal contact stiffness k n , contact stiffness ratio k n/ k s and inter-particle friction µ are to be determined.Usually, the normal contact stiffness k n can be determined by the initial stiffness of stress-strain curve, the inter-particle friction is related to the peak strength of samples, and the stiffness ratio is affected by the macro Poisson's ratio.Here, we took triaxial compression test of Fraser River sand under fully undrained condition with confining pressure of 200 kPa (Sivathayalan and Logeswaran, 2007) to calibrate the parameter.
Latin American Journal of Solids and Structures, 2018, 15(11), e134 3/15 We tried five k n values to conduct a series of simulations with empirical values of =0.5 µ and k n/ k s =1.0.The simulation results were compared with the experimental data, as shown in Figure 2. The initial slopes of the stress-strain curve (initial stiffness) were compared with the simulating results from DEM.One can see from Figure 2 that the numerical result of k n =1.5×10 7 N/m is best fit for the initial stiffness of stress-strain curve, which is then be chosen as the simulation parameters.
After that, we tried four friction coefficient values to conduct simulations with k n =1.5e7 and k n/ k s =1.0 .The simulation results were compared with the experimental data again, as shown in Figure 3.One can see that the numerical result of =0.75 µ is most consistent with the peak strength of sample.Finally, we adjusted the values of contact stiffness ratio k n/ k s to fit experimental results, as shown in Figure 4, and the value of k n/ k s =1.5 is then be determined.

Figure 2: Calibration of model parameters
Simulation was then carried out using the selected parameters calibrated from fully undrained test.A typical experiment under partial drainage condition ( α =0.1) were selected to assess preliminarily the ability of the developed DEM model and acceptability of chosen parameters.Figure 3 shows the comparison of simulation results with experimental results.One can see from Figure 3 the simulated macro-response from DEM has a good agreement with the results of laboratory tests.
One may discuss a little higher contact stiffness ratio are selected in simulations, however, the predicted macroresponse from DEM is roughly consistent with the results of laboratory tests.In this sense, the used parameters are acceptable.
Numerical investigation of diffuse instability in sandy soil using discrete element method under proportional strain path loading Latin American Journal of Solids and Structures, 2018, 15(11), e134 4/15

NUMERICAL SIMULATION SCHEME
The loading program is described by the following conditions for the imposed proportional strain test: where 1 dε , 3 dε is the increment of axial strain and radial strain respectively; v dε is the increment of volumetric strain.Using the common sign convention used in geomechanics, axial compression and volume outflow are considered as positive.Therefore, a positive α value corresponds to volumetric contraction and a negative α value corresponds to volumetric expansion.Depending on the value and sign of the parameter α the test can be classified into three types: It should be noticed that 1 α = corresponds to a one-dimensional test.Although the in-situ conditions are clearly dif- ferent, a constant strain ratio has been chosen for the whole test to study the influence of drainage condition on the instability of sands.
A total of 14 samples with different initial void ratio were selected.The void ratio range is: 0.618~0.98.The minimum and maximum strain increment ratio α are -0.5 and 0.5, respectively, with interval value of 0.05.Therefore 21× 14 proportional strain tests were simulated in all.The simulation scheme is shown in table 1.
It should be pointed out that this paper is aimed not to fit experimental data, but to investigate the micro-and macro-mechanism of diffuse instability in granular material.So, the simulation scheme ( α value and initial void ratio) is not coincident with the Fraser River sand test under generalized drainage boundary conditions (Sivathayalan and Logeswaran, 2007).

DIFFUSE INSTABILITY IN UNDRAINED TESTS ( 0 α = )
One of the key points to study the diffuse instability is how to determine the instability point.According to Chang et al. ( 2010), the instability conditions can be divided into two modes: Mode 1: If the determinant of stiffness matrix is zero.In this case, the sample bifurcates with the occurrence of unlimited displacement rates, which is considered as the onset of localization.
Mode 2: Loss of positiveness of second-order work.In this case, the second-order work could become zero while the determinant of stiffness matrix remains positive.
The diffuse instability belongs to Mode 2 failure.Hill's second-order work criterion has been widely applied to judge the stability of granular materials (Sibille et al., 2015;Hadda et al., 2015;Darve and Laouafa, 2015;Darve et al., 2004;Nicot et al., 2010;Daouadji et al., 2010;Sawicki and Świdziński, 2010).which states that a material, progressing from one stress state to another, is stable if the second-order work is strictly positive.For triaxial conditions, it can be written as following: Where ' 1 σ = axial effective stress; ' 3 σ = radial effective stress; 1 ε = axial strain; 3 ε = radial strain; q = deviatoric stress; and v ε = volumetric strain.The results of undrained tests with various initial void ratios are presented in Fig. 4. Corresponding test number is K2, K3, K4, K6, K8.As indicated in Fig. 4(a), For dense sand (e=0.658,e=0.715, e=0.748), the deviatoric stress continuously increases, and no maximum peak strength can be observed.For medium-dense sand (e=0.767), the deviatoric stress first reaches a peak, then decreases slightly but finally the strength is recovered.For loose sand (e=0.818),The peak value of deviatoric stress is reached rapidly, which is followed by a rapid decrease of the deviatoric stress down to a minimum strength.The minimum strength can be almost zero for loose sands, which represents the phenomenon called static liquefaction.
during the undrained test 0 v dε = , the second-order work is therefore given by: As the axial strain is increased, equation ( 5) vanishes only if 0 dq = .The deviatoric stress q becomes the control parameter and the peak value of deviatoric stress corresponds to the instability point.The second-order work is computed by equation ( 5) and plotted in Fig. 5 against the axial strain.
It is seen that, for the three dense samples of e=0.658, e=0.715 and e=0.748, the value of second-order work always keeps positive, so it is always in a stable state.The second-order work of the e=0.767 sample drops to a negative value near the strain of 1.4% and then rises to a positive value near the strain of 7.2%.This phenomenon is called partial failure and the stress state of 1 ε =7.2% is called the phase transition point.For the loose sample of e=0.818, the value of second-order work is always negative after the instability point, so there is no phase transition point, and it is always in unstable state after the instability point.From Fig. 4(b) we can see that the instability lines are within the critical state line, which indicate that the type of instability occurred in these three samples belongs to the diffuse instability.
In summary, larger initial void ratio sample is more easily to manifest instability, which is in accordance with the experimental research.With the increase of the initial void ratio, the diffuse instability can be categorized as two types: (1) Partial failure with the transition point; (2) Always in unstable state after the instability point.

DIFFUSE INSTABILITY VARIED WITH THE INITIAL VOID RATIO
To study the influence of initial void ratio on soil instability in other drainage conditions ( 0 α > or 0 α < ).We  respectively.It is clearly shown that, with the increase of the void ratio, the stress-strain curve gradually decreases, and the larger the void ratio is, the more prone to instability the sample becomes.Obviously, this rule is established not only in undrained condition, but also in other drainage conditions.What's different is that: when 0.5 α = all the five samples have properties of strain hardening, including e=0.767 and e=0.818 samples, which showed strain softening phenomenon in undrained tests.By contrast, when 0.5 α = − the dense samples of e=0.715 and e=0.748 showed strain softening phenomenon, which showed strain hardening phenomenon in undrained tests, and the strain hardening phenomenon for e=0.658 sample is very weak.Obviously, the value of α , representing the drainage condition, affects the stability of the soil in another dimension, which will be discussed in the following section.condition, but it can be foreseen that the sample will eventually lose its stability with the further reduction of α value.Obviously, under the same confining pressure, whatever the void ratio is, the sample with the smaller α value is more prone to instability.Therefore, the decrease of α value and the increase of the initial void ratio have the same effect on soil instability, which is driving soil from stable to unstable.
Latin American Journal of Solids and Structures, 2018, 15(11), e134 8/15 α ≠ , the deviatoric stress q will not be the control parameter any more.Substituting equation (3) into equation ( 4) yields the following equation: Therefore, the control parameter changes from deviatoric stress q to ' 3 q ασ + .However, 0 dq < is still a necessary condition for the onset of diffuse instability, for the convenience of discussion, we use 0 dq < as collapse criteria for all proportional tests.Since the strength will finally be recovered in partial failure, only the second type of diffuse instability are discussed in this paper.
The stable/unstable regional map of all proportional tests in e α − plane is drawn as shown in Figure 9.The graph is divided into two parts by undrained tests: the upper part is volume contraction area and the lower part is volume dilatation area.We take the simulation results from samples of e=0.767 as example to explain the stable/unstable regional map.There is no phase transition point after instability point when α = -0.4,-0.2 and 0 for sample of e=0.767.Therefore, the coordinate pairs of (0.767, -0.4), (0.767, -0.2) and (0.767, 0) related to unstable responses are contained in unstable region.However, the deviatoric stress continuously increases and no peak strength can be observed when α = 0.2 and 0.4.Therefore, the coordinate pairs of (0.767, 0.2), (0.767, 0.4) related to stable responses will be contained in stable region.
One can see from Figure 9 that some dense samples keep stable under undrained condition but become unstable in the volume dilatation area.On the contrary, some loose samples are unstable under undrained condition, but remain stable in the volume contraction area.The left side of the division line is close to negative infinity and the right side goes to positive infinity, which represents very dense sand (rock) and extremely loose sand(mud).This is only an idealization, which already has no properties of sand.
In general, the value of α and initial void ratio affect the instability of sand in two aspects under the same confining pressure.Each proportional test is a specific combination of α and e, and whatever the void ratio is, there is always a α value to make it unstable.Therefore, no matter what density of sand is, it is possible to lose stability if the drainage condition is suitable, and the looser sand is more prone to instability regardless of the drainage condition.

PORE WATER PRESSURE FORMULA
The calculation of pore water pressure in the discrete element model adopts the equivalent principle, that is, the total confining pressure is assumed to be constant.According to the effective stress theory, the sum of the increment of pore water pressure and the increment of effective confining pressure equals 0. This implies: = increment of effective confining pressure; u ∆ = increment of pore water pressure; 3 σ ∆ = increment of total confining pressure.
According to equation ( 7), the calculation formula of pore water pressure can now be expressed as: Where 0 u = initial pore water pressure (equals 0); ' 30 σ = initial effective confining pressure.According to the formula (8), the pore water pressure development curves of each sample in the proportional strain tests can be obtained.

DEVELOPMENT OF PORE WATER PRESSURE
The essential cause of instability in sandy soil is the decrease of effective stress.In the proportional strain tests, the decrease of effective stress is caused by the increase of pore water pressure, and the development of pore water pressure is closely related to the development of contraction or dilatation.
It is well known that the pore water pressure always equals zero in fully drained tests, so the volumetric strain of fully drained tests can be taken as the reference value to analyze the development of pore water pressure in proportional strain tests.The difference between the volumetric strain value of fully drained tests and proportional strain tests is defined as: ), which implies that point C correspond to the zero point of pore pressure.In other words, when the accumulative volumetric strain of the proportional strain test is equal to the accumulative volumetric strain of fully drained test, the accumulative pore water pressure will be equal to zero.Taking the BC line as a boundary we can see that  Where d α = strain increment ratio in fully drained test; p α = strain increment ratio in proportional strain test.According to the previous analysis, the value of pore water pressure in proportional strain test depends on the value of 12/15 drained test will not move, the volumetric strain curve of proportional strain test will decrease, resulting in the increase of v ε ∆ , so the pore water pressure will increase and the effective stress will decrease.As can be seen in Fig. 12 (b), the pore water pressure rises with the decrease of α value.It is the same thing if we fix α value and raise the initial void ratio e.
To sum up, the initial void e and the strain increment ratio α influence the development of dilatancy of the soil together, and then affect the development of the excess pore water pressure and the effective stress, and ultimately affect the stability of the soil.

CONCLUSIONS
In this paper, a series of proportional strain tests under different initial void ratio and different strain increment ratio were simulated using the discrete element method.The influence mechanism of the strain increment ratio and the initial void ratio on the instability of sandy soil was analyzed.By comparing with the volumetric strain curves of the fully drained tests, the development mechanism of the pore water pressure in proportional strain tests were discussed in detail.The preliminary conclusions of the study are as follows: (1) In the proportional strain test, the strain increment ratio α and the initial void ratio e work together to affect the instability of specimen.
No matter what the initial density of the sample is, it may be unstable if the drainage condition is suitable.Meanwhile, under the same drainage conditions, loose samples are more prone to instability.At a given confining pressure, there is a destabilizing boundary line in (3) The actual in situ drainage conditions will be neither fully drained nor fully undrained, and the imposed constant increment strain ratio is also impossible.However, the value of

Figure 3 :
Figure 3: Results verification with the experimental results.under partial drainage condition α=0.1
choose a series of proportional strain tests under different α value with the same 5 initial void ratio samples to analyze.For the sake of clarity, only the results of 0analyzed here, corresponding test numbers are A2, A3, A4, A6, A8 and U2, U3, U4, U6, U8.The simulation results of proportional tests under other α values are similar.

Figure 8 :
Figure 8: Simulated results under different α value for samples with the same void ratio

Figure 9 :
Figure 9: Stable / unstable region obtained from simulations results strain of fully drained test; p v ε = volumetric strain of proportional strain test.We chose the simulation results of e=0.838 sample with 0.4 α = as an example to analyze, and the results of other α values are similar.Fig. 10 shows correlation between pore-water pressure curve and volumetric strain curve for e=0.838 sample.The point A represents peak value of pore water pressure.The point C is the intersection point of volumetric strain curve for fully drained test and proportional strain test ( The point B is the zero point of pore-water pressure curve.Before the point B 0 u > , after the point B 0 u < .As shown in Figure 10 (a), point B and point C are approximately at the same strain level ( 0 Obviously, the value of the accumulative pore water pressure of proportional strain test depends on the value of v ε ∆ .To analyze the results of other α values also have the same rule, here are not to be discussed in detail due to limited space.It can be noticed from Fig. 10 (b) that the value of v ε ∆ reaches maximum at point D, which indicate that the strain increment ratio of two tests are equal (=0.5).The strain increment ratio (equation 4) essentially implies the development of the dilatancy during the loading of the soil.To further analyze the relationship between pore water pressure and the dilatancy of soil, we define α ∆ as follows:

Figure 10 :
Figure 10: Correlation between pore-water pressure and volumetric strain curve for sample with e=0.838.

Figure 11 :
Figure 11: Correlation between pore-water pressure and strain increment ratio during loading for sample with e=0.838.

Figure 12 :
Figure 12: (a) Volumetric strain curve of proportional strain tests with different α value; (b) Pore water pressure curve for proportional strain tests with different α value (e=0.838sample).
v ε ∆ , and the value of v ε ∆ depends on the relative position of the volumetric strain curve under the fully drained test and proportional strain test.If we fix the initial void ratio e and reduce the α value, the volumetric strain curve of fully Latin American Journal of Solids and Structures, 2018, 15(11), e134 The occurrence of diffuse instability is the result of effective stress reduction due to the development of pore water pressure.The accumulated pore water pressure depends on the value ofv ε ∆and the increment of pore water pressure is determined by the value of α ∆ .
be used as an index to reflect the liquefaction potential of sandy soil.

Table 1 :
Numerical simulation scheme