# The Application of MsPSO in the Rockfill Parameter Inversion of CFRD

Authors: Ye Zhu; Shichun Chi
Journal: Mathematical Problems in Engineering
 (2016)

---

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


---

## Body

## 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 Figure1, where the solid lines represent the material partition and the dotted lines represent the construction partition.Figure 1
Rockfill material zone and construction progress of the Shuibuya CFRD.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 Figure2 [14].Figure 2
Layout of in situ instrument measurement scheme in the maximum cross section of the Shuibuya CFRD.


## 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 DuncanE-B model [18], and the time-dependent creep deformation was modeled using the power function creep curve models [8].


### 3.1. Outline ofE-B Model

TheE-B 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 Et and bulk modulus B.The tangent modulusEt can be expressed as follows:(1)Et=Ei1-Rf·SL2,where Rf is failure ratio, Ei is initial deformation modulus, which expresses the compression-hardening property of rockfill materials, and SL is stress level. These parameters can be expressed as follows:(2)Ei=kPaσ3Pan,S=σ1-σ31-sin⁡φ2c1·cos⁡φ+2σ3sin⁡φ,where Pa is the standard atmospheric pressure, k is modulus number, n is exponent for defining the influence of the confining pressure on the initial modulus, c1 is cohesive strength, φ is friction angle, and σ1, σ3 are major and minor principal stresses, respectively.For rockfill materials, the cohesive strengthc1 is generally taken as zero. However, the angle of internal friction φ changes with compression stress:(3)φ=φ0-Δφlg⁡σ3Pa,where φ0 and Δφ are two constants for rockfill.The bulk modulusB is given by(4)B=kbPaσ3Pam1,where kb is bulk modulus number and m1 is bulk modulus exponent.TheE-B model parameters of the main and secondary rockfill of Shuibuya CFRD were introduced by Zhou et al. [14] and are shown in Table 1.Table 1
TheE-B model parameters of the Shuibuya CFRD materials.
Dam fill
ρ
d (g/cm3)
φ
0 (°)
Δ
φ (°)
k
n
R
f
k
b
m
1
Main rockfill
2.30
52
8.5
1100
0.35
0.81
600
0.1
Secondary rockfill
2.10
50
8.4
850
0.25
0.84
400
0.05


