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

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 1𝑐²𝜕²𝑝𝜕𝑡²=Δ𝑝,(2.1) where Δ is the Laplace operator and 𝑐 is the compression wave velocity𝑐=𝐾𝜌,(2.2) where 𝐾 is the bulk modulus of water and 𝜌 is the mass density. On the dam-reservoir interface (S1), 𝑝 satisfies the following boundary condition 𝜕𝑝𝜕𝑛=−𝜌𝑢̈^𝑛,(2.3) where 𝑛 stands for the outward normal of the boundary. At the reservoir bottom (S2), the following boundary condition 𝜕𝑝𝜕𝑛=0or𝑢̈^𝑛=0(2.4) applies. Neglecting the effect of surface waves, on the free surface (S3), 𝑝=0(2.5) applies. At infinity (S4), A Sommerfeld-type radiation boundary condition should be satisfied; namely, 𝜕𝑝𝜕𝑛=−𝑝𝑐̇.(2.6)

Without considering the material damping, the finite element formulation for dam-reservoir system can be partitioned as ⎡⎢⎢⎢⎢⎣−𝑀𝑄_𝑠00_𝑓𝑠0𝑀𝑀_{𝑓𝑓𝑏𝑓}𝑀_𝑓𝑏𝑀_𝑏𝑏⎤⎥⎥⎥⎥⎦⎧⎪⎪⎨⎪⎪⎩𝑢̈^𝑠𝑝̈_𝑓𝑝̈_𝑏⎫⎪⎪⎬⎪⎪⎭+⎡⎢⎢⎢⎢⎣𝐾^𝑠0𝐾𝑄_{𝑠𝑓𝑓𝑓}00𝐾_𝑏𝑓𝐾_𝑓𝑏𝐾_𝑏𝑏⎤⎥⎥⎥⎥⎦⎧⎪⎪⎨⎪⎪⎩𝑢^𝑠𝑝_𝑓𝑝_𝑏⎫⎪⎪⎬⎪⎪⎭=⎧⎪⎪⎨⎪⎪⎩𝑓^𝑠𝑓_𝑓−{𝑟}⎫⎪⎪⎬⎪⎪⎭,(2.7) 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 𝑀_𝑓=_𝑉𝑓1𝐾[𝑁]^𝑇𝐾[𝑁]𝑑𝑉,_𝑓=_𝑉𝑓1𝜌𝜕[𝜕𝑁𝑥]𝑓^𝑇𝜕𝜕[𝑁𝑥]+𝜕[𝜕𝑁𝑦]^𝑇𝜕𝜕[𝑁𝑦]+𝜕[𝜕𝑁𝑧]^𝑇𝜕𝜕[𝑁𝑧]𝑑𝑉,_𝑓=_𝑆𝑓1𝜌[𝑁]𝜕𝜕[𝑁𝑛]𝑑𝑆,(2.8) 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 {𝐿}𝑝+𝜌{𝑢}̈=0,(3.1) where {𝐿}={𝜕/𝜕𝑥𝜕/𝜕𝑦𝜕/𝜕𝑧} is the differential operator. The equation of continuity without considering the volumetric stress-strain relationship of compressible water is written as {𝐿}^𝑇{𝑢}̈=−1𝐾𝜕²𝑝𝜕𝑡².(3.2)

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 𝑦_𝑏(𝜂,𝜁)=[𝑁(𝜂,𝜁)]𝑦_𝑏,𝑧_𝑏(𝜂,𝜁)=[𝑁(𝜂,𝜁)]𝑧_𝑏.(3.3)

