**S ** **GLOBAL RESEARCH PUBLICATION Vol.1, No. 1, 2007** EISMIC ANALYSIS OF EARTH WALL GRAVITY DAMS USING DECOUPLED MODAL APPROACH
Adedeji, A. A. Department of Civil Engineering, University of Ilorin, PMB 1515, Ilorin, Nigeria. email: __amadeji@yahoo.com__
## ## ## ## **ABSTRACT**
This study aimed at the assessment of the effects of seismic loads on the earth dam for design purposes. The Unilorin concrete gravity dam is used in comparison with the earth wall gravity dam in seismic zone 1 with a peak ground acceleration of 0.05g. Finite element (FE) method of analysis was used by employing Lagrangian-Eulerian formulation of 4-node plain quadrilateral elements, with modal analysis used to decouple system dynamic equations, thereby saving computational time. The loadings were determined based on EM 1110-2-22000(1995), while the FE model is being implemented using the MATLAB programming tool. Though, the collapse and response of the earth gravity dam is higher than that of concrete gravity dam, the dam satisfies the stability and stress requirements in design. **KEY WORDS** *Sismic analysis, earth wall gravity dam, Finite element method*
**INTRODUCTION**
The structural response of a material to different loads determines how it will be economically utilized in the design process. The loads that structures are subjected to are either static or dynamic, or a combination of both. Static loads are time-independent while dynamic loads vary with time. Dynamic loads can be classified as cyclic, impact or moving but can occur simultaneously. Earthquake is a natural disaster that has claimed so many lives and destroyed lots of property. Earthquake hazards had caused the collapse and damage to continual functioning of essential services such as communication and transportation facilities, buildings, dams, electric installations, ports, pipelines, water and waste water systems, electric and nuclear power plants with severe economic losses. Earthquake is a major source of seismic forces that impinge on structures others are Tsunami, seethe etc. Earth wall is chosen as a material for the dam since its major constituent- earth is abundantly available and provides a sustainable solution. This necessitates the seismic analysis of concrete gravity dam. Finite element has been widely used in seismic analysis of concrete gravity dams (Waltz 1997, Lotfi 2003) with a defined approach as presented in this programme, using the most natural method based on the Lagrangian–Eulerian formulation. Smith (1985) claimed that dams are critical structures that should be made earthquake-resistant, since earthquakes cause severe damages, consequent huge economic and life losses. Hence, the need to prevent the occurrence of these earthquake hazards by carrying out seismic analysis of dams (Polyakov, 1985). Over the years, a lot of work has been done in making earth-fill and concrete dams earthquake-resistant, with advances in structural vibration and finite element method which have aid in the seismic
analysis of theses dams. However, not much attention has been paid to the seismic analysis of an earth wall dam. More importantly, most recent methods in seismic analysis of concrete gravity dams have not been employed for the seismic analysis of the earth wall dams. Therefore, the purpose of this is to apply the decoupled modal approach in the seismic analysis of Earth wall dams. Earthquakes had caused severe damages and consequently huge economic losses including losses of lives. On 27th December, 2004, The Punch newspaper reported that Tsunami, a seismic sea wave claimed a lot of lives, rendered many homeless and destitute, and destroyed essential services with Sri Lanka,, Indonesia, Somalia, India, etc been the major casualty countries. Nigeria is not left out from the occurrence of ground motions caused by earth movement as earth tremors had been recorded to occur some recent years back in Abeokuta and Ibadan. Most structures and infrastructure in the country are not earth-resistant. Also, earth wall, as a sustainable material is quite economical in construction of dams since there is abundant availability of good earth or soil, especially in localities closed to road construction and excavation sites. This study involved the use of finite element method for the dynamic analysis of both the dam and reservoir bodies, using a modal approach as the Lagrangian-Eulerian formulation forming the basis. Also, static loads (weights, hydrostatic pressures) are each visualized as being applied in one separate increment of time. The analytical computation of the modal approach procedure has been carried out and implemented using MATLAB programming tool. The pseudo-static seismic coefficient method was adopted in computing the seismic loads on the dam. The dam used as a case study was assumed to be in seismic zone 1 with seismic coefficient ranging between 0.0 and 0.05. The dam was analysed seismically using the decoupled modal approach and the results were compared with that of the concrete gravity dam.
**2. Seismic Response of Earth and Concrete gravity dams**
In earth dams, seismic forces or shaking can induce destabilising deformation or outright failure if not made earthquake resistant. A permanent simplified procedure can be adopted to estimate permanent horizontal displacements of the dams using finite element method that account for non-linear material behaviour and strength reduction due to liquefaction or strain softening. It has been shown ((Hatami, 2001) that the seismic performance of earth dams has been related to the nature and state of compaction of the fill material. Concrete dam’s structural safety and stability are jeopardized due to the hydrodynamic load of the reservoir that is subjected to ground motion
**2.1 Fluid –Structure System**__ __
During earthquake occurrence, the dam and reservoir body respond differently, as a result of hydrodynamic forces impinging on the fluid body and solid structure. As a result of this, interaction will occur between the fluid–solid structure interfaces as particles move relatively to the mesh points whereas, the meshes moves with the material particles (Bathe, 1996, Qixiang et al. 2000). Much research work has been carried out for the dynamic response of the fluid-solid structure systems. Several methods of analysis for the fluid-structure systems use finite element idealization in the non-linear dynamic response of the system (Fenves and Vargas- Loli, 1988).
In Fig.1(a, b) the following terms are defined as:
**3. LOADINGS **
**3.1 Static loads.**
The static loads are due to (i) The weight of the dam: the unit weight is assumed to be 19.62kN/m^{3} until an exact unit weight is determined from materials investigation., (ii) Hydrostatic pressure of the water in the reservoir and (iii) The uplift forces caused by hydrostatic pressure on the foundation at the interface of the dam and the foundation. Uplift forces are usually considered in stability and stress analysis to ensure structural adequacy and are assumed to be unchanged by earthquake forces.
**3.2 Dynamic loads**
Earthquake or seismic loads are the major dynamic loads (Major 1980, Schoeber 1981, Polyakov 1985, Wyatt (1989). being considered in the analysis and design of dams especially in earthquake prone areas. They are of kinematic origin and owe their existence to vibration caused in the structure by the movement of the earth’s surface during an earthquake. They have random characteristics and are regarded as deterministic in practical calculations to simplify the design model. According to EM1110- 2-2200(1995), the earthquake loading used in the design of concrete gravity dams are based on design earthquakes and site specific motions determined from seismological evaluation. The seismic coefficient method is used in determining the resultant location and sliding stability of dams. Seismic analysis of dams is performed for the most unfavourable direction, despite the fact that earthquake acceleration might take place in any direction. Seismic coefficient method of analysis is commonly known, as pseudostatic analysis, and is the ratio of the earthquake acceleration to the acceleration due to gravity. Fig.2 shows the dynamic loads on a gravity dam. There are different ways of computing earthquake loads on dams. The deterministic approach will be employed where the ground acceleration in terms of **g** (acceleration due to gravity) is specified for the region where the dam will be constructed. Hence, the exciting force on the structure is,
P_{(t)} = Ma_{x } (1)
and a_{x} = αg (2)
where **a**_{x}**, α, g **are the ground acceleration, seismic coefficient and acceleration due to gravity respectively.
Fig. 2: Seismically loaded gravity dam
From Fig. 2 and equation ( ) therefore, the equilibrium system is expressed as:
Pe_{x } = Ma_{x } = W_{α}g /g = W_{α} (3a)
In which a_{x} = αg and W = Mg (3b)
along vertical direction Pe_{w } =(2 * C_{e }* α * y * √ (h *y)) / 3 (4a)
and C_{e} = 51 / √ (1 – 0.72 * (h / (1000t_{e}))^{2}) (4b)
where Pe_{x}, M, a_{x}, W, α, g** **are the horizontal earthquake force on the dam, mass horizontal earthquake acceleration, weight, acceleration due to gravity and seismic coefficient respectively. Also Pe_{w}, h, t_{e} are the additional total water load down to depth y, total height of reservoir, and period of vibration respectively.
**3.3 Forced – damped Vibration.**
F
orced vibration is the vibration caused by a time- dependent disturbing forced. The governing equation for forced-damped vibration is
MÜ + CŮ + KU = P_{(t) } (5)
where M, C, K, Ü, Ů, U, and P_{(t)} are the mass, damping , stiffness, acceleration, velocity, displacement and exciting force on the body.
**4. Finite Element Method (FEM).**
The finite element method is the best method used in analysing the internal forces of element(s) by solving differential equation obtained through their discretisation in space dimensions (Adedeji, 2004). Lotfi (2003) and Iroko (2001) stated that finite element method has been the most widely used method for seismic analysis of concrete gravity dams. Finite element (FE) models are used for linear static and dynamic analyses and for non-linear analyses that account for interaction of the dam and foundation. The FEM provides the capability of modeling complex geometries and wide vibrations in material properties. Two-dimensional finite element analysis is generally appropriate for concrete gravity dams. Lotfi (2003) argued that the Lagrangian-Eulerian formulation is the most natural method, and employs nodal displacements and pressure degrees of freedom for the dam and reservoir region respectively. Bathe (1996) had previously put down this idea stating that in this type of formulation, the mesh points move but not necessarily with the material particles, due to fluid-solid coupling and interaction (Fenves and Vargas, 1988) that allows the approach to be used in modelling the fluid flow-structures interaction. The direct integration methods are used to solve the linear dynamic response of a system of finite elements There are different methods used in the direct integration process and are the Central difference, the Houbolt, the Wilson, the coupling of different integration operators and Newmark methods (Bathe, 1996). Lotfi (2003) claimed, that the direct integration process is usually carried out by applying Newmark’s algorithm or method, with the following assumptions used by (Bathe, 1996), are
^{t +Δt}Ů = ^{t}Ů + [(1-δ)^{t}Ü + δ^{t +Δt}Ü] Δt (6)
^{t +Δt}U =^{ t}U + ^{t}ŮΔt + [(0.5-α)^{t}Ü + α^{t +Δt}Ü] Δt^{2} (7)
where** **t, Δt, ^{t}Ü, ^{t}Ů, ^{t}U, δ, α** **are time, time step of integration, acceleration at time **t**, velocity at time **t**, displacement at time **t**, and integration parameters respectively. However, Lotfi (2001) reported that a non- symmetric system of equations to be solved at each time step will be encountered and not direct method by Newmark’s algorithm will not be computationally efficient.
**4.1 Decoupled Modal Technique **
Bathe (1996) stated that if the integration will be carried out for many time steps, it will be more effective to first transform the equilibrium equations of equation (8) into a form in which the step-by-step solution is less costly. Based on this, Lotfi (2003) applied the modal approach to obtain the natural frequencies and mode shape of the fluid–structure system by rewriting the equilibrium equation (5) as:
KΦ = ω^{2}MΦ (8)
where **ω,** **Φ, K, M **are** **the** **eigenvalues, eigenvectors, stiffness and mass of the system. The above equation is the generalised eigenproblem from which Φ (mode shapes or eigenvectors) and **ω** (frequencies or eigenvalues) are determined (Bathe, 1996). Thereafter, the solution at different time steps can be estimated based on the combination of the modes at the different time steps (Lotfi, 2003).
**4.2 Structural vibration**
It is well known that, the linear dynamic equation for the solid structure (dam), is given as in equal and are reproduced here.
MÜ + CŮ + KU = P_{(t) } (9) ##### Therefore, the linear dynamic equation for the fluid region (reservoir) is written as,
GÜ + LŮ + HU = R_{(t) } (10)
where M, C, K, Ü, Ů, U, and P_{(t)} are the mass, damping , stiffness, acceleration, velocity, displacement and exciting force on the dam. G, L, H, Ü, Ů, U, and P_{(t)} are the mass, damping , stiffness, acceleration, velocity, displacement and exciting force on the reservoir.
**4.3 Lagrangian-Eulerian formulation**
For the Lagrangian-Eulerian formulation for the fluid-structure system, both the Lagrangian and Eulerian formulations are coupled together. The Lagrangian formulation as stated before is used for formulating the solid domain as given in equation (11)
∫o_{V} ^{t}S_{ij} δ ^{t}ε_{ij} d^{0}V = ^{t}R (11)
and the Eulerian formulation is employed for the fluid domain,
^{t}τ_{ij} = -^{t}*p* *δ_{ij} + 2μ* ^{t}e_{ij }(12a)
^{t}e_{ij} = (^{t}ν_{i,j} + ^{t}ν_{j,i})/2 (12b)
where ^{t}S_{ij}, ^{t}e_{ij}, ^{0}V, ^{t}R, ^{t}τ_{ij}, ^{t}*p, *δ_{ij}, μ** **are the Piola-Kirchoff stress tensor at time **t, **Green-Lagrange strain tensor at time **t**, volume at time **t** = 0, Cauchy stress tensor at time **t**, pressure displacement at time **t**, Kronecker delta and dynamic viscosity respectively.
**Finite Element Model of a Dam-Reservoir System**
Discretisation of the dam-reservoir system is shown in Fig. 3.1 while in Fig. 3.2 shows two-dimensional plane solid and fluid quadrilateral continuum elements used in the discretisation of the coupled dam-reservoir system. Linear Dynamic Analysis for the Coupled Dam-Reservoir System is shown in details by Lotfi (2003), while the damping matrix is totally symmetric and the following relation holds.
K_{u} = -M_{u}^{T}_{ }(13)
It will be computationally inefficient to use the direct integration process (usually carried out by Newmark’s algorithm) since a non-symmetric system of equations to be solved at each time step will be encountered. On the other hand, the modal approach depends on obtaining the natural frequencies and mode shapes of the system.
**4.5 Decoupled Modal Approach**
The eigenvalue problem corresponding to equation 3.7 is
**K X**_{j} = λ_{j}**M** **X**_{j}_{ } (14)
where λ_{j} = eigenvalue of the system and eigenvector matrix It is clear that the eigenvalues of this system are real and free vibration modes exist. However, it is noted from the form of matrices **K** and M that the system is not symmetric and standard eigenvalue computation method are not directly applicable. It would be computationally costly since additional variables are required in using the available techniques to arrive at a symmetric form and reduce the problem to a standard eigenvalue one. On the contrary, eigenvalues and vectors extracted from the following eigenproblem will be used as:
**K**_{s}**X**_{j}** **= λ_{j}**M**_{s}**X**_{j}_{ }(15)
The eigenvectors of the (15) are not the true mode shapes of the coupled system but can be presumed as Ritz’s vectors and be combined to estimate the true solution. Standard methods are used to obtain the solutions of this substitute eigenproblem. Using the orthogonal condition and normalizing the modal matrix with respect to mass matrix, one could have:
**X**^{T}**M**_{s}**X ** = **I** (16a)
**X**^{T}**K**_{s}**X** = **Λ** (16b)
where **I** is the identity matrix and **Λ** is a diagonal matrix containing the eigenvalues of the symmetric substitute system. The following relations are derived from relations in equations (15) and (16).
**X**^{T}**MX** = I +**X**^{T}**M**_{u}**X ** (17a) **X**^{T }**KX = Λ – X**^{T}**M**_{u}^{T}**X** (17b)
As in modal techniques, the solution is written as a combination of different modes.
**U = XY ** (18)
The vector **Y** contains participation factors of the modes.** **Applying the Newmark’s technique for each new time step,
**KY**_{n+1} = **F**_{n+1 }(19)
**K** and **F**_{n+1} are denoted as the generalized effective stiffness matrix and the generalized effective force vector of the system at Time step n+1 respectively. They are defined as:
**K **= a_{0}(I + **X**^{T}**M**_{u}**X**) + a_{1}**C*** + (**Λ - X**^{T}**M**_{u}^{T}**X**) (20)
And, **F**_{n+1} + **F***_{n+1 } + ( I + **X**^{T}**M**_{u}**X**)(a_{0}**Y**_{n} + a_{2}**Y**_{n} + a_{3}**Y**_{n}) + **C***(a_{1}**Y**_{n} + a_{4}**Y**_{n} + a_{5}**Y**_{n}) (21)
In general, the vector of participation factors can be solved through relation in equation (19). Thereafter, the unknown vector is obtained by equation (15) as usual in the modal procedure, it must be also mentioned that the generalized effective stiffness matrix employed in equation (19) is inherently unsymmetrical. However, it may be easily transformed to a symmetric matrix by multiplying certain rows of the matrix in equation (19) by an appropriate factor. The symmetric matrices were utilized in the substitute eigenproblem that corresponds to the decoupled dam-reservoir system. Therefore, the natural frequencies and eigenvectors are actually related to this decoupled system. In actual programming, one can modify the usual subspace iteration routines to converge to the desired lowest modes of the dam first, and similarly for the finite reservoir region afterwards by an appropriate selection of initial vectors. Meanwhile, they can also be obtained as two separate problems. Assuming for simplicity, without loss of generality that the mode shapes for the dam are ordered first and the ones for the finite reservoir considered subsequently in the modal matrix. All similar arrangements for the eigenvalues in the diagonal matrix Λ and the relative matrix formations are expressed by Lotfi (2003). It is assumed that the damping matrix of the dam is of viscous type for the analysis carried out by modal approach in time domain. Therefore, this leads to following relationship, so that the structural damping matrix for the stress-strain relationship.
**C**_{1}* = 2**β**_{d}**Λ**^{1/2 } (22)
where **β**_{d} is the equivalent damping factor, which is assumed constant for all modes. Substituting equation (21) into equation (20) and (21), the expanded form of equation (19) is now concluded as ,
a_{0}**I**_{1}** + **a_{1}**C**_{1}*** + Λ**_{1}** -X**^{T}**B**^{T}**X**_{2}** Y**_{1}** F**_{1}
** = ** (23) a_{0}**X**_{2}^{T}**BX**_{1}** **a_{0}**I**_{2}** + **a_{1}**X**_{2}^{T}**LX**_{2}** + Λ**_{2}** Y**_{2}** **_{n+1}** F**_{2 n+1}
The vector of the participation factors is also assumed to be partitioned into two parts in this relation, and as before the indices 1, 2 correspond to dam and reservoir modes respectively. Equation (23) is equivalent to (19), which is initially an unsymmetrical system of equation. This matrix relation could become symmetric by multiplying the lower matrix equation by a factor of –**a**_{0}^{-1}, which yields the following equation.
a_{0}**I**_{1}**+ **a_{1}**C**_{1}***+ Λ**_{1}** -X**^{T}**B**^{T}**X**_{2}** Y**_{1}** F**_{1}
** = ** (24)
- **X**_{2}^{T}**BX**_{1}** **–a_{0}^{-1} (a_{0}**I**_{2}**+**a_{1}**X**_{2}^{T}**LX**_{2}**+Λ**_{2}**) Y**_{2}** **_{n+1}** **–a_{0}^{-1}**F**_{2 }_{ n+1}
The above equation (24) can now be readily solved for the vector of participation factors at each time step, and the original unknown vector is obtained from equation (15).
**5.** **ANALYSIS AND DISCUSSION OF RESULTS**
**Case Study: A Typical Analytical Example**
To verify the capability of this numerical procedure, the on-going construction of University of Ilorin (Unilorin) dam is considered. And since the country is generally less prone to earthquakes or tremors, the seismic zone of this area is taken as Zone 1 whose seismic coefficient ranges from 0 to 0.05g. The dam was previously designed as a concrete gravity dam whose data is utilized in this study. Fig. 3 shows the dam geometry which is to retain 1.8 million cubic meters of water with a reservoir length and height of 200m and 8m respectively.
**5. 2 Material Properties and Basic Data**
The concrete and earth wall is assumed to be homogenous and isotropic with the following basic properties in Table 1. Table 1. Properties of materials -
**Property of material** |
**Concrete ** |
**Earth** | Elastic modulus, E_{c} (GPa) Poisson ratio, ν_{c} Unit weight, γ_{c} (kN/m^{3}) Dry density γ_{d} Damping coefficient β_{d} | 22.75 0.20 24.80 - 0.05 | 21.2 0.15 19.62 18.20 0.05 |
**Finite Element Analysis and Formulation for the System**
In Fig. 4 a typical numbering of the dam system indicates the fluid- and structural wall domain, while an isoparameter of an element showing fluid-structural members domain is shown in Figs. 5 to 6.
Fig. 4. Numbering of the elements in both the fluid and structure domains
Fig. 5: Isoperimetric element definition
For the reservoir, the two element types are: ** **
Fig. 6: Element types in the fluid domain
For the dam body, the element types are as shown in Fig. 7. Fig. 7: Element types of the dam domain
**5.4. Programming the FEM model**
The steps carried out above are followed in the programming process. Figs. 8 and 9 below shows the MATLAB environment in which all the steps were performed.
Fig. 8: Function M-file for implementing the fluid domain model
Fig. 9: Function M-file for implementing the solid domain model
**5.5 Seismic Response of the Earth Wall and Concrete Gravity Dams **
The analysis is carried out by modal approach. The eigen-problem is solved in the first step based on the symmetric parts of the coupled equation of the system. The first two mode-numbers were considered to obtain the first two natural frequencies of both the earth wall and concrete gravity dam as shown in Table 2 |