# Seismic Damage Analysis of Concrete Gravity Dam Based on Wavelet Transform

Authors: Dunben Sun; Qingwen Ren
Journal: Shock and Vibration
 (2016)

---

## Abstract

The key to the dam damage assessment is analyzing the remaining seismic carrying capacity after an earthquake occurs. In this paper, taking Koyna concrete gravity dam as the object of study, the dynamic response and damage distribution of the dam are obtained based on the concrete damage plastic constitutive model. By using time-frequency localization performance of wavelet transform, the distribution characteristics of wavelet energy for gravity dam dynamic response signal are revealed under the action of different amplitude earthquakes. It is concluded by numerical study that the wavelet energy is concentrated in low-frequency range with the improving of seismic amplitude. The ultimate peak seismic acceleration is obtained according to the concentration degree of low-frequency energy. The earthquake damage of the dam under the moderate-intensity earthquake is simulated and its residual seismic bearing capacity is further analyzed. The new global damage index of the dam is proposed and the overall damage degree of the dam can be distinguished using defined formula under given earthquake actions. The seismic bearing capacity of the intact Koyna dam is 591 gal considering the dam-water interaction and its residual seismic bearing capacity after simulating earthquake can be calculated.


---

## Body

## 1. Introduction

Under earthquakes with the intensity greater than moderate level, the concrete dam may crack and its strength as well as stiffness will decrease. The continuum damage mechanics is an effective way to accurately model the degradation in the mechanical properties of concrete dam. Lee and Fenves [1, 2] proposed the plastic damage constitutive model of concrete and simulated failure processes of concrete structures as well as the Koyna dam under earthquake and acquired the consistent results with the actual earthquake damage. A lot of results have been achieved in analyzing the damage of concrete dams under the action of earthquake [3–5], in which the plastic damage models were used. To assess the structural damage state simulated using numerical methods, many kinds of damage indexes have been proposed [6–9]. Usually the damage index uses a set of structural response parameters that can be used to evaluate the damage to structures. These response parameters, which are generally used for damage assessment of buildings, can be portrayed as deformations, stiffness, energy dissipations, or their combination. For damage assessment of concrete dam, there are limited works in the field of global damage indexes [10, 11]. However, these indexes cannot reflect the residual bearing capacity of the structure.Structural damage usually causes a reduction in the stiffness of structure and will directly influence the dynamic response of the structure, in which the frequency components of vibration signal change with the stiffness reduction. Wavelet analysis has an excellent localized performance in time-frequency domain, which may be used to analyze a signal in an arbitrary resolution in time-frequency domain. Wavelet analysis is generally utilized in damage detection of structures [12–15] but is rarely used in seismic response analysis of structure [16]. So far, almost no article has been found to research the dynamic loading capacity of concrete dam by using wavelet analysis theory.In this paper, the damage distribution characteristics of an intact concrete gravity dam, the Koyna dam, are analyzed under the action of the same earthquake with different amplitudes, and the corresponding vibration signals of dam crest are extracted. By means of wavelet analysis, the energy change rule of vibration signals in low-frequency band is obtained with the improvement of earthquake amplitudes, from which the seismic bearing capacity of the intact dam is acquired. Subsequently, the damage distribution of the dam is simulated under moderate-intensity earthquakes and the residual seismic bearing capacity of an earthquake damaged dam is further calculated. The global damage index of the dam is proposed and the overall seismic damage degree of the dam is calculated.


## 2. Wavelet Analysis Theory

### 2.1. The Concept of Wavelet Transforms

Supposingf
(
t
)
∈
L
2
(
R
), L
2
(
R
) is square-integrable real-valued space, the continuous wavelet transform is defined as follows: (1)
W
f
a
,
b
=
a
∫
-
∞
∞
f
t
ψ
¯
a
t
-
b
d
t
,where ψ
(
t
) is wavelet mother function or basic wavelet, a is scale factor, and b is translation factor, and the value of b can be positive or negative.For the structures engineering, discrete time series signals are usually attained. To facilitate the application, it is necessary to use the discrete wavelet transform. For continuous waveletψ
(
a
t
-
b
), scale factor a and translation factor b can be expressed in discrete form by setting a
=
2
j, b
=
k, j
,
k
∈
Z, Z is integer domain. The discrete wavelet base function can be expressed as (2)
ψ
j
,
k
t
=
2
j
/
2
ψ
2
j
t
-
k
.The discrete wavelet transform is expressed as(3)
W
f
j
,
k
=
f
t
,
ψ
j
,
k
=
∫
R
f
t
ψ
¯
j
,
k
d
t
.


