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

Authors: Wang Feng; Chi Shichun; Li Shijie; Jia Yufeng
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

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. Compared with Duncan-Chang’s E-u and E-B models, the improved model could not only better describe the deformation properties of NRFM under the stress path loading, but also present the volumetric strain changing from dilatancy to contractancy with increasing transitional confining pressures. All simulations have demonstrated that the proposed M-A model is capable of modelling the mechanical behaviour of NRFM in the dam.


---

## Body

## 1. Introduction

Rockfill materials (RFMs) are widely used in the constructions of rockfill dams because of their high strength and low cost [1–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–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.Figure 1
Stress path of RFMs in the dam.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 stress-dilatancy 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 one-dimensional 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.


## 2. 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 Table1. 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.Table 1
Project and site for NRFM.
Project detail Depiction Project name Nuozhadu Hydropower Station Dam type Rockfill dam with clay core Site Puer City, Yunnan Province, China Type of NRFM Blasted rockTable 2
Physical and mechanical properties of rockfill.
Material Nuozhadu rockfill Particle shape Angular Specific densityρs (g/cm3) 2.63 Initial void ratioe0 0.32 Compressive strength (Mpa) 118.6Figure 2
Particle size distribution of NRFM.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 CO2 and air-free water (B-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.Figure 3
Stress paths of constant stress ratio.Figure 4
Stress paths of constant stress ratio with a transition.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 (Figure3), 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 (Figure4), the stress ratio was maintained constant (q/p=dq/dp=1.00) during the initial stage. Then, at each confining pressure, σ3t=300, 800, 1400, and 2000 kPa, the confining pressure was held constant and the axial load was increased (dq/dp=3.00) until the samples reached failure. Four samples were required according to the four transitional points.


## 3. 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 (Figure5). 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].Figure 5
Local coordinate at interparticle contact.


### 3.1. Local Coordinates and Definitions

As shown in Figure5, 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 a contact plane, the local forces FIα and the local movements δIα can be denoted as follows: FIα=Fnα,Fsα,Ftα and δIα=δnα,δsα,δtα, where the subscripts n,s, and t represent three orthogonal unit vectors that form the local coordinate system. The vector n is inward normal to the contact plane. Vectorss andt are on the contact plane.


### 3.2. Force-Displacement Relations of the Contacted Grains

#### 3.2.1. Elastic Part

Suppose that the relationship of the local forcesFIα and the local movements δIα is nonlinear elasticity. The contact stiffness of the plane includes normal stiffness Knα and shear stiffness Ksα and Ktα. The elastic stiffness is given by(1)dFnα=KnαdΔnαe,dFsα=KsαdΔsαe,dFtα=KtαdΔtαe.The tensor forms of (1) are as follows:(2)dFIα=KIJαdΔJαe,I,J=n,s,t.It is clear thatKIJα is a diagonal matrix, which is as follows:(3)KIJα=Knα000Ksα000Ktα.In the global coordinate system, (2) can be expressed as(4)dfiα=kijαdδjαe.Then(5)kijα=BiIαK′IJαBjJα,Bα=-cosγ-sin⁡γcos⁡β-sin⁡γsin⁡β-sinγcos⁡γcos⁡βcos⁡γsin⁡β0-sin⁡βcos⁡β,where, dfiα and dδjαe are the increments of the local forces and the local movements in the global coordinate system. kijα is the stiffness matrix in the global coordinate system. Bα 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](6)Knα=Kn0Fnαl2n,Ksα=Ktα=Kr0Fnαl2n,where Fnα is the contact force in the normal direction. Kn0, Kr0, and n are material constants. l is the branch length between two grains. For two spherical grains, l is the same as the grain size d. According to Hertz-Mindlin’s theory [36], n=1/3.


#### 3.2.2. Plastic Part

Stress dilatancy is a well-known phenomenon in RFMs [37–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(7)FnαdSnαp+TαdSαp=Fnαtan⁡φdSαp.Then(8)dSnαpdSαp=tan⁡φ-TαFnα,where dSnαp 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 Tα and the increase of the share plastic movement dSαp can be given by(9)Tα=Fsα2+Fnα2,dSαp=dSsαp2+dStαp2.The yield function is assumed to be of Mohr-Coulomb type(10)FFnα,κ=Tα-FnακSαp=0,where κ is the hardening function, which is defined by a hyperbolic curve (Figure 6) as follows:(11)κ=Kp0tan⁡ϕu·SαpFnαtan⁡ϕu+Kp0Sαp,where ϕu is the interparticle friction angle and Kp0 is the initial slope of the hyperbolic curve, which can be expressed as follows:(12)Kp0=kn0.0125+10fn.Figure 6
The curve of the hardening function in theκ-Sap plane.Figure 7
Calculation flowchart of the M-A model.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 dδiαp can be expressed by(13)dδxαpdδyαpdδzαp=dSαptan⁡φ-TαFnαFsαTαFtαTαBα.


### 3.3. Macro-Micro Relationship

The stress-strain relationship for the sample can be determined by integrating the behaviour of NRFM interparticle contacts in all orientations. 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–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:(14)dεij=Aik-1∑α=1Ndδjαlkα,Aik=∑α=1Nliαlkα,where dδjα is the incremental relative displacement between two contacted particles. lkα is the branch vector joining the centres of two contacted particles. N is the total number of contact orientations, according to Oda’s research [41], which can be expressed as(15)NV=12πd31+ee.Using this hypothesis, the mean force on the contact plane of a given orientation is(16)dfjα=dσijAik-1lkαV.The stress increment can be obtained by the contact forces and branch vectors for all contacts as follows [42, 43]: (17)dσij=1V∑α=0Ndfjαliα.When the contact numberN is sufficiently large in an isotropic packing, the summation of flexibility tensor in (14) can be written in integral form [44], given by (18)dεij=Aik-1N∫02π∫0πdδjαγ,βlkαγ,βξγ,βsin⁡γdγdβ,(19)Aik=N∫02π∫0πliαγ,βlkαγ,βξγ,βsin⁡γdγdβ,where ξ is a distribution density function, and in this paper, its value is 1/4π, which means that the sample is isotropic.


## 3.1. Local Coordinates and Definitions

As shown in Figure5, 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 a contact plane, the local forces FIα and the local movements δIα can be denoted as follows: FIα=Fnα,Fsα,Ftα and δIα=δnα,δsα,δtα, where the subscripts n,s, and t represent three orthogonal unit vectors that form the local coordinate system. The vector n is inward normal to the contact plane. Vectorss andt are on the contact plane.


## 3.2. Force-Displacement Relations of the Contacted Grains

### 3.2.1. Elastic Part

Suppose that the relationship of the local forcesFIα and the local movements δIα is nonlinear elasticity. The contact stiffness of the plane includes normal stiffness Knα and shear stiffness Ksα and Ktα. The elastic stiffness is given by(1)dFnα=KnαdΔnαe,dFsα=KsαdΔsαe,dFtα=KtαdΔtαe.The tensor forms of (1) are as follows:(2)dFIα=KIJαdΔJαe,I,J=n,s,t.It is clear thatKIJα is a diagonal matrix, which is as follows:(3)KIJα=Knα000Ksα000Ktα.In the global coordinate system, (2) can be expressed as(4)dfiα=kijαdδjαe.Then(5)kijα=BiIαK′IJαBjJα,Bα=-cosγ-sin⁡γcos⁡β-sin⁡γsin⁡β-sinγcos⁡γcos⁡βcos⁡γsin⁡β0-sin⁡βcos⁡β,where, dfiα and dδjαe are the increments of the local forces and the local movements in the global coordinate system. kijα is the stiffness matrix in the global coordinate system. Bα 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](6)Knα=Kn0Fnαl2n,Ksα=Ktα=Kr0Fnαl2n,where Fnα is the contact force in the normal direction. Kn0, Kr0, and n are material constants. l is the branch length between two grains. For two spherical grains, l is the same as the grain size d. According to Hertz-Mindlin’s theory [36], n=1/3.


### 3.2.2. Plastic Part

Stress dilatancy is a well-known phenomenon in RFMs [37–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(7)FnαdSnαp+TαdSαp=Fnαtan⁡φdSαp.Then(8)dSnαpdSαp=tan⁡φ-TαFnα,where dSnαp 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 Tα and the increase of the share plastic movement dSαp can be given by(9)Tα=Fsα2+Fnα2,dSαp=dSsαp2+dStαp2.The yield function is assumed to be of Mohr-Coulomb type(10)FFnα,κ=Tα-FnακSαp=0,where κ is the hardening function, which is defined by a hyperbolic curve (Figure 6) as follows:(11)κ=Kp0tan⁡ϕu·SαpFnαtan⁡ϕu+Kp0Sαp,where ϕu is the interparticle friction angle and Kp0 is the initial slope of the hyperbolic curve, which can be expressed as follows:(12)Kp0=kn0.0125+10fn.Figure 6
The curve of the hardening function in theκ-Sap plane.Figure 7
Calculation flowchart of the M-A model.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 dδiαp can be expressed by(13)dδxαpdδyαpdδzαp=dSαptan⁡φ-TαFnαFsαTαFtαTαBα.


## 3.2.1. Elastic Part

Suppose that the relationship of the local forcesFIα and the local movements δIα is nonlinear elasticity. The contact stiffness of the plane includes normal stiffness Knα and shear stiffness Ksα and Ktα. The elastic stiffness is given by(1)dFnα=KnαdΔnαe,dFsα=KsαdΔsαe,dFtα=KtαdΔtαe.The tensor forms of (1) are as follows:(2)dFIα=KIJαdΔJαe,I,J=n,s,t.It is clear thatKIJα is a diagonal matrix, which is as follows:(3)KIJα=Knα000Ksα000Ktα.In the global coordinate system, (2) can be expressed as(4)dfiα=kijαdδjαe.Then(5)kijα=BiIαK′IJαBjJα,Bα=-cosγ-sin⁡γcos⁡β-sin⁡γsin⁡β-sinγcos⁡γcos⁡βcos⁡γsin⁡β0-sin⁡βcos⁡β,where, dfiα and dδjαe are the increments of the local forces and the local movements in the global coordinate system. kijα is the stiffness matrix in the global coordinate system. Bα 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](6)Knα=Kn0Fnαl2n,Ksα=Ktα=Kr0Fnαl2n,where Fnα is the contact force in the normal direction. Kn0, Kr0, and n are material constants. l is the branch length between two grains. For two spherical grains, l is the same as the grain size d. According to Hertz-Mindlin’s theory [36], n=1/3.


## 3.2.2. Plastic Part

Stress dilatancy is a well-known phenomenon in RFMs [37–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(7)FnαdSnαp+TαdSαp=Fnαtan⁡φdSαp.Then(8)dSnαpdSαp=tan⁡φ-TαFnα,where dSnαp 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 Tα and the increase of the share plastic movement dSαp can be given by(9)Tα=Fsα2+Fnα2,dSαp=dSsαp2+dStαp2.The yield function is assumed to be of Mohr-Coulomb type(10)FFnα,κ=Tα-FnακSαp=0,where κ is the hardening function, which is defined by a hyperbolic curve (Figure 6) as follows:(11)κ=Kp0tan⁡ϕu·SαpFnαtan⁡ϕu+Kp0Sαp,where ϕu is the interparticle friction angle and Kp0 is the initial slope of the hyperbolic curve, which can be expressed as follows:(12)Kp0=kn0.0125+10fn.Figure 6
The curve of the hardening function in theκ-Sap plane.Figure 7
Calculation flowchart of the M-A model.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 dδiαp can be expressed by(13)dδxαpdδyαpdδzαp=dSαptan⁡φ-TαFnαFsαTαFtαTαBα.


## 3.3. Macro-Micro Relationship

The stress-strain relationship for the sample can be determined by integrating the behaviour of NRFM interparticle contacts in all orientations. 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–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:(14)dεij=Aik-1∑α=1Ndδjαlkα,Aik=∑α=1Nliαlkα,where dδjα is the incremental relative displacement between two contacted particles. lkα is the branch vector joining the centres of two contacted particles. N is the total number of contact orientations, according to Oda’s research [41], which can be expressed as(15)NV=12πd31+ee.Using this hypothesis, the mean force on the contact plane of a given orientation is(16)dfjα=dσijAik-1lkαV.The stress increment can be obtained by the contact forces and branch vectors for all contacts as follows [42, 43]: (17)dσij=1V∑α=0Ndfjαliα.When the contact numberN is sufficiently large in an isotropic packing, the summation of flexibility tensor in (14) can be written in integral form [44], given by (18)dεij=Aik-1N∫02π∫0πdδjαγ,βlkαγ,βξγ,βsin⁡γdγdβ,(19)Aik=N∫02π∫0πliαγ,βlkαγ,βξγ,βsin⁡γdγdβ,where ξ is a distribution density function, and in this paper, its value is 1/4π, which means that the sample is isotropic.


## 4. 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 (Figure8) was used to determine the interparticle parameters.Figure 8
Flowchart of the PSO algorithm.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 pbest, and the global best position of the entire particle swarm is denoted gbest. For each iteration, the velocity and particle location are expressed by(20)Vidk+1=ωVidk+c1r1Pidk-Xidk+c2r2Pgdk-Xidk,Xidk+1=Xidk+Vk+1,id,where Vidk and Xidk are the velocity and the position of the ith particle in the previous iteration, respectively. r1 and r2 are random numbers ranging from 0 to 1. c1 and c2 are acceleration coefficients, which are the positive parameters; here, they are equal to 2. Pidk and Pgdk are pbest and gbest, 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(21)Er=1n∑i=1nx-xt2,wherex is the calculated value, which represents the axial strain and volume strain, respectively; xt is the test data; and n is the total number of test data. As the number of search steps increases, Er will gradually become smaller and satisfy the end conditions.Table3 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 Er 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.Table 3
Parameter-search range of the M-A model.
Material K n 0 (N/m) Φ u (°) K n a (N/m) φ (°)  Nuozhadu rockfill 2.0 × 10 8 – 3.0 × 10 8 39-40 1.5× 107–2.5 × 107 35.6–36.6Under unloading-reloading conditions, the stress-strain curve is nearly a straight line during the test (Figure14). 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 Kna, a significant phenomenon was found: parameter Kna varies regularly with the mean principal stress P1 at the unloading point. After plotting lg(Kna/Pa)-lg(P1/Pa) 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:(22)Knα=Kn1PaP1Paw,where Kn1 is the intercept of a straight line and w is the slope. The normal stiffness Kna 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)).Figure 9
Relationship between normal stiffness and the principal stress at the unloading point.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 Figure15). 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–48]. From the analysis of the fitting parameter φ, we found that the dilatancy angle φ decreased as the transitional confining pressure σ3t increased. As shown in Figure 10, the dilatancy angle φ and the transitional confining pressure σ3t exhibit a good linear relationship under the logarithmic coordinates. Here, the dilatancy angle of the contacted particles φ can be determined as follows:(23)φ=φ0-Δφlg⁡σ3tPa,where φ0 and Δφ are parameters whose values are equal to the intercept and slope of the line, respectively.Figure 10
Relationship between the dilation angles and the transition confining pressures.