### 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:(5)ε1t=ε1f1-t-λ1,where ε1f is the limit of axial creep strain and λ1 is a coefficient related to the axial creep strain rate. The time t is calculated in hours. The limit of axial creep strain ε1f that combines both the confining pressure and the stress level viscous responses is given by the following relationship:(6)ε1f=β·σ3,(7)β=c2·SL1-dSL,(8)ε1f=c2·SL1-dSL·σ3,where σ3 is confining pressure, ε1f is proportional to σ3, the coefficient β is the hyperbolic function of SL, and c2 and d are the creep parameters.According to the creep test fitting curve betweenλ1 under various SL and σ3, λ1 is independent of SL, and the relation between λ1 and σ3 can be expressed by the power function as follows:(9)λ1=η·σ3-m2.Thus, the axial creep model of rockfill can be described by (5)–(9) and the creep parameters c2, d, η, and m2.The relation of the volume creep strainεvt can also be expressed by the power function of time as follows:(10)εvt=εvf1-t-λv,where εvt is the limit of volume creep strain and λv is a coefficient related to the volume creep strain rate. The limit of axial creep strain εvt is defined as(11)εvf=αv+βv·σ3,where αv and βv are both functions of the stress level, defined by the following relationship:(12)αv=cαSLdα,βv=cβSLdβ.The limit of volume creep strain can be written as(13)εvf=cα·SLdα+cβ·SLdβ·σ3,where λv is assumed as a constant: λv=const.Thus, the volume creep model of rockfill can be described by (10)–(13) and the creep parameters cα, dα, cβ, dβ, and λv.The creep test parameters of the main and the secondary rockfill of Shuibuya CFRD are shown in Table2.Table 2
Parameters of the power function creep curve model of the Shuibuya CFRD materials [14].
Dam fill
c
2
d
η
m
2
c
α
d
α
c
β
d
β
λ
v
Main rockfill
0.2892
0.8465
0.0831
0.3899
0.4445
2.0827
0.4360
1.6383
0.0678
Secondary rockfill
1.4296
0.2852
0.1116
1.0384
0.3808
0.6521
0.1936
1.0061
0.0851In three-dimensional calculations, the confining pressureσ3 is replaced by p, which equals σ1+σ2+σ3/3.According to (5) and (10), the creep strain ratio ε˙t has the following forms:(14)ε˙1t=λ1·ε1f·t-λ1+1=λ1·ε1f1-ε1tε1f1+1/λ1,ε˙vt=λv·εvf·t-λv+1=λv·εvf1-εvtεvf1+1/λv.The time history is divided into several incrementsΔt, and ε1t and εvt are obtained through summation:(15)ε1t=∑ε˙1Δt,εvt=∑ε˙vΔt.Therefore, the increments of shear creep and volume creep are given by the following, respectively:(16)Δε1c=ε˙1Δt,Δεvc=ε˙vΔt.The relation among shear creep strainεs, axial creep strain ε1, and volume creep strain εv under ordinary conditions of triaxial test is given as follows: (17)εs=233ε1-εv.The above equation is also suitable for ε˙st, εst, and εsf.If the Prandtl-Reuss flow law is adopted on theπ plane, the strain velocity tensor can be expressed as(18)ε˙=13ε˙vI+12ε˙sqS,where ε˙ is a strain velocity tensor, I is a unit tensor, S is the deviatoric stress tensor, and q is deviator stress, which is given by(19)q=12σ1-σ22+σ2-σ32+σ3-σ121/2.Then the creep increment tensorΔεc can be written by(20)Δεc=13ΔεvcI+12ΔεscqS.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:(21)ΔRc=∭BTDΔεcdxdydz.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 Figure1. 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.Figure 3
The load steps in the FEM analysis.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].


## 3.1. Outline ofE-B Model

TheE-B 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 Et and bulk modulus B.The tangent modulusEt can be expressed as follows:(1)Et=Ei1-Rf·SL2,where Rf is failure ratio, Ei is initial deformation modulus, which expresses the compression-hardening property of rockfill materials, and SL is stress level. These parameters can be expressed as follows:(2)Ei=kPaσ3Pan,S=σ1-σ31-sin⁡φ2c1·cos⁡φ+2σ3sin⁡φ,where Pa is the standard atmospheric pressure, k is modulus number, n is exponent for defining the influence of the confining pressure on the initial modulus, c1 is cohesive strength, φ is friction angle, and σ1, σ3 are major and minor principal stresses, respectively.For rockfill materials, the cohesive strengthc1 is generally taken as zero. However, the angle of internal friction φ changes with compression stress:(3)φ=φ0-Δφlg⁡σ3Pa,where φ0 and Δφ are two constants for rockfill.The bulk modulusB is given by(4)B=kbPaσ3Pam1,where kb is bulk modulus number and m1 is bulk modulus exponent.TheE-B model parameters of the main and secondary rockfill of Shuibuya CFRD were introduced by Zhou et al. [14] and are shown in Table 1.Table 1
TheE-B model parameters of the Shuibuya CFRD materials.
Dam fill
ρ
d (g/cm3)
φ
0 (°)
Δ
φ (°)
k
n
R
f
k
b
m
1
Main rockfill
2.30
52
8.5
1100
0.35
0.81
600
0.1
Secondary rockfill
2.10
50
8.4
850
0.25
0.84
400
0.05


