Numerical Uncoupling of Domains in Dam-Reservoir Problem

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. )e 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 fluidmay 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.


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.
en obtained pressure will be applied on the dam, and new configuration of dam is calculated.e 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].erefore, only symmetric equations for each subsystem is calculated [3][4][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.is approach needs iteration in each time step even in linear problems, so this approach is involved with number of solutions.e 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-reservoirfoundation interaction in frequency domain covering several aspects of the problem [10][11][12][13].Frequency domain approach is computationally very efficient.However, concepts of the frequency domain formulation are more difficult than those of time domain formulation.e 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].ese analytical solutions are really efficient in time but they have geometrical and boundary condition limitation and also they cannot be used in nonlinear problems.
e 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.e principal aim of this study is to employ frontal solution for solving the damreservoir system.
e 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.en 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.
e 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) is scheme is capable to model the nonlinear problems if time domain formulation is adopted

Governing Equation of the Coupled Dam-Reservoir Problem
e dam-reservoir system is classified as a coupled problem where two fields each governed by a second-order differential equation, interact at their interface.e equation of motion for dam structure due to earthquake motion is given as where In above integration on the interface of the two domain, vector n { } is normal vector on interface and defined to be positive, going from the solid domain into the fluid domain.
[N s ] and [N f ] are shape function of structure and fluid, respectively.
[M] and [K] are defined as below 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 so [C] is taken as where α s and β s are constants that can be calculated from the below relation where ζ i is the damping ratio and ω i is ith natural frequency of the structure.e equation of motion for reservoir domain due to earthquake motion is given as [4] [ where where β b is the acoustic impedance ratio of the reservoir bottom material, c w is velocity of sound in the water, h is depth of the reservoir, and Ω f denotes the fluid domain.s f , s b , and s t 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: 2 Shock and Vibration M 0 Equation ( 8) can be rewritten for dynamic loads as follows: M 0 By using step by step integration scheme of Newmark [19], the displacement and hydrodynamic pressure and their first derivative can be discretized as follows: 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 4 where Multiplying the reservoir motion equation by −βΔt 2 /ρ f renders the asymmetric coupled matrix of Equation ( 12) to a symmetric one as follows: 4 where

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]. is 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 pressure p   is substituted into Equation (1)  is approach reduces coupled problem to subsystems.erefore, 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 Shock and Vibration convergence of the solution of the subsystems at each time step.
is approach considers the effect of interaction, but it cannot be generalized as an exact solution for interaction problem.e 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 phenomenon Although 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.

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.erefore, the element matrices, such as stiffness, are computed again any time needed.e 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].erefore 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.e number of the variables in the front is named front width.e 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.
e 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.
e procedure of the assembly for structure and reservoir by frontal method is illustrated in Figure 1. 4 Shock and Vibration e solution basis for frontal method is Gaussian elimination.e interaction equations of dam and reservoir contain solid DOFs ({u}) and liquid DOFs ({p}).Equation ( 13) can be rewritten as follows: where X { } � u p  , F { } � F s F f  , and [K] is effective dynamic matrix as defined in Equation (13).
{p} can be divided to reservoir pressure DOFs p r   and the interface pressure DOFs p f  : Proposed solution procedure in this study is based on the frontal method; however, the order of elimination is somehow different.e solution procedure is illustrated in Figure 2. 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 p r   and solid DOFs {u}, not the interface pressure DOFs p f  , appear for the last time, the corresponding equations can be eliminated.In other words, the pressure variables on interface p f   are kept in front.Gaussian elimination is used for elimination of the variable {x s } (sth equation in front).It must be noted that variable {x s } is p r   or {u}, not p f  , at this step.e equation coefficients and right hand side term are extracted corresponding to {x s } by using Gaussian elimination as follows: Shock and Vibration 5 where 1 < i < front width and i < j < front width.
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 Figure 1(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.en, 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.

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.e 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.32 m 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]   Shock and Vibration are (case A) EI 9.8437 × 10 9 ton•m 2 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: I tw 3 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.
e 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/m 3 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 gures, good agreement is found between the responses.For comparison purpose, the model was analyzed by using two di erent programs: rst program solved the equations by frontal method and second program solved the equations by iterative scheme.e results are presented in Figure 8 and Figure 9.As shown in the gures, excellent

Shock and Vibration
agreement between the results is achieved from proposed method and iterative solution.
For second example, Pine Flat concrete gravity dam with 121.92 m height and reservoir with 116.12 m depth and 580 m length (Figure 10) has been subjected to Taft Lincoln acceleration (Figure 11).e 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/m 3 , 0.2, and 22.75 GPa, respectively.For the water, unit weight and sound speed in the water were taken as 1000 kg/m 3 and 1440 m/s, respectively.
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].
As shown in Figure 12, good agreement is found between the results.
In Figure 13 and Figure 14, the response of the Pine Flat is investigated by using two di erent schemes, rst frontal method and second iterative method.
As shown in the gures, excellent agreement between the results is achieved from proposed method and iterative solution.
Table 1 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.

Conclusions
Frontal solution is proposed for solving the dynamic analysis of concrete gravity dams.Coupled equations of damreservoir interaction in time domain are obtained and solved by frontal method.A Fortran program is developed and its accuracy is veri ed by comparing results with published results.Two examples of dam-reservoir model were analyzed.e results of frontal solution were compared with the results of iterative method.e major advantages of the frontal solution are as follows: (1) e 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 di erent order variables is eliminated.(4) e 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) e presented scheme does not have the geometrical and boundary condition limitation.(6) e presented scheme is capable of considering the nonlinear phenomenon.
[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:

Figure 1 :
Figure 1: Frontal assembly procedure for two coupled domains.(a) e schematic model including fluid elements, structure elements, and nodes before beginning the assembly.(b) e first element of each media at the beginning of assembly.(c) e second element of each media for assembling.(d) Active nodes and deactivated nodes at the final step of assembly.

Figure 9 :
Figure 9: Hydrodynamic pressure at bottom of the reservoir.