Dynamic responses of concrete-faced rockfill dam due to different seismic motion input methods

Abstract

The concrete-faced rockfill dam valley foundation was considered as an open energy system and a reasonable non-uniform seismic motion input method was applied to the dynamic analysis of a concrete-faced rockfill dam based on the generalized plastic model. First, the corresponding program was validated by means of the scattering question of an idealized semicircle valley. Subsequently, the seismic elasto-plastic finite element analyses were performed to compare and investigate the performances of a concrete-faced rockfill dam under different seismic motion input methods. The results show that the dynamic responses of the concrete-faced rockfill dam are decreased by 10%–30% approximately with the use of non-uniform seismic motion input method. As a result, the traditional uniform seismic motion input method would overestimate the responses of the dam. From the perspective of seismic safety evaluation, the overestimations would disturb the reasonable assessment of the aseismic capacity of the dam. Moreover, the slope stability analysis results might be conservative and unreasonable due to overestimating the accelerations during the earthquake.

Introduction

Concrete-faced rockfill dam (CFRD) is a type of embankment dam that is being widely used and rapidly developed. Several CFRDs over 200 m in height are being built or to be built in China’s strong seismic regions, such as Houziyan, Gushui, and Lawa. Since the dam failures due to strong earthquakes will cause catastrophic consequences, the safety of these CFRDs during earthquake is an important engineering and social problem and has been well recognized.^1,2 Numerical analysis of CFRDs during strong motion is an effective method to evaluate the seismic safety and has been investigated by many researchers.^3–6 Although some sophisticated numerical methods have been developed for seismic analysis of high dams, it is still immature in studying realistic earthquake behavior of a complete CFRD valley foundation system especially involving nonlinearities of rockfill and foundation materials. Generally, the fixed boundary and uniform seismic motion input (USMI) method proposed by Clough ⁷ were adopted in most of the studies, which means CFRDs were fixed on the shaking table (the same ground motion acts simultaneously). The foundation flexibility can be considered in this method but the interaction between CFRD and foundation was not appropriate. In addition, with the increase in dam height, the dam axial length even over 1000 m, the wave propagation effect is obvious in the dam body. It is well known that the CFRD valley foundation is an open energy system, so the different interaction between them and the wave radiation effects due to infinite foundation during earthquake should be considered. Therefore, the reasonable safety evaluation of CFRDs has to be researched through non-uniform seismic motion input (NUSMI) method. ⁸

To cover the limit of the USMI method, various artificial boundaries and corresponding seismic motion input method^9–16 have been investigated. The essence of the artificial boundary is to allow scattered wave from the generalized structure through the truncated foundation boundaries toward infinite domain without reflecting. The typical artificial boundaries include transmitting boundary, ¹¹ absorbing boundary, ¹² viscous boundary, ¹³ and viscous-spring boundary.^14,15 Viscous-spring artificial boundary has been applied in dynamic response of underground structure, ¹⁶ arch dam,^17–19 earth and rockfill dam, ²⁰ and so on. Zhao et al. ¹⁶ found that the result of the method has high precision compared with analytical result, and this method was used in dynamic analysis of a cross-river shield tunnel. Zhang et al. ¹⁹ developed the artificial boundary and investigated the arch dam–canyon interaction, and the results showed that the use of viscous-spring boundary model leads to about 25%–40% reduction in dynamic response of the arch dam. Zhou et al. ²⁰ investigated the effects of the method on the interaction between earth-core rockfill dam and foundation based on equivalent linear model (Hardin–Drnevich model). However, this method has rarely been used in CFRDs especially based on the generalized plastic model. Reasonable dynamic analysis results are particularly important for assessing the dam safety, including the slope stability, the foundation liquefaction, the slab damage, and the aseismic measures design. So, in this study, on the basis of the generalized plastic model ²¹ for the rockfill materials, the generalized plastic interface model ²² and element with asymmetric nodes ²³ between slab and cushion layer contact surface established by the authors, the plastic damage constitutive model for concrete proposed by Lee and Fenves, ²⁴ the viscous-spring boundary, and the corresponding NUSMI method ¹⁵ are used by the authors. Numerical verification of scattering question of an idealized semicircle valley is performed. Subsequently, for comparing effects of different seismic motion input methods on the CFRDs’ dynamic responses, including the dam acceleration, concrete face slab damage development, and distribution and the permanent deformation, a series of seismic analyses of the CFRDs constructed on the bedrock are performed through two-dimensional (2D) finite element numerical analysis. Finally, some conclusions are given.

