APPLICATION OF A NONLINEAR DILATANCY MODEL IN THE STRESS AND DEFORMATION ANALYSIS OF A HIGH CONCRETE-FACED ROCKFILL DAM

Original scientific paper The concrete-faced rockfill dam (CFRD) is a widely applied geotechnical structure. The stress and deformation conditions of CFRD are absolutely critical for the safety of dam. However, only a few constitutive models have been reported to appropriately simulate and predict the specific stress and deformation of high CFRDs. In this paper, a new nonlinear dilatancy constitutive model was introduced to analyze CFRDs. Based on the generalized potential theory, the stress-strain matrix [Dpq] for the nonlinear dilatancy mode was also demonstrated. The numerical simulation result by the nonlinear dilatancy model was compared with that by a large-scale tri-axial compression test. It is concluded that the behaviours of rockfill materials can be well simulated by the nonlinear dilatancy model. In addition, the nonlinear dilatancy model was applied to simulate the step-by-step construction and impounding of the reservoir process in the Chinese Shuibuya CFRD. From comparative studies of the numerical results and in-situ monitoring data, a good agreement can be explored between the measured and computed settlements. This result also indicates that a 3D finite element analysis based on the nonlinear dilatancy model can be applied to predict the stress and deformation of CFRDs.


Introduction
The concrete-faced rockfill dam (CFRD) is a form of rockfill dam consisting of cushion, transition, main rockfill, and downstream rockfill zones [1].In recent decades, CFRDs have been widely applied for their multiple advantages, such as good safety, better adaptability to topography and geology, application of local materials, cost-effectiveness, simple construction and short construction periods [2].Especially in recent years, many high CFRDs have been constructed, such as the Shuibuya CFRD with a height of 233m, and the Hongjiadu CFRD with a height of 179m in China [3][4][5].In the construction of high CFRDs, the deformation of the rockfill zone and the stress of the concrete slab in the dam are key factors for the dam's safety.Thus, engineers and researchers pay particular attention to these issues during the design of high CFRDs.various numerical methods are involved to analyze and solve problems in the construction and operation of the CFRD.The finite element method (FEM) is a powerful tool to predict the deformations of dam and slab joints and analyze the stress distribution of the slab, and have been applied in many studies of CFRDs [6 ÷ 8].However, the reliability of FEM results significantly depends on the adaptability of the specific constitutive model for certain rockfill materials.
Many constitutive models have been proposed to predict the mechanical behaviours of rockfill materials.
Duncan et al. [9] demonstrated a hyperbolic model for the nonlinear stress-strain relationship of soil, which is often applied for stress and deformation analysis of rockfill dams.Naylor et al. [10] introduced a nonlinear K-G model for FEM analysis in geotechnical engineering and the construction process of the Beliche dam was predicted with this method.Shen [11] proposed a double yield surface model for stress-strain analysis of rockfill dams in China.Lollino et al. [12] demonstrated a Lade model of sand, with which the behaviour of rockfill was mimicked in their plane-strain FEM analysis.Xu et al. [14,15] proposed a generalized plasticity model for rockfills and simulated the construction process of the Zipingpu CFRD.Varadarajan et al. [16] attempted to apply the disturbed state concept model to capture the stress-strain behaviours of rockfills.Xiao et al. [17,18] proposed both a 3D bounding surface model and a state-dependent mode to explore the behaviours of rockfills.A new double yield surface model was introduced to characterize the strength and deformation behaviours of rockfills [19].Many studies are derived from practical tests, thus the basic deformation properties of rockfills can be described and reflected.Due to the relatively simple parameters and clear physical meanings, the Duncan-Chang E-B model is one of the most widely applied models in CFRD construction simulation.However, in this particular model the dilatancy and plastic deformation of rockfill materials are not taken into account.Additional details should be considered to describe the constitutive models for rockfills.Furthermore, it is necessary to verify the feasibility of the constitutive model in terms of the description of the behaviour of rockfill materials and reproduction of the slab deformation and stress through numerical simulation of high CFRDs construction using model parameters of its rockfills.
In order to predict the deformation and stabilization of high CFRDs better, a new nonlinear dilatancy constitutive model for rockfill is introduced in this paper.Based on the generalized potential theory, the stress-strain matrix [D pq ] for the nonlinear dilatancy mode is also demonstrated.The nonlinear dilatancy model is employed to predict the stress-strain behaviour of rockfill materials in stress space.The numerical simulation results are compared with data from a laboratory large-scale tri-axial compression test.Based on the above nonlinear dilatancy model, taking Shuibuya CFRD as an example, numerical simulation is conducted.Finally, the numerical results are compared with in-situ monitoring data.