The Cartesian coordinates of a point (𝑥,𝑦,𝑧) and inside the semi-infinite reservoir are expressed as𝑦(𝜉,𝜂,𝜁)=𝑦𝑥_𝑏(𝜉)=𝑥_𝑏+𝜉,𝑧(𝜉,𝜂,𝜁)=𝑧_𝑏(𝜂,𝜁)=[𝑁(𝜂,𝜁)]𝑦_𝑏,(𝜂,𝜁)=[𝑁(𝜂,𝜁)]𝑧_𝑏,(3.4) where the coordinate 𝜉 is equal to 0 on the vertical boundary. The Jacobian matrix for coordinate transformation from (𝑥,𝑦,𝑧) to (𝜉,𝜂,𝜁) is expressed as [𝐽(𝜂,𝜁)]=⎡⎢⎢⎢⎢⎣𝑥^,𝜉𝑦_,𝜉𝑧_,𝜉𝑥_,𝜂𝑦_,𝜂𝑧_,𝜂𝑥_,𝜁𝑦_,𝜁𝑧_,𝜁⎤⎥⎥⎥⎥⎦=⎡⎢⎢⎢⎢⎣10𝑦0_𝑏,𝜂0𝑧_𝑏,𝜂0𝑦_𝑏,𝜁𝑧_𝑏,𝜁⎤⎥⎥⎥⎥⎦.(3.5)

For a three-dimensional problem, 𝑑𝑉=||𝐽||𝑑𝜉𝑑𝜂𝑑𝜁,(3.6) where |𝐽| is the determinant of the Jacobian matrix. The partial differential operator defined in (3.1) is expressed as {𝐿}=[𝐽(𝜂,𝜁)]⁻¹𝜕𝜕𝜉𝜕𝜕𝜂𝜕𝜕𝜁^𝑇=𝑏¹𝜕𝜕𝜉+𝑏²𝜕𝜕𝜂+𝑏³𝜕𝜕𝜁,(3.7) with𝑏¹=⎧⎪⎪⎨⎪⎪⎩100⎫⎪⎪⎬⎪⎪⎭,𝑏²=||𝐽1||⎧⎪⎪⎨⎪⎪⎩𝑧0^𝑏,𝜁−𝑦_𝑏,𝜁⎫⎪⎪⎬⎪⎪⎭,𝑏₃=||𝐽1||⎧⎪⎪⎨⎪⎪⎩−𝑧0^𝑏,𝜂𝑦_𝑏,𝜂⎫⎪⎪⎬⎪⎪⎭.(3.8)

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 {𝑝_𝑏}={𝑝(𝜉=0,𝑡)}. One isoparametric element 𝑆^𝑒 is shown in Figure 2. The hydrodynamic pressure field 𝑝=𝑝(𝜉,𝜂,𝜁,𝑡) is obtained by interpolating the nodal hydrodynamic pressure functions𝑝=[𝑁(𝜂,𝜁)]{𝑝}.(3.9)

Substituting (3.9) and (3.7) into (3.1), the fluid particle acceleration {𝑢}̈={𝑢(̈𝜉,𝜂,𝜁)} is expressed as {𝑢}̈=−1𝜌𝐵¹{𝑝}_,𝜉+𝐵²{𝑝},(3.10) with 𝐵¹=𝑏¹[𝑁(𝜂,𝜁)],𝐵²=𝑏²[𝑁(𝜂,𝜁)]_,𝜂+𝑏³[𝑁(𝜂,𝜁)]_,𝜁.(3.11)

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 𝐸⁰{𝑝}_,𝜉𝜉−𝐸²{𝑝}−1𝑐²𝐸⁰{𝑝}̈=0,(3.12) where [𝐸⁰], [𝐸²], and [𝑀⁰] are coefficient matrices𝐸⁰=_𝑆𝜉𝐵¹^𝑇1𝜌𝐵¹||𝐽||𝑑𝜂𝑑𝜁=_𝑆𝜉[𝑁(𝜂,𝜁)]^𝑇1𝐸𝜌[𝑁(𝜂,𝜁)]||𝐽||𝑑𝜂𝑑𝜁,₂=_𝑆𝜉𝐵²^𝑇1𝑀𝜌𝐵²||𝐽||𝑑𝜂𝑑𝜁,₀=_𝑆𝜉[𝑁(𝜂,𝜁)]^𝑇1𝐾[𝑁(𝜂,𝜁)]||𝐽||𝑑𝜂𝑑𝜁=𝑐1²𝐸⁰.(3.13)

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{𝑟}=𝐸⁰{𝑃}_,𝜉.(3.14)

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𝐸⁰{𝑃}_,𝜉𝜉−𝐸²{𝑃}+𝜔²𝑐²𝐸⁰{𝑃}=0,(3.15) and the amplitudes of the acoustic nodal load {𝑅}={𝑅(𝜉,𝜔)} are expressed as {𝑅}=−𝐸⁰{𝑃}_,𝜉.(3.16)

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):𝐸²[Φ]=𝐸⁰[Φ]𝜆^2𝑗ℎ²,(4.1) where ⟨𝜆^2𝑗⟩ 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[Φ]^𝑇𝐸⁰[Φ]=[𝐼].(4.2)

