Testing and Micromechanical Modelling of Rockfill Materials Considering the Effect of Stress Path

We have extended the micromechanics-based analytical (M-A) model to make it capable of simulating Nuozhadu rockfill material (NRFM) under different stress paths. Two types of drained triaxial tests on NRFM were conducted, namely, the stress paths of constant stress ratio (CSR) and the complex stress paths with transitional features. The model was improved by considering the interparticle parameter variation with the unloading-reloading cycles and the effect of the stress transition path. The evolution of local dilatancy at interparticle planes due to an externally applied load is also discussed. ComparedwithDuncan-Chang’s E-u and EBmodels, the improvedmodel could not only better describe the deformation properties ofNRFMunder the stress path loading, but also present the volumetric strain changing fromdilatancy to contractancywith increasing transitional confining pressures. All simulations have demonstrated that the proposed M-Amodel is capable of modelling the mechanical behaviour of NRFM in the dam.


Introduction
Rockfill materials (RFMs) are widely used in the constructions of rockfill dams because of their high strength and low cost [1][2][3][4][5][6][7].Normally, the large-scale conventional triaxial compression (LCTC) test (i.e., the confining pressure is constant) is used to study the mechanical properties of RFMs [8][9][10][11][12][13][14].However, during such tests, the loading path is not consistent with that in actual engineering applications.Observation data of dams and the results of finite element calculation analyses have shown that the stress path of RFMs can be approximated as the stress ratio remaining constant during dam construction, and it will become a transitional stress path during reservoir filling (Figure 1) [15].Moreover, the mechanical properties of RFMs are stress-path dependent [16].In other words, the stress-strain relationships of the RFMs vary greatly under different loading paths.Therefore, it is necessary to consider the influence of the stress path in the constitutive modelling of RFM.
Many LCTC tests have been conducted for examining the behaviour of RFMs.Marsal [17,18] conducted a series of super LCTC tests (the sample was 1130 mm in height and 180 mm in diameter).He found that the peak state friction angle of the RFMs decreases with increasing confining pressures.Indraratna and Salim [19] investigated the stressdilatancy phenomenon of RFMs during the LCTC test.Then, they proposed a relative constitutive equation to improve the stress-dilatancy relationship by incorporating particle breakage.Araei et al. [20] conducted the LCTC tests of RFMs under cyclic loading conditions and studied the effect of loading frequency on the sample.Xiao et al. [21,22] performed several LCTC tests of RFMs with different initial void ratios and investigated the effects of density and pressure on the strength and deformation of RFMs.However, the confining pressure of the LCTC test was constant, and under the stress path, the changing confining pressure had a significant impact on the mechanical behaviour (such as strength, stress-strain relationship, and dilatancy) of the RFMs.Although many stress path tests of granular material had been conducted, these tests primarily focused on sand.Few stress path tests have been reported for RFMs.Liu et al. [23] conducted several large-scale triaxial tests of Yixing RFM under different stress paths way.Subsequently, they determined that the stress path has little influence on the shear strength of the RFMs.Yang et al. [24] performed large-scale triaxial tests of RFM under three types of loading conditions, in which the confining pressure, mean principal stress, and the principal stress ratio were constant.They observed that the plastic work could well reflect the breakage of coarse-grained material under different loading conditions.However, the loading conditions of these tests were not consistent with those in an actual dam.Stress path tests, which can reflect the real stress state of the RFMs in a dam, have not been conducted.
Based on large-scale triaxial tests, many constitutive models have been proposed for RFMs.Schofield and Wroth [25] proposed the critical state theory, which represents the beginning of modern soil mechanics.Many scholars have employed the concept of critical state line (CSL) in models for RFMs.However, the critical state theory does not apply to RFMs under low stress conditions (where the confining pressure is less than 500 kPa), because the critical state phenomenon of the sample is often not observed.Duncan and Chang [26] proposed Duncan-Chang's E-u and E-B models.The parameters of these two models can easily be selected from the test and have explicit physical meanings.Therefore, these models are widely used for simulating the deformation behaviour of RFMs.Finite element analyses using these two models have shown that the calculation of the settlement displacement is acceptable, but the calculated horizontal displacement considerably differs from the observation data, which indicates that the model based on the LCTC tests is still incomplete.Several models based on stress path tests have been proposed for RFMs by fitting the test curve [27], but these models are too simple and the theoretical studies are not sufficiently thorough.Although many constitutive models have been proposed for RFMs, the model based on the stress path test has not been fully investigated.With the development of the high rockfill dam, a thorough study of the stress path model for RFMs is urgently needed.
In recent years, increasingly more researchers have begun to model granular materials from the particle length scale.The discrete element method is a commonly used approach for studying the effect of particle properties on macroparticle assemblies.Liu et al. [28] investigated breakage of the rockfill particles in a dam model using Particle Flow Code (PFC) software.Discrete element modelling of sand under onedimensional compression and creep conditions was studied by McDowell and De Bono [29] using the PFC software; however, the calculation speed of PFC is slow; it is difficult to analyse the actual engineering.The M-A model not only appropriately considers the microstructure of soils but also is easier to be applied in engineering practice.Based on this model, Chang and Hicher [30] simulated the conventional triaxial test of Hostun sand under drained and undrained conditions.The simulation results showed that the model could well describe the stress-strain behaviour of the sand during the test.Chang et al. [31,32] introduced a new slip mechanism of the contact plane in the M-A model; consequently, the model could well reflect the dilatancy phenomenon of Hostun sand under compression and extension loading conditions.Yin et al. [33] proposed a new formulation that accounts for the stress reversal on a contact plane and the density state-dependent dilatancy.Subsequently, the M-A model was extended to be capable of simulating sand under cyclic loading.Previous studies of the M-A model primarily focused on the mechanical properties of sand.The models for RFMs have not been investigated as intensively as the models for sand.Li proposed the M-A model for RFMs based on LCTC tests, and the simulation results indicated that the model could well describe the deformation properties of RFMs during the test [34].However, modelling of RFMs under the stress path has not been conducted.
In this paper, the M-A model is extended to make it capable of simulating the triaxial tests of NRFM under the stress path.The model is improved by considering the parameter variations under unloading-reloading conditions and the stress transition path.The local dilatancy relations at interparticle planes due to applied load are also discussed.The features of the M-A model are described in comparison with Duncan-Chang's models.Finally, the overall applicability of the proposed model is evaluated through comparisons of the prediction and test results.