### 2.2. The Concept of Relative Wavelet Energy

Wavelet transform (WT) is a time-frequency transform operation, which converts a time-domain signal into the time-scale plane according to (1). It provides us with a tool to inspect the relatively narrow frequency bands over a relatively short time window. The WT algorithm can be expressed by the tower-style decomposition as Figure 1. The analysis phase of the wavelet transform decomposes a signal into frequency bands that are localized in time and scale. The window size (sale) used in wavelet transform is chosen to be shot at high frequencies and long at low frequencies, providing good time resolution at high-frequency and good frequency resolution at low frequencies. Because of this localization property, wavelets are very good in isolating singularities and irregular structures in signals.Figure 1
Orthogonal decomposition relationship of original signalV
0.This decomposition relationship can be expressed as(4)
V
0
=
V
N
⨁
j
=
1
N
W
j
,
V
j
⊕
W
j
=
V
j
-
1
,where V
j represents the low-frequency part of the original signal and W
j represents the high-frequency detail component. It can be seen from the tower-style decomposition that any signal can be decomposed into the high-frequency part and low-frequency part. The low-frequency part can be decomposed to the next level again. For the binary wavelet, the frequency range of high-frequency band is half of the upper grade frequency band. The development of the damage can be detailedly observed by analyzing structural vibration signals of different frequency domain.Decomposing the original signal into the low-frequency coefficientA
N and high-frequency coefficients D
i
(
i
=
1,2
,
…
,
N
), then the high-frequency band energy for each level can be expressed as E
i
=
D
i
2, where D
i denotes second-order norms of one-dimensional array D
i. In order to facilitate writing, low-frequency band energy is expressed in E
N
+
1, such that E
1
,
E
2
,
…
,
E
N
+
1 constitutes the wavelet energy spectrum of signal f
(
t
) in the N scales. The relative wavelet energy of each frequency band can be expressed as follows: (5)
γ
i
=
E
i
∑
i
=
1
N
+
1
E
i
.


## 2.1. The Concept of Wavelet Transforms

Supposingf
(
t
)
∈
L
2
(
R
), L
2
(
R
) is square-integrable real-valued space, the continuous wavelet transform is defined as follows: (1)
W
f
a
,
b
=
a
∫
-
∞
∞
f
t
ψ
¯
a
t
-
b
d
t
,where ψ
(
t
) is wavelet mother function or basic wavelet, a is scale factor, and b is translation factor, and the value of b can be positive or negative.For the structures engineering, discrete time series signals are usually attained. To facilitate the application, it is necessary to use the discrete wavelet transform. For continuous waveletψ
(
a
t
-
b
), scale factor a and translation factor b can be expressed in discrete form by setting a
=
2
j, b
=
k, j
,
k
∈
Z, Z is integer domain. The discrete wavelet base function can be expressed as (2)
ψ
j
,
k
t
=
2
j
/
2
ψ
2
j
t
-
k
.The discrete wavelet transform is expressed as(3)
W
f
j
,
k
=
f
t
,
ψ
j
,
k
=
∫
R
f
t
ψ
¯
j
,
k
d
t
.


## 2.2. The Concept of Relative Wavelet Energy

