Mathematical Problems in EngineeringVolume 2011 (2011), Article ID 210624, 23 pagesdoi:10.1155/2011/210624

Research Article

Finite Element Analysis of Dam-Reservoir Interaction Using High-Order Doubly Asymptotic Open Boundary

Yichao Gao,¹ Feng Jin,¹ Xiang Wang,² and Jinting Wang¹

¹State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China²Changjiang Institute of Survey, Planning, Design, and Research, Wuhan 430010, China

Received 30 May 2011; Accepted 23 August 2011

Academic Editor: Otto von Estorff

Copyright © 2011 Yichao Gao 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.

Abstract.

The dam-reservoir system is divided into the near field modeled by the finite element method, and the far field modeled by the excellent high-order doubly asymptotic open boundary (DAOB). Direct and partitioned coupled methods are developed for the analysis of dam-reservoir system. In the direct coupled method, a symmetric monolithic governing equation is formulated by incorporating the DAOB with the finite element equation and solved using the standard time-integration methods. In contrast, the near-field finite element equation and the far-field DAOB condition are separately solved in the partitioned coupled methodm, and coupling is achieved by applying the interaction force on the truncated boundary. To improve its numerical stability and accuracy, an iteration strategy is employed to obtain the solution of each step. Both coupled methods are implemented on the open-source finite element code OpenSees. Numerical examples are employed to demonstrate the performance of these two proposed methods.

1. Introduction

The coupled analysis of dam-reservoir interaction has great significance for the design and safety evaluation of concrete dams under earthquakes. The finite element method and substructure method are often applied for the analysis of dam-reservoir system (Figure 1). The dam structure as well as the near-field reservoir with irregular geometry is discretized with finite elements. The remaining part of the reservoir, called the far field, is simplified as a semi-infinite layer with constant depth. On the truncated boundary, which separates the near and far field, the equations of motion and continuity should be satisfied simultaneously. In the early studies, frequency-domain analysis methods [1, 2] are often used for linear problems. However, for a nonlinear analysis of the dam-reservoir system, it is necessary to develop a direct time domain analysis. Zienkiewicz and Bettess [3] as well as Küçükarslan et al. [4] studied fluid-structure interaction in the time domain by imposing the Sommerfeld radiation condition [5]. Tsai and his coworkers [6–8] established the time-domain models for the dynamic interaction analysis of dam-reservoir system by using the transmitting boundary. This approach is temporally global; that is, it requires the evaluation of convolution integrals. Boundary element method is undoubtedly a powerful numerical tool for analysis of problems involving unbounded domain. When the boundary element method [9–13] is applied to the direct time-domain analysis of dam-reservoir interaction, the formulation is temporally and spatially global. As a result, great numerical effort is required to solve the transient problems. This hinders its application to long-time analysis of large-scale engineering problems. The scaled boundary finite element method is a semianalytical method which is efficient for the analysis of dam-reservoir interaction. Details of this approach are given, for example, by Li et al. [14] and Lin et al. [15].

Figure 1: A typical dam-reservoir system.

The high-order open boundaries are promising alternatives for the simulation of semi-infinite reservoir in the analysis of dam-reservoir system. The high-order open boundaries [16, 17] are of increasing accuracy as the order of approximation increases. Moreover, the formulations are temporally local so that they are computationally efficient. As demonstrated by Prempramote et al. [18], these high-order open boundaries are singly asymptotic at the high-frequency limit and are appropriate for the radiative fields, where virtually all of the field energy is propagating out to infinity [19]. However, in some classes of application, such as a semi-infinite reservoir with constant depth (also known as a wave guide), a cutoff frequency exists. When the excitation frequency is close to or below the cutoff frequency, the wave field is largely nonradiative. In such cases, the high-order transmitting boundaries break down at late times in a time domain analysis [18]. To model an unbounded domain with the presence of nonradiative wave fields, one advance is the introduction of the doubly asymptotic boundaries [19–22]. Thus, the dynamic stiffness is exact at both the high-frequency and the low-frequency limit (i.e., statics), with its formulation spatially global. However, the highest order denoting the accuracy reported in the literature is three [23].