NUSMI method

Viscous-spring artificial boundary

Based on the assumption that the scattering wave could be regarded as 2D cylindrical wave and three-dimensional (3D) spherical wave, the 2D and 3D viscous-spring artificial boundary was proposed by Liu and Lu. ¹⁵ As illustrated in Figure 1, the viscous-spring artificial boundary can be approximately equivalent as parallel spring–dashpot system installing on the cut-off model boundary. This boundary is efficient and convenient for finite element procedure with satisfied accuracy and is used herein. The spring coefficient K and viscous dashpots coefficient C can be calculated as follows ¹⁴

$$K=\alpha \frac{G}{r}\sum _{i=1}^I{A}_i$$

(1)

$$C=\rho c\sum _{i=1}^I{A}_i$$

(2)

where ρ and G are the medium mass density and shear modulus, respectively. r is the distance from the scattering source to the artificial boundary. For normal direction and tangential direction boundary, c, respectively, denotes the velocity of P wave and S wave. $\alpha$ is the coefficient of boundary parameter in different direction. $\alpha$ is equal to 2.0 in normal direction and 1.5 in tangential direction according to the reference literature by J Liu and Li. ²⁵ $\sum _{i=1}^I{A}_i$ is the area controlled by the boundary node.

Figure 1.

Viscous-spring artificial boundary.

On the basis of viscous-spring artificial boundary, the authors developed a type of element (viscous-spring boundary element (VSBE)) to substitute the traditional concentrated viscous-spring artificial boundary. The VSBE is an interface element with spring and damping as illustrated in Figure 2. The spring and dashpot’s parameters are calculated by equations (1) and (2). Compared with the traditional concentrated viscous-spring artificial boundary, in which two or three boundary elements should be added to perform 2D or 3D analysis, only one boundary element is needed using VSBE.

Figure 2.

Viscous-spring interface element.

The equivalent nodal load

Different from the USMI method, amplitudes and phases of earthquake are different along boundary. ¹⁵ Viscoelastic artificial boundary can be easily combined with finite element mesh in NUSMI method. The VSBE is used by spring and damping belongs to boundary of the model. By the equivalent nodal load on the boundary nodes, it can realize seismic motion input with VSBE. The equivalent nodal load was used to simulate the stress boundary of the actual wave field and it can be calculated as follows ¹⁴

$${\mathit{F}}_b={\mathit{R}}_b^{\mathrm{e}f}+{\mathit{C}}_b{\stackrel{·}{\mathit{u}}}_b^{\mathrm{e}f}+{\mathit{K}}_b{\mathit{u}}_b^{\mathrm{e}f}$$

(3)

where ${\mathit{K}}_b$ and ${\mathit{C}}_b$ are the additional matrix of element stiffness and damping caused by viscous-spring artificial boundary, respectively. ${\mathit{F}}_b$ is the equivalent nodal load vector on the boundary nodes and realizes seismic motion input with VSBE. ${\mathit{u}}_b^{\mathrm{ef}}$ , ${\stackrel{·}{\mathit{u}}}_b^{\mathrm{e}f}$ , and ${\mathit{R}}_b^{\mathrm{ef}}$ are the displacement vector, velocity vector, and stress vector at boundary in the free field, respectively.

Scattering question of semicircular valley

The scattering question of semicircular valley under vertical incidence of SV wave is analyzed to validate the accuracy of NUSMI method ¹⁵ and the program was developed by the authors.

As illustrated in Figure 3, the incidence angle θ is equal to 0° under wave vertical incidence. The valley radius r₀ is equal to 210 m. The mass density of material ρ = 2.65 g/cm³, and the shear modulus G = 5.292 GPa, and Poisson’s ratio υ = 0.33. The height of model y_b = 525 m and the width x_b = 1050 m. The finite element grid is shown in Figure 4, in which height and width of the outer element both are equal to 17.5 m. The expression of non-dimensional frequency is defined as η = 2r₀/λ = ωr₀/πc_s (0.5 and 2.0 are selected), in which the wavelength λ = 1 m, ω is the circular frequency, and c_s is the shear wave velocity.

Figure 3.

Semicircle valley under wave incidence.