## 5. 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 thex-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 x-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)).Figure 11
Stress-dilatancy relations for selected contact plane orientations for NRFM.
(a) (b) (c) (d) (e) (f) (g) (h)  (i)In order to study the local dilatancy behaviour of the particles, the triaxial test of the type two with the transition confining pressureσ3t=300 kPa (see Figure 4) was chosen to be calculated. According to the employed local dilatancy equation (see (8)), the local dilatancy rate dSnap/dSap 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 Figures11(b)–11(i), respectively. During the CSR stage with q/p=1, as the load increased, the local dilatancy rate also increased. All the values are positive, which indicate that the contact planes are contracted. The 65° contact plane reaches the highest mobilized ratio (close to π/4+Φu/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.Figure12 shows the relationship of the local dilatancy rate dSnap/dSap and the global stress ratio q/p. 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].Figure 12
Local plastic movement increment versus global stress ratio for selected contact plane orientations.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).
(a) (b) (c) (d)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.
(a) (b)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.
(a) (b)


## 6. Comparison and Verification of Constitutive Models

### 6.1. 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 σ3t=800 kPa (Figure 13). The comparisons of the predicted results are as follows.As shown in Figures13(a) and 13(c), the shapes of the q~ε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 q~ε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 Figures13(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.


### 6.2.  M-A Model Predictions