As a result, the inverse of the eigenvector matrix can be obtained by the matrix multiplication[Φ]⁻¹=[Φ]^𝑇𝐸⁰.(4.3)

Premultiplying (4.1) with [Φ]^𝑇 results in[Φ]^𝑇𝐸²[Φ]=𝜆^2𝑗ℎ².(4.4)

The relationship between amplitudes of the hydrodynamic pressure and amplitudes of the modal hydrodynamic pressure {𝑃}={𝑃(𝜉,𝜔)} is defined as {𝑃}=[Φ]𝑃.(4.5)

Substituting (4.5) into (3.15) premultiplied with [Φ]^𝑇 and using (4.2) and (4.3) lead to a system of decoupled equations 𝑃_{𝑗,𝜉𝜉}+1ℎ²𝑎²⁰−𝜆^2𝑗𝑃_𝑗=0,(4.6) with the dimensionless frequency𝑎₀=𝜔ℎ𝑐,(4.7) where index 𝑗 indicates the modal number. Substituting (4.5) into (3.16), the acoustic nodal force vector is expressed as {𝑅}=−𝐸⁰[Φ]𝑃_,𝜉.(4.8)

The amplitude of the modal nodal force vector {𝑅}={𝑅(𝜉,𝜔)} is defined as 𝑅=−ℎ𝑃_,𝜉or𝑅_𝑗=−ℎ𝑃_𝑗,𝜉.(4.9)

Premultiplying (4.8) and using (4.2) and (4.9) yield𝑅=ℎ[Φ]^𝑇{𝑅}.(4.10)

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 𝑅_𝑗=𝑆_𝑗𝑎₀𝑃_𝑗.(4.11)

By eliminating 𝑅_𝑗 and 𝑃_𝑗 from (4.6), (4.9), and (4.11), an equation for the modal dynamic stiffness coefficient is derived as 𝑆^𝑗𝑎₀²+𝑎²⁰−𝜆^2𝑗=0.(4.12)

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 𝑆_𝑗𝑎₀=𝑖𝑎₀−𝜆^2𝑗𝑌_1(𝑖)𝑌^𝑗(1)(𝑎₀)⁻¹,𝑎₀=(−1)^𝑖2𝑖𝑎₀𝑌^(𝑖)1,𝑗−𝜆^2𝑗𝑌^𝑗(𝑖+1)(𝑎₀)⁻¹𝑖=1,2,…,𝑀_𝐻.(4.13)

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 𝑌^{(𝑀𝐻𝑗+1)}(𝑎₀). Denoting the residual term for mode 𝑗 as 𝑌_𝐿,𝑗𝑎₀=𝑌^{(𝑀𝐻𝑗+1)}𝑎₀.(4.14) The continued fraction solution for 𝑌_𝐿,𝑗(𝑎₀) at the low frequency limit is expressed as 𝑌_𝐿,𝑗𝑎₀=(−1)^𝑀𝐻+1𝜆_𝑗+(−1)^𝑀𝐻+1𝑖𝑎₀−𝑖𝑎₀²𝑌_{(𝑖)𝐿,𝑗}𝑌^(1)𝐿,𝑗(𝑎₀)⁻¹,𝑎₀=(−1)^{𝑀𝐻+𝑖+1}2𝜆_𝑗−𝑖𝑎₀²𝑌^{(𝑖+1)𝐿,𝑗}(𝑎₀)⁻¹𝑖=1,2,…,𝑀_𝐿.(4.15)

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 1ℎ𝑟^𝑗̃=1𝑐𝑝̇^𝑗̃−1ℎ𝜆^𝑗0=−1𝑝̃^𝑗(1),(4.16)ℎ𝜆^𝑗𝑝_𝑗̃−2𝑐𝑝̇̃^𝑗(1)−1ℎ𝜆^𝑗0=−1𝑝̃^𝑗(2),(4.17)ℎ𝜆^𝑗𝑝̃^{𝑗(𝑖−1)}+(−1)^𝑖2𝑐𝑝̇̃^𝑗(𝑖)−1ℎ𝜆^𝑗0=−1𝑝̃^𝑗(𝑖+1)𝑖=2,3,…,𝑀_𝐻,(4.18)ℎ𝜆^𝑗𝑝^{(𝑀𝐻)𝑗}̃+(−1)^𝑀𝐻+11ℎ𝜆^𝑗𝑝̃^(0)𝐿,𝑗+(−1)^𝑀𝐻+11𝑐𝑝̇̃^(0)𝐿,𝑗−1𝑐𝑝̇̃^(1)𝐿,𝑗0=−1,(4.19)𝑐𝑝̇̃^{(𝑖−1)𝐿,𝑗}+(−1)^{𝑀𝐻+𝑖+1}2ℎ𝜆^𝑗𝑝̃^{(𝑖)𝐿,𝑗}−1𝑐𝑝̇̃^{(𝑖+1)𝐿,𝑗}𝑖=1,2,…,𝑀_𝐿,(4.20) where {𝑝̃^(𝑖)}(𝑖=1,2,…,𝑀_𝐻) and {𝑝̃^𝐿(𝑖)}(𝑖=0,1,2,…,𝑀_𝐿) are the auxiliary variables defined in modal space.