Figure 4.

Finite element mesh of semicircle valley.

Due to no analytical solution for benchmark problem, the calculated results are calibrated with the solution using indirect boundary integral method by Wong. ²⁶ The comparisons are illustrated in Figures 5 and 6; it can be seen that the NUSMI input method could be used to simulate the responses of the canyon with the satisfied accuracy.

Figure 5.

Surface displacement amplitude under SV wave vertical incidence (η = 0.5): (a) horizontal and (b) vertical.

Figure 6.

Surface displacement amplitude under SV wave vertical incidence (η = 2.0): (a) horizontal and (b) vertical.

Interface element with asymmetric nodes

To study the dynamic damage of the concrete face slab, the mesh size of the slab should be controlled less than 0.3 m.^6,23,27 However, the transition of the finite element mesh from the slab to the cushion layer and rockfill materials would be complex and some irregular mesh will influence the calculation accuracy. Therefore, interface element with asymmetric nodes was applied by the authors to resolve this problem. ²³ The interface element contains arbitrary number of nodes and can solve the problem of mesh transition between the slab and the cushion layer. Meanwhile, it has the advantages in modeling and the adaptation of the change of dam zoning. In addition, on the basis of ensuring the accuracy of the calculation, the element can greatly reduce the number of elements and improve the computational efficiency.

The typical interface element with asymmetric nodes is shown in Figure 7. The fine mesh connected to the nodes at the top surface and the coarse mesh connected to the nodes at the bottom surface. In simulation, the interface element is divided into n − 1 sub-elements using the n − 1 virtual nodes placed on the bottom surface that will be finally condensed out, and each sub-element can be regarded as a traditional slip element. The relationship between the displacement vectors $\{\delta \}$ and $\{\delta \}\prime$ with and without virtual nodes is

$$\{\delta \}\prime =[T]\{\delta \}$$

(4)

$$[T]=\left[\begin{array}{c}I\hfill \\ \multicolumn{1}{c}{t}\end{array}\right]$$

(5)

$$\mathit{t}=\left[\begin{array}{ccccccc}{N}_n^n\hfill & 0\hfill & N_1^n\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_n^n\hfill & 0\hfill & N_1^n\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{N_n^i}& 0\hfill & N_1^i\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_n^i\hfill & 0\hfill & N_1^i\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{N_n^2}& 0\hfill & N_1^2\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_n^2\hfill & 0\hfill & N_1^2\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]$$

(6)

$$N_1^i=\frac{x_i-x_n}{L},N_n^i=\frac{x_1-x_i}{L}$$

(7)

where [T] is a matrix of (4n) × (2n + 4), L is the length of the interface, and i is the node number.

Figure 7.

The interface element with asymmetric nodes.

The stiffness matrix of the interface element, $[K]\prime$ , that contains the degrees of freedom (DOFs) of the virtual nodes can be obtained from the stiffness matrix of each sub-element according to its number of DOFs, and $[K]\prime$ is a 4n × 4n matrix. The stiffness matrix, $[K]$ , of the interface element after the virtual nodes have been condensed out can be obtained from the following equation

$$[K]=[T{]}^T[K]\prime [T]$$

(8)

Compared with the interface element with asymmetric nodes with 2 Gauss points and the traditional Goodman element, the results show a good agreement and the maximum difference is less than 1%. It indicated that the accuracy of proposed interface element is acceptable. ²³

Finite element analysis of CFRD

The GEODYNA program, ²⁸ which was developed by Institute of Earthquake Engineering, Dalian University of Technology, was adopted for static and dynamic finite element analysis in this study. Using GEODYNA, the dynamic responses were investigated by the authors.^{6,20,21,23,27}

Finite element model

A 2D CFRD with a height of 200 m is used as the finite element model. The upstream slope is 1:1.4, and the downstream slope is 1:1.5. The dam was filled in 114 layers, and the concrete face slab was constructed in three stages (the slab top heights of each stage are 60, 130, and 200 m, respectively). The water level is 190 m, and impoundment begins after the completion of second stage of face slab pouring.

