A constitutive model for cyclic densification of coarse-grained soil filler for the high-speed railway subgrade considering particle breakage

Purpose – This study aims to propose a semiempirical and semitheoretical cyclic compaction constitutive model of coarse-grained soil filler for the high-speed railway (HSR) subgrade under cyclic load. Design/methodology/approach – According to the basic framework of critical state soil mechanics and in view of the characteristics of the coarse-grained soil filler for the HSR subgrade to bear the train vibration load repeatedly for a long time, the hyperbolic empirical relationship between particle breakage and plastic work wasderived.Consideringtheinfluenceofcyclicvibrationtimeandstressratio,theparticlebreakagecorrection functionofcoarse-grainedsoilfillerfortheHSRsubgradeundercyclicloadwasproposed.Accordingtotheclassicaltheoryofplasticmechanics,theshearingdilatationequationofthecoarse-grainedsoilfillerforthe HSRsubgradeconsideringparticlebreakagewasmodifiedandobtained.Asemiempiricalandsemitheoreticalcycliccompactionconstitutivemodelofcoarse-grainedsoilfillerfortheHSRsubgradeundercyclicloadwas furtherestablished.ThebackwardEulermethodwasusedtodiscretizetheconstitutiveequation,builda numerical algorithm of “ elastic prediction and plastic modification ” and make a secondary development of the program to solve the cyclic compaction model. Findings – Through the comparison with the result of laboratory triaxial test under the cyclic loading of coarse-grained soil filler for the HSR subgrade, the accuracy and applicability of the cyclic compaction model wereverified.Resultsshowthatthemodelcanaccuratelypredictthecumulativedeformationcharacteristicsof coarse-grainedsoilfillerfortheHSRsubgradeunderthetrainvibrationloadingrepeatedlyforalongtime.Itconsiderstheeffectsofparticlebreakageandstressratio,whichcanbeusedtocalculateandanalyzethestress anddeformationevolutionlawofthesubgradestructureforHSR.


