Abstract
An intelligent algorithm that simultaneously analyzes multiple rockfill parameters is proposed. First, the paper introduces the operation and monitoring condition of the Shuibuya concrete-faced rockfill dam (CFRD). Then the constitutive rockfill models and the FEM analysis procedure are introduced in this paper. Third, the MsPSO intelligent algorithm was adopted to inverse the rockfill parameters. The recalculated displacement of Shuibuya CFRD using the inversed rockfill parameters is presented, and a satisfactory result was obtained, indicating that the inversion method is correct and effective. The method developed in this paper can be adopted in any geotechnical engineering parameter inversion.
1. Introduction
Concrete-faced rockfill dams (CFRDs) have rapidly developed in recent decades due to their strong adaptability to topography, geology, and climate and easily available construction materials. CFRDs were constructed with increasingly tall heights in recent years [1, 2]. Rockfill is used as the support structure and maintains the stability of CFRDs. Creep deformation accounts for the majority of rockfill dam postdeformations. Large rockfill creep deformation causes separation between the slabs and the cushion layer or even cracks in concrete slab [3], thus impairing the integrity of seepage control systems and weakening the structure or even threatening the safety of the dam [4]. To avoid this potential risk in tall CFRDs, rockfill creep problem must be taken seriously; it is necessary to understand the deformation characteristics of rockfill to make the rockfill compatible with the concrete face in order to improve the design of tall CFRDs.
Currently, the common method for studying the character of rockfill creep references laboratory creep tests [5–8]. At first, the laboratory creep tests were conducted under low confining pressure. Early in 1985, Parkin [5] presented rate methods and performed laboratory rockfill creep tests on the odometer. In China, the earliest rockfill creep test was completed with stress-controlled triaxial equipment by Shen [6] in 1991, and the exponent creep curve model was the first creep model type. As rockfill creep under low confining pressure laboratory tests cannot reflect the behavior of 200 m high CFRDs, Cheng and Ding [8] proposed a power function rockfill creep curve model under high confining pressure. In their study, many groups of laboratory triaxial creep tests with different confining pressures were performed; the creep deformation over time can be expressed by a power function. Until recently, the power function creep curve model has been applied widely for the creep calculation of tall CFRDs.
It is known that the in situ geomechanical properties of rockfill may differ from those determined by laboratory tests; one of the most obvious differences is called the “scale effect” [9, 10]. Therefore, to understand the in situ geomechanical properties of rockfill, it is better to perform back-analysis based on constructed CFRDs. Displacement inversion is an effective method to check and modify the parameters of rockfill models. Moreover, the back-analysis results are beneficial to the safety assessment of the project itself and to improving the design of subsequent structures of the same type.
Intelligent algorithm such as the artificial neural network (ANN) and the particle swarm optimization (PSO) algorithm were proposed to solve the geotechnical engineering inversion problem in recent years and achieved good results [11–16]. Generally, the optimal values of the parameters to be determined are progressively approximated through iteration by minimizing the error function. Because of the complexity of problems in geotechnical engineering, numerical calculation is often used; the time-consuming finite element method (FEM) calculation is performed frequently in the inversion, so the rate of convergence is slow and sometimes the inversion fails due to nonlinear problems on a large scale. Furthermore, the result is affected by the initial values, and a local minimum or premature convergence is likely to be obtained; as a result, the solution is sometimes unstable.
Zhou et al. reported the analyses of actual measured deformations resulting from continuous monitoring of the Shuibuya CFRD, which is currently the tallest CFRD in the world [14]. The actual settlement records and laboratory test deformation parameters of the main and secondary rockfill were introduced. In this paper, a displacement inversion for rockfill parameters is performed using the multiswarm particle swarm optimization (MsPSO) algorithm [17]. This work focused on vertical displacements, and the available monitoring data span nearly 5 years, covering the construction phase, the first filling stage of the reservoir, and 2 years of its operational time [14].
2. The Shuibuya CFRD
The Shuibuya CFRD was constructed in Badong City in Hubei Province, China, which is 117 km from Enshi City, located upstream, and 92 km from the Geheyan water-power plant downstream. The normal water level of the reservoir is 400 m, the maximum dam height is 233 m, the upstream dam slope is 1 : 1.4, and the downstream integrated dam slope is 1 : 1.4. The dam comprises five rockfill materials, including bedding material (IIA), transition material (IIB), primary rockfill (IIC), secondary rockfill (IID), and downstream rockfill (IIE). The boundary between the primary and secondary rockfill zones begins at an elevation of 380 m in the axes of the dam and ends at the downstream bottom of the dam with a slope of 1 : 0.2. The construction of the dam began in 2002 and was completed in 2007 and the filling of the reservoir was conducted in several sequential steps over the period from October 2006 to November 2008. The concrete dam face was constructed over three periods: first stage (January 2005 to March 2005), second stage (January 2006 to April 2006), and third stage (January 2007 to March 2007). The maximum cross section and the construction process of the Shuibuya CFRD are shown in Figure 1, where the solid lines represent the material partition and the dotted lines represent the construction partition.
The vertical displacements inside the dam body were measured by hydraulic overflow settlement gauges. They are distributed in three important cross sections: 0+132, 0+212, and 0+356. More specifically, a total of 11 monitoring lines were placed in the three cross sections, five in section 0+212 (at elevations of 235 m, 265 m, 300 m, 335 m, and 370 m) and three each in section 0+132 and section 0+356 (at elevations of 300 m, 335 m, and 370 m). The layout of the in situ instrument measurement scheme of the 0+212 cross section is shown in Figure 2 [14].
3. Displacement Calculation
The displacement of Shuibuya CFRD was calculated using the FEM. The instantaneous elastic-plastic deformation of rockfill materials was modeled using the Duncan - model [18], and the time-dependent creep deformation was modeled using the power function creep curve models [8].
3.1. Outline of - Model
The - model adopts an incremental Hooke’s law to describe the nonlinearity of rockfill materials and Mohr-Coulomb’s law to model the shear failure. This model has two moduli, deformation modulus and bulk modulus .
The tangent modulus can be expressed as follows:where is failure ratio, is initial deformation modulus, which expresses the compression-hardening property of rockfill materials, and is stress level. These parameters can be expressed as follows:where is the standard atmospheric pressure, is modulus number, is exponent for defining the influence of the confining pressure on the initial modulus, is cohesive strength, is friction angle, and , are major and minor principal stresses, respectively.
For rockfill materials, the cohesive strength is generally taken as zero. However, the angle of internal friction changes with compression stress:where and are two constants for rockfill.
The bulk modulus is given bywhere is bulk modulus number and is bulk modulus exponent.
The - model parameters of the main and secondary rockfill of Shuibuya CFRD were introduced by Zhou et al. [14] and are shown in Table 1.
3.2. The Power Function Creep Curve Model
The power function creep curve model in [8] is employed in this paper. This model was proposed through the creep test results which showed that the correlation of the axial creep and time can be fitted by the power function in the double logarithmic:where is the limit of axial creep strain and is a coefficient related to the axial creep strain rate. The time is calculated in hours. The limit of axial creep strain that combines both the confining pressure and the stress level viscous responses is given by the following relationship:where is confining pressure, is proportional to , the coefficient is the hyperbolic function of , and and are the creep parameters.
According to the creep test fitting curve between under various and , is independent of , and the relation between and can be expressed by the power function as follows:
Thus, the axial creep model of rockfill can be described by (5)–(9) and the creep parameters , , , and .
The relation of the volume creep strain can also be expressed by the power function of time as follows:where is the limit of volume creep strain and is a coefficient related to the volume creep strain rate. The limit of axial creep strain is defined aswhere and are both functions of the stress level, defined by the following relationship:
The limit of volume creep strain can be written aswhere is assumed as a constant: .
Thus, the volume creep model of rockfill can be described by (10)–(13) and the creep parameters , , , , and .
The creep test parameters of the main and the secondary rockfill of Shuibuya CFRD are shown in Table 2.
In three-dimensional calculations, the confining pressure is replaced by , which equals .
According to (5) and (10), the creep strain ratio has the following forms:
The time history is divided into several increments , and and are obtained through summation:
Therefore, the increments of shear creep and volume creep are given by the following, respectively:
The relation among shear creep strain , axial creep strain , and volume creep strain under ordinary conditions of triaxial test is given as follows: The above equation is also suitable for , , and .
If the Prandtl-Reuss flow law is adopted on the plane, the strain velocity tensor can be expressed aswhere is a strain velocity tensor, is a unit tensor, is the deviatoric stress tensor, and is deviator stress, which is given by
Then the creep increment tensor can be written by
So far, the numerical simulation method has been completely derived for CFRDs [15, 19]. The initial strain technique [20] is used to implement this creep model. The equivalent nodal force caused by rockfill creep can be calculated using the following equation:
By loading the equivalent force on the corresponding nodes, the displacement caused by rockfill creep can be calculated.
3.3. The Finite Element Calculation
In this paper the discretized mesh for the FEM analysis is composed of 8318 elements and 9138 nodes.
The interface between the slab and the cushion was simulated using Goodman contact elements, which were also applied to simulate the slab joints and the peripheral joints. The concrete face slab was simulated using a linear elastic model with an elastic modulus of 25.5 GPa and a Poisson ratio of 0.167 [21].
The process of construction and water storage is divided into 30 loading steps based on the real construction schedule, as shown in Figure 1. Because the creep deformation of the rockfill occurs while the load is applied, the creep is analysis for the whole period of the CFRD construction and operation process. Basically, each loading step is followed by an equivalent loading step for creep deformation calculation. As a result, the creep process consisted of 29 steps. Thus, there are 59 steps used during FEM analysis, as shown in Figure 3.
This paper used the static finite element calculation program developed by Ganchen Gu et al. with Fortran software and added with independent development of a creep deformation equivalent loading calculation module. The use of the program to calculate dam creep deformation has been previously published [22].
4. Parameter Sensitivity Analysis
The variability of strength parameters , and damage ratio are small for many engineering experiences because many methods have been developed to calculate them [23] and the method is essentially perfect; therefore, their experimental values can be used directly. It is known that the main structure of CFRD is composed of the main and the secondary rockfill. It is necessary to inverse the two types of material parameters separately; without the three parameters listed above, there are parameters that need to be inversed. The calculation work of inversing all 26 parameters is extremely involved; as a result, it is necessary to analyze the sensitivity of these parameters and eliminate those parameters with low sensitivity. Sensitivity Indicators (SIs) of the parameters can be calculated using the modified Morris method [24] as follows:where represents the SI of the rockfill model parameters to the th measurement point’s settlement, is the quantity of calculation times, is the rockfill model parameters value at the th calculation, is the experimental value of the rockfill model parameters, represents the calculated settlement at the th measurement point at the th calculation, and represents the calculated settlement at the th measurement point using the experimental rockfill model parameters.
The measurement points at the elevation of 300 m were analyzed as an example here. The SIs of the main and secondary rockfill parameters are shown in Figure 4.
(a) SI of the - model parameters
(b) SI of the creep model parameters
Several characteristics are shown by Figure 4.(1)The SI of - model parameter is much larger than that of the creep model parameter, and there are significant quantitative differences between the two. As a result, the range of - model parameters value should be smaller, while the range of the creep model parameters value should be relatively larger.(2)The experimental values of the parameters that have relatively lower SIs, such as , , , and , can be used directly.(3)The measurement points buried on the upstream side are more effected by the main rockfill parameters, as shown by the SI of the secondary rockfill parameters becoming smaller as the measurement points get closer to the upstream side. Thus, the measurement points buried on the far upstream side can be used to inverse the main rockfill model parameters, and the points buried on the far downstream side can be used to inverse the secondary rockfill model parameters; the remaining points can be used to verify the rationality of this method.
5. The Objective Function
The difference between the calculated and observed displacements at the measuring points could be taken as an objective function. To investigate the vertical displacements inside the dam body, the settlements of the measurement points in cross section 0+212 are analyzed here. As shown by Figure 2, the measurement points in cross section 0+212 were placed in 5 lines (at elevations of 370 m, 335 m, 300 m, 265 m, and 235 m). The weight of the measurement points at the same elevation needs to be averaged to reflect the effect of height. Because the rockfill creep deformation is time-dependent, several time periods should be analyzed; in this paper, there are 9 time periods, according to Figure 1. In the absence of specific instructions, the weight of the measurement points at different time can be averaged. The mathematical expression of the objective function can be expressed aswhere denotes time, denotes elevation, is the quantity of measuring points at the th elevation, denotes a measurement point, is a group of rockfill model parameters to be inversed, is the calculated value of displacement at th measuring point at time , and is the corresponding observed value of displacement at th measuring point at time .
The optimization problem includes the objective function and the scope of variables; the variable here was , and its scope can be described as follows:where is the quantity of the rockfill model parameters to be inversed and and are the up and down limits of the inversed parameters, respectively.
It should be mentioned that the settlement value of each measurement point at the beginning was not zero, and the meaning of this value was not clear. To overcome this problem during the process of inversion, the paper used the settlement difference between measurement points at different time periods. This is shown in the last part of the paper.
6. MsPSO Displacement Inversion Method
6.1. MsPSO Algorithm
The MsPSO algorithm that was proposed in [17] is a heterogeneous search approach based on four subswarms. Their definition and update rules are presented as follows.
(1) Define the basic subswarms and for exploitation search, in which particles implement the classic velocity and position update rules [25]:where th particle is represented by in -dimensional space and its velocity is . Each particle will dynamically adjust its trajectory based on its historical best position () and the best position () discovered by the whole swarm during the search process. The superscript “” denotes the basic subswarms and and are both positive constants that are known as cognitive and social parameters, respectively. and are vectors of random numbers generated by the uniform distribution between 0 and 1.
(2) Define the adaptive subswarm , in which the particles can adjust the flight directions adaptively by learning from better particles in the basic subswarms. The velocity of th particle, which employs the fitness values from the basic subswarms to update, is updated by the following equation:where and are the fitness values of the particles in the subswarm and , respectively. is the sum of the fitness value: that is, . Similar to and , and are also vectors of random numbers.
The position update rule in is the same as that in (26).
(3) Define the exploration subswarm , in which the velocity and position update rules of th particle are defined aswhere , , and are called the impact factors and .
The cooperation and communication model for the subswarms is shown in Figure 5. In the graph the subswarms and use the information from the basic subswarms and to update. As can be seen, subswarm is the only subswarm that shares all the useful information from the others. The four subswarms cooperate with each other directly or indirectly and all contribute to finding the global optimum .
The steps of inversion by MsPSO are as follows.(1)Initialization: set the maximum number of the iterations and the swarm size . The inertia weight decreased linearly from 0.9 to 0.4, , , , and .(2)Defining subswarms: specify the size of each subswarm .(3)Randomly initialize the position and the velocity of the particles.(4)Update the position of the particles according to (25)~(28) and set , .(5)Termination check: if , move to step (6); otherwise, set , and move to step (4).(6)End the procedure and obtain the inversed rockfill model parameters.
It should be mentioned that the fitness of each particle depends on the objective function. Because the FEM calculated values cannot completely be consistent with the actual measured values, the value of will never reach 0. It is difficult to determine the value of the ; if the = 0, the procedure will iterate for = 500 times. However, as early as approximately 30 iterations, the optimal value no longer changes. Therefore, to save procedure time, the termination check of step (5) was changed to set a certain number of iterations. Once optimal value has been kept for the times, the procedure will stop, and the optimal value is the global optimum. In this case the maximum number and the will no longer need to be set.
The particles in this study are the rockfill parameters, and their numbers and ranges should be determined before inversion. This paper combined the value of experimental parameters shown in Tables 1 and 2 and the sensitivity analysis to perform the work. According to (8), to avoid the denominator becoming negative, the value of parameter cannot be too large. The number and range of the inversed parameters are shown in Table 3.
6.2. Procedures of the Method
Because the FEM analysis needs half an hour to run once and because FEM analysis is performed frequently during the inversion process, considerable computational resources are needed to perform the inversion, and the overall time cost is unacceptable. Therefore, in this study, an RBF neural network was trained to establish the mapping relationship between the rockfill model parameters and the calculated displacement [22].
Figure 6 shows the flowchart of the displacement inversion method proposed in this study.
6.3. Results of the Inversion
The values of the inversed rockfill model parameters and the corresponding objective values are shown in Table 4.
The recalculated displacement values of Shuibuya CFRD using the inversed rockfill model parameters are shown in Figures 7–9.
It should be mentioned that some of the measurement data are not ideal, such as the data at 3-1#. As 3-1# is closer to the upstream side compared with 3-2#, its settlement should be smaller in theory; however, the data show the opposite. This type of measurement data cannot be used for inversion. Among the measurement points in Figures 7–9, the measurement data at 3-2# and 3-3# were used to inverse the main rockfill model parameters, and the measurement data at 3-7# and 3-8# were used to inverse the secondary rockfill model parameters; the other data were used to verify the rationality of the method in this study. By comparing the results, it can be seen that the plots depict a similar trajectory of calculated settlements and measured settlements, indicating that the method of inversion in this study is correct and effective.
As shown in Table 4, the inversed values of - model parameters are approximately equal to their experimental values, while there is a certain gap between the inversed values of creep model parameters and their experimental values.
The creep model proposed by Cheng and Ding in [8] was based on the laboratory creep test, which lasted for approximately 1000 hours. According to (10) and the corresponding experimental values of rockfill creep model parameters, the following can be obtained. For the main rockfill, For the secondary rockfill,As a result, when the laboratory creep test treats the volume creep strain after 1000 hours as the limit volume creep, a contradiction occurred.
In laboratory test, rockfill commonly takes a few hours to complete most of the creep deformation. However, the observed creep deformation from actual engineering lasts for years or even decades. The most important factor that causes the difference is restricting the size of rockfill samples in the laboratory test. Rockfill creep deformation is a process of particle rolling, slipping, rearranging, and the sharp angle crushing. All of these processes end very soon in laboratory tests; furthermore, because the size of the rockfill sample is limited, it is difficult to measure the post creep deformation. Therefore, the observed limit creep deformation is incomplete, which is the reason for the contradiction mentioned above.
This paper concludes that the major problem of laboratory rockfill creep test is the measurement of the limit creep deformation, which causes the relevant parameters to have problems during the calculation of the actual dam creep deformation. As shown in (8) and (13), , , , , , and are positively related to the limit creep deformation. Because the measured limit creep deformation by laboratory test is incomplete, the obtained value of these parameters will be smaller than their actual value. Thus, it is reasonable that the values of the inversed parameters are larger than the experimental values.
7. Conclusion
Some conclusions were obtained in this study.(1)The MsPSO algorithm used in this paper has 36 evolution times, and the recalculated values of FEM using the inversed parameters are in good agreement with the actual measured values. Therefore, the algorithm has a good effect in terms of convergence accuracy and speed. The main reason for this result is that MsPSO with four related subswarms can search extensively and effectively, which overcomes the problem of falling into a local optimum in the search process. As long as those engineering problems can be simplified down to optimization problems, this method can be used efficiently and accurately.(2)There is a certain gap between the inversed parameters and the experimental parameters. It was found that the measurement of the limit creep deformation is incomplete in laboratory tests, so the corresponding parameters do not apply to actual engineering. Thus, the method of obtaining rockfill creep model parameters through back-analysis according to the actual measured data appears to be more reasonable.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This research was supported by the National Science Foundation of China (Grants no. 51179024 and no. 51379029). This financial support is gratefully acknowledged.