Recently, a novel high-order doubly asymptotic open boundary for one-dimensional scalar wave equation is proposed by Prempramote et al. [18] by extending the work in [24]. This high-order doubly asymptotic open boundary is capable of accurately mimicking the unbounded domain over the entire frequency (i.e., from zero to infinity). This open boundary condition is constructed by using the continued fraction solution of dynamic stiffness matrix without explicitly evaluating its solution at discrete frequencies. When applied for a semi-infinite layer with a constant depth, the constants of the continued fraction solutions with any order are determined explicitly and recursively. Excellent accuracy and stability for long-time transient analysis are reported.

Wang et al. [25] extend the high-order doubly asymptotic open boundary condition by Prempramote at al. [18] to the analysis of the hydrodynamic pressure of a semi-infinite reservoir with constant depth. By applying the sequential staggered implicit-implicit partition algorithm, the high-order doubly asymptotic open boundary is coupled with the general purpose finite element software ABAQUS to analyze the gravity dam-reservoir interaction system. Numerical examples demonstrate the high accuracy and long-time stability of this proposed technique.

Givoli et al. [26] proposed a new finite element scheme for the solution of time-dependent semi-infinite wave-guide problems by incorporating with a high-order open boundary. Two versions of finite element formulation, namely, the augmented version and the split version, are proposed. Good performance of this scheme is demonstrated in numerical examples, but the global mass and stiffness matrices of the augmented formulation are nonsymmetric.

In this paper, two coupled numerical methods for dam-reservoir interaction analysis based on the excellent high-order doubly asymptotic open boundary [25] are proposed. In the direct coupled method, the high-order doubly asymptotic open boundary is directly incorporated with the near-field finite element equation. A monolithic governing equation for the whole dam-reservoir system is formulated with sparse and symmetric coefficients matrices, which can be solved using the standard time-integration methods. In the partitioned coupled method, the near-field finite element equation and the far-field high-order doubly asymptotic boundary condition are separately solved. The high-order doubly asymptotic open boundary is split into two parts. The first part is proved to be the Sommerfeld radiation boundary, which can be included in the damping matrix of the near-field finite element equation. The second part includes all the high-order terms and is governed by a system of first-order ordinary differential equations. These two sets of equations are solved by a sequential staggered implicit-implicit partition algorithm. To improve the stability and accuracy of this partitioned coupled method, an iteration strategy is employed to obtain the solution of each step. Both of these two coupled methods are numerically implemented on the open-source finite element code OpenSees to analyze the gravity dam-reservoir and arch dam-reservoir interaction.

This paper is organized as follows. In Section 2, the finite element formulation of dam-reservoir system is addressed. In Section 3, the scaled boundary finite element method is applied to 3-dimensional semi-infinite layer with constant depth. The governing equation on the truncated boundary is obtained. In Section 4, the scaled boundary finite element equation is decoupled by modal decomposition. Based on the continued fraction solutions of the dynamic stiffness, a high-order doubly asymptotic open boundary is constructed by introducing auxiliary variables. In Section 5, two coupled numerical methods are presented: the direct coupled method and the partitioned coupled method. Both numerical methods are implemented on the open source finite element code OpenSees. In Section 6, numerical examples of a gravity dam and an arch dam are presented. In the final section, conclusions are stated.

2. Modeling of Dam-Reservoir System

A typical dam-reservoir system is shown in Figure 1. The reservoir is divided into two parts: the near field with irregular geometry and the far field extending to the infinity with constant depth. The dam structure and near field reservoir are dicretized with finite elements, and the far field reservoir is modeled with high-order doubly asymptotic boundary.

The hydrodynamic pressure is denoted as , and the acceleration of water particle is denoted as . Assuming the water in the reservoir to be compressible, inviscid, and irrotational with small amplitude movements, the hydrodynamic pressure in the reservoir satisfies the scalar wave equation where is the Laplace operator and is the compression wave velocity where is the bulk modulus of water and is the mass density. On the dam-reservoir interface (S1), satisfies the following boundary condition where stands for the outward normal of the boundary. At the reservoir bottom (S2), the following boundary condition applies. Neglecting the effect of surface waves, on the free surface (S3), applies. At infinity (S4), A Sommerfeld-type radiation boundary condition should be satisfied; namely,