Introduction
In recent years, the construction of high-speed railway (HSR) projects in China has ushered in vigorous development. By the end of 2019, China's HSR total mileage has registered 35,000 km and accounts for more than 60% of the HSR length throughout the world, ranking first worldwide. According to the Medium-and Long-term Railway Network Planning (Mid-term Adjustment in 2018), China will continue to build HSR on a certain scale in the future. However, the development of HSR in China has stepped into a long-term safe and stable operation stage from large-scale construction. The stable railway subgrade structural layer serves as a vital guarantee for the safe and smooth operation of HSR (Cai, Ye, Yan, Wei, & Yao, 2020). Figure 1 shows a typical profile of ballastless track subgrade structure of HSR in China, which is mainly divided into surface layer of subgrade bed, bottom layer of subgrade bed and embankment body below subgrade bed (Tian, 2014b;Ye, Cai, Chen, Yang, & Chen, 2020a). Generally, the surface layer of subgrade bed is filled with graded crushed stones, featuring good engineering characteristics. The bottom layer of subgrade bed is filled with Groups A and B fillers, and the thickness of the bottom layer is 2.3 m, which generally accounts for more than half of the total height of the subgrade. Therefore, it is an important factor causing the deformation of the subgrade structure. In order to quantitatively describe the deformation characteristics of this subgrade structure under train load and evaluate the safety of subgrade operation in the future decades, it is necessary to establish a numerical calculation model, and choosing a suitable dynamic constitutive model is an important guarantee for the accuracy and reliability of numerical simulation results.
For laboratory tests, with the rapid development of large-scale indoor dynamic triaxial test equipment, many scholars have carried out a lot of research work on the dynamic characteristics of coarse-grained subgrade filler for HSR under cyclic loading (Hu, Wang, & Zhang, 2001;Tian, 2014a;Ning, 2018;Liu, Xie, Gao, & Zhan, 2013;Zhou et al., 2016;Mei, Leng, Liu, Nie, & Xu, 2017;Zhang, Zhao, & Zhai, 2017;Bian, Jiang, Shen, & Chen, 2011;Wang, 2001). Physical and mechanical parameters related to the state of subgrade bed filler are tested and described in all aspects in these studies. The effects of different dynamic stress amplitude, confining pressure and water content on the cumulative deformation of coarsegrained soil samples are studied. The evolution law of cumulative deformation with the cyclic loading vibration time is revealed, and the empirical expressions of the relationships between cumulative deformation and the cyclic vibration time, dynamic stress amplitude and confining pressure are established. However, these empirical models established based on the test results can only capture the cumulative deformation generated under each cyclic loading and cannot simulate the cyclic loading and unloading process during actual dynamic effect.
In terms of constitutive model, several dynamic constitutive models of soil have been developed so far. Equivalent linear model, a typical representative of viscoelastic model, was first proposed by Seed in 1968, which approximately considered the dynamic nonlinear property of soil by equivalent linearization method (Zhang & Ling, 2016). Hardin-Drnevich model (Hardin & Drnevich, 1972), Ramberg-Osgood model (Ramberg & Osgood, 1943), bilinear model and some combined curve models all belong to different forms of equivalent linearization models. Mroz first put forward the theory of plastic modulus field in 1967, which represented the beginning of the study on dynamic elastic-plastic model of soil (Morz, 1967). Iwan, Provest, Zienkie-wicz, Yang, Elgamal and others have established their own multi-yield surface models (Iwan, 1967;Provest, 1978;Zienkiewicz & Morz, 1984;Yang, Elgamal, & Parra, 2003). In China, Wang Jianhua, Xu Gancheng, Chen Shengshui and Zhuang Haiyang have established their own multi-yield surface models (Wang & Yao, 1994;Xu, Xie, & Zheng, 1995;Qiu & Zhang, 2011;Zhuang, Chen, & Zhu, 2006). Then, Dafalias and Popov (1975) proposed a more simplified boundary surface model based on Cambridge model. Tabbaa and Wood (1989) promoted the modified Cambridge model to a two-sided dynamic hardening model capable of describing the hysteretic response characteristics of clay under cyclic loading. However, the above-mentioned dynamic constitutive models generally have two problems: first, it can only simulate the situation where the number of cyclic loading is small (e.g. seismic load); second, grain crushing effect under cyclic loading is not considered.
Based on the constitutive theory of Indraratna, Thakur, Vinod and Salim (2012), the theory of critical state soil mechanics and classical plastic mechanics, in combination with the cyclic loading triaxial test, a constitutive model study on cyclic densification of coarse-grained soil filler in the HSR subgrade considering grain crushing is considered in this paper.
2. Deformation characteristics of coarse-grained soil filler for the HSR subgrade 2.1 Compactness Under cyclic loading, the coarse-grained soil filler of HSR subgrade generally features nonlinearity, stress dependence, etc. Even when the stress level is very close to the static strength, the coarse-grained soil filler still shows strong compactness (Indraratna et al., 2012;Ye et al., 2020b). Under the initial loading conditions, the axial strain develops rapidly and partially rebounds during unloading, and each subsequent cyclic loading produces a continuous plastic strain increment. However, the amplitude of the plastic strain increment gradually decreases as the vibration load time increases. This is because the stiffness of coarse-grained soil filler becomes larger due to the previous stress under cyclic loading, and then the cumulative strain decreases under subsequent load. Under cyclic loading, the development of cumulative deformation of coarse-grained fillers of HSR subgrade is a gradual development process, and each single load will produce a small strain increment. A large number of dynamic triaxial tests show that the strain generated under the initial loading conditions is large, but the subsequent strain increments are continuously reduced under the premise that the stress amplitude and confining pressure remain unchanged, which is the compactness of coarse-grained fillers under cyclic loading. Deeply understanding and revealing this characteristic is conducive to guiding the research on the impact of vibration A constitutive model for the HSR subgrade load of roller during construction and high-speed train load during operation on compaction and vibration effect of coarse-grained subgrade fillers.

Crushing
Under static and dynamic loads, different degrees of grain crushing are generated for HSR coarse-grained subgrade filler. Hardin (1985) proposed to describe the crushing effect of coarse-grained soil filler under cyclic loading with relative grain crushing rate Rr, namely where S MPQ is the total grain crushing degree, which is expressed by the area surrounded by the initial grading curve, the gradation curve after loading and the straight line with particle size equal to 0.075 mm in the gradation curve of coarse-grained soil filler for HSR subgrade before and after loading shown in Figure 2; S MPN is the degree of potential grain crushing, which is expressed by the area surrounded by the initial gradation curve, 100% straight line and straight line with particle size equal to 0.075 mm in Figure 2.
3. Establishment of the cyclic densification model of coarse-grained soil filler for the HSR subgrade 3.1 Model assumption The stress path of the coarse-grained soil filler for the HSR subgrade caused by cyclic triaxial loads under stress control conditions is shown in Figure 3. First, samples are subject to confining pressure and isotropically consolidated to the initial mean stress σ p,0 , then reaches the maximum mean stress σ p,max (point B) under initial stress and then unloads to the minimum mean stress σ p,min (point A). The subsequent loading and unloading events take place between points A and B along the path AB. In order to establish the dynamic constitutive model of coarse-grained soil filler under cyclic loading in a simple and convenient manner, the following two assumptions are proposed.
(1) Under the initial loading conditions, the coarse-grained soil filler produces elasticplastic deformation. The maximum stress state reached by the initial loading (2) There is one elastic surface that varies with the vibration time. The deformation outside the elastic surface is plastic deformation, while the deformation inside the elastic surface is regarded as elastic deformation.