For an order 𝑀_𝐻=𝑀_𝐿 doubly asymptotic continued fraction solution, the residual term 𝑝̃^{(𝑀𝐿+1)𝐿,𝑗}=0 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 𝑝̃^𝑗(1),…,𝑝^{(𝑀𝐻)𝑗}̃,𝑝̃^(0)𝐿,𝑗,𝑝̃^(1)𝐿,𝑗,…,𝑝̃^{(𝑀𝐿)𝐿,𝑗} 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𝑝𝑝^(𝑖)_𝐿(𝑖)=[Φ]𝑝̃^(𝑖)𝑖=1,2,…,𝑀_𝐻,=[Φ]𝑝̃^𝐿(𝑖)𝑖=0,1,2,…,𝑀^𝐿,(5.1) 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 {𝑟}=1{0}=−[𝐴]{𝑝𝑐}−𝐸⁰2{𝑝}̇−[𝐴]𝑝⁽¹⁾,(5.2){0}=−[𝐴]𝑝_(𝑖−1)𝑐𝐸⁰𝑝̇⁽¹⁾−[𝐴]𝑝⁽²⁾,(5.3)−(−1)^𝑖2{0}=−[𝐴]𝑝_(𝑀𝐻)𝑐𝐸⁰𝑝̇^(𝑖)−[𝐴]𝑝^(𝑖+1)𝑖=2,3,…,𝑀_𝐻,(5.4)+(−1)^𝑀𝐻+1[𝐴]𝑝^𝐿(0)+(−1)^𝑀𝐻+11{0}=−1𝑐𝐸⁰𝑝̇^𝐿(0)−1𝑐𝐸⁰𝑝̇^𝐿(1),(5.5)𝑐𝐸⁰𝑝̇^{𝐿(𝑖−1)}+(−1)^{𝑀𝐻+𝑖+1}2[𝐴]𝑝^𝐿(𝑖)−1𝑐𝐸⁰𝑝̇^𝐿(𝑖+1)𝑖=1,2,…,𝑀^𝐿(5.6) with the symmetric and positive definite matrix [𝐴]=1ℎ[Φ]^−𝑇𝜆_𝑗[Φ]⁻¹.(5.7)

The residual term [𝐸⁰]{𝑝̇^{(𝑀𝐿𝐿+1)}}/𝑐 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𝑀_𝑔𝑢̈_𝑔+𝐶_𝑔𝑢̇_𝑔+𝐾_𝑔𝑢_𝑔=𝑓_𝑔.(5.8)