## 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:(5)ε1t=ε1f1-t-λ1,where ε1f is the limit of axial creep strain and λ1 is a coefficient related to the axial creep strain rate. The time t is calculated in hours. The limit of axial creep strain ε1f that combines both the confining pressure and the stress level viscous responses is given by the following relationship:(6)ε1f=β·σ3,(7)β=c2·SL1-dSL,(8)ε1f=c2·SL1-dSL·σ3,where σ3 is confining pressure, ε1f is proportional to σ3, the coefficient β is the hyperbolic function of SL, and c2 and d are the creep parameters.According to the creep test fitting curve betweenλ1 under various SL and σ3, λ1 is independent of SL, and the relation between λ1 and σ3 can be expressed by the power function as follows:(9)λ1=η·σ3-m2.Thus, the axial creep model of rockfill can be described by (5)–(9) and the creep parameters c2, d, η, and m2.The relation of the volume creep strainεvt can also be expressed by the power function of time as follows:(10)εvt=εvf1-t-λv,where εvt is the limit of volume creep strain and λv is a coefficient related to the volume creep strain rate. The limit of axial creep strain εvt is defined as(11)εvf=αv+βv·σ3,where αv and βv are both functions of the stress level, defined by the following relationship:(12)αv=cαSLdα,βv=cβSLdβ.The limit of volume creep strain can be written as(13)εvf=cα·SLdα+cβ·SLdβ·σ3,where λv is assumed as a constant: λv=const.Thus, the volume creep model of rockfill can be described by (10)–(13) and the creep parameters cα, dα, cβ, dβ, and λv.The creep test parameters of the main and the secondary rockfill of Shuibuya CFRD are shown in Table2.Table 2
Parameters of the power function creep curve model of the Shuibuya CFRD materials [14].
Dam fill
c
2
d
η
m
2
c
α
d
α
c
β
d
β
λ
v
Main rockfill
0.2892
0.8465
0.0831
0.3899
0.4445
2.0827
0.4360
1.6383
0.0678
Secondary rockfill
1.4296
0.2852
0.1116
1.0384
0.3808
0.6521
0.1936
1.0061
0.0851In three-dimensional calculations, the confining pressureσ3 is replaced by p, which equals σ1+σ2+σ3/3.According to (5) and (10), the creep strain ratio ε˙t has the following forms:(14)ε˙1t=λ1·ε1f·t-λ1+1=λ1·ε1f1-ε1tε1f1+1/λ1,ε˙vt=λv·εvf·t-λv+1=λv·εvf1-εvtεvf1+1/λv.The time history is divided into several incrementsΔt, and ε1t and εvt are obtained through summation:(15)ε1t=∑ε˙1Δt,εvt=∑ε˙vΔt.Therefore, the increments of shear creep and volume creep are given by the following, respectively:(16)Δε1c=ε˙1Δt,Δεvc=ε˙vΔt.The relation among shear creep strainεs, axial creep strain ε1, and volume creep strain εv under ordinary conditions of triaxial test is given as follows: (17)εs=233ε1-εv.The above equation is also suitable for ε˙st, εst, and εsf.If the Prandtl-Reuss flow law is adopted on theπ plane, the strain velocity tensor can be expressed as(18)ε˙=13ε˙vI+12ε˙sqS,where ε˙ is a strain velocity tensor, I is a unit tensor, S is the deviatoric stress tensor, and q is deviator stress, which is given by(19)q=12σ1-σ22+σ2-σ32+σ3-σ121/2.Then the creep increment tensorΔεc can be written by(20)Δεc=13ΔεvcI+12ΔεscqS.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:(21)ΔRc=∭BTDΔεcdxdydz.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 Figure1. 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.Figure 3
The load steps in the FEM analysis.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 Rf 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 7+9-3×2=26 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:(22)Sky=1n-1∑i=1n-1yki+1-yki/yk0di+1-di/d0,where Sky represents the SI of the rockfill model parameters to the kth measurement point’s settlement, n is the quantity of calculation times, d(i) is the rockfill model parameters value at the ith calculation, d(0) is the experimental value of the rockfill model parameters, yki represents the calculated settlement at the kth measurement point at the ith calculation, and yk0 represents the calculated settlement at the kth 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 Figure4.Figure 4
SI of the parameters to the measurement points settlement at the elevation of 300 m.
(a)
SI of theE-B model parameters
(b)
SI of the creep model parametersSeveral characteristics are shown by Figure4.(1)
The SI ofE-B 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 E-B 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 asm1, η, m2, and λv, 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 Figure2, 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 as(23)FX=∑i=19∑j=151nij∑k=1nijukXuk∗-12,where i denotes time, j denotes elevation, nij is the quantity of measuring points at the jth elevation, k denotes a measurement point, X is a group of rockfill model parameters to be inversed, ukX is the calculated value of displacement at kth measuring point at time i, and uk∗ is the corresponding observed value of displacement at kth measuring point at time i.The optimization problem includes the objective function and the scope of variables; the variable here wasX, and its scope can be described as follows:(24)X=x1,x2,…,xd,…,xD,xdmin≤xd≤xdmaxd=1,2,…,D,where D is the quantity of the rockfill model parameters to be inversed and xdmax and xdmin 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 subswarmsS1 and S2 for exploitation search, in which particles implement the classic velocity and position update rules [25]:(25)Vi1/2t+1=ωVi1/2t+c1r1Pi1/2-Xi1/2t+c2r2gbest-Xi1/2t,(26)Xit+1=Xit+Vit+1,where ith particle is represented by Xi=Xi1,Xi2,…,XiD in D-dimensional space and its velocity is Vi=Vi1,Vi2,…,ViD. Each particle will dynamically adjust its trajectory based on its historical best position (Pi) and the best position (gbest) discovered by the whole swarm during the search process. The superscript “(1)/(2)” denotes the basic subswarms and c1 and c2 are both positive constants that are known as cognitive and social parameters, respectively. r1 and r2 are vectors of random numbers generated by the uniform distribution between 0 and 1.(2) Define the adaptive subswarmS3, in which the particles can adjust the flight directions adaptively by learning from better particles in the basic subswarms. The velocity of ith particle, which employs the fitness values from the basic subswarms to update, is updated by the following equation:(27)Vi3t+1=ωγγ1Vi1t+1+γγ2Vi2t+1+Vi3t+c1r3Pi3-Xi3t+c2r4gbest-Xi3t,where γ1 and γ2 are the fitness values of the particles in the subswarm S1 and S2, respectively. γ is the sum of the fitness value: that is, γ=γ1+γ2. Similar to r1 and r2, r3 and r4 are also vectors of random numbers.The position update rule inS3 is the same as that in (26).(3) Define the exploration subswarmS4, in which the velocity and position update rules of ith particle are defined as(28)Vi4t+1=Vi1t+1+Vi2t+1-Vi3t+1,Xi4t+1=α1Xi4t+α2Pi4+α3gbest+Vi4t+1,where α1, α2, and α3 are called the impact factors and α1+α2+α3=1.The cooperation and communication model for the subswarms is shown in Figure5. In the graph the subswarms S3 and S4 use the information from the basic subswarms S1 and S2 to update. As can be seen, subswarm S4 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 gbest.Figure 5
Information sharing model of MsPSO.The steps of inversion by MsPSO are as follows.(1)
Initialization: set the maximum number of the iterationstmax=500 and the swarm size n=40. The inertia weight ω decreased linearly from 0.9 to 0.4, c1=c2=1.49, α1=1/6, α2=1/3, and α3=1/2.(2)
Defining subswarms: specify the sizen/4 of each subswarm SSS=1,2,3,4.(3)
Randomly initialize the position and the velocity of the particles.(4)
Update the position of the particles according to (25)~(28) and set gbestS=arg⁡minfXiS, gbest=arg⁡minfgbestS.(5)
Termination check: ifgbest≤goal, move to step (6); otherwise, set t=t+1, 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 ofgbest will never reach 0. It is difficult to determine the value of the goal; if the goal = 0, the procedure will iterate for tmax = 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 goal 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 Tables1 and 2 and the sensitivity analysis to perform the work. According to (8), to avoid the denominator becoming negative, the value of parameter d cannot be too large. The number and range of the inversed parameters are shown in Table 3.Table 3
Ranges of the parameters to be inversed.
Material
K
n
K
b
c
2
d
c
α
d
α
c
β
d
β
Main rockfill
1000~1200
0.3~0.4
500~700
0.2~2.6
0.8~1.6
0.4~2.8
2~4.8
0.4~2
1.5~3.1
Secondary rockfill
750~950
0.2~0.3
300~500
1.5~3.5
0.2~1
0.3~2.7
0.5~3.3
0.2~2.2
1~2.6