Considering the shear dilatancy equation of grain crushing
Under triaxial stress, based on the principle of energy conservation, the following relationships of stress, strain and grain crushing can be obtained (Indraratna et al., 2012).
where σ q and σ p are deviatoric stress and mean stress, respectively; dε v and dε 1 are the increments of axial strain and volume strain, respectively; w f is the basic friction angle; E B is the energy consumed per unit volume of grain crushing. According to the study in Jia, Xu, Chi, Xiang and Zhou (2017), E B can be expressed by plastic work W p per unit volume, while R r has a hyperbolic relationship with W p , namely where a and b are empirical parameters.
For the indoor triaxial test of coarse-grained soil filler of HSR subgrade under cyclic loading, W p can be calculated by the area enclosed by the hysteretic curve, and then the R r under cyclic loading can be obtained by Formula (3). Then, W p and R r are corresponding to the cumulative deformation one by one, and finally the relative grain crushing rate dW p de 1 and the plastic work consumption rate dRr dε 1 per unit volume can be obtained, respectively, which can be expressed by linear relationship (Indraratna et al., 2012) as where β is the empirical parameter. Therefore, in combination with the expression forms of stress invariant and strain invariant, formula (2) can be rewritten as follows: where dε s,p and dε v,p are plastic body strain and shear strain increments respectively; ηM is the critical state stress ratio; η i is the current stress ratio.
In actual engineering, the coarse-grained soil filler of HSR subgrade is affected by many factors. The cyclic stress ratio η c and the vibration time N having a significant impact on the material performance are adopted as the influencing factors, and the influence on the grain crushing rate dRr dε 1 is considered. In order to describe the loading and unloading amplitude of vibration load under different working conditions, a multiplier ln σp; max σ p; min is introduced to establish the evolution function of grain crushing rate with strain development, namely where f SR and f N are the correction functions of the influence of cyclic stress ratio and cyclic vibration time on grain crushing rate, respectively.
In σ q Àσ p stress space, according to the relationship ε s;p ¼ 2 3 ðε 1 − ε 3 Þ between shear strain ε s,p and axial strain ε 1 and radial strain ε 3 under axisymmetric stress and taking the Poisson's ratio of 0.3, ε 1 5 1.15 ε s;p can be obtained, in combination with formula (6) and substituted into formula (5), Considering some constant parameters in formula (7), for the convenience of subsequent expression, formula (7) is rewritten as (8) is the shear dilatancy equation of coarse-grained soil filler under cyclic loading considering the grain crushing effect.

Incremental expression form of constitutive equation
In order to accurately describe the cumulative deformation characteristics and to facilitate numerical programming, the classical elastic-plastic theory often decomposes the total strain where ε e and ε p are elastic strain vectors and plastic strain vectors, respectively.
The increment form of the above formula can be written as where dε, dε e and dε p are the total strain increment vector, the elastic strain increment vector and the plastic strain increment vector, respectively.
Similarly, for the strain invariants corresponding to the σ q − σ p stress space, the shear strain increment dε s and volume strain increment dε v can be written as dε s ¼ dε s;e þ dε s;p (11) where dε s,e , dε s,p , dε v,e , and dε v,p are elastic shear strain increment, plastic shear strain increment, elastic volume strain increment and plastic volume strain increment, respectively. In the elastic range, dε s,e and dε v,e can be calculated by the following formulas, respectively.
where e i is the initial void ratio at the beginning of shear; G is the shear modulus and k is the slope of the rebound curve of the isotropic compression rebound test.
The strain vector for the plastic phase can be obtained from the following formula (Indraratna et al., 2012).
where f h , f g and f f are the hardening function, the plastic potential function and the yield function, respectively.