Without considering the material damping, the finite element formulation for dam-reservoir system can be partitioned as where , , and are the mass matrix, static stiffness matrix, the coupling matrix between structure and acoustic fluid, respectively, and is the external force vector. Subscript denotes the degrees of freedom of the dam structure; subscript denotes the degrees of freedom on the truncated boundary; subscript denotes the degrees of freedom of the near-field water except for those on the truncated boundary. The interaction load applied to the semi-infinite reservoir (far field) by the near-field reservoir is denoted as , so the external load applied to the near-field reservoir on the truncated boundary is equal to . The mass and stiffness matrices of water treated as acoustic fluid are expressed as where denotes the shape function of finite elements.

To solve (2.7) for the dam-reservoir system, the relationship between the interaction load and the hydrodynamic pressure of the semi-infinite reservoir is formulated in the following section.

3. Summary of the Scaled Boundary Finite Element Method for Semi-Infinite Reservoir with Constant Depth

The scaled boundary finite element method is a semianalytical method developed to model unbounded domains with arbitrary shape [27]. The scaled boundary finite element formulation for the two-dimensional semi-infinite reservoir with constant depth was described in detailed [25]. For the sake of completeness, a brief summary of the equations necessary for the interpretation of high-order doubly asymptotic open boundary for hydrodynamic pressure is presented in this section.

To facilitate the coupling with the finite elements of the near-field reservoir, the truncated boundary is discretized by elements that have the same nodes and shape function as the finite elements. The derivation is summarized for three-dimensional semi-infinite reservoir with a vertical boundary (Figure 2). Streamlined expressions are presented as follows.

(a)

(b)

Figure 2: Semidiscretization of the truncated boundary.

For an acoustic fluid, the relationship between acceleration and hydrodynamic pressure is equivalent to that between stress and displacement in stress analysis and is expressed as where is the differential operator. The equation of continuity without considering the volumetric stress-strain relationship of compressible water is written as

The vertical truncated boundary of the semi-infinite reservoir is specified by a constant coordinate . It is discretized by two-dimensional elements (Figure 2(a)). A typical element is shown in Figure 2(b). The geometry of an element is interpolated using the shape functions formulated in the local coordinates and as

The Cartesian coordinates of a point and inside the semi-infinite reservoir are expressed as where the coordinate is equal to 0 on the vertical boundary. The Jacobian matrix for coordinate transformation from to is expressed as

For a three-dimensional problem, where is the determinant of the Jacobian matrix. The partial differential operator defined in (3.1) is expressed as with

Along horizontal lines passing through a node on the boundary, the nodal hydrodynamic pressure functions are introduced. On the boundary, the nodal hydrodynamic pressure follows as . One isoparametric element is shown in Figure 2. The hydrodynamic pressure field is obtained by interpolating the nodal hydrodynamic pressure functions

Substituting (3.9) and (3.7) into (3.1), the fluid particle acceleration is expressed as with

Applying Galerkin’s weighted residual technique in the circumferential directions to (3.2), the scaled boundary finite element equation for the three-dimensional semi-infinite reservoir with constant depth is obtained as where , , and are coefficient matrices

The coefficient matrices of (3.12) are evaluated by standard numerical integration techniques in the finite element method. Similar to the static stiffness and mass matrices in the finite element method, both of the coefficients and are sparse and positive definite. This formulation is the same as that of two-dimensional case [25]. It is applicable for 3-dimensional case with vertical truncated boundary of arbitrary geometry, such as the arch dam-reservoir system.

The acoustic nodal load vector on a vertical boundary with a constant is obtained based on virtual work principle and is expressed as

Assuming the time-harmonic behavior and ( is the excitation frequency) with the amplitudes of the hydrodynamic pressure , (3.12) is transformed into the frequency domain as and the amplitudes of the acoustic nodal load are expressed as

4. High-Order Doubly Asymptotic Open Boundary for Hydrodynamic Pressure

The derivation of high-order doubly asymptotic open boundary for hydrodynamic pressure is implemented based on the modal expression of (3.15). The streamlined expressions in [25] are summarized in this section.

4.1. Modal Decomposition of Scaled Boundary Finite Element Equation

Following the procedure in detailed [25], the system of ordinary differential equations in (3.15) can be decoupled by a modal transformation. The modes are obtained from the following generalized eigenvalue problem ( stands for a diagonal matrix): where is the diagonal matrix of positive eigenvalues, is a characteristic length (e.g., the depth of the semi-infinite layer) to nondimensionlize the eigenvalues, and are the matrix of eigenvectors representing the modes, which are normalized as

As a result, the inverse of the eigenvector matrix can be obtained by the matrix multiplication

Premultiplying (4.1) with results in