### 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].Figure6 shows the flowchart of the displacement inversion method proposed in this study.Figure 6
Displacement inversion method based on an RBF neural network and MsPSO.


### 6.3. Results of the Inversion

The values of the inversed rockfill model parameters and the corresponding objective values are shown in Table4.Table 4
The values of the inversed parameters.
Material
K
n
K
b
c
2
d
c
α
d
α
c
β
d
β
Main rockfill
F
(
X
)
=
2.05
Experimental
1100
0.35
600
0.2892
0.8465
0.4445
2.0827
0.4360
1.6383
Inversion
1080
0.33
590
0.574
1.223
0.856
3.329
0.555
1.913
Secondary rockfill
F
(
X
)
=
1.73
Experimental
850
0.25
400
1.4296
0.2852
0.3808
0.6521
0.1936
1.0061
Inversion
848
0.22
379
1.85
0.580
0.410
1.510
0.314
2.056The recalculated displacement values of Shuibuya CFRD using the inversed rockfill model parameters are shown in Figures7–9.Figure 7
Cumulative settlements over time of the main rockfill area at an elevation of 300 m.Figure 8
Cumulative settlement over time of the secondary rockfill area at an elevation of 300 m.Figure 9
Cumulative settlement over time of the center point at elevations of 235 m, 265 m, and 335 m.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 Figures7–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 Table4, the inversed values of E-B 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,(29)εv1000=εvf1-1000-0.0678=0.374εvf.For the secondary rockfill,(30)εv1000=εvf1-1000-0.0851=0.444εvf.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), c2, d, cα, dα, cβ, and dβ 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.


