# Numerical Uncoupling of Domains in Dam-Reservoir Problem

Authors: Saba Golchin; Reza Attarnejad; Shahram Vahdani
Journal: Shock and Vibration
 (2018)

---

## Abstract

Flexibility of dam structure affects the hydrodynamic pressure acting on the dam. Several approaches have been proposed to consider this effect. Most of these approaches are involved with an iterative scheme. Of course solving the total numerical model including the dam and the reservoir is the most accurate method, but it has certain deficiencies. Using the frontal solution method of total model, dam structure, and fluid domain and keeping the interface degrees of freedom in the front is proposed in the current study. Having the solution of the interface degrees of freedom, the structure and fluid may be analyzed separately. The main advantage of the method lies in the fact that the accuracy of the results is the same as analysis of the total model, no iteration is necessary, combination of Lagrangian and Eulerian formulations for solid and fluid may be used, and the unknown variables are of the same order. Performing the analysis in time domain extends the method to nonlinear analysis if required.


---

## Body

## 1. Introduction

Although the dams are very stiff structures and experience small amount of deformations during earthquakes, these deformations are the source of drastic change in hydrodynamic pressure. Several approaches have been proposed to solve the coupled equations. One approach to include this interaction, uses iterative schemes, where first the dam geometry is assumed to be constant when the reservoir is solved. Then obtained pressure will be applied on the dam, and new configuration of dam is calculated. The iteration will be continued until the convergence is achieved. Fellippa and Park [1] described staggered solution for coupled problems in detail. Staggered approach reduced coupled problem to subsystems [2]. Therefore, only symmetric equations for each subsystem is calculated [3–5]. However, the stability of the method is conditional and must be cared about convergence of the solution of the subsystems at each time step. This approach needs iteration in each time step even in linear problems, so this approach is involved with number of solutions. The methodology may consider the effect of interaction but it cannot be exact solution for interaction problems.In the other approach, a complete model of the dam and reservoir is formed. Moreover, the two formulation methods of Eulerian (suitable for reservoir formulation) and Lagrangian (suitable for structure formulation) have to be combined or one of them has to be selected. Wilson and Khalvati [6], Akkose, Bayraktar, and Dumanoglu [7], and Akkose and Simsek [8] used the Lagrangian approach, and displacement was taken as a variable of the solid and fluid domain. In Eulerian approach, displacement was taken as a variable of the solid domain and pressure or velocity was taken as a variable of the fluid domain [3, 4, 9].Both of the abovementioned solution methods may be performed in frequency or time domain. A lot of work has been done by Chopra and coadjutor on dam-reservoir-foundation interaction in frequency domain covering several aspects of the problem [10–13]. Frequency domain approach is computationally very efficient. However, concepts of the frequency domain formulation are more difficult than those of time domain formulation. The main disadvantage of the frequency domain is that it cannot be used in nonlinear problems.Simulation method in time domain needs massive storage space because of solving all of equations including structure and fluid as well as far field variables for each time step. Furthermore, ill conditioning is dramaticaly conflictive and occurs during solving the equations of very different order variables.Many analytical research studies have been carried out about dam-reservoir system. Tsai, Lee, and Ketter [14] first presented a semianalytical method to express the far field radiation condition in reservoir. Lee and Tsai [15] also developed a closed form solution for fluid domain. Miquel and Bouaanani [16] first presented a simplified evaluation for seismic response of the dam-reservoir systems, and then they used this simplified approach to determine the modal dynamic response of flexible beam-fluid systems [17]. These analytical solutions are really efficient in time but they have geometrical and boundary condition limitation and also they cannot be used in nonlinear problems.The proposed solution method in this study is a branch of complete model solution with the exception that only the selected degrees of freedom on the interface of dam and reservoir will be solved in time domain. The principal aim of this study is to employ frontal solution for solving the dam-reservoir system. The frontal method for solving the equations was presented by Irons [18] for the first time. By using frontal solution, all equations are assembled but only a selection of the degrees of freedoms is solved. In the other words, by considering the interface pressure degrees of freedom and keeping them as the front, total model including the displacement of the structure and the pressure of the fluid is assembled by the frontal method. Then only the interface degrees of freedom in front, which is pressure acting on the interface, is solved. After finding the values of the pressure which act on the interface of the dam, the structure can be analyzed due to seismic force and dynamic pressure, and the fluid can be analyzed due to determined boundary conditions, separately.The achievements can be explained as follows:(1)
In the first step of the scheme, only the value of the pressure which acts on the interface is solved without losing any accuracy comparing with other solution methods which solve the entire system(2)
In the second step of the scheme, Lagrangian or Eulerian approach can be used, which was better, for the solid domain and the fluid domain(3)
This scheme is capable to model the nonlinear problems if time domain formulation is adopted