The finite element mesh of the CFRD is shown in Figure 8; the detailed mesh view of the concrete face slab, cushion layer, transition layer, and interface elements is shown in Figure 8. Seven “observation points,” as shown in Figure 9, were chosen to compare the dynamic responses due to different seismic motion input methods. The concrete face slab is divided into five layers along the normal direction; the largest element size is 20 cm, for the convenience of studying the damage development process. The slab element uses 4-node isoparametric elements and the contact surface between the slab and cushion layer uses asymmetric 19-node interface element. The model contains 18,636 elements and 20,472 nodes in total, including 8960 slab elements.

Figure 8.

The finite element mesh of the concrete-faced rockfill dam.

Figure 9.

The detailed mesh view of the concrete face slab, cushion layer, transition layer, and interface elements.

The dynamic water pressure on face slabs was considered by added mass element during earthquake. The added mass element was defined on the nodes outside of the face slabs. The added mass intensity can be calculated by “Specifications for seismic design of hydraulic structures” (DL 5073-2000). ²⁹

Material parameters

The generalized plasticity model parameters for rockfill materials ²¹ are listed in Table 1. The contact surface between the concrete face slab and the cushion layer uses a generalized plasticity interface model, and the computational parameters ²³ are listed in Table 2. The plastic damage model proposed by Lee and Fenves ²⁴ was used to simulate the concrete face slab. The property parameters used in the analysis are provided in Table 3. Because the compressive strength of concrete is considerably higher than the tensile strength, only the tensile damage of the slab concrete is considered in this study. The bedrock material has the following parameters: elastic modulus E = 10 GPa, Poisson’s ratio ν = 0.25, mass density ρ = 2450 kg/m³, and damping value = 0.05.

Table 1.

Parameters of the generalized plastic model for rockfill materials.

G 0	K 0	Mg	Mf	αf	αg	H 0	HU 0	ms
1000	1400	1.8	1.38	0.45	0.4	1800	3000	0.5
mv	ml	mu	rd	γDM	γu	β 0	β 1
0.5	0.2	0.2	180	50	4	35	0.022

Table 2.

Parameters of the generalized plastic interface model for contact surface between the concrete face slab and the cushion layer.

Ds 0 (kPa)	Dn 0 (kPa)	Mc	er	λ	a (kPa0.5)	b	c
1000	1500	0.88	0.4	0.091	224	0.06	3.0
a	rd	km	Mf	k	H 0 (kPa)	fh	t/m
0.65	0.2	0.6	0.65	0.5	8500	2	0.1

Table 3.

Parameters of the plastic damage model for concrete face slab.

ρ (kg/m3)	E (GPa)	ν	ft (MPa)	fc (MPa)	lc (m)	Gt (N m−1)
2450	31	0.167	3.48	27.6	0.2	325

Input of seismic motion

The seismic motion input to the mechanical model is based on “Specifications for seismic design of hydraulic structures” (DL 5073-2000) ²⁸ with a specified artificial seismic wave spectrum. The seismic wave history (the acceleration response of outcrop) is shown in Figure 10. The horizontal peak acceleration of seismic waves along the river direction is 0.3 g, and the vertical peak acceleration takes two-thirds of that along the river direction. Acceleration amplification response spectra of the standard design and the input of seismic motion are shown in Figure 11.

Figure 10.

The acceleration history of outcrop: (a) horizontal direction and (b) vertical direction.

Figure 11.

The acceleration amplification response spectrum.

Analysis of computational results

The step-by-step construction of the dam was simulated first in the analysis. The initial stresses induced in the concrete face slab and dam body during water impounding are also considered in the analysis and are subsequently applied for the seismic analysis. The dynamic responses of the CFRD were analyzed by two different seismic motion input methods. In case 1, the USMI method was used with massless foundation and the NUSMI method was used with mass of foundation in case 2. In case 1, the massless foundation model is assumed that the bedrock is linear elastic and massless. The boundary of the bedrock is fixed and the seismic inertia force is applied to the dam body, so that seismic loads along the boundary are the same. The same finite element mesh was used for the two cases except boundary element. In this way, the calculated results by different seismic motion input methods could be compared reasonably.

The dam body acceleration responses

Figure 12 illustrates the acceleration response history of the points A, B, and C at the points on dam axis of different elevation during 10–15 s. The maximum acceleration responses along the dam axis are shown in Figure 13. It can be seen that the acceleration histories are similar tendency by different input methods under wave vertical incidence, and the peak value by NUSMI method decreased about 10%–25% and postponed, which indicated that the NUSMI method could reflect the radiation damping and the wave propagation effect in the bedrock and dam body.