Here, subscript 𝑔 denotes global. [𝑀_𝑔], [𝐶_𝑔], and [𝐾_𝑔] are the global mass matrix, global damping matrix, and global stiffness matrix, which are expressed as 𝑀_𝑔=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐌^𝑠−𝐐_𝑓𝑠00𝐌_𝑓𝑓𝐌_𝑓𝑏𝐶0𝐌_𝑏𝑓𝐌0_𝑏𝑏⋱00⋱⋱00_𝑔0⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=1𝑐⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣000000𝐂00^∞0𝐘0₁₍₁₎⋱⋱⋱0𝐘_(𝑀𝐻)100𝐘0_(0)𝐿1−𝐂_∞−𝐂_∞⋱0⋱⋱−𝐂_∞𝐾_𝑔−𝐂_∞0⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐊^𝑠0𝐊𝐐_{𝑠𝑓𝑓𝑓}0𝐊_𝑓𝑏0𝐊_𝑏𝑓𝐊−𝐀_𝑏𝑏−⋱0𝐀−⋱⋱𝐀−−0𝐀𝐀𝐘−𝐀_(0)𝐿00𝐘0_(1)𝐿0⋱⋱⋱0𝐘0_{(𝑀𝐿)𝐿0}⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(5.9) with the block matrices𝐶_∞=𝐸⁰,𝑌_1(𝑖)𝑌^(0)𝐿0=(−1)^𝑀𝐻+1[𝐴],𝑌^(0)𝐿1=(−1)^𝑀𝐻+1𝐸⁰,𝑌_(𝑖)𝐿0=(−1)^𝑖2𝐸⁰𝑖=1,2,…,𝑀_𝐻,=(−1)^{𝑀𝐻+𝑖+1}2[𝐴]𝑖=1,2,…,𝑀_𝐿.(5.10)

The global load vector [𝑓_𝑔] and unknowns [𝑢_𝑔] are expressed as (the semicolon “;” stands for the vertical concatenation of vectors) 𝑢_𝑔=𝑓_𝑔𝑢_𝑠^𝑇;𝑝_𝑓^𝑇;𝑝_𝑏^𝑇;𝑝⁽¹⁾^𝑇;⋯;𝑝^{(𝑀𝐿)𝐿}^𝑇^𝑇,=𝑟_𝑠^𝑇;𝑟_𝑓^𝑇;{0}^𝑇;{0}^𝑇;⋯;{0}^𝑇^𝑇.(5.11)

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 𝑁_𝑠+𝑁_𝑓+𝑁_𝑓𝑏+(𝑀_𝐻+𝑀_𝐿+2)𝑁_𝑓𝑏, 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 (𝑀_𝐻+𝑀_𝐿+1)𝑁_𝑓𝑏 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{𝑟}=1𝑐𝐸⁰{𝑝}̇−ℎ1[Φ]^−𝑇𝜆_𝑗𝑝̃⁽¹⁾.(5.12)

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 ⎡⎢⎢⎢⎢⎣−𝑀𝑄_𝑠00_𝑓𝑠0𝑀𝑀_{𝑓𝑓𝑏𝑓}𝑀_𝑓𝑏𝑀_𝑏𝑏⎤⎥⎥⎥⎥⎦⎧⎪⎪⎨⎪⎪⎩𝑢̈^𝑠𝑝̈_𝑓𝑝̈_𝑏⎫⎪⎪⎬⎪⎪⎭+1𝑐⎡⎢⎢⎢⎢⎣000000𝐸00₀⎤⎥⎥⎥⎥⎦⎧⎪⎪⎨⎪⎪⎩𝑢̇^𝑠𝑝̇_𝑓𝑝̇_𝑏+⎡⎢⎢⎢⎢⎣𝐾_𝑠⎫⎪⎪⎬⎪⎪⎭0𝐾𝑄_{𝑠𝑓𝑓𝑓}00𝐾_𝑏𝑓𝐾_𝑓𝑏𝐾_𝑏𝑏⎤⎥⎥⎥⎥⎦⎧⎪⎪⎨⎪⎪⎩𝑢^𝑠𝑝_𝑓𝑝_𝑏⎫⎪⎪⎬⎪⎪⎭=⎧⎪⎪⎪⎨⎪⎪⎪⎩𝑓^𝑠𝑓_𝑓[Φ]_−𝑇𝜆_𝑗𝑝̃⁽¹⁾ℎ⎫⎪⎪⎪⎬⎪⎪⎪⎭.(5.13)

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 𝐾_𝐴𝑧_𝐴,𝑗(𝑡)+𝑐ℎ𝜆_𝑗𝐶_𝐴𝑧̇_𝐴,𝑗(𝑡)=𝑓_𝐴,𝑗(𝑡),(5.14) where the coefficient matrices are expressed as 𝐾_𝐴=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣01⋱⋱1⋱01(−110)^𝑀𝐻(−1)0_𝑀𝐻+1⋱2⋱⋱0(−1)0_{𝑀𝐻+𝑀𝐿}𝐶_𝐴2⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣(−1)⁰02⋱⋱0(−1)⋱_𝑀𝐻−102(−10)_𝑀𝐻1⋱01⋱⋱110⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦.(5.15)