The relationship between amplitudes of the hydrodynamic pressure and amplitudes of the modal hydrodynamic pressure is defined as

Substituting (4.5) into (3.15) premultiplied with and using (4.2) and (4.3) lead to a system of decoupled equations with the dimensionless frequency where index indicates the modal number. Substituting (4.5) into (3.16), the acoustic nodal force vector is expressed as

The amplitude of the modal nodal force vector is defined as

Premultiplying (4.8) and using (4.2) and (4.9) yield

This equation transforms the amplitude of the acoustic nodal force vector to the amplitude of the modal force vector. The modal dynamic stiffness coefficient is defined as

By eliminating and from (4.6), (4.9), and (4.11), an equation for the modal dynamic stiffness coefficient is derived as

4.2. Doubly Asymptotic Continued Fraction Solution for Modal Dynamic Stiffness

Based on a doubly asymptotic solution of the modal dynamic stiffness coefficient , a temporally local open boundary [18] is constructed for a single mode of wave propagation. The solution of (4.12) is expressed as a doubly asymptotic continued fraction. An order high-frequency continued fraction is constructed first as

It is demonstrated by Prempramote et al. [18] that the high-frequency continued fraction solution does not converge when the excitation frequency is below the cutoff frequency. To determine a valid solution over the whole frequency range, an order low-frequency continued fraction solution is constructed for the residual term . Denoting the residual term for mode as The continued fraction solution for at the low frequency limit is expressed as

The doubly asymptotic continued fraction solution is constructed by combining the high-frequency continued fraction solution in (4.13) with the low-frequency solution in (4.15) using (4.14).

4.3. High-Order Doubly Asymptotic Open Boundary

Following the procedure developed for the modal space [18], the acoustic force-pressure relationship in the time domain is formulated by using the continued fraction solution of the dynamic stiffness coefficient and introducing auxiliary variables. A system of first-order ordinary differential equations with symmetric coefficient matrices is obtained, which represents a temporally local open boundary condition.

Substituting the continued fraction solution of the modal dynamic stiffness coefficient (4.13)–(4.15) into (4.11) and applying the inverse Fourier transform, the time-domain high-order doubly asymptotic open boundary in the modal space is expressed as where and are the auxiliary variables defined in modal space.

For an order doubly asymptotic continued fraction solution, the residual term applies. Equations (4.16)–(4.20) constitute a system of first order ordinary differential equations relating the interaction load , hydrodynamic pressure and the modal auxiliary variables in the modal space. This system of ordinary differential equations is a temporally local high-order doubly asymptotic open boundary condition for the semi-infinite reservoir with constant depth, which is directly established on the nodes of a vertical boundary. This boundary condition can be coupled seamlessly with finite element method.

5. Coupled Numerical Methods for Dam-Reservoir Interaction Analysis

Based on the high-order doubly asymptotic open boundary, two coupled numerical methods will be presented in this section: the direct coupled method and the partitioned coupled method. Both coupled numerical methods are implemented on the open-source object-oriented finite element code OpenSees, that is, the open system for earthquake engineering simulation.

5.1. The Direct Coupled Method

In the direct coupled method, by incorporating the high-order doubly asymptotic open boundary with the near-field finite element equation, a monolithic governing equation for the near-field structure and auxiliary variables representing the far-field reservoir is formulated.

To derive a symmetric monolithic formulation, modal transform is applied to the auxiliary variables defined in modal space and as where and are the auxiliary variables defined in real space.

Using (4.5), (4.10), (5.1), and left-multiplying (4.16)–(4.20) by yield with the symmetric and positive definite matrix

The residual term in (5.6) is setting to zero. Substituting (5.2)–(5.6) into (2.7) leads to a global linear system of ordinary differential equations in the time-domain

Here, subscript denotes global. , , and are the global mass matrix, global damping matrix, and global stiffness matrix, which are expressed as with the block matrices

The global load vector and unknowns are expressed as (the semicolon “;” stands for the vertical concatenation of vectors)

This system of linear equations describes the complete dam-reservoir system taking account of the influence of the semi-infinite reservoir. Similar to the near-field finite element formulation (2.7), the coefficient matrices of this formulation are sparse and symmetric. This dynamic system can be solved using a standard time-integration method, such as the Newmark family schemes.