Nonlinear shear dilatancy model 2.1 The stress-strain relationship in the dilatancy model
Considering the effects upon volume variation from shear stress, the incremental form of the stress-strain relationship could be expressed as [20,21]: Eqs. ( 1) and ( 2) can be inverted as follows: (3) (4) where K p is the bulk modulus, K q is the modulus of dilatancy, G is the shear modulus, p is the mean stress, q is the generalized shear stress, and ε v and ε s are bulk strain and generalized shear stain, respectively.
With the tests results of triaxial compression tests, the following equations can be obtained: (5) where σ 1 is major principle stress in the vertical direction, σ 3 is minor principle stress in the horizontal direction, ε 1 is axial strain and ε 3 is horizontal strain.
In order to clarify the relationship between the parameters of K p , K q , G and the strain state, it is assumed that the bulk strain in soil is composed of elastic strain and dilatancy strain.They are marked with the superscript e and d respectively, as follows: (6) (7) Generalized Hooke's Law defines the relationship between elastic stress and strain, while dilatancy strain obeys the Rowe Dilatancy Equation [13].(8) where K f is the minimum bulk energy ratio of the soil, which was measured in the test.
The elastic Poisson's ratio μ is presumed to be a constant.Based on conventional triaxial compression tests, the derivation of relations among the parameters of K p , K q , G and the stress state will be illustrated in the sections that follow.

Bulk modulus Kp
The elastic bulk strain and stress p are linked by bulk modulus K p , in the form of Hooke's Law: ( ) where μ is elastic Poisson's ratio, and E ur is Young's modulus of soils, just as in the Duncan-Chang model, where parameters K ur and n are measured in a cyclic loading test, P a is atmospheric pressure.

Dilatancy modulus K q
Referring to the Rowe dilatancy equation, under conditions of triaxial stress, the dilatancy strain can be described as: (11) Similarly, the vertical strain of the triaxial test can be divided into vertical elastic strains and vertical dilatancy strains, written as: (12) where E t is the tangent modulus, under triaxial conditions, it can be described as the following.(13) Based on Eqs.(11), (12) and ( 13), the dilatancy modulus K q can be described as: (14)

Tangent modulus E t
Tangent modulus E t is included in the expression of both dilatancy modulus Eq. ( 14) and shear modulus Eq. ( 17), just as it is shown in the Duncan-Chang model as follows: ( ) where S is the stress level, it can be described as: Based on the results of various tests, at certain times, Eq. ( 18) needs to be modified when it is applied to coarse grain soils.The modified equation can be expressed as: where K is modulus number, n is the exponent for defining the effects of confining pressure on the initial modulus, R f is failure ratio, c is cohesion, φ is internal friction angle, P a is atmospheric pressure, and α and β are material constants, compared with the Duncan-Chang model, which in the expression of E t , α=1 and β=2 are determined.By fitting the whole E t −S curves under different confining pressures, the parameters α and β can be obtained.
For coarse grain soils, cohesion c is generally taken to have the value of 0. However, the friction angle φ is changed with confining stress, where 0 ϕ and ϕ ∆ are constant values.There are 11 parameters involved in the nonlinear shear dilatancy model in total, and they can be listed as c, φ, φ 0 , K, n, R f , α, β, K ur , K f , and μ.All of these parameters can be achieved by conventional triaxial compression tests.