Two Types of Large Triaxial Tests
RFMs were taken from the site of Nuozadu rockfill dam, which is the highest rockfill dam in China (ranked the third highest in the world).Specific details of the study site are presented in Table 1.These materials are used as the main rockfill of the core-dam body.Table 2 presents the basic properties of the NRFM, including particle shape, specific density, initial void ratio, and compressive strength.As shown in Figure 2, the particle size distribution of the prototype NRFM is reduced using the parallel gradation technique [35], with a maximum particle size of 60 mm.The uniformity coefficient Cu and the curvature coefficient Cc of the NRFM are 8.8 and 1.47, respectively.The test equipment is a large triaxial apparatus.The specimen is 302 mm in diameter and 655 mm in height.The maximum confining pressure of the device is 2500 kPa, the maximum axial force is 1000 kN (compression), and the maximum axial displacement is 200 mm.The specimen was compacted and then saturated with CO 2 and air-free water (-value was above 0.96 after saturation).The specimen was consolidated at an initial confining pressure of 40 kPa.Then, it was loaded at the set value of stress ratio (as shown in Figures 3 and 4) under the drained condition, which was consistent with the practical working conditions of NRFM in the dam.Because the pore water pressure was zero, all the stresses in this paper refer to the effective stresses.
Ten samples were prepared for the tests.According to the loading method, the samples could be divided into two types.For the first type of test (Figure 3), the sample was loaded under the paths of a constant stress ratio (CSR).When the confining pressure reached 1600 kPa, a cycle of unloading and reloading was applied under the same CSR.After the cyclic loading, the loading was continued until the confining pressure  3 reached 2500 kPa.Because there were six different stress ratios in the tests, six samples were required.
For the second type of test (Figure 4), the stress ratio was maintained constant (/ = / = 1.00) during the initial stage.Then, at each confining pressure,  3 = 300, 800, 1400, and 2000 kPa, the confining pressure was held constant and the axial load was increased (/ = 3.00) until the samples reached failure.Four samples were required according to the four transitional points.