Plastic potential function
In σ q Àσ p stress space, the direction of the plastic strain increment is perpendicular to the plastic potential surface, so Put formula (16) into formula (8), we can get dσ q dσ p þ It is not difficult to find out that the solution of formula (17) is the plastic potential function. However, it can be seen from formula (15) that precise plastic potential function is not required, but only partial differential forms of plastic potential function for σ q and σ p are required. With reference to the simplified calculation method in Indraratna et al. (2012), partial derivatives of plastic potential function for partial stress σ q and mean stress σ p are obtained, respectively: A constitutive model for the HSR subgrade

Yield function
In this paper, the model only considers the plastic deformation of coarse-grained soil filler of HSR subgrade caused by shear stress, while the influence of hydrostatic pressure on plastic deformation of fillers is ignored. In σ q Àσ p stress space, the yield function f f can be expressed as follows: When the loading path of cyclic loading coincides with the yield function, i.e. the yield function is 0, the material yield occurs, which is characterized by plastic deformation.

Plastic strain calculation under initial loading conditions
Using formula (15)-(20) and in combination with the method proposed by Indraratna and Salim to derive the hardening function from undrained stress conditions (Salim & Indraratna, 2004), the plastic strain under initial loading conditions can be obtained as follows: where a is the parameter related to the initial stiffness of coarse-grained subgrade filler; σ p,cu0 and σ p,cs are the stresses at the intersection of the starting point of the undrained stress path and the critical state line in the current stress state; σ p,cu0(i) and σ p,cs(i) are the initial values of σ p,cu0 and σ p,cs at the moment of loading under stress.
Since the current stress ratio may approach the critical stress ratio under low confining pressure, which will lead to a larger shear strain increment value, according to Indraratna instead of η i , so that the value of stress ratio in numerical calculation can be greater than the critical stress ratio, thus ensuring the stability and robustness of numerical calculation. According to Indraratna et al. (2012), formulas (21) and (8)