Stress-strain matrix [D pq ]
Based on mathematical principles, Yang proposed the generalized potential theory for the soil constructive model and the corresponding multi-potential surface constructive model [22].The traditional soil constitutive models are established in the p−q space, and the effects of lode angle and rotation of principle stress are ignored.The general relationship between plastic-strain increment and stress increment can be expressed as: (22) where A, B, C and D are plastic coefficients.When there is only one total plastic strain direction, the relationship among plastic coefficients must meet the following condition: 0 = − BC AD .Based on multi-potential surface theory, the following relation can be established: (23) where k φ (k = 1, 2, 3) is three mathematical potential functions which are independent of the linear gradient vector, only in the (p, q) space, it can be treated as a twodimensional issue, and in this case, k = 1, 2, the simplest potential functions can be expressed as: φ and 2 φ are substituted into Eq. ( 23) and the following Eq.( 24) can be achieved: Then, the definitions of p and q are substituted into Eq.( 24): (25) From Eqs. ( 22), ( 24) and ( 25), the constitutive relationship can be expressed: (26) In addition, the incremental stress matrix is: (27) [D e ] is the elastic matrix, and spherical stress increment dp and deviation stress increment dq can be expressed as: (28) Based on Eqs. ( 26), ( 27) and (28), Eq. ( 29) can be obtained as the following: The description of Eq. ( 29) is then changed into the following equation: (30) where: 11 22 12 21 A a a a a = − Based on Eq. ( 26), ( 27) and (30), the following Eq.(31) can be achieved as: (31) where: Referring to the above decomposition principle, the volumetric strain could also be divided into the elastic part and plastic part.
The volumetric strain increment induced by elastic strain can be expressed as: The volumetric strain increment induced by plastic strain can be expressed as: (33) So far, the volumetric strain induced by the elastic strain can be obtained by the generalized Hooke's Law, and for volumetric strain induced by plastic strain, it is only related to dilatancy strain, comparing Eq. (33) with Eq. ( 22): ≠ is enough to meet the condition of 0 AD BC − = , and the uniqueness of plastic strain increment can also be fulfilled.Now, the matrix is improved to become a nonsymmetric matrix.
When Eq. ( 34) is substituted into Eqs.( 30) and (31), the new expression for stress increment can be obtained as follows: (35) where: (36) Thus, the calculating matrix [D pq ] can be obtained for the nonlinear dilatancy constitutive model, showed as follows: (37) According to Eq. (37), it is indicated that, when there is no observed dilatancy and K q = ∞ , under this condition, [D pq ] would be regressed back to a symmetric matrix in a dual-modulus model.

Numerical simulation test
The 3D constitutive relationship and its matrix form under confining pressure are validated with numerical simulation.ABAQUS software was re-developed by applying the nonlinear shear dilatancy (NSD) model.The calculation results from numerical simulations with the nonlinear shear dilatancy (NSD) model and Duncan-Chang E-B model (E-B) are compared with that from the triaxial draining compression test for a certain rockfill material.
This rockfill material is mainly composed of sandstone, and the largest particle size in the test is 60 mm.The experimental dry density is 2,27 g/m 3 , while the relative density is 0,9.
The values of model parameters are obtained from the above triaxial draining tests and presented in Tab. 1.These parameters are then applied to simulate the triaxial tests.A plan view and typical cross section of the dam are shown in Fig. 3 and Fig. 4 respectively.The upstream dam slope is 1:1,4 and the downstream dam slope is 1:1,25.There are four access berms at the downstream section, with a width of 4,5 m.The dam is comprised of six rockfill materials, including IA, IIA, IIIA, IIIB, IIIC and IIID.The boundary of the main and sub rockfill zone begins at an altitude of 380 m in the axes of the dam while it ends at the downstream bottom of the dam with a slope of 1:0,2.Zoning and construction details are listed in Tab.The concrete face is divided into 55 slabs and there are 27 vertical joints with 16 m spacing and 28 vertical joints with 8 m spacing.The face slab is constructed with C25 concrete.The reinforcement is placed in the centre of thickness of the slab.The steel ratio is 0,4 % along the slope and 0,35 % in the horizontal direction.
The construction phases and time of each phase are provided in Tab. 3

Settlement monitoring system
The Shuibuya CFRD is one of the most significant engineering projects in China, and consequently an improved and detailed settlement monitoring system has been established to monitor the deformation of the Shuibuya CFRD.Settlement and horizontal displacement inside the dam body are measured by water tube settlement gauges and alignment displacement sensors.They are set up in three important cross-sections, 0 + 132, 0 + 212 and 0 + 356, in which section 0 + 212 is the maximum cross-section of the Shuibuya CFRD.During the filling of the dam body, these gauges and meters were set up from bottom to top.Five monitoring lines are placed in section 0 + 212 and the specific locations are at elevations of 235 m, 265 m, 300 m, 340 m and 370 m, respectively.Three monitoring lines are placed in section 0 + 132 and 0 + 356, with the specific locations being at elevations of 300 m, 340 m, and 370 m, respectively.Therefore, there are 11 monitoring lines and 72 monitoring points in total.The locations of all devices at the maximum crosssection 0 + 212 are shown in Fig. 6 and the settlement gauge positions are listed in Tab. 6.In addition, the stress and displacement of the slabs and the displacement of joints are also measured in this system.34 strain gauges and 74 reinforcement sensors are installed in the slabs on sections 0 + 132, 0 + 212, 0 + 356 and 0 + 448.The stress of the slabs can be observed.46 joint sensors are set up to observe the displacement of vertical joints, and 13 joint meters are installed to observe the displacement of the peripheral joints.