Constitutive Model
In the proposed constitutive model, the rockfill sample is considered to be a collection of grains.The deformation of the sample is generated by moving contacted grains in all orientations (Figure 5).Thus, the stress-strain relationship of the NRFM can be obtained from the average of the mobilized deformation of all contact planes between the rockfill grains.The forces and movements at all contact planes are used to obtain the macroscopic stress and strain tensors (Figure 7).The macroscopic stiffness tensor is obtained by integrating the stiffness properties at interparticle contacts.A statically  constrained microstructure is assumed for the micro-macro links, which means that the forces on each contact plane are assumed to be equal to the resolved components of the macroscopic stress tensor [30].

Local Coordinates and Definitions.
As shown in Figure 5, an auxiliary local coordinate is established on each contact plane.The orientation of a contact plane between two grains is defined by the vector perpendicular to the plane.For the  contact plane, the local forces    and the local movements    can be denoted as follows: , where the subscripts , , and  represent three orthogonal unit vectors that form the local coordinate system.The vector  is inward normal to the contact plane.Vectors s and t are on the contact plane.

Force-Displacement Relations of the
The tensor forms of (1) are as follows: It is clear that    is a diagonal matrix, which is as follows: In the global coordinate system, (2) can be expressed as Then where,    and    are the increments of the local forces and the local movements in the global coordinate system.   is the stiffness matrix in the global coordinate system.  is the coordinate transformation matrix. and  are the angles between the vectors and coordinate system (Figure 5).
The value of the elastic stiffness can be estimated from Hertz-Mindlin's formulation [36] where    is the contact force in the normal direction. 0 ,  0 , and  are material constants. is the branch length between two grains.For two spherical grains,  is the same as the grain size .According to Hertz-Mindlin's theory [36],  = 1/3.

Plastic Part.
Stress dilatancy is a well-known phenomenon in RFMs [37][38][39].During triaxial tests, the plastic sliding of the particles made the contact plane move, and then shear dilation occurs [31,32].Assuming that that plastic work for a contact plane due to both normal and shear movements is equal to the energy loss due to friction at the contact, the dilatancy effect can be described by Then where    is the increment of the normal plastic movement. is the dilatancy angle.tan  represents the obliquity when the plastic normal movement is zero, which is related to the phase transformation line of the rockfill assembly.The shear force   and the increase of the share plastic movement   can be given by The yield function is assumed to be of Mohr-Coulomb type where  is the hardening function, which is defined by a hyperbolic curve (Figure 6) as follows: where   is the interparticle friction angle and  0 is the initial slope of the hyperbolic curve, which can be expressed as follows: The plastic flow in the direction normal to the contact plane is governed by the stress-dilatancy equation in (8).Thus the flow rule is nonassociated.In the global coordinate system the incremental plastic movements    can be expressed by  13) Equation ( 16) Equation ( 4) Equation ( 18) Equations ( 10)-( 11) During the integration process, a relationship is required to link the macro and micro variables.Liao et al. [40] proposed that the stress-strain relationship for a granular system could be derived based on the hypothesis that the mean field of displacement is the best fit of actual particle displacements.This relationship has been introduced in the M-A model by Chang et al. and could well describe the stress-strain behaviour of the sand [30][31][32][33].As a kind of granular materials, this hypothesis is also applied to the NRFM.The relation between the macrostrain and interparticle displacement is as follows: where    is the incremental relative displacement between two contacted particles.   is the branch vector joining the centres of two contacted particles. is the total number of contact orientations, according to Oda's research [41], which can be expressed as Using this hypothesis, the mean force on the contact plane of a given orientation is The stress increment can be obtained by the contact forces and branch vectors for all contacts as follows [42,43]: When the contact number N is sufficiently large in an isotropic packing, the summation of flexibility tensor in ( 14) can be written in integral form [44], given by where  is a distribution density function, and in this paper, its value is 1/4, which means that the sample is isotropic.