## 2. Governing Equation of the Coupled Dam-Reservoir Problem

The dam-reservoir system is classified as a coupled problem where two fields each governed by a second-order differential equation, interact at their interface. The equation of motion for dam structure due to earthquake motion is given as(1)Mu¨+Cu˙+Ku=FI−Mu¨g+Qp,where M, C, and K are mass, damping, and stiffness matrices of the dam structure, respectively; FI is the static load vector; u¨g is the ground acceleration; u, u˙, and u¨ are vectors of displacement, velocity, and acceleration of structure body, respectively; p is the vector of hydrodynamic pressure; and Q is coupling matrix defined as below(2)Q=∫NsnNfTds.In above integration on the interface of the two domain, vectorn is normal vector on interface and defined to be positive, going from the solid domain into the fluid domain. Ns and Nf are shape function of structure and fluid, respectively.M and K are defined as below(3)M=∑∫ΩsNsTρsNsdΩ,K=∑∫ΩsBTD′BdΩ,where Ωs denotes the structure domain, ρs is mass density of the structure, B is strain matrix, and D' is elastic constants matrix.In this study, Rayleigh damping is assumed, and soC is taken as(4)C=αsM+βsK,where αs and βs are constants that can be calculated from the below relation(5)αs+βsωi2=2ωiζi,where ζi is the damping ratio and ωi is ith natural frequency of the structure.The equation of motion for reservoir domain due to earthquake motion is given as [4](6)Gp¨+Dp˙+Hp=−ρfQTu¨+u¨g,where G, D, and H are quasimass, damping, and stiffness matrices of the reservoir domain, respectively; ρf is the density of the fluid; and p˙ and p¨ are vectors of first and second derivation of the hydrodynamic pressure. G, D, and H are defined as follows:(7)G=∑1cw2∫ΩfNfTNfdΩ+∑1g∫sfNfTNfds,D=∑1cw∫stNfTNfds+∑1βbcw∫sbNfTNfds,H=∑∫Ωf∂∂xNfT∂∂xNf+∂∂yNfT∂∂yNfdΩ+∑π2h∫stNfTNfds,where βb is the acoustic impedance ratio of the reservoir bottom material, cw is velocity of sound in the water, h is depth of the reservoir, and Ωf denotes the fluid domain. sf, sb, and st are surface of the reservoir, bottom of the reservoir, and truncation boundary, respectively.Equations (1) and (6) can be written in matrix form as follows:(8)M0ρfQTGu¨p¨+C00Du˙p˙+K−Q0Hup=FI−Mu¨g−ρfQTu¨g.Equation (8) can be rewritten for dynamic loads as follows:(9)M0ρfQTGu¨p¨+C00Du˙p˙+K−Q0Hup=−Mu¨g−ρfQTu¨g.By using step by step integration scheme of Newmark [19], the displacement and hydrodynamic pressure and their first derivative can be discretized as follows:(10)un+1=un+Δtu˙n+Δt221−2βu¨n+2βu¨n+1,u˙n+1=u˙n+1−αu¨n+αu¨n+1Δt,(11)pn+1=pn+Δtp˙n+Δt221−2βp¨n+2βp¨n+1,p˙n+1=p˙n+1−αp¨n+αp¨n+1Δt,where n indicates the number of the time steps and α and β are integration parameters which determined stability and accuracy of the Newmark method. In this study, α and β are chosen as 0.5 and 0.25, respectively. Substitution of Equation (10) and (11) in Equation (9) results(12)4Δt2M+2ΔtC+K−Q4Δt2ρfQT4Δt2G+2ΔtD+Hun+1pn+1=FsFff,where(12a)Fs=−Mu¨g+M4Δt2un+4Δtu˙n+u¨n+C2Δtun+u˙n,(12b)Fff=−ρfQTu¨g+ρfQT4Δt2un+4Δtu˙n+u¨n+G4Δt2pn+4Δtp˙n+p¨n+D2Δtpn+p˙n.Multiplying the reservoir motion equation by−βΔt2/ρf renders the asymmetric coupled matrix of Equation (12) to a symmetric one as follows:(13)4Δt2M+2ΔtC+K−Q−QT−1ρfG−Δt2ρfD−Δt24ρfHun+1pn+1=FsFf,where(13a)Ff=Δt24QTu¨g−QTun+Δtu˙n+Δt24u¨n−G1ρfpn+Δtρfp˙n+Δt24ρfp¨n−DΔt2ρfpn+Δt24ρfp˙n.