FEM analysis of stress and deformation of Shuibuya CFRD 4.1 Computational model and parameters
A three-dimensional finite element analysis of Shuibuya CFRD was conducted by ABAQUS software.In this finite element model, the dam body and concrete face are composed of 101,870 elements and 11,646 elements respectively.Most elements involved are eight nodes of hexahedral elements while some prism elements are applied in the boundary zone.The three-dimensional finite element mesh of the Shuibuya CFRD is shown in Fig. 7.The results obtained from this computational procedure are exactly conformed to the actual construction schedule.Firstly, blocks (1) to ( 5) of the dam body up to EL. 319,2 m are simulated.In each block, layer-by-layer elements are reactivated to simulate the process of construction.The thickness of each layer is less than 6m.Before the placement of stage-I concrete slab, an algorithm for slope paving [23] is applied for the upstream slope of EL. 319,2 m.This allows the freshly cast concrete slab to be kept in strict contact with the special cushion zone.Then the element of the stage-I concrete slab ( 6) is activated to simulate the dam construction.The procedure is repeated until the whole dam body is completed.The impoundment process is simulated by increasing water levels and water pressure is simulated as a surface force on the slabs.The bottom boundary of the dam is fixed in three directions (X, Y, Z).A nonlinear shear dilatancy model is applied to analyze the Shuibuya rockfill materials.The maximum particle size of the Shuibuya rockfill materials is 800 mm, and it is impossible to determine the model parameters by triaxial tests on the actual materials.The actual main rockfill and sub rockfill have a dry density of 2,18 g/m 3 and 2,15 g/m 3 after compaction, respectively.In this study, the largest particle size of the rockfill materials in the test is 60 mm.However, the studies have shown that smaller rockfill particle size results in higher shear strength and higher rockfill stiffness [26].On the basis of this consideration, a smaller density of rockfill materials is chosen to determine the model parameters.The dry density of main rockfill and sub rockfill tested is 2,08 g/m 3 and 2,05 g/m 3 , respectively.The lower dry density results in a lower modulus and shear strength of the rockfill [27].The model parameters are given in Tab. 4. Goodman contact elements [24], a kind of nonthickness element, are involved between the concrete face slab and cushion material, the parameters for the interface are also determined by laboratory tests and presented in Tab. 5, with modulus factor (k 1 ), modulus exp (n 1 ), failure ratio (R f1 ), friction angle (φ 1 ) and cohesion (c 1 ) also included.The Goodman contact elements with certain thickness are applied to simulate slab joints and peripheral joints.In Xu and Zou [14], it was reported that the shear stiffness of slab joints had little effect on the slab stresses.Therefore, the stiffness of the joints in two shear directions is assumed to be 10 MPa/m.The compressive stiffness of the joints is assumed to be 25000 MPa/m (the compressive stiffness of the wood in the joints).The concrete face slab is simulated by the linear elastic model, in which elastic modulus E = 25500 MPa, Poisson's ratio υ = 0,167 and density ρ = 2,40 g/cm 3 .

Stress and deformation of the rockfill
The simulated settlements of sections 0 + 212, 0 + 132 and 0 + 356 in the construction stage are compared with that of in-situ monitoring results, as shown in Fig. 8(a) ÷ 8(c).The calculated settlements and in-situ observed settlements at elevations of 235 m, 265 m, 300 m, 340 m and 37 0m of cross section 0 + 212 are demonstrated in Tab.6(a) ÷ 6(e).According to the comparison, the numerical results of the settlement are in good agreement with observed results.From Tab. 6(a) ÷ 6(e), the predicted settlements are smaller than in-situ monitoring results.This difference may result from the creep deformation of rockfill materials and the settlement of bedrock foundation.
The settlement and horizontal deformation contours of maximum cross section 0+212 under normal water level are demonstrated in Fig. 9.It is clear that the settlement of the dam during the normal water level period is 217 cm, which is approximately 0,93 % of the dam height.The maximum settlement is located at the middle height of dam body.The maximum horizontal displacement values towards upstream and downstream are 23,2 cm and 86,7 cm, respectively.Compared with the settlement, the calculated values of horizontal displacement are relatively small.Similar results are also reported by Zhou [3,4]