The unknown vector {𝑧_𝐴,𝑗(𝑡)} consists of the modal auxiliary variables of mode 𝑗 (the semicolon “;” stands for the vertical concatenation of vectors)𝑧_𝐴,𝑗(𝑡)=𝑝̃^𝑗(1);⋯;𝑝^{(𝑀𝐻)𝑗}̃;𝑝̃^(0)𝐿,𝑗;𝑝̃^(1)𝐿,𝑗;⋯;𝑝̃^{(𝑀𝐿)𝐿,𝑗}.(5.16)

The only nonzero entry on the right-hand side is the modal hydrodynamic pressure 𝑝_𝑗̃ obtained from (4.5)𝑓_𝐴,𝑗(𝑡)=−𝑝_𝑗̃;0;⋯;0;0;0;⋯;0.(5.17)

It is demonstrated by Wang et al. [25] that (5.14) is stable up to the order 𝑀_𝐻=𝑀_𝐿=100 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 𝑡_𝑛+1 is determined by the last-solution extrapolation predictor [28, 29]𝑝̃⁽¹⁾^𝑝𝑛+1=𝑝̃⁽¹⁾_𝑛.(5.18)

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 {𝑧_𝐴,𝑗}₀=0 for each mode in (5.14).

Step 2. Extract {𝑝̃⁽¹⁾}_𝑛 from {𝑧_𝐴,𝑗}_𝑛=0 of each mode and assigned to {𝑝̃⁽¹⁾}^𝑝𝑛+1 in (5.18).

Step 3. Form the right-hand term of (5.13), compute {𝑢}_𝑛+1 and {𝑝}_𝑛+1 by using an implicit method.

Step 4. Calculate modal hydrodynamic pressure {𝑝}̃_𝑛+1 and form the right-hand term of (5.14).

Step 5. Compute {𝑧_𝐴,𝑗}_𝑛+1 for each mode by using an implicit method.

Step 6. Increment 𝑛 to 𝑛+1 and go to Step 2.

Since the predicted vector {𝑝̃⁽¹⁾}^𝑝𝑛+1 has been used rather than {𝑝̃⁽¹⁾}_𝑛+1 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, {𝑝̃⁽¹⁾}^𝑝𝑛+1 in Step 2 is extracted by the last computed {𝑧_𝐴,𝑗}_𝑛+1 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 𝑀_𝐻+𝑀_𝐿+1, 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 𝑀_𝐻=𝑀_𝐿=10 high-order doubly asymptotic open boundary condition is of excellent accuracy. So, the order 𝑀_𝐻=𝑀_𝐿=10 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 𝐸=35 GPa, Posisson’s ration 𝑣=0.2, and mass density 𝜌=2400 kg/m³. The water in the reservoir has a pressure wave velocity 𝑐=1438.7 m/s and mass density 𝜌=1000 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 2×7200/1438.7≈10s.

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.

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.