Research Article  Open Access
Numerical Uncoupling of Domains in DamReservoir Problem
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.
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 damreservoirfoundation 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 damreservoir 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 damreservoir systems, and then they used this simplified approach to determine the modal dynamic response of flexible beamfluid 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 damreservoir 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 DamReservoir Problem
The damreservoir system is classified as a coupled problem where two fields each governed by a secondorder differential equation, interact at their interface. The equation of motion for dam structure due to earthquake motion is given aswhere , , and are mass, damping, and stiffness matrices of the dam structure, respectively; is the static load vector; is the ground acceleration; , , and are vectors of displacement, velocity, and acceleration of structure body, respectively; is the vector of hydrodynamic pressure; and is coupling matrix defined as below
In above integration on the interface of the two domain, vector is normal vector on interface and defined to be positive, going from the solid domain into the fluid domain. and are shape function of structure and fluid, respectively.
and are defined as belowwhere denotes the structure domain, is mass density of the structure, is strain matrix, and is elastic constants matrix.
In this study, Rayleigh damping is assumed, and so is taken aswhere and are constants that can be calculated from the below relationwhere is the damping ratio and is ith natural frequency of the structure.
The equation of motion for reservoir domain due to earthquake motion is given as [4]where , , and are quasimass, damping, and stiffness matrices of the reservoir domain, respectively; is the density of the fluid; and and are vectors of first and second derivation of the hydrodynamic pressure. , , and are defined as follows:where is the acoustic impedance ratio of the reservoir bottom material, is velocity of sound in the water, h is depth of the reservoir, and denotes the fluid domain. , , and 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:
Equation (8) can be rewritten for dynamic loads as follows:
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 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) resultswhere
Multiplying the reservoir motion equation by renders the asymmetric coupled matrix of Equation (12) to a symmetric one as follows:where
3. Most Common Solutions for Coupled DamReservoir 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 pressure is substituted into Equation (1) to obtain , , and . It is better to select the calculated pressure from previous time step as arbitrary pressure in order to achieve rapid convergence.(2)Calculated is substituted in Equation (6) to obtain , , and .(3)Calculated is substituted in Equation (1) to obtain , , and .(4)Calculated is substituted in Equation (6) to obtain , , and .(5)Check for convergence: if and 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 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.
4. Frontal Solution for Coupled DamReservoir 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 Figure 1.
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:where , and is effective dynamic matrix as defined in Equation (13).
{p} can be divided to reservoir pressure DOFs and the interface pressure DOFs :
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 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 and solid DOFs {u}, not the interface pressure DOFs , appear for the last time, the corresponding equations can be eliminated. In other words, the pressure variables on interface 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 or {u}, not , at this step. The equation coefficients and right hand side term are extracted corresponding to {x_{s}} by using Gaussian elimination as follows:where 1 < < front width and < < 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. 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 structurefluid 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 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 × 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 () and modulus of elasticity (E) as follows:where 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.
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/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 figures, good agreement is found between the responses.
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 Figure 8 and Figure 9. As shown in the figures, excellent 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). 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/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 VargasLoli [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 different schemes, first frontal method and second iterative method.
As shown in the figures, 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 timeconsuming than iterative method, in example one and two, respectively.

6. 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 verified by comparing results with published results. Two examples of damreservoir 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
References
 C. A. Felippa and K. C. Park, “Staggered transient analysis procedures for coupled mechanical systems: formulation,” Computer Methods in Applied Mechanics and Engineering, vol. 24, no. 1, pp. 61–111, 1980. View at: Publisher Site  Google Scholar
 J. C. Pedro and P. Sibanda, “An algorithm for the strongcoupling of the fluidstructure interaction using a staggered approach,” ISRN Applied Mathmatics, vol. 2012, Article ID 391974, p. 14, 2012. View at: Publisher Site  Google Scholar
 M. Ghaemian and A. Ghobarah, “Staggered solution schemes for dam reservoir interaction,” Journal of Fluids and Structures, vol. 12, no. 7, pp. 933–948, 1998. View at: Publisher Site  Google Scholar
 F. Kalateh and R. Attarnejad, “Finite element simulation of acoustic cavitation in the reservoir and effects on dynamic response of concrete dams,” Finite Element in Analysis and Design, vol. 47, no. 5, pp. 543–558, 2011. View at: Publisher Site  Google Scholar
 M. A. HaririArdabili, S. M. SeyedKolbadi, and M. R. Kianoush, “FEM based parametric analysis of a typical gravity dam considering input excitation mechanism,” Soil Dynamics and Earthquake Engineering, vol. 84, pp. 22–43, 2016. View at: Publisher Site  Google Scholar
 E. I. Wilson and M. Khalvati, “Finite element for dynamic analysis of fluid solid system,” International Journal for Numerical Methods in Engineering, vol. 19, no. 11, pp. 1657–1668, 1983. View at: Publisher Site  Google Scholar
 M. Akkose, A. Bayraktar, and A. A. Dumanoglu, “Reservoir water level effects on nonlinear dynamic response of arch dams,” Journal of Fluids and Structures, vol. 24, no. 3, pp. 418–435, 2008. View at: Publisher Site  Google Scholar
 M. Akkose and E. Simsek, “Nonlinear seismic response of concrete gravity dams to near fault ground motions including damwatersedimentfoundation interaction,” Applied Mathematical Modelling, vol. 34, no. 11, pp. 3685–3700, 2010. View at: Publisher Site  Google Scholar
 F. Han, J. Yao, C. Wang, and H. Zhu, “Bow flare water entry impact prediction and simulation based on moving particle semiimplicit turbulence method,” Shock and Vibration, vol. 2018, Article ID 7890892, 16 pages, 2018. View at: Publisher Site  Google Scholar
 J. F. Hall and A. K. Chopra, “Two dimensional dynamic analysis of concrete gravity and embankment dams including hydrodynamic effects,” Earthquake Engineering and Structural Dynamics, vol. 10, no. 2, pp. 305–332, 1982. View at: Publisher Site  Google Scholar
 G. L. Fenves and A. K. Chopra, “Effects of reservoir bottom absorption on earthquake response of concrete gravity dams,” Earthquake Engineering and Structural Dynamics, vol. 11, no. 6, pp. 809–829, 1983. View at: Publisher Site  Google Scholar
 K. B. Fok and A. K. Chopra, “Frequency response function for arch arch dams: hydrodynamic and foundation flexibility effect,” Earthquake Engineering and Structural Dynamics, vol. 14, no. 5, pp. 769–795, 1986. View at: Publisher Site  Google Scholar
 K. B. Fok and A. K. Chopra, “Hydrodynamic and foundation flexibility effects in earthquake response of arch dams,” Journal of Engineering Mechanics, vol. 112, no. 8, pp. 1810–1828, 1987. View at: Publisher Site  Google Scholar
 C. S. Tsai, G. C. Lee, and R. L. Ketter, “A semi analytical method for time domain analyses of dam reservoir interaction,” International Journal for Numerical Methods in Engineering, vol. 29, no. 5, pp. 913–933, 1990. View at: Publisher Site  Google Scholar
 G. C. Lee and C. S. Tsai, “Time domain analysis of damreservoir system. I: exact solution,” Journal of Engineering Mechanics, vol. 117, no. 9, pp. 1990–2006, 1991. View at: Publisher Site  Google Scholar
 B. Miqual and N. Bouaanani, “Simplified evaluation of the vibration period and seismic response of gravity dam water systems,” Engineering Structures, vol. 32, no. 8, pp. 2488–2502, 2010. View at: Publisher Site  Google Scholar
 B. Miqual and N. Bouaanani, “Efficient modal dynamic analysis of flexible beam fluid systems,” Applied Mathematical Modeling, vol. 39, no. 1, pp. 99–116, 2014. View at: Publisher Site  Google Scholar
 B. M. Irons, “A frontal solution program for finite element analysis,” International Journal for Numerical Methods in Engineering, vol. 2, no. 1, pp. 5–32, 1970. View at: Publisher Site  Google Scholar
 O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: The Basic Fifth Edition, vol. 1, ButterworthHeinemann, Oxford, UK, 2000.
 G. Fenves and L. M. VargasLoli, “Nonlinear dynamic analysis of fluid structure systems,” Journal of Engineering Mechanics, vol. 114, no. 2, pp. 219–240, 1988. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Saba Golchin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.