## 3. Most Common Solutions for Coupled Dam-Reservoir Problem

Different approaches have been proposed to solve the coupled equations. Most popular approach is iterative schemes [4]. This approach needs iteration in each time step even in linear problems. A concept of iterative method at any time step can be summarized as follows, where i indicates iteration step:(1)
An arbitrary pressurep is substituted into Equation (1) to obtain ui, u˙i, and u¨i. It is better to select the calculated pressure from previous time step as arbitrary pressure in order to achieve rapid convergence.(2)
Calculatedu¨i is substituted in Equation (6) to obtain pi, p˙i, and p¨i.(3)
Calculatedpi is substituted in Equation (1) to obtain ui+1, u˙i+1, and u¨i+1.(4)
Calculatedu¨i+1 is substituted in Equation (6) to obtain pi+1, p˙i+1, and p¨i+1.(5)
Check for convergence: ifpi+1−pi and ui+1−ui are not lower than a reasonable tolerance, cycle 3 and 4 are repeated until convergence is achieved.This approach reduces coupled problem to subsystems. Therefore, only symmetric equations for each subsystem is calculated (Equations (1) and (6)). However, the stability of the method is conditional and care must be taken about convergence of the solution of the subsystems at each time step. This approach considers the effect of interaction, but it cannot be generalized as an exact solution for interaction problem.The direct numerical solution of dam and reservoir together in time domain may be considered as the best interaction solution for the following reasons:(1)
Except the theoretical solution, it is the most accurate one(2)
It does not have the geometrical and boundary condition limitations(3)
It is capable of considering the nonlinear phenomenonAlthough it has its own deficiencies,(1)
the solution is extremely expensive,(2)
it has to use uniform formulation; Lagrangian or Eulerian for both the dam and reservoir which Lagrangian and Eulerian formulation are suitable for the dam structure and the reservoir, respectively,(3)
solving the equations of different order of variables has its own concerns, although it is solvable.


## 4. Frontal Solution for Coupled Dam-Reservoir Problem