Equation (5.8) is of order , where is the number of degrees of freedom of dam structure, is the number of degrees of freedom of near-field reservoir except for those of the truncated boundary, and is the number of degrees freedom of truncated boundary. Compared with the near-field finite element equation (2.7), additional auxiliary degrees of freedom are introduced in (5.8). From the point view of numerical implementation, (5.8) can be implemented on existing finite element code without any special treatment and solved by existing finite element solver. The block coefficient matrices and in (5.10) are evaluated only once in the analysis. However, the formulation of the direct coupled method involves additional auxiliary variables, which requires more computational effort and memory to solve this dynamic system.

5.2. The Partitioned Coupled Method

In the partitioned coupled method, the near-field finite element equation and the far-field high-order doubly asymptotic boundary condition are separately solved. They are coupled by the interaction force on the truncated boundary. The deviation of the partitioned coupled method is detailed in [25]. Streamlined expressions are presented as follows.

Using (4.3) and (4.10), (4.16) left-multiplied by is rewritten as

Substituting (5.12) into (2.7), the finite element formulation of the dam-reservoir system considering the interaction between the near-field and the semi-infinite reservoir is expressed as

The high-order doubly asymptotic open boundary is split into two parts. The first part is the additional damping term on left-hand side of (5.13), which is equivalent to the Sommerfeld radiation boundary as demonstrated in [25]. The second part is the coupling term on the right-hand side of (5.13) representing the contribution of the high-order terms of the doubly asymptotic boundary. It can be regard as an external load acted on the truncated boundary. For efficiency consideration in the numerical implementation, the coupling term is evaluated in the modal space.

Assembling (4.17)–(4.20) multiplied by , a system of ordinary differential equation for the modal auxiliary variables is formulated as where the coefficient matrices are expressed as

The unknown vector consists of the modal auxiliary variables of mode (the semicolon “;” stands for the vertical concatenation of vectors)

The only nonzero entry on the right-hand side is the modal hydrodynamic pressure obtained from (4.5)

It is demonstrated by Wang et al. [25] that (5.14) is stable up to the order at least.

Equations (5.13) for the near field and (5.14) for the far field are coupled by the auxiliary variables defined in the modal space. Using the same time integration scheme, such as the trapezoidal rule, these two sets of equations are solved by a sequential staggered implicit-implicit partitioned procedure proposed in [28, 29]. The value of the modal auxiliary variables at time station is determined by the last-solution extrapolation predictor [28, 29]

The auxiliary variables are obtained by integrating (5.14) for prescribed hydrodynamic pressure . The algorithm proposed in [25] to solve (5.13) and (5.14) proceeds as follows.

Step 1. Initialize variables and in (5.13) and for each mode in (5.14).

Step 2. Extract from of each mode and assigned to in (5.18).

Step 3. Form the right-hand term of (5.13), compute and by using an implicit method.

Step 4. Calculate modal hydrodynamic pressure and form the right-hand term of (5.14).

Step 5. Compute for each mode by using an implicit method.

Step 6. Increment to and go to Step 2.

Since the predicted vector has been used rather than in solving (5.13) and (5.14), this algorithm may lead to a numerical instability or poor accuracy. To avoid these, the solution algorithm for the partitioned coupled method in this paper is modified. The solution process within one step given by Step 2 to Step 5 is repeated a number of times in an iterative manner [26]. In each additional cycle, in Step 2 is extracted by the last computed in Step 5. Numerical experiments demonstrate that stability and accuracy are improved by performing a few additional iterations.

The order of (5.14) is only , and little computational effort is required to solve (5.14). Consequently, additional memory required for the solution of partitioned method is less than that of the direct coupled method.

This partitioned coupled method and the sequential staggered implicit-implicit partition algorithm [25] are both implemented in the open-source finite element code OpenSees rather than ABAQUS; as a result, the computational efficiency is greatly improved without any time costing restart analysis in ABAQUS.

6. Numerical Examples

Two numerical examples are analyzed to evaluate the accuracy and efficiency of the two present coupled numerical methods. The first example is a gravity dam with an irregular near field reservoir. The open boundary is employed to represent the regular semi-infinite reservoir of a constant depth. The time integration for these two coupled numerical methods is performed by the trapezoidal integration.

It is demonstrated by Wang et al. [25] that an order high-order doubly asymptotic open boundary condition is of excellent accuracy. So, the order open boundary condition is chosen for these two numerical examples.

6.1. Gravity Dam