## 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 subswarmsS1 and S2 for exploitation search, in which particles implement the classic velocity and position update rules [25]:(25)Vi1/2t+1=ωVi1/2t+c1r1Pi1/2-Xi1/2t+c2r2gbest-Xi1/2t,(26)Xit+1=Xit+Vit+1,where ith particle is represented by Xi=Xi1,Xi2,…,XiD in D-dimensional space and its velocity is Vi=Vi1,Vi2,…,ViD. Each particle will dynamically adjust its trajectory based on its historical best position (Pi) and the best position (gbest) discovered by the whole swarm during the search process. The superscript “(1)/(2)” denotes the basic subswarms and c1 and c2 are both positive constants that are known as cognitive and social parameters, respectively. r1 and r2 are vectors of random numbers generated by the uniform distribution between 0 and 1.(2) Define the adaptive subswarmS3, in which the particles can adjust the flight directions adaptively by learning from better particles in the basic subswarms. The velocity of ith particle, which employs the fitness values from the basic subswarms to update, is updated by the following equation:(27)Vi3t+1=ωγγ1Vi1t+1+γγ2Vi2t+1+Vi3t+c1r3Pi3-Xi3t+c2r4gbest-Xi3t,where γ1 and γ2 are the fitness values of the particles in the subswarm S1 and S2, respectively. γ is the sum of the fitness value: that is, γ=γ1+γ2. Similar to r1 and r2, r3 and r4 are also vectors of random numbers.The position update rule inS3 is the same as that in (26).(3) Define the exploration subswarmS4, in which the velocity and position update rules of ith particle are defined as(28)Vi4t+1=Vi1t+1+Vi2t+1-Vi3t+1,Xi4t+1=α1Xi4t+α2Pi4+α3gbest+Vi4t+1,where α1, α2, and α3 are called the impact factors and α1+α2+α3=1.The cooperation and communication model for the subswarms is shown in Figure5. In the graph the subswarms S3 and S4 use the information from the basic subswarms S1 and S2 to update. As can be seen, subswarm S4 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 gbest.Figure 5
Information sharing model of MsPSO.The steps of inversion by MsPSO are as follows.(1)
Initialization: set the maximum number of the iterationstmax=500 and the swarm size n=40. The inertia weight ω decreased linearly from 0.9 to 0.4, c1=c2=1.49, α1=1/6, α2=1/3, and α3=1/2.(2)
Defining subswarms: specify the sizen/4 of each subswarm SSS=1,2,3,4.(3)
Randomly initialize the position and the velocity of the particles.(4)
Update the position of the particles according to (25)~(28) and set gbestS=arg⁡minfXiS, gbest=arg⁡minfgbestS.(5)
Termination check: ifgbest≤goal, move to step (6); otherwise, set t=t+1, 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 ofgbest will never reach 0. It is difficult to determine the value of the goal; if the goal = 0, the procedure will iterate for tmax = 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 goal 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 Tables1 and 2 and the sensitivity analysis to perform the work. According to (8), to avoid the denominator becoming negative, the value of parameter d cannot be too large. The number and range of the inversed parameters are shown in Table 3.Table 3
Ranges of the parameters to be inversed.
Material
K
n
K
b
c
2
d
c
α
d
α
c
β
d
β
Main rockfill
1000~1200
0.3~0.4
500~700
0.2~2.6
0.8~1.6
0.4~2.8
2~4.8
0.4~2
1.5~3.1
Secondary rockfill
750~950
0.2~0.3
300~500
1.5~3.5
0.2~1
0.3~2.7
0.5~3.3
0.2~2.2
1~2.6