In common methods of matrix solution, the element’s information is gathered in element records and is kept on peripheral. Therefore, the element matrices, such as stiffness, are computed again any time needed. The advantage lies on the fact that almost all memory is free for assemblage of global stiffness.Frontal solution of large matrices is based on special assembly. As soon as the coefficients of an equation are assembled from the contributions of all relevant elements, the corresponding variable can be eliminated [18]. Therefore in frontal method, the complete matrices are neither generated nor stored. Of course, the results of the two methods are the same in linear analysis. Since there is no superposition in frontal method, it can be extended to nonlinear analysis.In frontal method, the equations which are being formed at any given instant, their corresponding nodes, and degrees of freedom are named the front. The number of the variables in the front is named front width. The equations, nodes, and degrees of freedom (DOFs) belonging to the front are named active; those which have passed through the front and have been eliminated are named inactive. Active nodes can be deactivated after the last appearance. By eliminating variables as soon as their assembly have been completed, core storage is made available for variables yet to be assembled. The advantage lies on the fact that required memory for the solution will be reduced to the front.Taking advantage of frontal technique to keep the pressure DOFs on the interface of dam and reservoir in the front can highly help the solution. The procedure of the assembly for structure and reservoir by frontal method is illustrated in Figure1.Figure 1
Frontal assembly procedure for two coupled domains. (a) The schematic model including fluid elements, structure elements, and nodes before beginning the assembly. (b) The first element of each media at the beginning of assembly. (c) The second element of each media for assembling. (d) Active nodes and deactivated nodes at the final step of assembly.The solution basis for frontal method is Gaussian elimination. The interaction equations of dam and reservoir contain solid DOFs ({u}) and liquid DOFs ({p}). Equation (13) can be rewritten as follows:(14)K¯X=F,where X=up,F=FsFf, and K¯ is effective dynamic matrix as defined in Equation (13).{p} can be divided to reservoir pressure DOFs pr and the interface pressure DOFs pf:(15)p=prpf.Proposed solution procedure in this study is based on the frontal method; however, the order of elimination is somehow different. The solution procedure is illustrated in Figure2. In assembly procedure, the elements are considered each in turn according to prescribed order, which caused the front width to be reduced. Whenever a new element is introduced, coefficient of dynamic stiffness and dynamic load of the element are summed either into existing equations, if the nodes are already active, or into new equations which have to be included in the front if the nodes are being active for the first time. If some reservoir pressure DOFs pr and solid DOFs {u}, not the interface pressure DOFs pf, appear for the last time, the corresponding equations can be eliminated. In other words, the pressure variables on interface pf are kept in front. Gaussian elimination is used for elimination of the variable {xs} (sth equation in front). It must be noted that variable {xs} is pr or {u}, not pf, at this step. The equation coefficients and right hand side term are extracted corresponding to {xs} by using Gaussian elimination as follows:(16a)K¯ij∗=K¯ij−K¯isK¯sjK¯ss,(16b)Fi∗=Fi−K¯isFsK¯ss,where 1 < i  < front width and i < j < front width.Figure 2
Frontal solution procedure for dam-reservoir problem.When assembly for all elements is finished, and all the variables are eliminated except the interface pressure variable. All the interface pressure DOFs remain in the front as shown in Figure1(d). Afterwards, the interface pressure DOFs can be eliminated by using Gaussian elimination and the coefficients of the interface pressure equation stored away.It was shown that in frontal solution, the pressure variables on interface can be kept in front and can be solved without continuing the solution for all equations. Then, the achieved results can be used as the boundary condition for solving the reservoir and as loading profile for solving the dam. In this way, no iteration is needed, there is no loss of accuracy, and each media can be solved by using its suitable formulation. Obviously, all deficiencies are removed without losing the accuracy.


## 5. Numerical Results