A typical gravity dam-reservoir system with an irregular near field is shown in Figure 3(a), which is the same as the flexible dam example in [25]. The dam body has a modulus of elasticity  GPa, Posisson’s ration , and mass density  kg/m³. The water in the reservoir has a pressure wave velocity  m/s and mass density  kg/m³.

(a)

(b)

Figure 3: A gravity dam-reservoir system with irregular near field: (a) Geometry; (b) Mesh.

The finite element mesh is shown in Figure 3(b). The system is divided into three parts: the dam body, the near-field reservoir, and the far-field semi-infinite reservoir with constant depth. The dam body is discretized with eight-node solid elements, the near-field reservoir with 156 eight-node acoustic fluid elements, and the dam-reservoir interface with 13 three-node interface elements. The far-field reservoir is modeled by 13 three-node quadratic line elements on the truncated boundary, which share the same nodes and are compatible with those of the near-field acoustic fluid elements. The total number of nodes of the whole model is 653.

To verify the results, an extended mesh covering a far-field reservoir region of 7200 m is analyzed. The size of extended mesh is sufficiently large to avoid the pollution of the dam response by the waves reflected on the truncated boundary for a time duration of .

The El Centro earthquake ground motion (see Figure 4) is imposed as the horizontal acceleration at the base of the dam. The time step is chosen as 0.002 s during which the pressure wave travels about one third of the distance between two adjacent nodes. The partitioned coupled method is performed without any additional iteration; that is, the solution process is the same as that in [25]. The responses of hydrodynamic pressure at heel of the gravity dam of the first 20 s are plotted in Figure 5. Excellent agreement between the solutions of both coupled numerical methods and the extended mesh solution is observed during the first 10 s.

Figure 4: Time history of El Centro earthquake.

Figure 5: Hydrodynamic pressure at heel of the gravity dam under El Centro ground motion with a time step of 0.002 s.

To demonstrate the improvement of the modified solution algorithm for partitioned coupled method, the time step is increased to 0.02 s. The responses of the first 20 s are plotted in Figure 6. As it is shown, the results of the solution algorithm proposed by Wang et al. [25] tend to be divergent and inaccurate. However, both the solutions of direct coupled method and partitioned coupled method agree with the solution of extended mesh very well. The number of additional iterations within each step for partitioned coupled method recorded is usually one and no more than four in this example.

Figure 6: Hydrodynamic pressure at heel of the gravity dam under El Centro ground motion with a time step of 0.02 s.

The computer time for the gravity dam example list in Table 1 is recorded on a PC with a 2.93 GHz Intel Core i7 CPU. High computational efficiency of both coupled methods is observed. The computer time of the direct coupled method is about one tenth of that of the extended mesh. The computational effort of the partitioned coupled method is directly associated with the number of additional iterations. When there is no additional iteration, the computer time of the partitioned coupled method is nearly equal to that of the direct coupled method.

Table 1: Computer time for the gravity dam example.

Time step (sec)
Extended mesh
Direct coupled method
Partitioned coupled method
Algorithm in [25]

0.002
2980 s
242 s
241 s

0.02
219 s
25 s
46 s
25 s

6.2. Arch Dam

An arch dam-reservoir system is shown in Figure 7. The physical properties of dam body and water are the same as that in the example of gravity dam. The arch dam is of a height of 22 m and the near-field reservoir covers a region of the dam height. The dam body is discretized with 272 twenty-node hexahedron solid elements, the near-field reservoir with 1088 twenty-node hexahedron acoustic fluid elements, and the dam-reservoir interface with 136 eight-node interface elements. The far-field reservoir is modeled by 136 eight-node quadratic elements on the truncated boundary, which share the same nodes and are compatible with those of the near-field acoustic fluid elements. The total number of nodes of the whole model is 6652. To verify the results, a similar extended mesh covering the region of 7200 m is analyzed.

Figure 7: An arch dam-reservoir system.

Similar to the gravity dam example, the El Centro earthquake ground motion is imposed as the horizontal ( direction) acceleration at the base of the arch dam. The time step is chosen as 0.02 s. The responses of the hydrodynamic pressure at the heel of the arch dam are plotted in Figure 8. The solution of the direct coupled method agrees with solution of the extended mesh and is long time stable. However, as for the partitioned coupled method, numerical divergence is observed at the early time of the analysis. Numerical instability of such coupling strategy is also reported in [12]. It can be expected as the partitioned coupled method is conditional stable; that is, the integration time step is limited by stability limits. When the time step is greater than the stability limit, numerical instability may occur, for example, the results of the arch dam example. As the dam-reservoir system is quite complicated, it is difficult to determine the stability limits of different application cases. When the partitioned coupled method is applied, smaller time step should be used.