Stress and deformation of the concrete face slab
The response of the CFRD during reservoir impoundment was investigated by the proposed numerical procedure.The deflection and along-slope displacement of the concrete face slab at the normal water level stage are plotted in Fig. 10.The maximum bending deflection is 62,2 cm at EL. 233,23 m.In the slope direction, downward displacement is observed at the upper layer of the face slab and upward displacement is observed at the bottom of the face slab.There is a compressive state in the whole face slab.Both the stress along the slope direction and dam axes are explored in Fig. 11.The maximum compressive stress along the slope direction is 15,58 MPa, which is located at the bottom of the face slab near the river valley.The maximum axial compressive stress of the slab is 19,78 MPa.This is observed near the middle part of the riverbed at the height of about 1/4 of the dam.The tension stress zone only distributes near the bottom of the face slab at the two banks, with the maximum value of 0,85 MPa.There is a good consistency between the predicted stress distributions and the in-situ observed results [25].model.The application of a tensor form and matrix form for a nonlinear dilatancy model was also demonstrated, with a finite element analysis.Numerical analysis of the triaxial compression test was carried out for the rockfill of the Shuibuya CFRD.The predicted results of the nonlinear dilatancy model are able to achieve good agreement with laboratory tests results.This indicates that the proposed model could be applicable to elucidate the stress and strain variations of rockfill materials in general.In particular, dilation/contraction features can be captured precisely.Thus, the matrix form [D pq ] for the nonlinear dilatancy model is proved to be reasonable.
Based on the above nonlinear dilatancy model, and taking the Shuibuya CFRD as an example, a numerical simulation was conducted.The stress and deformation during the construction and reservoir filling processes were simulated.Compared with the available in-situ strain monitoring data, the simulated deformation trends and values were in good agreement with the in-situ data.The slab stress during reservoir filling were also analyzed and discussed.
From the above results, it can be concluded that this nonlinear dilatancy model for rockfill materials is valid for predicting the deformation of CFRD during the construction and reservoir filling stage.
The numerical analysis does not consider the creep deformation of rockfill, which may underestimate the deformation of the dam and the stresses of the slab, especially during the reservoir filling stage.Therefore, the creep deformation of Shuibuya CFRD will be subsequently investigated by numerical procedure.

Figure 1
Fig.1presents the test results of the stress-strain behavior and volume change of rockfills, together with the predicted results.The NSD model predictions closely match the test data, and those of Duncan-Chang E-B model cannot predict mechanical behaviours such as the volumetric strain expansion.This means that the NSD model can capture the experimental trend, and the obvious deviation occurs between the test data and Duncan-Chang E-B model.They show that volumetric dilation/contraction of rockfill materials can be appropriately described by the NSD model, and it indicates that the 3D constitutive relationship and its matrix form for the NSD model are reasonable.

2 .
The thickness of the concrete face is 0,3 m at the top with an elevation of 405 m and 1,1 m at the bottom with an elevation of 177 m.The thickness of the concrete face obeys the relation: d = 0,3 + 0,0035•H, where d is the thickness of the concrete face and H is the vertical distance below the top elevation of 405 m.

Figure 2 A
Figure 2 A general view of the Shuibuya CFRD

Figure 3 Figure 4
Figure 3 Plan view of the Shuibuya CFRD and Fig. 5.The dam's construction began in January 2003, with the filling of the dam and pouring of the face slab completed in October 2006 and March 2007 respectively.The concrete slab was poured in three phases.The top elevations for the three phases are EL.278 m, EL. 340 m and EL.405 m respectively.The water impoundment level reached EL. 400 m in March 2008.

Figure 5
Figure 5 Construction phase of the Shuibuya CFRD

Figure 6
Figure 6 Locations of devices at maximum cross-section 0 + 212

Figure 7
Figure 7 Three-dimensional FE mesh of the Shuibuya CFRD

Figure 9
(a) Settlement (unit: cm) (b) Horizontal displacement (unit: cm) Displacement contours of cross section 0+212 at the full storage stage (a) Normal (Normal face upwards indicated positive) (b) Slope direction (Sloping upwards indicated positive) Figure 10 Contour lines of slab displacement at the full storage stage (a) Horizontal (minus value indicates compression) (b) Slope direction (minus value indicates compression) Figure 11 Contour lines of slab displacement at the full storage stage 5 Conclusion A new nonlinear dilatancy constitutive model for rockfill of a CFRD was introduced in this study.The advantages of the hyperbolic E-B model and Rowe Dilatancy Equation were incorporated in the proposed

Table 1
Material parameters of the tested material

Table 2
Materials and construction details of the Shuibuya CFRD

Table 3
Construction phase and time

Table 4
Material parameters in NSD model for the Shuibuya CFRD.

Table 5
Goodman element parameters of the concrete rockfill for the . Additionally, the major principle stress and minor principle stress are 4,21 MPa and 1,86 MPa, respectively.

Table 6 (
c) Comparison of settlements at EL. 300 m (in October 2006)

Table 6 (d)
Comparison of settlements at EL. 340 m (in October 2006)

Table 6 (
e) Comparison of settlements at EL. 370 m (in October 2006)