Figure 12.

The acceleration responses history of observation points in dam body: (a) point A, (b) point B, and (c) point C.

Figure 13.

The maximum acceleration responses along the dam axis.

The dam deformation

In Figure 14, the settlement histories of “observation point” E located at the dam crest calculated by different seismic motion input methods are shown. The maximum settlement of point E calculated by the NUSMI method decreases 18% compared with that by the USMI method. Figures 15 and 16 illustrate the contours of the horizontal and vertical permanent deformation in the dam body at the end of earthquake. It can be seen that the horizontal and vertical permanent deformation to the downstream direction by the NUSMI method is significant lower (about 20%) than that by the USMI method. It can be concluded that due to the radiation damping of infinite foundation, the lower dam deformation in both horizontal and vertical directions will be expected in the dam body.

Figure 14.

The settlement history of point E by different seismic motion input methods.

Figure 15.

The horizontal permanent deformation in the dam (unit: m): (a) USMI and (b) NUSMI.

Figure 16.

The vertical permanent deformation in the dam (unit: m): (a) USMI and (b) NUSMI.

The stress and damage in the slab

Two typical slab elements (damage occurred and without damage) and stress histories along the slope direction during the earthquake are depicted in Figure 17. The maximum stresses along the slope direction in the middle layer slab by different seismic motion input methods are compared in Figure 18. It can be seen that the tensile and compressive stress tendencies by the two different seismic motion input methods are similar under wave vertical incidence, but the tensile damage range and the maximum compressive stress are reduced due to the radiation damping when the NUSMI method is implemented. Also, because of the plastic damage model used for the concrete face slab, the maximum tensile stress did not exceed the tensile strength of the concrete.

Figure 17.

The slope direction stress history of typical concrete face slab elements: (a) without tensile damage and (b) with tensile damage.

Figure 18.

The maximum slope direction stress of concrete face slab elements: (a) tensile stress and (b) compressive stress.

Moreover, from the tensile damage distribution as illustrated in Figures 19 and 20, penetrating damage mainly occurred in the slab at the location of 0.75 H (H is the dam height) when the NUSMI method is used, which is different from that by the USMI method. During Wenchuan earthquake in 2008, the tensile damage occurred in the slab at the elevation of 845 m (0.8 H).

Figure 19.

The tensile damage distribution of the concrete face slab at the end of earthquake: (a) USMI and (b) NUSMI.

Figure 20.

The tensile damage distribution of the concrete face slab from 0.6 to 0.9 H at the end of earthquake: (a) USMI and (b) NUSMI.

From the dynamic responses study, it is concluded that the USMI method overestimates the dam responses because of ignoring the radiation damping; in this study, the reduction reaches 10%–30% in terms of peak accelerations, compressive stresses, and crest displacements.

Conclusion

The main purpose of this study is to evaluate the effects of the seismic motion input methods on the dynamic responses of CFRDs based on the generalized plastic model. The VSBE was proposed and used for the NUSMI. The NUSMI method and the elasto-plastic finite element procedure for both construction and dynamic analysis were developed and verified. A generalized plastic model for rockfill material has been used and the interface element with asymmetrical nodes is applied to describe the behavior of the interface between concrete face slab and cushion layer. The seismic responses of the CFRD calculated by different seismic motion input methods are compared and investigated. Compared with USMI method, the results show that the acceleration peak value of dam body by NUSMI method decreased about 10%–25%; the horizontal and vertical permanent deformation is about 20% lower and the compressive stress of face slab decreased about 25% at three-fourths of the dam height. It means that the use of NUSMI method leads to about 10%–30% reduction in dynamic response of the CFRD. It indicates that the traditional USMI method would overestimate the responses of the dam. From the perspective of seismic safety evaluation, the overestimations would disturb the reasonable assessment of the aseismic capacity of the dam. Moreover, the slope stability analysis results might be conservative and unreasonable due to overestimating the accelerations during the earthquake.

The 3D dynamic response of CFRDs, particularly which of the Zipingpu Dam during the 2008 Wenchuan earthquake, will be subsequently investigated by extending this numerical procedure, and the propagation effect of the seismic waves on the slab damage and dam deformation mechanism will be discussed. Also, the NUSMI method could be used to study the dynamic responses of the dam constructed on the deep overburden layers considering the nonlinear properties of overburden materials.