Plastic strain calculation under subsequent loading conditions
As above mentioned in this paper, under cyclic loading, the first loading of coarse-grained soil filler of HSR subgrade generates the maximum strain, and the cumulative deformation gradually tends to be stable with the increase of vibration times, namely, showing compactness. From the perspective of mathematical expression of constitutive equation, parameters representing material rigidity can be modified to simulate the compactness of subgrade filler. Therefore, the stiffness in the plastic strain formula under initial loading condition is modified, and a new plastic shear strain increment formula under the subsequent loading condition is obtained: where α 1 is the parameter considering kinematic hardening effect; α 2 is a parameter considering the influence of subsequent loading stress ratio; α 3 is a parameter considering the influence of vibration times; μ and δ are model empirical constants related to cyclic vibration times; jjABjj, jjACjj and jjADjj are the distances between two points in the stress space of σ q Àσ p (see Figure 3); σ p,e is the elastic average stress (corresponding to point D in Figure 3). α 1 is introduced to accurately simulate the Bauschinger effect. Since the material is elastic after a certain number of vibrations, an elastic surface is introduced to describe the effect of the loading history. The elastic surface hardens with the increase of vibration. However, σ p,e change in the given stress path AB, and the relationship between σ p,e and vibration time N given (Indraratna et al., 2012) is In Figure 3, point D is used to control the deformation of the material under subsequent loading conditions, and point C represents the current stress position. If C is higher than D, the material shows elastoplasticity, otherwise, the material shows elasticity. α 2 is introduced to simulate the effect of the magnitude of stress ratio. When the filler is subjected to greater stress, namely, the closer the stress ratio is to the critical stress ratio, the greater the stiffness of the filler and the smaller the deformation. α 3 is introduced to simulate the influence of cyclic vibration times. It is a function related to vibration times N. With the increase of loading times, the parameter α 3 gradually increases, which further reflects the increase of stiffness of the filler and the decrease of cumulative deformation.
After the shear strain is calculated by the above formulas, the plastic volume strain increment can be calculated by formula (24). Finally, the total shear strain and volume strain increments can be calculated by formulas (11) and (12).
4. Implementation of the cyclic densification model of coarse-grained soil filler for the HSR subgrade in the ABAQUS procedure 4.1 Numerical implementation process of cyclic densification model The ABAQUS is a piece of finite element software, with a wide application range and powerful functions. It can simulate various problems from simple linear to complex nonlinear, which only provides standard finite element analysis program, but also has good openness. Users can develop user subprogram interface and application program interface with the help of the secondary development platform so as to expand the function of the master program and widely apply in the actual engineering (Cao, 2014). This paper adopts FORTRAN language to prepare a UMAT subprogram of cyclic densification model of coarse-grained soil filler for the HSR subgrade based on the secondary development interface of ABAQUS. The main process of UMAT subprogram is shown in Figure 4. A constitutive model for the HSR subgrade 4.2 Numerical implementation and verification of the cyclic densification model Based on the above description of the cyclic densification model and finite element numerical implementation method, two-dimensional symmetry method is adopted, and the finite element subprogram is prepared to predict the triaxial test results of coarse-grained subgrade filler. Figure 5 is a schematic diagram of numerical simulation of the cyclic triaxial test. For convenience, the two-dimensional plane of the cylindrical triaxial sample is simulated, i.e. the dash area of the triaxial sample in Figure 5(a), and the two-dimensional eight-node in ABAQUS is used for quadratic-reduced integration unit; Figure 5(b) shows the mesh division and boundary conditions of the ABAQUS model. The left boundary of the model is the central axis, without lateral deformation. The lower boundary of the model restrains vertical and radial deformation. The top of the model is fixed radially, but free deformation is allowed in the vertical direction. The triaxial test process of cyclic loading is simulated with two analysis steps: first, defining the model boundary and confining pressure; second, applying vertical cyclic loading to the model.
In this paper, numerical calculation is carried out for triaxial test of coarse-grained soil under cyclic loading in Bian et al. (2016) and Tian et al. (2019). In numerical calculation, the model parameters are taken as follows.
According to the indoor vibration and compaction test results, the relation curve between the relative grain crushing rate and the plastic work per unit volume is drawn, as shown in Figure 6. Formula (3) is used to fit the test results, and the empirical parameters a and b of the coarse-grained soil filler of the HSR subgrade are equal to 1,855 and 0.083, respectively.  The hysteresis loop area enclosed by the stress-strain curve obtained by triaxial test under cyclic loading is used to calculate W p , and the parameters a and b determined above are used to calculate R r under cyclic loading. Then, W p and R r are corresponding to the cumulative deformation, and finally, the relative grain crushing rate dR r de 1 and the energy consumption rate per unit volume dW p de 1 are obtained, respectively, as shown in Figure 7. It can be seen from Figure 7 that linear relation is shown between the two, and the empirical parameter β obtained by fitting equals to 1,855.
According to the triaxial test results in Tian et al. (2019), the relationship between the f SR and the cycle stress ratio and between the f N and the cyclic vibration times can be obtained, as shown in Figure 8. Based on the least square method principle, the empirical expressions of f SR and f N are obtained by fitting . σ d is the vertical dynamic stress amplitude of the triaxial test; σ 3 is the confining pressure of triaxial test; C SR1 , C SR2 , C SR3 , C N1 , C N2 and C N3 are the fitting parameters.   Table 1 shows the specific values of the fitting coefficients of formulas (26) and (27) for the correction function of relative grain crushing rate. Other parameters involved in the simulation of the triaxial test of the coarse-grained soil filler of HSR subgrade under cyclic loading can be obtained according to the "trial and error" method. The specific values are shown in Table 2. Figure 9 shows the cumulative deformation test results and numerical calculation results of the coarse-grained soil filler of HSR subgrade under cyclic loading. It can be seen  Figure 9 that the calculation results of the cyclic densification model proposed in this paper can more accurately predict different confining pressures and accumulative deformation under stress amplitude but also can accurately predict that the cyclic densification characteristics that the cumulative strain will continue to stabilize as the cumulative vibration time continues to increase. This fully shows that the theoretical formula of the cyclic densification model and the finite element calculation method proposed in this paper are accurate and reasonable.

Conclusion
(1) Based on the stress characteristics of long-term repeated low-amplitude cyclic loading of the HSR subgrade, for the coarse-grained soil filler of HSR subgrade, according to the critical state soil mechanics theory, the empirical expressions for grain crushing and plastic work of filler are given, and the shear dilatancy equation of coarse-grained soil filler that can consider the influence of grain crushing and stress ratio is corrected and proposed.
(2) Based on the classical elastic-plastic theory, the correction function of stress ratio and cyclic vibration times is introduced. Through elastic-plastic decoupling, the cyclic densification constitutive model of subgrade coarse-grained soil filler based on the modified shear dilatancy equation is derived.
(3) The numerical analysis and solution program of "elastic prediction-plastic correction" of the proposed model is established. language, and the programmed calculation of the constitutive model is realized through the secondary development program interface provided by the finite element software ABAQUS.
(4) Compared with the triaxial test of coarse-grained soil of HSR subgrade under cyclic loading done previously, it shows that the cyclic densification model considering grain crushing proposed in this paper can more accurately simulate the long-term cumulative deformation and cyclic densification characteristics of fillers.