## 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].Figure6 shows the flowchart of the displacement inversion method proposed in this study.Figure 6
Displacement inversion method based on an RBF neural network and MsPSO.


## 6.3. Results of the Inversion

The values of the inversed rockfill model parameters and the corresponding objective values are shown in Table4.Table 4
The values of the inversed parameters.
Material
K
n
K
b
c
2
d
c
α
d
α
c
β
d
β
Main rockfill
F
(
X
)
=
2.05
Experimental
1100
0.35
600
0.2892
0.8465
0.4445
2.0827
0.4360
1.6383
Inversion
1080
0.33
590
0.574
1.223
0.856
3.329
0.555
1.913
Secondary rockfill
F
(
X
)
=
1.73
Experimental
850
0.25
400
1.4296
0.2852
0.3808
0.6521
0.1936
1.0061
Inversion
848
0.22
379
1.85
0.580
0.410
1.510
0.314
2.056The recalculated displacement values of Shuibuya CFRD using the inversed rockfill model parameters are shown in Figures7–9.Figure 7
Cumulative settlements over time of the main rockfill area at an elevation of 300 m.Figure 8
Cumulative settlement over time of the secondary rockfill area at an elevation of 300 m.Figure 9
Cumulative settlement over time of the center point at elevations of 235 m, 265 m, and 335 m.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 Figures7–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 Table4, the inversed values of E-B 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,(29)εv1000=εvf1-1000-0.0678=0.374εvf.For the secondary rockfill,(30)εv1000=εvf1-1000-0.0851=0.444εvf.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), c2, d, cα, dα, cβ, and dβ 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.



---

*Source: 1096967-2016-05-26.xml*