A special computer code was generated to determine the dynamic response of structure-fluid systems by using frontal method in basis of abovementioned algorithm in two dimensional cases. The generated code uses the eight node plain strain elements for discretization of the structure domain and four node elements for discretization of the fluid domain. In this study, fluid is considered as compressible, inviscid, and irrotational, and concrete is assumed to be homogenous and isotropic.As a first example, the dam reservoir system analyzed by Lee and Tsai [15] is considered. Dam with vertical upstream with 180 m height and reservoir has been subjected to ramp acceleration as shown in Figure 3. In this analysis, the width of the dam was w=32.32m and reservoir length was 5H = 900 m (H is structure height) as shown in Figure 4. Parameters assumed to describe the structural system in [15] are (case A) EI = 9.8437 × 109 ton·m2 and weight per unit length of the structure = 36 ton/m. Lee and Tsai [15] have studied this example analytically. For analysis of dam as plain strain, we need to convert these quantities to unit weight (ρs) and modulus of elasticity (E) as follows:(17)I=tw312=1∗32.32312=2813.4m3⇒E=rigiditisI=35GPa,ρs=36000w∗t=3600032.32∗1=1113.86kg/m3,where I is area moment of inertia related to the z axis. z axis is the third axis in a Cartesian coordinate system. For 2D modeling, element thickness t = 1 was assumed.Figure 3
Ramp acceleration.Figure 4
Geometrical detail of structure-reservoir system.The material of the dam was assumed to be linearly elastic. For the concrete of the dam, Poisson’s ratio was taken as 0.2. For the water, unit weight and sound speed in the water were taken as 1000 kg/m3 and 1440 m/s, respectively.In order to evaluate the accuracy and capability of the proposed solution, results are presented and compared with those presented by Lee and Tsai [15]. Crest displacement of the dam without and with interaction is presented in Figure 5 and Figure 6, respectively. Hydrodynamic pressure at bottom of the reservoir is presented in Figure 7. As shown in the figures, good agreement is found between the responses.Figure 5
Crest displacement without interaction.Figure 6
Crest displacement with interaction.Figure 7
Hydrodynamic pressure at bottom of the reservoir.For comparison purpose, the model was analyzed by using two different programs: first program solved the equations by frontal method and second program solved the equations by iterative scheme. The results are presented in Figure8 and Figure 9. As shown in the figures, excellent agreement between the results is achieved from proposed method and iterative solution.Figure 8
Crest displacement with interaction.Figure 9
Hydrodynamic pressure at bottom of the reservoir.For second example, Pine Flat concrete gravity dam with 121.92 m height and reservoir with 116.12 m depth and 580 m length (Figure10) has been subjected to Taft Lincoln acceleration (Figure 11). The ground motion has a peak acceleration of 0.18 g. For the concrete, unit weight, Poisson’s ratio, and modulus of elasticity were taken as 2500 kg/m3, 0.2, and 22.75 GPa, respectively. For the water, unit weight and sound speed in the water were taken as 1000 kg/m3 and 1440 m/s, respectively.Figure 10
Geometrical detail of dam-reservoir system.Figure 11
The Taft Lincoln earthquake record (1952).Crest displacement obtained from present study and response obtained by Fenves and Vargas-Loli [20] previously (case: without cavitation) are illustrated in Figure 12. It should be considered that ground motion is scaled to a peak acceleration of 1.0 g in [20].Figure 12
Pine Flat crest displacement, peak acceleration scaled to 1.0 g.As shown in Figure12, good agreement is found between the results.In Figure13 and Figure 14, the response of the Pine Flat is investigated by using two different schemes, first frontal method and second iterative method.Figure 13
Pine Flat crest displacement, peak acceleration scaled to 1.0 g.Figure 14
Hydrodynamic pressure at Pine Flat heel, peak acceleration scaled to 1.0 g.As shown in the figures, excellent agreement between the results is achieved from proposed method and iterative solution.Table1 presents the calculation time of two solution scheme for two examples. As shown, frontal solution is found 1.37 and 1.65 times less time-consuming than iterative method, in example one and two, respectively.Table 1
Execution time for two solution scheme.
Frontal method
Iterative method
Example 1
20
27.36
Example 2
28.2
46.53


## 6. Conclusions

Frontal solution is proposed for solving the dynamic analysis of concrete gravity dams. Coupled equations of dam-reservoir interaction in time domain are obtained and solved by frontal method. A Fortran program is developed and its accuracy is verified by comparing results with published results. Two examples of dam-reservoir model were analyzed. The results of frontal solution were compared with the results of iterative method. The major advantages of the frontal solution are as follows:(1)
The accuracy of the solution is the same as solution of the total equations in each time step.(2)
Separate solution of the structure and reservoir causes the suitable formulation for each one to be used, namely, Lagrangian and Eulerian, respectively.(3)
Solving very different order variables is eliminated.(4)
The presented scheme takes less system memory and calculation time than iterative method. For two models examined in this paper, execution time was reduced to 37 and 65 percent for example 1 and example 2, respectively.(5)
The presented scheme does not have the geometrical and boundary condition limitation.(6)
The presented scheme is capable of considering the nonlinear phenomenon.



---

*Source: 4369240-2018-10-08.xml*