Figure 8: Hydrodynamic pressure at heel of the arch dam under El Centro ground motion with a time step of 0.02 s.

The computer time of the direct coupled method recorded is 1463 s, which is about one ninth of that of the extended mesh, that is, 12547 s. It also demonstrates the high computational efficiency of the direct coupled method.

7. Conclusions

Two coupled numerical methods were developed for the dam-reservoir interaction analysis by incorporating the finite element method with the excellent high-order doubly asymptotic open boundary. The dam-reservoir system is divided into the near-field with irregular geometry and the far-field by the truncated boundary. In the direct coupled method, a global monolithic equation for the whole dam-reservoir system is formulated with sparse and symmetric coefficient matrices, which can be solved by the standard finite element solver. In the partitioned coupled method, the near-field finite element equation and the high-order open boundary condition are separately solved. They are coupled by the interaction force applied on the truncated boundary. The partitioned coupled method is achieved by using a sequential staggered implicit-implicit procedure. To improve the numerical stability and accuracy of the algorithm, an iteration strategy was employed to obtain the solution of each step.

Both of the two coupled numerical methods are implemented on the open-source finite element code OpenSees. Numerical experiments demonstrated the high efficiency and accuracy of both coupled numerical methods. The memory required for the solution of the partitioned method is less than that of the direct coupled method. Although the numerical stability and accuracy of the partitioned coupled method can be improved by additional iterations within each step, the partitioned coupled method is conditionally stable yet. Its stability is also related to the predictor. Further research should be carried out to improve the stability of the partitioned coupled method. In contrast, the direct coupled method is unconditionally stable if only an unconditionally stable time integration algorithm such as trapezoidal integration is chosen. Consequently, larger time steps can be used in the direct coupled method than that in the partitioned coupled method.

Acknowledgments

This research was supported by the State Key Laboratory of Hydroscience and Engineering of Tsinghua University under Grant no. 2008-TC-2, the National Science Foundation of China under Grants nos. 90510018, and 90715041, and by the National Key Technology R&D Program under Grant no. 2008BAB29B05. The authors are also grateful to Professor Chongmin Song’s instructive suggestions on this research.

References

A. K. Chopra and P. Chakrabarti, “Earthquake analysis of concrete gravity dams including dam water foundation rock interaction,” Earthquake Engineering and Structural Dynamics, vol. 9, no. 4, pp. 363–383, 1981.

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.

O. C. Zienkiewicz and P. Bettess, “Fluid-structure dynamic interaction and wave forces. An introduction to numerical treatment,” International Journal for Numerical Methods in Engineering, vol. 13, no. 1, pp. 1–16, 1978.

S. Küçükarslan, S. B. Coşkun, and B. Taşkin, “Transient analysis of dam-reservoir interaction including the reservoir bottom effects,” Journal of Fluids and Structures, vol. 20, no. 8, pp. 1073–1084, 2005.

A. Sommerfeld, Partial Differential Equations in Physics, Academic Press, New York, NY, USA, 1949.

C. S. Tsai, G. C. Lee, and R. L. Ketter, “Semi-analytical method for time-domain analyses of dam-reservoir interactions,” International Journal for Numerical Methods in Engineering, vol. 29, no. 5, pp. 913–933, 1990.

C. S. Tsai and G. C. Lee, “Time-domain analyses of dam-reservoir system. II: substructure method,” Journal of Engineering Mechanics, vol. 117, no. 9, pp. 2007–2026, 1991.

R. Yang, C. S. Tsai, and G. C. Lee, “Explicit time-domain transmitting boundary for dam-reservoir interaction analysis,” International Journal for Numerical Methods in Engineering, vol. 36, no. 11, pp. 1789–1804, 1993.

T. Touhei and T. Ohmachi, “A FE-BE method for dynamic analysis of dam-foundation-reservoir systems in the time domain,” Earthquake Engineering and Structural Dynamics, vol. 22, no. 3, pp. 195–209, 1993.