Wavelet transform (WT) is a time-frequency transform operation, which converts a time-domain signal into the time-scale plane according to (1). It provides us with a tool to inspect the relatively narrow frequency bands over a relatively short time window. The WT algorithm can be expressed by the tower-style decomposition as Figure 1. The analysis phase of the wavelet transform decomposes a signal into frequency bands that are localized in time and scale. The window size (sale) used in wavelet transform is chosen to be shot at high frequencies and long at low frequencies, providing good time resolution at high-frequency and good frequency resolution at low frequencies. Because of this localization property, wavelets are very good in isolating singularities and irregular structures in signals.Figure 1
Orthogonal decomposition relationship of original signalV
0.This decomposition relationship can be expressed as(4)
V
0
=
V
N
⨁
j
=
1
N
W
j
,
V
j
⊕
W
j
=
V
j
-
1
,where V
j represents the low-frequency part of the original signal and W
j represents the high-frequency detail component. It can be seen from the tower-style decomposition that any signal can be decomposed into the high-frequency part and low-frequency part. The low-frequency part can be decomposed to the next level again. For the binary wavelet, the frequency range of high-frequency band is half of the upper grade frequency band. The development of the damage can be detailedly observed by analyzing structural vibration signals of different frequency domain.Decomposing the original signal into the low-frequency coefficientA
N and high-frequency coefficients D
i
(
i
=
1,2
,
…
,
N
), then the high-frequency band energy for each level can be expressed as E
i
=
D
i
2, where D
i denotes second-order norms of one-dimensional array D
i. In order to facilitate writing, low-frequency band energy is expressed in E
N
+
1, such that E
1
,
E
2
,
…
,
E
N
+
1 constitutes the wavelet energy spectrum of signal f
(
t
) in the N scales. The relative wavelet energy of each frequency band can be expressed as follows: (5)
γ
i
=
E
i
∑
i
=
1
N
+
1
E
i
.


## 3. Plastic Damage Constitutive Mode

Concrete structure may crack under an earthquake and the structural strength as well as stiffness will decrease. The concrete damage plastic (CDP) model provided by the software ABAQUS can simulate concrete damage and failure behavior. CDP model assumes that the concrete material would be damaged due to tensile cracking and compression crushing. The evolution of the corresponding yield or failure surface is controlled by the tensile equivalent plastic strainε
~
t
p
l and the compression equivalent plastic strain ε
~
c
p
l. The stress-strain relationship is expressed by the following formula:(6)
σ
=
1
-
d
σ
-
=
1
-
d
D
0
e
l
:
ε
-
ε
p
l
,where D
0
e
l is the initial elastic stiffness matrix of the material; σ is the total stress tensor; σ
- is the effective stress tensor; d is the scalar damage variable, 0
≤
d
≤
1
,
d
=
0 for undamaged state, d
=
1 for completely damaged state; ε
p
l is plastic strain. Plastic-strain rate can be expressed as follows by adopting the nonassociated flow rule:(7)
ε
˙
p
l
=
λ
∂
G
σ
-
∂
σ
-
,where λ is the plastic consistent parameters, G is the plastic flow potential function, and the Drucker-Prager hyperbolic function has been chosen here which can be expressed as (8)
G
=
ζ
σ
t
0
tan
⁡
φ
2
+
q
-
2
-
p
-
tan
⁡
φ
,where φ is the expansion angle measured in p
-
q  plane at high confining pressure; σ
t
0 is the tensile stress when the uniaxial tension failure occurs; ζ is eccentricity ratio.The yield function expressed in effective stress modified by Lee and Fenves has been taken as the yield condition of the concrete, which can be written in effective stress:(9)
F
=
1
1
-
α
q
-
-
3
α
p
-
+
β
ε
-
p
l
σ
-
m
a
x
-
γ
-
σ
-
m
a
x
-
σ
-
c
ε
c
p
l
,where α, γ are the dimensionless constant; p
- is the effective hydrostatic pressure; q
- is Mises equivalent effective stress; σ
- is effective stress tensor; σ
-
m
a
x is the maximum eigenvalue of effective stress tensor σ
-; symbol  denotes the ramp function as x
=
(
1
/
2
)
(
x
+
x
); β
(
ε
-
p
l
) is a dimensionless variable which can be written as a function of the two effective plastic-strain parameters, ε
-
t
p
l and ε
-
c
p
l:(10)
β
ε
-
p
l
=
σ
-
c
ε
-
c
p
l
σ
-
t
ε
-
t
p
l
1
-
α
-
1
+
α
.The elastic stiffness is restored as the cracks change from opening to closing and the load changes from tension state to compression state under cyclic loading. Thus the damage variable can be expressed as follows:(11)
d
=
1
-
1
-
r
d
c
1
-
t
d
t
,where r and t are the stiffness recovery factors when the internal stress appears opposite in signs under cyclic loading.


## 4. The Seismic Damage Analysis of Gravity Dam