Calibrations of Model Parameters
The model parameters are divided into two parts: the global parameters, which can be chosen from the test, and the interparticle parameters, which cannot be obtained from the test directly.Here, a particle swarm optimization (PSO) algorithm (Figure 8) was used to determine the interparticle parameters.
PSO is a population-based heuristic search technique.Its original concept came from the movements of birds or particles [45].This algorithm possesses a fast convergence rate and high estimate precision of sensor bias.Therefore, it is widely used in geotechnical engineering.In this algorithm (Figure 8), each potential solution is called a "particle."The solution is searched by a particle swarm flying in the hyperspace.Based on the optimization problem, the objective function of the optimization problem is evaluated for each particles position in the search space.The fitness values of each position are calculated.The best-known positions of each particle are denoted best, and the global best position of the entire particle swarm is denoted gbest.For each iteration, the velocity and particle location are expressed by where    and    are the velocity and the position of the th particle in the previous iteration, respectively. 1 and  2 are random numbers ranging from 0 to 1.  1 and  2 are acceleration coefficients, which are the positive parameters; here, they are equal to 2.    and    are best and best, respectively.In this paper, the fitness value of the objective function, which determines the POS performance, is represented by the mean squared error, and it is expressed as where x is the calculated value, which represents the axial strain and volume strain, respectively;   is the test data; and  is the total number of test data.As the number of search steps increases,   will gradually become smaller and satisfy the end conditions.Table 3 shows the search range of the parameters.According to the initial stress increments and the constitutive model, the optimizer calculates the strain increments.Then, the error   can be calculated through comparison with test data.The candidate solution is improved throughout the iteration process.When the solution reaches the best fitness, the model parameters are obtained.Considerable efforts were devoted to determining the interparticle parameters of the model under the stress paths of different CSRs.In particular, under unloading-reloading conditions and the stress transition paths, the interparticle parameters might be different; specific instructions are as follows.
Under unloading-reloading conditions, the stress-strain curve is nearly a straight line during the test (Figure 14).Therefore, we assumed that only linear elastic deformation of the particle occurred during this stage (plastic deformation could be neglected).From the analysis of the fitting parameter     , a significant phenomenon was found: parameter    varies regularly with the mean principal stress  1 at the unloading point.After plotting lg(   /  ) − lg( 1 /  ) in Figure 9, a good linear relationship is observed.Therefore, the value of the normal stiffness for two elastic spheres could be expressed as follows: where  1 is the intercept of a straight line and  is the slope.The normal stiffness    is similar to the loading and unloading modulus of the conventional triaxial test [26], because the two types of tests have similarly shaped deformation curves (see Figure 13(a)).
Under the stress transition path, as the transition confining pressures decrease, the volumetric strain changes from continuous contractancy to final regular dilatancy (as shown in Figure 15).The employed local dilatancy equation (shown as ( 8)) could well describe the volume changes of sand [31,32] under compression and extension loading conditions.However, for the NRFM, the large rockfill grains are easily broken as the pressure increases.The abrasion of the rock might reduce the friction coefficient of the contacted particles [46][47][48].From the analysis of the fitting parameter , we found that the dilatancy angle  decreased as the transitional confining pressure  3 increased.As shown in Figure 10, the dilatancy angle  and the transitional confining pressure  3 exhibit a good linear relationship under the logarithmic coordinates.Here, the dilatancy angle of the contacted particles  can be determined as follows: where  0 and Δ are parameters whose values are equal to the intercept and slope of the line, respectively.