R. J. Câmara, “A method for coupled arch dam-foundation-reservoir seismic behaviour analysis,” Earthquake Engineering and Structural Dynamics, vol. 29, no. 4, pp. 441–460, 2000.

O. Czygan and O. von Estorff, “Fluid-structure interaction by coupling BEM and nonlinear FEM,” Engineering Analysis with Boundary Elements, vol. 26, no. 9, pp. 773–779, 2002.

D. Soares Jr., O. von Estorff, and W. J. Mansur, “Efficient non-linear solid-fluid interaction analysis by an iterative BEM/FEM coupling,” International Journal for Numerical Methods in Engineering, vol. 64, no. 11, pp. 1416–1431, 2005.

A. Seghir, A. Tahakourt, and G. Bonnet, “Coupling FEM and symmetric BEM for dynamic interaction of dam-reservoir systems,” Engineering Analysis with Boundary Elements, vol. 33, no. 10, pp. 1201–1210, 2009.

S. M. Li, H. Liang, and A. M. Li, “A semi-analytical solution for characteristics of a dam-reservoir system with absorptive reservoir bottom,” Journal of Hydrodynamics, vol. 20, no. 6, pp. 727–734, 2008.

G. Lin, J. Du, and Z. Hu, “Dynamic dam-reservoir interaction analysis including effect of reservoir boundary absorption,” Science in China Series E: Technological Sciences, vol. 50, no. 1, pp. 1–10, 2007.

S. V. Tsynkov, “Numerical solution of problems on unbounded domains. A review,” Applied Numerical Mathematics, vol. 27, no. 4, pp. 465–532, 1998.

D. Givoli, “High-order local non-reflecting boundary conditions: a review,” Wave Motion, vol. 39, no. 4, pp. 319–326, 2004, New computational methods for wave propagatio

S. Prempramote, C. Song, F. Tin-Loi, and G. Lin, “High-order doubly asymptotic open boundaries for scalar wave equation,” International Journal for Numerical Methods in Engineering, vol. 79, no. 3, pp. 340–374, 2009.

T. Geers, “Singly and doubly asymptotic computational boundaries,” in Proceedings of the IUTAM Symposium on Computational Methods for Unbounded Domains, pp. 135–141, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1998.

T. L. Geers, “Doubly asymptotic approximations for transient motions of submerged structures,” Journal of the Acoustical Society of America, vol. 64, no. 5, pp. 1500–1508, 1978.

P. Underwood and T. L. Geers, “Doubly asymptotic, boundary-element analysis of dynamic soil-structure interaction,” International Journal of Solids and Structures, vol. 17, no. 7, pp. 687–697, 1981.

T. L. Geers and B. A. Lewis, “Doubly asymptotic approximations for transient elastodynamics,” International Journal of Solids and Structures, vol. 34, no. 11, pp. 1293–1305, 1997.

T. L. Geers and B. J. Toothaker, “Third-order doubly asymptotic approximations for computational acoustics,” Journal of Computational Acoustics, vol. 8, no. 1, pp. 101–120, 2000.

M. H. Bazyar and C. Song, “A continued-fraction-based high-order transmitting boundary for wave propagation in unbounded domains of arbitrary geometry,” International Journal for Numerical Methods in Engineering, vol. 74, no. 2, pp. 209–237, 2008.

X. Wang, F. Jin, S. Prempramote, and C. Song, “Time-domain analysis of gravity dam-reservoir interaction using high-order doubly asymptotic open boundary,” Computers and Structures, vol. 89, no. 7-8, pp. 668–680, 2011.

D. Givoli, B. Neta, and I. Patlashenko, “Finite element analysis of time-dependent semi-infinite wave-guides with high-order boundary treatment,” International Journal for Numerical Methods in Engineering, vol. 58, no. 13, pp. 1955–1983, 2003.

C. Song and J. P. Wolf, “The scaled boundary finite-element method—alias consistent infinitesimal finite-element cell method—for elastodynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 147, no. 3-4, pp. 329–355, 1997.

K.-C. Park, “Partitioned transient analysis procedures for coupled-field problems: stability analysis,” Journal of Applied Mechanics, vol. 47, no. 2, pp. 370–376, 1980.

K. C. Park and C. A. Felippa, “Partitioned transient analysis procedures for coupled-field problems: accuracy analysis,” Journal of Applied Mechanics, vol. 47, no. 4, pp. 919–926, 1980.