Concrete structure may appear to crack under an earthquake and structural strength as well as stiffness is also reduced. The concrete damage plastic model provided by the software ABAQUS can simulate concrete damage and failure behavior of dam under earthquakes. In this research, taking Koyna gravity dam as the example of study, the wavelet energy distribution characteristics for the seismic response are analyzed based on the plastic damage model. The Koyna gravity dam is located on the Koyna River in the west of the Indian Peninsula. The concrete gravity dam is 103 m high and 70 m wide at its base; section size is shown in Figure2. Upstream water head was 91.7 m when a 6.5 magnitude earthquake shook the region in 1967. The Koyna earthquake with maximum acceleration measured at the foundation gallery of 0.47 g in horizontal direction normal to the dam axis was recoded. The dynamic excitation is represented by accelerogram whose horizontal component is time histories shown in Figure 3. The calculation parameters of materials are selected to be as follows: elasticity modulus E
=
3.0
×
1
0
4 Mpa, Poisson ratio ν
=
0.2, concrete density ρ
=
2630 kg/m3, concrete compressive strength 24.1 Mpa, and concrete tensile strength 2.9 Mpa.Figure 2
Dam profile.Figure 3
Input seismic waves.The wavelet energy distribution and variation characteristics for dam crest displacement are analyzed by mainly considering the following two cases: (1) regardless of the dam-water interaction (case A); (2) considering the dam-water interaction (case B). Studies have shown that the response of the dam is not sensitive to the foundation stiffness [17]. So the rigid foundation has been considered for the model in all of the analyses. The base of the dam has been tied to the foundation, and there is no sliding at the dam-foundation interface. This model was adopted for the seismic damage analysis of concrete dam by [2, 18]. The dam-reservoir dynamic interactions can be modeled in a simple form using the Westergaard added mass technique. The added mass per unit area of the upstream wall is given in approximate form by the expression [19] 7
/
8
γ
w
h
(
h
-
y
), where γ
w denotes the density of water and h the height of reservoir water. Rayleigh proportional damping is assumed in dam dynamic analysis, taking damping ratio ζ
=
0.05. Concrete material damage is defined by the relationship of damage factor and cracking displacement [20], as shown in Table 1. The 4-node bilinear plane strain quadrilateral reduced integration elements are used for discretizing finite element meshes of the dam.Table 1
Concrete plasticity calculation parameters.
Damage factor
0.00
0.381
0.617
0.763
0.853
0.909
0.944
0.965
0.979
0.987
0.992
Cracking displacement (10−4 m)
0.00
0.66
1.23
1.73
2.20
2.65
3.08
3.51
3.94
4.38
4.82The dynamic responses of the dam are calculated under the same earthquake acceleration waves with different Acceleration Scale Factors (ASF) and displacement signals of the dam crest are extracted. The signals are decomposed into several different frequency bands by using wavelet analysis. The following four factors should be considered in order to select reasonable wavelet: tight branch, symmetry, regularity, and vanishing moments. In order to calculate conveniently, the db5 wavelet is selected here. Because sampling frequency of the signals is 50 Hz, according to the sampling theorem, the highest frequency of the signal is 25 Hz. By a preliminary FFT transform, it is known that the amplitudes of the Koyna wave are mainly concentrated in frequency range 0~26 Hz. The first natural frequency of the dam is 3 Hz or so after the preliminary calculation. The signals are decomposed into 3 levels by using the db5 wavelet of Daubechies wavelets. The frequency band range decomposed is shown in Table2, from which it can be seen that the lowest frequency band is in 0~3.125 Hz containing the first natural frequency of the dam. By decomposing displacement signals into 3 levels, the relative wavelet energy of all frequency bands can be calculated using formula (5) as shown in Table 3. It can be seen from Table 3 that, for case A, most of the relative wavelet energy is concentrated in the frequency range 0~6.25 Hz, namely, in the third decomposition level. When the ASF are 0.1~0.4, no crack appears to the dam and the relative wavelet energy is no change. This shows that relative wavelet energy and the damage of the dam have the very good consistency. Along with seismic amplitude increase, the damage of the dam increases and stiffness decreases, leading to the decrease of the natural frequency of the dam and the main reaction energy drifting to low-frequency range which can be proved by the data in Table 3. When the ASF arrives to 1.4, the low-frequency relative wavelet energy γ
4 increases a little. This means that the absorption energy in low-frequency range γ
4 becomes saturated, namely, the dam is very close to (or arrives at) the limit state of seismic bearing capacity. And this point can be verified from the corresponding damage cloud picture of the dam (Figure 4). In Figure 4 there is a table indicating DAMAGET. It can be seen that the tensile damage variable in the black areas is close to 1. Damage almost passes through the dam neck, causing the final collapse of the dam.Table 2
Natural vibration frequency of Koyna dam.
Frequency order
1
2
3
4
5
6
7
8
Frequencyω
i (Hz)
3.062
8.0735
11.071
15.827
24.096
24.346
24.525
29.825Table 3
Relations of the acceleration scale factor and relative wavelet energy (intact dam).
ASF
Relative wavelet energy
γ
1
γ
2
γ
3
γ
4
Case A
Case B
Case A
Case B
Case A
Case B
Case A
Case B
0.1
1e − 4
2.04e − 4
0.0044
8.51e − 4
0.4677
0.0072
0.5278
0.9917
0.2
1e − 4
2.53e − 4
0.0044
0.0026
0.4677
0.0252
0.5278
0.9719
0.4
1e − 4
3.92e − 4
0.0044
0.0078
0.4677
0.0769
0.5278
0.9149
0.6
2.51e − 4
4.81e − 4
0.0055
0.0074
0.2851
0.1026
0.7091
0.8895
0.8
3.02e − 4
2.2e − 4
0.0037
0.0040
0.1561
0.0561
0.8399
0.9397
1.0
1.89e − 4
1.06e − 4
0.0020
0.0018
0.0671
0.0275
0.9307
0.9706
1.05
1.57e − 4
1.81e − 4
0.0017
0.0018
0.0394
0.0259
0.9587
0.9721
1.2
1.50e − 4
1.51e − 4
0.0017
0.0020
0.0272
0.0161
0.9709
0.9817
1.25
1.10e − 4
1.48e − 4
0.0013
0.0018
0.0240
0.0151
0.9746
0.9829
1.30
1.31e − 4
1.45e − 4
0.0015
0.0014
0.0191
0.0150
0.9794
0.9832
1.40
1.03e − 4
0.0015
0.018
0.9804
1.45
6.85e − 5
0.0013
0.0178
0.9807Figure 4
Damage distribution at ASF 1.4 (case A).From Table3, it can be seen that, for case B, when the ASF are 0.1~0.4, relative wavelet energy is mainly concentrated in the low-frequency range. This shows that the hydrostatic pressure plays a main role. When the ASF increases from 0.6 to 1.30, the relative wavelet energy is concentrated in the low-frequency range and becomes saturated at ASF 1.25. By comparing the two cases in Table 3, it can be seen that the seismic bearing capacity of the dam decreases slightly, considering the interaction of dam-reservoir. Figure 5 shows the damage cloud picture of limit state of seismic bearing capacity that is almost the same as Figure 4. Figure 6 gives the maximum principal stress distribution of the dam body under ASF 1.25.Figure 5
Damage distribution at ASF 1.25 (case B).Figure 6
The maximum principal stresses distribution at ASF 1.25 (case B).To check the correctness of the approaches presented in this paper, the results of cracking development of Koyna concrete gravity dam subjected to the same earthquake wave shown in Figure3 are compared in Figure 7. In the figure, the result in this paper is consistent with other researcher’s [21]. It can be seen from Figure 7 that the crack does not penetrate the whole section of the dam, indicating that the dam does not reach the limit state of bearing capacity. The data in Table 3 shows that, for case B, the dam reaches the limit state of bearing capacity when ASF is 1.25.Figure 7
Comparison of damage distribution and maximum crest displacement.
(a)
The result of this paper. Maximum crest displacement: 0.0507 m
(b)
Other researcher’s result. Maximum crest displacement: 0.047 mFor case A, Figure8 shows the input acceleration scale factor versus peak displacements. It can be seen in Figure 8 that the dam crest displacement amplitudes have no obvious increase with the increase of ASF. Figure 9 is the horizontal displacement time-history curve of the crest under the earthquake at ASF 1.4, from which it can be seen that displacement is convergent. Figures 8 and 9 indicate that the seismic bearing capacity of the dam cannot be determined by B-R theory [22].Figure 8
Relationship between the displacement amplitudes of dam crest and ASF.Figure 9
Horizontal displacement time history of the dam crest (case A).The dynamic response of the intact dam is analyzed and its seismic bearing capacity is acquired as in the previous one (shown in Table3). From the seismic code it can be known that for the dam withstanding the moderate-intensity earthquake, local damage may occur, but it has residual seismic capacity and can be used by performing a simple repair. Assuming that the dam body is subjected to seismic action of amplitude 280 gal (ASF 0.6 in Figure 3), the damage cloud picture of the dam (shown in Figure 10) can be calculated using ABAQUS Software. Figure 10 shows that the local damage occurs in the dam neck, showing nearly 6 m deep crack from the downstream surface at the macroscopic level. Using the ABAQUS restart function, the dynamic response of the dam with 6 m deep initial crack in its neck can be calculated under earthquake with different ASF. The relative wavelet energy γ
1~γ
4 (shown in Table 4) in 4 frequency bands can be obtained using the same wavelet processing methods for the dam crest displacement signals. It can be seen from Table 4 that, for case A, the dam body reaches the limit state of seismic bearing capacity when ASF is 1.15 (i.e., 543.95 gal). Based on the seismic bearing capacity of the dam, the global damage index can be defined as follows: (12)
D
I
=
1
-
R
R
0
,where R is the residual seismic bearing capacity of earthquake damaged dam and R
0 is the seismic bearing capacity of intact dam. For case A, it can be calculated that the global damage index of the dam is 0.18 using formula (12). For case B, the residual seismic bearing capacity of earthquake damaged dam is 1.1 multiples of seismic wave and the corresponding damage cloud picture is shown in Figure 11. By comparing the seismic bearing capacity of the intact dam (in Table 3) with that of an earthquake damaged dam (in Table 4), it can be calculated that the global damage index of the dam is 0.12 using formula (12). Therefore, the global damage index of the gravity dam subjected to the known seismic action is over 10% whether considering the dam-water interaction or not. This indicates that the gravity dam can be used by performing simple repairs.Table 4
Relations of ASF and relative wavelet energy (earthquake damaged dam).
Seismic multiples
Relative wavelet energy
γ
1
γ
2
γ
3
γ
4
Case A
Case B
Case A
Case B
Case A
Case B
Case A
Case B
0.4
8.16e − 5
1.07e − 4
0.0059
0.0036
0.3465
0.3494
0.6475
0.6469
0.6
9.57e − 5
1.96e − 4
0.0052
0.0035
0.2299
0.3046
0.7648
0.6916
0.8
2.05e − 4
2.59e − 4
0.0046
0.0076
0.1269
0.1974
0.8683
0.7947
1
1.64e − 4
2.29e − 4
0.0028
0.007
0.0569
0.1291
0.9402
0.8636
1.1
1.73e − 4
2.85e − 4
0.0027
0.0045
0.0466
0.1007
0.9506
0.8945
1.15
1.37e − 4
3.04e − 4
0.0028
0.0062
0.0293
0.0988
0.9677
0.8948
1.25
1.58e − 4
0.002
0.0301
0.9678Figure 10
Damage distribution at ASF 0.6.Figure 11
Damage distribution at ASF 1.1 for earthquake damaged dam.


## 5. Conclusion

The following conclusions can be drawn from the above theoretical and numerical analyses.(1)
The vibration signals of the dam are decomposed into several different frequency bands by using wavelet analysis; thus energy distribution of the subband can be seen clearly. With the increase of the earthquake wave magnitude, the damage of the dam increases, leading to the energy concentration in the low-frequency range. The dam body gets to the limit state of seismic bearing capacity when the low-frequency energy reaches saturation. Therefore, the corresponding numerical value of dam seismic capacity can be obtained.
(2)
The damage distribution of the dam under the action of earthquake can be simulated and predicted by means of a concrete damaged plasticity model. Considering the damage as the initial damage of the dam, the numerical value of residual seismic bearing capacity for the earthquake damaged dam can be calculated by using the wavelet analysis method. The overall damage degree of the dam with initial damage can be evaluated by comparing the changes of seismic bearing capacity before and after earthquake damage. Considering the dam-water interaction, the global damage index slightly reduces and the seismic bearing capacity of concrete dam lows a little.
(3)
In the paper, we proposed a simple and feasible way for the earthquake damage assessment.



---

*Source: 6841836-2016-03-08.xml*