The Evolution of Local Dilatancy at Interparticle Planes
The stress-dilatancy behaviour of the NRFM calculated using the M-A model is based on the mobilization of contacted particles.Therefore, it is necessary to study the local dilatancy behaviour on individual contact planes between two particles.In the large triaxial test, the applied loads are axisymmetric about the -axis.Therefore, the orientation of the contact plane can be represented by an inclined angle , which can be measured by the branch vector and the -axis, as shown in Figure 11(a).Eight contact planes were selected for the investigation: Υ = 15 ∘ , 25 ∘ , 35 ∘ , 45 ∘ , 55 ∘ , 65 ∘ , 75 ∘ , and 85 ∘ (shown in Figure 11(a)).
In order to study the local dilatancy behaviour of the particles, the triaxial test of the type two with the transition confining pressure  3 = 300 kPa (see Figure 4) was chosen to be calculated.According to the employed local dilatancy equation (see (8)), the local dilatancy rate    /  was defined here.When the value of the local dilatancy rate is negative, the contact plane moves upward and shear-induced dilation occurs; conversely, when the value is positive, the contact plane moves downward and the contraction occurs (Figures 11(g) and 11(h)).
The stress-dilatancy relations for the eight selected contact orientations are plotted in Figures 11(b  are contracted.The 65 ∘ contact plane reaches the highest mobilized ratio (close to /4 + Φ  /2 = 64.6 ∘ ).When stress transition occurred, a clear turning point was observed, as shown in Figures 11(b)-11(i).As the axial force increased, the local dilatancy rate decreased.The contact planes are contracted/dilated at different degree for the given applied stress ratio.The 65 ∘ and 75 ∘ contact planes initially exhibited contraction.Then, they mobilized to higher stress levels and reached the dilation region (Figures 11(g) and 11(h)).The other planes only exhibited contraction.
Figure 12 shows the relationship of the local dilatancy rate    /  and the global stress ratio /.The contact planes are contracted/dilated to different degrees for a given applied stress ratio.Note that the local dilatancy relations, with the planes greater than 55 ∘ , are no longer linear when plotted against the global stress ratio.Similar nonlinear relationships were also found during the compression and extension test of Hostun sand simulated by Chang et al. [31,32].

Comparison of Predicted Results with Three Models.
Duncan-Chang's E-u and E-B models are widely used for calculating the deformation of RFMs.Here, we verify the validity of the M-A models by comparing the predicted results with those of Duncan-Chang's models.The parameters of the models could be obtained from the conventional triaxial tests ( 3 = constant) of the NRFM.Detailed information is available in the literature [49].The values of the parameters are listed in Table 5.Because there are too many samples to be tested, two samples were chosen to be simulated.One was loaded under a principal stress ratio of  1 / 3 = 4.0.The other sample had a transition confining pressure of  3 = 800 kPa (Figure 13).The comparisons of the predicted results are as follows.
As shown in Figures 13(a) and 13(c), the shapes of the  ∼  1 curves simulated using the E-u and E-B models are the same.The reason for this result is that, under the same stress conditions, the values of the tangent modulus are equal for the two models.The shapes of the  ∼  1 curves are similar to that of the test result.However, the calculated values are clearly lower than the test results.The reason might be that the parameters are chosen from the conventional triaxial test; thus, the model cannot reflect the effect of the stress path on the deformation of the samples.The parameters of the M-A model were selected from the stress path test; thus, the simulation results and the experimental results are well fitted.
As shown in Figures 13(b) and 13(d), the  V ∼  1 curves simulated using the E-u and E-B models are both significantly different from the test results.One reason for this result is that the models cannot reflect the effect of the stress path.Another reason is that they cannot describe the dilatancy behaviours of the NRFM due to the limit of elasticity theory.The dilatancy equation (see (8)) of the M-A model could well describe the dilatancy behaviours of the NRFM.Therefore, better simulation results can be obtained.was chosen to determine the normal stiffness    under the unloading-reloading conditions.The dilatancy angle  under the stress transition path could be obtained from (23).The diameter of the particle was chosen as the characteristic particle size  50 = 18 mm suggested by Chang and Yin [31].All the parameters are listed in Tables 2 and 4. For several tests, without being used for the parameter determination, they were calculated; the calculations could be considered as genuine model predictions.For the unextended M-A model, the details information could be seen in the literature [30].

M-A Model
The comparisons between model predictions and type one test are shown in Figure 14 in terms of the stressstrain-volume relationship on NRFM.The unextended M-A model overestimates the deviatoric stress especially under    relationships of NRFM (Figure 15(b)), while the model without considering the effect of stress path cannot well represent the stress-strain relationship (overestimating the dilation of the NRFM).Under the CSR path, the change in the behaviour of the sample is similar to that of the type one test at  1 / 3 = 2.5.However, at the transition confining pressure, the stress path moves towards the failure line and the axial strain rapidly increases.The volume of the sample changes from contraction to dilatancy as the transition confining pressure decreases.Because the load path of the test is similar to the stress variation of the actual RFM in the dam, this model could well describe the constitutive relationship of the NRFM under the paths of a constant stress ratio during dam construction and the transitional stress paths upon reservoir filling.

Conclusions and Discussion
A micromechanical model considering the effect of stress path was proposed for NRFM.The main conclusions and discussion are summarized below.
The stress path of RFMs during dam construction can be approximated as the CSR load, and it will become a transitional feature upon reservoir filling.Corresponding to these engineering conditions, two types of triaxial drained tests were conducted to simulate the actual stress variations of the NRFM.In the first type of test, no failure occurred, and the NRFM exhibited contraction during the entire loading period.However, in the second type of test, the samples reached the failure point.Dilatancy became increasingly more serious as the transition confining pressure decreased.
At the particle length scale, the movements of the local contacted planes can well reflect the deformation behaviour of the NRFM sample.The calculation speed of the model is fast, for the mean force and the movements of the contact plane were used to obtain the macro stress and strains.However, the model cannot reflect the deformation characteristics of specific particles.
The M-A model has been extended.For the NRFM under unloading-reloading conditions, the new stiffness of the particles can well predict the linear relationship of the stressstrain curves; for the tests under the stress transition path, the local dilatancy angle can well describe the volume of the sample changing from contractancy to dilatancy.Compared with Duncan-Chang's models, the M-A model could better describe the deformation behaviour of the NRFM under the stress path.Model predictions for the drained triaxial tests have demonstrated that the present micromechanical approach is capable of modelling the deformation behaviours of NRFM under actual stress paths in the dam.

Figure 1 :
Figure 1: Stress path of RFMs in the dam.

Figure 3 :
Figure 3: Stress paths of constant stress ratio.

Figure 4 :Figure 5 :
Figure 4: Stress paths of constant stress ratio with a transition.

Figure 6 :
Figure 6: The curve of the hardening function in the -  plane.

Figure 7 :
Figure 7: Calculation flowchart of the M-A model.
Initialize the particle swam Calculate the fitness of all particles Search for pbest and gbest Update the particles' location and velocity Calculate the fitness of all particles Update pbest and gbest Satisfy the end conditions No Yes Obtain the optimum parameters of the model

Figure 9 :
Figure 9: Relationship between normal stiffness and the principal stress at the unloading point.

Figure 12 :
Figure 12: Local plastic movement increment versus global stress ratio for selected contact plane orientations.

Figure 13 :
Figure 13: Comparisons between three types of model simulations and test results on NRFM: type one test, (a) and (c): type two test, (b) and (d).

Figure 14 :
Figure 14: Model simulations of type one test.(a) Comparisons between unextended M-A model simulations and test results on stress-strain relationships of NRFM; (b) comparisons between extended M-A model simulations and test results on stress-strain relationships of NRFM.

Figure 15 :
Figure 15: Model simulations of type two test.(a) Comparisons between unextended M-A model simulations and test results on stress-strain relationships of NRFM; (b) comparisons between extended M-A model simulations and test results on stress-strain relationships of NRFM.

Table 1 :
Project and site for NRFM.

Table 2 :
Physical and mechanical properties of rockfill.
Suppose that the relationship of the local forces    and the local movements    is nonlinear elasticity.The contact stiffness of the plane includes normal stiffness    and shear stiffness    and    .The elastic stiffness is given by Contacted Grains 3.2.1.Elastic Part.

Table 3 :
Parameter-search range of the M-A model.

Table 4 :
Parameters of the M-A model.

Table 5 :
Parameters of Duncan-Chang's E-u and E-B models.Figure15shows the comparisons between the model predictions and type two test.The extended model considering the effect of stress path can well capture the stress-strain (Figure 14(b)).Under CSR loading, all the axial strains are less than 15%.This result indicates that no failure occurred throughout the entire range of tests.Under the unloadingreloading condition, the elastic contact of the particles (proposed in the extended model) could well reflect the linear relationship of the stress-strain curves of the sample.The calculated volumetric-strain curves appear to be linear relationships and match the test data well (Figure14(b)).Because little changes in the volume occurred during this stage, the calculated volume changes are neglected.