Under CSR loading, the parameters of the extended M-A model were calibrated by using the test data ofσ1/σ3=1.5, 2.5, and 4.0. Equation (22) was chosen to determine the normal stiffness Knα 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 D50=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].Table 4
Parameters of the M-A model.
(a) K n 0 (N/m) n K r 0 / K n 0 Φ u (°) 2.5× 108 0.33 0.5 39.17°  (b) K n 1 (N/m) w φ 0 (°) Δ φ (°) 3.98 × 10 4 0.52 36.91 ° - 0.99Table 5
Parameters of Duncan-Chang’s E-u and E-B models.
Material Φ (°) ΔΦ (°) R f K N D G F K b m Nuozhadu rockfill 54.37° 10.47° 0.72 1491 0.24 7.65 0.41 0.2 683 0.10The comparisons between model predictions and type one test are shown in Figure14 in terms of the stress-strain-volume relationship on NRFM. The unextended M-A model overestimates the deviatoric stress especially under unloading-reloading conditions (Figure 14(a)), while the extended model proposed in this paper could well predict the stress-strain-volume behaviour of NRFM under CSR loading (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 unloading-reloading 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 (Figure 14(b)). Because little changes in the volume occurred during this stage, the calculated volume changes are neglected.Figure15 shows 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 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.


## 6.1. 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 σ3t=800 kPa (Figure 13). The comparisons of the predicted results are as follows.As shown in Figures13(a) and 13(c), the shapes of the q~ε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 q~ε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 Figures13(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.


## 6.2.  M-A Model Predictions

Under CSR loading, the parameters of the extended M-A model were calibrated by using the test data ofσ1/σ3=1.5, 2.5, and 4.0. Equation (22) was chosen to determine the normal stiffness Knα 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 D50=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].Table 4
Parameters of the M-A model.
(a) K n 0 (N/m) n K r 0 / K n 0 Φ u (°) 2.5× 108 0.33 0.5 39.17°  (b) K n 1 (N/m) w φ 0 (°) Δ φ (°) 3.98 × 10 4 0.52 36.91 ° - 0.99Table 5
Parameters of Duncan-Chang’s E-u and E-B models.
Material Φ (°) ΔΦ (°) R f K N D G F K b m Nuozhadu rockfill 54.37° 10.47° 0.72 1491 0.24 7.65 0.41 0.2 683 0.10The comparisons between model predictions and type one test are shown in Figure14 in terms of the stress-strain-volume relationship on NRFM. The unextended M-A model overestimates the deviatoric stress especially under unloading-reloading conditions (Figure 14(a)), while the extended model proposed in this paper could well predict the stress-strain-volume behaviour of NRFM under CSR loading (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 unloading-reloading 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 (Figure 14(b)). Because little changes in the volume occurred during this stage, the calculated volume changes are neglected.Figure15 shows 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 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.


## 7. 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 stress-strain 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.



---

*Source: 7630541-2016-11-22.xml*