A non-reflecting boundary condition for the finite element modeling of infinite reservoir with layered sediment

Advances in Water Resources 29 (2006) 1515–1527 www.elsevier.com/locate/advwatres

A non-reflecting boundary condition for the finite element modeling of infinite reservoir with layered sediment Indrani Gogoi *, Damodar Maity Department of Civil Engineering, Indian Institute of Technology Guwahati, North Guwahati 781 039, Assam, India Received 4 May 2005; received in revised form 24 September 2005; accepted 9 November 2005 Available online 27 December 2005

Abstract The design of seismic resistant concrete gravity dam necessitates accurate determination of hydrodynamic pressure developed in the adjacent reservoir. The hydrodynamic pressure developed on structure is dependent on the physical characteristics of the boundaries surrounding the reservoir including reservoir bottom. The sedimentary material in the reservoir bottom absorbs energy at the bottom, which will affect the hydrodynamic pressure at the upstream face of the dam. The fundamental parameter characterizing the effect of absorption of hydrodynamic pressure waves at the reservoir bottom due to sediment is the reflection coefficient. The wave reflection coefficient is determined from parameters based on sediment layer thickness, its material properties and excitation frequencies. An analytical or a closed-form solution cannot account for the arbitrary geometry of the dam or reservoir bed profile. This problem can be efficiently tackled with finite element technique. The need for an accurate truncation boundary is felt to reduce the computational domain of the unbounded reservoir system. An efficient truncation boundary condition (TBC) which accounts for the reservoir bottom effect is proposed for the finite element analysis of infinite reservoir. The results show the efficiency of the proposed truncation boundary condition.  2005 Elsevier Ltd. All rights reserved. Keywords: Truncation boundary condition; Hydrodynamic pressure; Reservoir bottom absorption; Reflection coefficient; Layered sediment

1. Introduction A concrete gravity dam impounds water and entails high risk, as failure of a dam particularly during an earthquake may be disastrous for human life and property at the downstream. Hence, design and analysis of such a structure should be carried out with high level of accuracy. The dynamic behavior of a dam is complicated by the proximity of the reservoir and the flexible strata consisting of porous sediment on which it rests. The developed hydrodynamic pressure on structure is dependent on the physical characteristics of the bound-

*

Corresponding author. Tel.: +91 361 2300415; fax: +91 361 2690762. E-mail address: [email protected] (I. Gogoi). 0309-1708/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.11.004

aries surrounding the reservoir including reservoir bottom. For simplification of the analytical procedures, the reservoir bottom is generally considered to be rigid, which does not represent the actual behavior of the system. Such a representation results in erroneous estimation of the hydrodynamic pressure as no attempt is made to model the absorption of the pressure waves at the bottom surface. The hydrodynamic forces in the reservoir are affected by radiation of waves towards infinity and wave absorption at the reservoir bottom. When the reservoir bottom is considered to be rigid, the pressure waves are reflected from the reservoir bed and thus hydrodynamic pressure is over-estimated. Due to the absorption at the reservoir bottom, the magnitude of hydrodynamic pressure due to ground motion will be reduced. Therefore, the hydrodynamic pressure exerted on the adjacent structure will be less and consequently

1516

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

the displacement and stress field in the structure will be affected. Thus, an accurate evaluation of hydrodynamic pressure on the dam or offshore structure must consider the effect of sediments of the reservoir bottom. Studies pertaining to development of hydrodynamic pressure indicate that the presence of sediments at the reservoir bottom have significant effects on the response of dam [3,10,6] as they absorb part of energy of the hydrodynamic waves. To incorporate the reservoir bottom absorption effect in the evaluation of hydrodynamic pressure, the sediment layer has been modeled either as viscoelastic or porous media. The reflection coefficient is the fundamental parameter that characterizes the absorption effect that can be either considered to be frequency independent [9,13] or frequency dependent [10,5]. Use of such frequency dependent reflection coefficient would affect the hydrodynamic pressure determined using a dynamic solution procedure in the frequency domain. Studies conducted by Chuhan et al. [5] reveal that consideration of a viscoelastic or a porous model of the sediment layer gives similar results. The reservoir can be analyzed using various techniques and boundary conditions. For reservoir having simple boundary, classical solution treating reservoir as continuum may be used effectively. However, if the reservoir boundary is not regular, numerical methods such as finite element method may be used. Due to varied geometric forms of dam and irregular geometry of the reservoir, the finite element method (FEM) is recognized to be one of the most powerful numerical tools for solving such problems. However, the use of FEM for discretization of reservoir, a computational difficulty exists if the length of the reservoir is infinite. In such practical situation the infinite reservoir is truncated by imposing an artificial boundary at a certain distance away from the structure. A suitable boundary condition must be imposed at truncation surface to evaluate the effect of radiation damping, so that no spurious wave reflection occurs from the boundary. Such a boundary condition is called a non-reflecting boundary condition (NRBC). It is evident from the review on NRBC [8] that in the last 35 years NRBCs have been developed from zero-order Sommerfeld boundary condition to high order local boundaries. However, it is important that the NRBC is practically implementable and computationally efficient for a complex problem such as a dam–reservoir system considering the effect of reservoir bottom absorption. To evaluate the hydrodynamic response of a dam– reservoir system, the formulation of the radiation condition of the unbounded reservoir is the key issue in the analysis of fluid–structure interaction problems. The far boundary condition proposed by Maity and Bhattacharyya [11] is found to be effective in modeling the effects of radiation damping of the infinite reservoir. It was derived from the equation [15] that prescribes

hydrodynamic pressure at any point in the reservoir. The principle advantages of this technique are that it is quite simple to implement in the finite element formulation, which does not disturb the symmetrical and banded form of the matrix coefficients in the fluid equation. However, in case of absorptive reservoir bottom beyond the truncated boundary, it is necessary to include this condition at the truncation surface for accurate representation of the infinite reservoir for finite element analysis. Therefore, it did not include the effect of energy dissipation of the reservoir due to absorption of pressure waves at its bottom. Further, this boundary condition could not account for excitation frequency greater than the resonant frequency. Thus, modification of the truncation boundary condition proposed by Maity and Bhattacharyya [11] will lead to a more appropriate boundary condition representing the realistic behavior of a reservoir. In the present analysis, the infinite water domain has been analyzed by imposing truncation boundary condition at a certain distance away from the structure. The analysis is restricted to two-dimensional system and ground excitation is considered to be horizontal and in the direction perpendicular to the height of the structure. In the present work, the water is considered linearly compressible, inviscid and of small amplitude of motion. The governing equation of water is expressed in term of pressure only to reduce the number of degrees of freedom of the system. The proposed numerical model for dynamic analysis of a reservoir utilizes the wave reflection coefficient approach to account for the sediment properties at the reservoir bottom. The effect of depth of sediment layer resting on a rock on reflection coefficient is presented. A new TBC is proposed for the finite element analysis of infinite reservoir. This TBC eliminates the inadequacies of the TBC proposed by Maity and Bhattacharyya [11] through (i) incorporating the effect of reservoir bottom absorption of the infinite reservoir and (ii) evaluating the TBC in complex form which can account for excitation frequencies beyond the fundamental frequency of the reservoir. The effect of reservoir bottom absorption and the gravity waves on the development of hydrodynamic pressure due to horizontal ground motion is investigated. The numerical studies show that the proposed finite element algorithm in time domain could effectively evaluate the hydrodynamic pressure.

2. Mathematical modeling 2.1. Governing equation for water The governing equation to determine the magnitude of hydrodynamic pressure generated due to small amplitude vibration of compressible but non-viscous water in

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

a two-dimensional reservoir can be expressed by the classical Helmholtz wave equation as r2 pðx; y; tÞ ¼

1 € pðx; y; tÞ c2

ð1Þ

Here p(x, y, t) is hydrodynamic pressure, t is the time variable, x and y are space variables and c is the acoustic wave speed in water. 2.2. Boundary conditions Boundary conditions of the reservoir are as shown in Fig. 1. The origin of the horizontal x-axis and vertical yaxis is considered at the bottom of the dam–reservoir interface. The boundary conditions are defined at any point (x, y) relative to the origin.

2.2.3. At the reservoir bed interface (Cr) The absorption of pressure waves at the bottom of the reservoir is modeled by using the technique proposed by Hall and Chopra [9]. Thus the condition may be expressed as op ðx; 0; tÞ ¼ qan þ ixqpðxÞ on

1 op €p þ ¼0 g oy

ð2Þ

2.2.2. At the dam–reservoir interface (Cs) At the dam–reservoir interface, the pressures should satisfy op ð0; y; tÞ ¼ qaeixt on

ð3Þ

where aeixt is the horizontal component of the ground acceleration in which, x is the circular frequency pffiffiffiffiffiffiffi of vibration, q is the density of water and i ¼ 1; n is the outwardly directed normal to the elemental surface along the interface.

Y

The frequency independent reflection coefficient, a is given as 1  qqc s cs

ð6Þ

1 þ qqc s cs

Here, qs is the mass density of sediment, cs ispthe comffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pression wave velocity in the reservoir bed = ðEs =qs Þ, where Es = elastic modulus of the sediment. The sediment deposit may consist of layers of varying thickness and material characteristics. This necessitates an evaluation of reflection coefficient that can account for the varying thickness and material characteristics. While considering the effect of sediment layers, the equivalent reflection coefficient (a) is determined from the average of the reflection coefficients a1 and a2 for sediment layers of thickness ds1 and ds2 respectively as a¼

ða1 d s1 þ a2 d s2 Þ ðd s1 þ d s2 Þ

ð7Þ

A more realistic value of frequency dependent reflection coefficient a(x) can be obtained by considering a layer of single-phase viscoelastic sediment layer of height ds over

1 .. + ∂p = 0 g p ∂y

y=H

Truncation Boundary

Infinite Reservoir

= -ρ a eiw t

p=0, x=∞

∂n

∂p

ð4Þ

where n is the outwardly directed normal to the elemental surface along the interface. The coefficient q is given as   1 1a q¼ ð5Þ c 1þa

a¼ 2.2.1. At the free surface (Cf) Considering the effects of surface waves of the water, the boundary condition of the free surface is taken as

1517

∂p / ∂n = iωqp

X x=0

x=L

Fig. 1. Boundary conditions in a dam–reservoir system.

1518

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

a rigid rock. The complex frequency dependent reflection coefficient, B can be determined from the following equations given in matrix form [5] as 2

1

6 6 ik 6 6 60 4

iE1 k 1

iE1 k 1

0

qf x2

qf x2

0

eik1 d s

eik1 d s

eik2 d s

E1 k 1 eik1 d s E1 k 1 eik1 d s 3 2 3 2 B 1 7 6 7 6 6 A1 7 6 ik 7 7 6 7 6 6 7¼6 7 6 B1 7 6 0 7 5 4 5 4 0 A2

0

E2 k 2 eik2 d s

3 7 7 7 7 7 5

ð8Þ

In Eq. (8) E1 and E2 are elastic modulus, k1 and k2 are wave numbers in sediment and rock respectively. A1, A2 and B1 are constants obtained from the solution of one-dimensional wave equation. The detailed derivation can be found in the above reference. 2.2.4. At the truncation boundary (Ct) The specification of the far boundary condition is one of the most important features in the FE analysis of a semi-infinite or infinite reservoir. This is due to the fact that the developed hydrodynamic pressure, which affects the response of the structure, is dependent on the truncation boundary condition. Application of Sommerfeld radiation condition [16] at the truncation boundary leads to op ðL; y; tÞ ¼ 0 ox

ð9Þ

L represents the distance between the structure and the truncation boundary. Incorporating the effect of reservoir bottom absorption, Sharan [13] has incorporated the following condition: op f ðL; y; tÞ ¼  pðL; y; tÞ ox H

ð10Þ

pðx; y; tÞ ¼ 2qaH

1 X k2m I m ðkm xÞ e ðWm Þeixt ðbm Þk m m¼1

ð11Þ

where, a is the applied acceleration at the fluid–structure interface in the normal direction and RH 9 I m ¼ H1 0 Wm dy > > > > Wm ¼ 2k1m ½ðkm þ xqÞeikm y þ ðkm  xqÞeikm y  > > > > >  2 > > ð2m1Þp 2 > > þ i2xq=H km ¼ = 2H qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð12Þ > > k m ¼ k2m  X2 > > >   > > 2 v > 2 2 > bm ¼ ðkm  x q Þ H  k2 þv2 þ ixq > > m > > ; 2 x X ¼ x=c and v ¼ g The value of v can be considered to be zero [2] if the effect of surface gravity waves is neglected. For a rigid reservoir bottom, q = 0.0 and a = 1.0, for which the pressure obtained is real valued. Generally, determination of the eigenvalues km requires use of advanced programming techniques or specialized finite element program like Earthquake Analysis of Concrete Gravity Dams, EAGD [7]. Here, an alternative strategy [1] have been used to determine approximate values of km. Differentiating the hydrodynamic pressure in Eq. (11) with respect to x, one may arrive in the following expression: 1 X op k2m I m ðkm xÞ ðx; y; tÞ ¼ 2qaH e ðWm Þeixt ox ðb Þ m m¼1

ð13Þ

To obtain a finite model of the infinite reservoir, it is truncated at a surface, where the normal pressure gradient may [14] be expressed as op op ¼ ¼ fm p_ on ox From Eqs. (11), (13) and (14), fm is obtained as P k2m I m ðk m xÞ i 1 ðWm Þ m¼1 bm e fm ¼  P 2 1 Xc m¼1 bkmm kI mm eðkm xÞ ðWm Þ

ð14Þ

ð15Þ

Here, f is a complex damping parameter imposed at the truncation boundary. The main drawback of the above condition is that it is non-local and does not produce reasonably accurate results at second and third cutoff frequencies.

In the above equation, at the truncated boundary, the value of x will be L. The condition at the truncated boundary of the fluid domain considering the effect of non-reflective waves can be obtained from the plane wave equation [16] and can be expressed as

2.3. Proposed truncation boundary condition (TBC)

op 1 ¼  p_ on c

The proposed condition along the truncation boundary for the compressible fluid domain is derived from the wave equation. The general solution of Eq. (1) satisfying Eqs. (2)–(4) and the radiation condition, can be solved by method of separation of variables [2] to obtain pressure at any point (x, y) as

To obtain a finite model of the infinite reservoir from the analytical solution and to incorporate the effects of nonreflecting wave at the truncation surface the Eqs. (14) and (16) are combined and expressed as follows:   op 1 ¼ fm  p_ ð17Þ on c

ð16Þ

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

2.4. Finite element implementation By the use of Galerkin method, and assuming pressure to be the nodal unknown the discretized form of Eq. (1) may be written in two-dimension as  X  Z 1 X 2 € N i pi  2 N i pi dX ¼ 0 Nj r ð18Þ c X where, Nj is the interpolation function for the weight function and X is the region under consideration. Using Green’s theorem Eq. (18) may be transformed to  Z  oN j X oN i oN j X oN i  pi þ pi dX ox ox oy oy X Z Z X X oN j 1 oC pi ¼ 0 N i dX € Nj pi þ N j  2 c X on C ð19Þ

polate the generalized acceleration at any point on the fluid–structure interface in terms of generalized nodal accelerations of an element. At the reservoir bed, {Br} may be expressed as fBr g ¼ ixq½Rr fpg where ½Rr  ¼

XZ

ð26Þ

½N T ½T ½N dC

From the condition specified by Eq. (17) at the truncation boundary, the following expression emerges:   1 _ fBt g ¼  ½Rtc  þ ½Rt  fpg ð28Þ c where, ) PR T ½Rtc  ¼ ½N  fm ½N  dC Ct PR T ½Rt  ¼ ½N  ½N  dC Ct Substitution of all terms in Eq. (20) gives

½Ef€ pg þ ½Gfpg ¼ fBg

Here,

in which, 9 PR T > ½E ¼ c12 ½N  ½N dX > X > =   R P T T o o o o ½N  ½N  þ ½N  ½N  dX ½G ¼ X ox ox oy oy > > > PR T op ; fBg ¼ ½N  dC ¼ fB g þ fB g þ fB g þ fB g f s r t C on ð21Þ Here the subscript f, s, r and t stand for the free surface, structure–fluid interface, and reservoir bottom–fluid interface and truncation surface respectively. According to the boundary conditions for the fluid domain, if linearised surface wave condition is adopted, the same may be written in finite element form as 1 fBf g ¼  ½Rf f€ pg g in which, XZ ½Rf  ¼ ½N T ½N dC

ð22Þ

ð23Þ

Cf

At the structure–fluid interface, if {a} is the vector of nodal accelerations of generalized coordinates, {Bs} may be expressed as fBs g ¼ q½Rs fag

ð24Þ

where ½Rs  ¼

XZ

T

½N  ½T ½N s dC

ð25Þ

Cs

Here, [T] is the transformation matrix for generalized accelerations of a point on the fluid structure interface and [Ns] is the matrix of shape functions used to inter-

ð27Þ

Cr

in which i varies from 1 to total number of nodes and C represents the boundaries of the fluid domain. The whole system of Eq. (19) may be written in a matrix form as ð20Þ

1519

_ þ ½Gfpg ¼ q½Rs fag ½Ef€pg þ ½Afpg ½E ¼ ½E þ g1 ½Rf  ½A ¼ 1c ½Rt  þ ½Rtc  ½G ¼ ½G  ixq½Rr 

9 > = > ;

ð29Þ

ð30Þ

ð31Þ

For any prescribed acceleration at the fluid–structure interface, Eq. (30) is solved to obtain the hydrodynamic pressure in the fluid domain. In the present formulation, the Newmark time integration scheme has been adopted which is unconditionally stable with the parameters b and d as d P 0.5 and b P (0.5 + d)2.

3. Numerical results The feasibility of the proposed non-reflecting boundary condition at the artificially truncated boundary of the reservoir with sediment layer is examined. A computer code incorporating the algorithm described above is developed to determine the hydrodynamic pressure distribution on a two dimensional structure exposed to infinite fluid media. The structure is considered to be rigid and subjected to harmonic excitation in the horizontal direction. The full depth of the reservoir is considered as 70 m. The mass density of water q is considered as 1000 kg/m3 and the acoustic velocity of water, c as 1440 m/s. The amplitude of the external sinusoidal excitation is assumed to be equal to the gravitational acceleration g. The pressure coefficient at the bottom of the upstream face (cb = pb/qaH) is determined for different values of Tc/H. At Tc/H P 4, the hydrodynamic pressure coefficient is generally maximum at the base. In the proposed boundary condition the

1520

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

effect of reservoir bottom absorption is included, considering different values of reflection coefficients arbitrarily as 0.95 and 0.50. The finite element model consists of four noded isoparametric elements. Nv and Nh represent the number of elements in vertical and horizontal directions respectively. A convergence study has been carried out to arrive on a suitable time step for solution of fluid domain by the Newmark’s integration technique. In the present analysis, the number of time steps per cycle of excitation, Nt is considered as 32 at which the results are found to converge sufficiently. Thus the time step for the analysis of the infinite reservoir is T/32. Fig. 3. Variation of coefficient fm along dam height for a = 0.5.

3.1. Convergence of coefficient, fm researchers have considered it as either constant [13] or zero (Sommerfeld). 3.2. Effectiveness of proposed TBC To ascertain the effectiveness of the present truncation boundary condition, the hydrodynamic pressure coefficient is evaluated considering various geometries of the dam–reservoir system. The coefficient, fm for finding various results have been calculated by considering the value of m as 10 for all cases unless it is mentioned. 3.2.1. Dam with vertical upstream face The hydrodynamic pressure at the vertical upstream face of a dam is computed subjected to harmonic excitation. The maximum absolute hydrodynamic pressure can be obtained at an interval of 0.25T. However, to avoid initial unsteadiness the maximum hydrodynamic pressure coefficient at a time of 1.25T has been considered. The results are compared with closed form and other methods [1,13] and Sommerfeld). The hydrodynamic pressure coefficient, cp is the ratio of the hydrodynamic pressure to maximum hydrostatic pressure (cp = p/qaH). A comparison of hydrodynamic pressure coefficient cp along the dam height is plotted at Tc/H = 100, a = 0.5 (Fig. 4); Tc/H = 4, a = 0.95 (Fig. 5) and Tc/H = 1, a = 0.5 (Fig. 6). It is evident from

y/H

To get the effect of the unbounded fluid domain at the truncation surface, fm is determined Eq. (15) numerically considering reservoir bottom absorption. The value of m is assumed to be a large number, which leads to precise determination of normal derivative of complex valued pressure at the truncation boundary. The coefficient, fm obtained by Sharan [13] considers only the first term (m = 1) of the infinite series. The value thus obtained is real and does not depend on the location (x, y) of the truncation boundary. In the present analysis, fm is determined considering first several terms of the infinite series till the value is sufficiently converged. Therefore more accuracy is expected in the finite element solution as the proposed coefficient, fm is complex and accounts for the location of truncation boundary. The efficiency of the proposed truncation boundary condition can be established from the convergent results that is obtained by summation of few terms of the infinite series in Eq. (15). The convergence of the absolute values of the coefficient, fm at a distance L = 0.5H, a = 0.5 and at different depths of the fluid domain are plotted in Fig. 2. It is observed that the magnitude of fm converges for m = 10 for all cases (i.e., at various distance L and reflection coefficient, a). The variation of the absolute value of the coefficient, fm along the reservoir depth is plotted in Fig. 3. It is clearly seen that the coefficient fm is not constant at all though the previous

Fig. 2. Convergence of coefficient, fm.

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

Exact (Bouaanani 2003) Proposed (L = 0.1H) Sharan (L = 0.1H) 0.0

0.1

0.2

0.3

0.4 cp

0.5

0.6

0.7

0.8

Fig. 4. Distribution of absolute hydrodynamic pressure along a vertical surface (Tc/H = 100, a = 0.5).

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527 1.0

1521

0.8 Exact (Bouaanani 2003)

0.8

Sharan (L = 0.1H)

0.75

Proposed (L = 0.1H) cb

y/H

0.6 0.4

0.7 Sharan (1992)

0.2

0.65

Present Exact

0.0 0.0

1.0

2.0

3.0

4.0

5.0

0.6

cp

Fig. 5. Absolute hydrodynamic pressure distribution along a vertical surface (Tc/H = 4, a = 0.95).

0

0.2

0.4

0.6

0.8

1

L/H

Fig. 8. Convergence of absolute hydrodynamic pressure at the base (Tc/H = 100, a = 0.95).

1.0 4.60 Sharan (1992) Present Exact

4.58

0.8

4.56 4.54

y/H

0.6

4.52 cb

Exact (Bouaanani 2003) Proposed (L = 0.2H)

0.4

Sharan (L = 0.2H)

4.50 4.48 4.46

0.2

4.44 0.0 0.00

4.42 0.05

0.10

0.15 cp

0.20

0.25

0.30

0.1

0.3

0.5

0.7

0.9

L/H

Fig. 6. Absolute hydrodynamic pressure distribution along a vertical surface (Tc/H = 1, a = 0.5).

Fig. 9. Convergence of absolute hydrodynamic pressure at the base (Tc/H = 4, a = 0. 95).

the results that the present TBC is accurate even if the boundary is truncated at a closer distance with a wide range of Tc/H. The total hydrodynamic pressure, ct is the ratio of total hydrodynamic pressure to total hydrostatic pressure. The plot in Fig. 7 shows a comparison of total hydrodynamic pressure coefficient (ct) at the dam interface for Tc/H = 100, L/H = 0.1 and a = 0.5, which further confirms the accuracy of the proposed boundary condition.

A convergence test with respect to the truncated length of the reservoir (L/H) of the finite element model is conducted for different values of Tc/H. The hydrodynamic pressure coefficients obtained for different values of Tc/H are plotted in Figs. 8–11 and tabulated in Table 1 for reflection coefficient values of 0.95 and 0.5. It is interesting to note that the results obtained by the present TBC are undistinguishable from the exact solution

0.8 1.5 Bouaanani (2003)

1.0

0.75

Proposed, L = 0.1H

ct

cb

0.5 0.0 0.00 -0.5

0.20

0.40

0.60

0.80

0.7 Sharan (1992)

1.00

0.65

Present Exact

-1.0

0.6 0

-1.5 Time (second)

Fig. 7. Time history of the absolute total hydrodynamic pressure coefficient on a vertical surface (Tc/H = 100, a = 0.5).

0.2

0.4

0.6

0.8

1

L/H

Fig. 10. Convergence of absolute hydrodynamic pressure at the base (Tc/H = 100, a = 0.5).

1522

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

be concluded that the present truncation boundary condition has effectively accounted for the effect of the reservoir bottom absorption for all the important frequency range that is significant in earthquake response analysis of dams. The proposed TBC condition has the advantage of producing more accurate results for very low to very high ranges of frequencies and even at resonant frequency for both the highly absorptive and reflective cases even by using a very coarse mesh. Since an earthquake excitation may be a composition of randomly different frequencies, the present TBC is found to be very efficient computationally.

1.50

cb

1.25

1.00

Sharan (1992) Present Exact

0.75 0.1

0.3

0.5 L/H

0.7

0.9

Fig. 11. Convergence of absolute hydrodynamic pressure at (Tc/ H = 4, a = 0.5).

at a = 0.5 and 0.95 for Tc/H = 100 (Figs. 8 and 10) even at a very closer distance. Very small error is observed at Tc/H = 10 (Table 1). At the resonant frequency, having a = 0.5 and 0.95, the present TBC is convergent at L = 0.2H; whereas the Sharan’s TBC is observed to be convergent at L = 0.5H (Fig. 11). The present TBC is also suitable at a distance of L = 0.2H (Table 1) for Tc/H = 1.0. It is further noted that the reflection coefficient of the reservoir bottom affects the accuracy of the solution at excitation frequencies higher than the fundamental frequency of the reservoir. These conclusions are in agreement with Sharan [13] carried out in the frequency domain analysis. At lower excitation frequencies (Tc/H = 10, 100), the reflection coefficient does not have effect on the accuracy of the results obtained. It can thus

3.2.2. Dam with inclined upstream face The inclination of the upstream face of the dam (Fig. 12) would have effect on the development of hydrodynamic pressure [4,13] in the dam–reservoir system. No classical solution considering reservoir bottom effect for inclined upstream face is observed in the open literature accessible to the present authors. The effectiveness of the proposed TBC is tested by comparing the hydrodynamic pressure evaluated at a large distance (L = 4.0H) using Sommerfeld boundary condition. The hydrodynamic pressure at the base of the inclined upstream dam face for different Tc/H ratio is presented in Tables 2 and 3. The results show that the infinite reservoir may be truncated at a relatively closer distance away from the structure in case of lower excitation frequency (i.e., Tc/ H = 100). Further, it is observed that at Tc/H = 4.0, the accuracy on the development of hydrodynamic pressure is less at a closer distance. Thus, a length of reservoir equal to its depth should be considered for a fully

Table 1 Comparison of truncation boundary condition for vertical surface (h = 90) Tc/H

L/H

Mesh size Nv · Nh

a

Exact

Truncation boundary condition

cm

cm

Sommerfeld

1

4

10

100

0.1 0.2 0.1 0.2

20 · 2 20 · 4 20 · 2 20 · 4

0.95

0.1073

0.5

0.1414

0.2 0.5 1.0 0.2 0.5 1.0

20 · 4 20 · 10 20 · 20 20 · 4 20 · 10 20 · 20

0.95

4.4364

0.5

1.2187

0.1 0.2 0.1 0.2

20 · 2 20 · 4 20 · 2 20 · 4

0.95

0.8155

0.50

0.8097

0.02 0.10 0.02 0.10

50 · 1 20 · 2 50 · 1 20 · 2

0.95

0.7431

0.50

0.7430

Sharan [13] % Error

cm

1.9892 1.2883 1.2343 0.2306

1853.87 1200.65 772.91 63.08

0.1552 0.1251 0.1397 0.143

60.7411 24.0546 12.3857 9.61 1.0873 1.3319

1270.00 442.00 179.00 688.6 10.74 9.34

1.7299 1.6351 1.7514 1.6493 8.7381 4.7213 8.741 4.7205

Proposed BC % Error

cm

% Error

44.64 16.59 1.20 1.13

0.1072 0.1073 0.1397 0.1414

0.09 0.00 1.20 0.00

4.4804 4.455 4.441 1.2572 1.2225 1.2214

0.99 0.41 0.09 3.16 0.31 0.22

4.4369 4.4365 4.4364 1.2206 1.2210 1.2204

0.01 0.00 0.00 0.15 0.19 0.14

112.15 101.21 116.30 103.69

0.684 0.8533 0.7178 0.7547

16.12 5.35 14.16 6.79

0.7971 0.8155 0.7964 0.8092

2.26 0.00 1.64 0.06

1075.96 535.39 1076.41 535.31

0.6559 0.7074 0.6559 0.7073

11.73 4.80 11.73 4.80

0.7417 0.7427 0.7417 0.7424

0.18 0.06 0.18 0.08

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

1523

ervoir is considered to be inclined (Lhi = 0.5H), beyond which the reservoir bed is assumed to be horizontal (Lh). To the authors’ knowledge, there is no classical solution available for the analysis of hydrodynamic pressures for a non-horizontal reservoir bed considering the effects of compressibility of water and reservoir bottom absorption. Hence, the parametric study has been conducted by comparing results obtained by applying Sommerfeld boundary condition or present TBC at L = 4.0H, which gave identical results. Since the effect of reservoir bottom is almost negligible for higher values of Tc/H, as is evident from the preceding cases, here only Tc/ H = 10, 4 and 1 is considered and results are shown in graphical form. The pressure distribution for Tc/ H = 10 and a = 0.5 are plotted in Figs. 14 and 15 for different direction of bed inclinations. The results show that the proposed TBC is also effective enough even when the reservoir bed is inclined. Similar observations have been reflected from Fig. 16 with a higher excitation frequency (Tc/H = 1). However, the accuracy could not be arrived at a closer distance if the reservoir is excited at its resonant frequency (Tc/H = 4) as seen from Fig. 17. A similar observation in time domain analysis was also reported by [12].

H

θ L Fig. 12. Typical finite element model of infinite reservoir with inclined upstream dam face.

inclined dam for accurate evaluation of hydrodynamic pressure at resonant frequency. At Tc/H = 1.0, convergence is obtained within a very short length of reservoir (L = 0.1H). However, since the proposed TBC is derived assuming upstream face of the reservoir to be vertical, the accuracy and effectiveness of the TBC decreases with the increase of dam inclination.

3.3. Effect of surface gravity waves on hydrodynamic pressure

3.2.3. Dam with vertical upstream face and inclined reservoir bed A typical finite element discretization is shown in Fig. 13 for the analysis of infinite reservoir with inclined bed. The length of reservoir up to half the depth of res-

The present algorithm accounts for the effect of surface gravity waves in the finite analysis of infinite reservoir. Fig. 18 shows the influence of surface gravity wave

Table 2 Comparison of truncation boundary condition for inclined surface (h = 60) Tc/H

L/H

a

Sommerfeld L = 4H

Truncation boundary condition Sommerfeld cm

1

4

10

100

0.10 0.50 0.10 0.50

0.95

0.1039

0.5

0.1385

0.10 0.50 1.00 0.10 0.50 1.00

0.95

1.1846

0.5

0.9087

0.95

0.5208

0.5

0.5153

0.95

0.4710

0.5

0.4710

0.1 0.5 1.0 0.1 0.5 1.0 0.02 0.2 0.02 0.2

Sharan [13] % Error

cm

Proposed BC % Error

cm

% Error

0.1187 0.1037 0.142 0.139

14.24 0.19 2.53 0.36

0.1056 0.1041 0.1398 0.1390

1.64 0.19 0.94 0.36

0.1040 0.1039 0.1391 0.1390

0.10 0.00 0.43 0.36

0.6439 0.7195 0.9942 0.6572 0.7091 0.8771

45.64 39.26 16.07 27.68 21.97 3.48

0.6488 0.8557 1.1914 0.4862 0.5103 0.9080

45.23 27.76 0.57 46.49 43.84 0.08

0.7431 0.9858 1.1835 0.5065 0.6300 0.9080

37.27 16.78 0.09 44.26 30.67 0.08

0.8093 0.6176 0.5455 0.8091 0.6172 0.5452

55.40 18.59 4.74 57.02 19.77 5.80

0.4765 0.4953 0.5417 0.4753 0.4947 0.5183

8.51 4.90 4.01 7.76 4.00 0.58

0.4788 0.4954 0.5186 0.4775 0.4948 0.5183

8.06 4.88 0.42 7.34 3.98 0.58

1.5045 0.8925 1.5045 0.8925

219.40 89.47 219.41 89.48

0.4456 0.4699 0.4456 0.4699

5.40 0.24 5.40 0.24

0.4517 0.4710 0.4517 0.4710

4.11 0.01 4.10 0.01

1524

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

Table 3 Comparison of truncation boundary condition for inclined surface (h = 30) Tc/H

L/H

Sommerfeld L = 4H

a

Truncation boundary condition Sommerfeld cm

1

4

10

0.1 0.5 0.1 0.5

0.95

0.1241

0.5

0.1344

0.1 0.5 1.0 0.1 0.5 1.0

0.95

0.8329

0.5

0.5432

0.95

0.2628

0.5

0.2553

0.95

0.2400

0.5

0.2400

0.1 0.5 1.0 0.1 0.5 1.0

100

0.02 0.20 0.02 0.20

Sharan [13] % Error

cm

% Error

0.81 12.09 4.17 9.08

0.1275 0.1240 0.1336 0.1336

2.74 0.08 0.60 0.60

0.1261 0.1241 0.1335 0.1339

1.61 0.00 0.67 0.37

0.5242 0.6133 0.809 0.3216 0.4572 0.5407

37.06 26.37 2.87 40.80 15.83 0.46

0.3310 0.4269 0.8332 0.1996 0.3050 0.5430

60.26 48.75 0.04 63.25 43.85 0.04

0.5164 0.5975 0.8332 0.2371 0.3658 0.5430

38.00 28.26 0.04 56.35 32.66 0.04

0.3622 0.3025 0.2802 0.3620 0.3024 0.2800

37.82 15.11 6.62 41.79 18.45 9.67

0.2291 0.2493 0.2605 0.2287 0.2490 0.2588

12.82 5.14 0.88 10.42 2.47 1.37

0.2292 0.2581 0.2605 0.2287 0.2550 0.2555

12.79 1.79 0.88 10.42 0.12 0.08

0.5311 0.3794 0.5311 0.3794

121.29 58.08 121.66 58.35

0.1714 0.2408 0.1714 0.2394

28.58 0.33 28.46 0.13

0.1715 0.2393 0.1715 0.2394

28.54 0.29 28.42 0.13

y/H

C

H

L

Proposed BC % Error

0.1231 0.1391 0.1288 0.1466

A

ϕ

cm

B

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

Sommerfeld, Lh/H = 4.0, ϕ = 10 Proposed, Lh/H = 0.2, ϕ = 10 0.0

0.1

0.2

0.3

0.4

Lhi L

0.5

0.6

0.7

0.8

0.9

|cp|

Fig. 14. Hydrodynamic pressure distribution at reservoir bed slope of / = +10, Tc/H = 10.

in the development of hydrodynamic pressure at the base of the dam–reservoir interface. It is interesting to note that a small peak occurs at an early stage (x/ xn < 0.1) due to the presence of gravity waves. However, the second peak coincides with the first mode when the effects of surface gravity waves are not considered. The results show that there is a negligible influence of surface gravity waves in the hydrodynamic pressure beyond the excitation frequency ratio (x/xn) of 0.1. Here, x is the excitation frequency and xn is the natural frequency of the reservoir with rigid bottom. It is also observed from Fig. 18 that the changes in reflection coefficient do not affect the amplitude of the first resonant peak consider-

y/H

Fig. 13. Finite element discretization of an infinite reservoir having an inclined bed. 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

Sommerfeld, Lh/H = 4.0, ϕ = -10 Proposed Lh/H = 0.1, ϕ = -10 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

|cp|

Fig. 15. Hydrodynamic pressure distribution at reservoir bed slope of / = 10, Tc/H = 10.

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527 Table 4 Effect of gravity waves on hydrodynamic pressure coefficient

1.0 0.9 0.8 0.7 0.6 y/H

1525

Tc/H

x/xn

200.0

0.02

v 0.00

Sommerfeld, Lh = 4.0H Present, Lh = 0.5H Present, Lh = 0.1H

0.5 0.4 0.3 0.2

0.04 100.0

0.1 0.0 0.00

a

0.04

0.00 0.17

0.05

0.10

0.15 |cp|

0.20

0.25

0.30

4.0

1.0

0.00

Fig. 16. Hydrodynamic pressure distribution for inclined reservoir bed slope at Tc/H = 1, a = 0.95.

106.42 1.0

4.0

0.00 1703.00

1.0

cp

1.0 0.5 1.0 0.5

0.7416 0.7416 1.0219 1.0219

1.0 0.5 1.0 0.5

0.74209 0.74206 0.80942 0.80937

1.0 0.5 1.0 0.5

54.801 1.2182 54.8090 1.2184

1.0 0.5 1.0 0.5

0.1143 0.1422 0.1143 0.1422

0.9 0.8

y/H

0.7 0.6

Table 5 Effect of reservoir bottom absorption on fundamental frequency of the reservoir

0.5 0.4

Sommerfeld, Lh = 4.0H Present, Lh = H Present, Lh = 0.5H Present, Lh = 0.2H

0.3 0.2 0.1 0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

|cp|

Fig. 17. Hydrodynamic pressure distribution for inclined reservoir bed at Tc/H = 4, a = 0.5.

Reflection coefficient, a

Fundamental frequency, xnf

Time period, Tn (s)

0.25 0.50 0.95 1.00

28.0 31.0 32.0 32.3

0.224 0.203 0.196 0.195

frequencies gravity waves have negligible effect and reservoir bottom absorption have more prominent effect on the development of hydrodynamic pressure.

1.4 α = 0.25, χ > 0.0 α = 0.50, χ > 0.0 α = 0.50, χ = 0.0

1.3

cb

1.2

3.4. Effect of reservoir bottom absorption on fundamental frequency of reservoir

1.1 1.0 0.9 0.8 0.7 0.6 0.0

0.2

0.4

0.6 ω /ω n

0.8

1.0

1.2

Fig. 18. Effect of gravity waves on hydrodynamic pressure coefficient.

ing the surface gravity waves. This is in agreement with findings reported by Bouaanani et al. [2]. The authors attribute the occurrence of this phenomenon to the presence of ice layer at the reservoir surface. However, it is observed from the present numerical results that the first peak may occur even in the absence of the ice layer. The results in Table 4 shows further that the effect of surface gravity waves is prominent in the low frequencies (i.e., at Tc/H = 200 and 100) where the absorptive reservoir bottom does not have much effect. For higher excitation

The fundamental frequency (xnf) of the reservoir with absorptive bottom has been determined and presented in Table 5. The results show that the fundamental frequency is reduced elongating the fundamental time period of the reservoir in case of lower magnitude of the reflection coefficient. This is because the energy is being dissipated as the flexibility of the reservoir bottom is increased. Similar observations are seen from Fig. 19 where the hydrodynamic pressure coefficient at the base of the upstream face of the dam is reduced substantially for lower magnitude of reflection coefficient. 3.5. Sediment effect on reflection coefficient The results presented in the previous sections emphasize that the proposed truncation boundary condition is effective in calculating the hydrodynamic pressure using the finite element technique. Assuming a higher reflection coefficient (a = 1.0) implies a highly reflective reservoir

1526

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

Fig. 19. Effect of reservoir bottom absorption on hydrodynamic pressure coefficient, without gravity waves.

0.07

0.9

0.06

Reflection Coefficient

1.0

Reflection Coefficient

0.8 0.7 0.6 0.5 0.4

Single Sediment Layer, E1 Two Sediment Layer, E1/E2 = 0.5 Two Sediment Layer, E1/E2 = 2.0

0.05 0.04 0.03 0.02

0.3

0.01

0.2

0

0.1

0

0.2

0.4

0.6

0.8

1

ds/H

0.0 0

500

1000

1500

2000

2500

3000

3500

4000

Fig. 21. Effect of sediment layers on reflection coefficient.

Density (kg/m3)

Fig. 20. Variation of reflection coefficient.

1.2 Reflection coefficient

bottom, which may yield exaggerated values of hydrodynamic pressure. Hence, the characteristics of the overlying sediment layer play a significant role in determining the reflection coefficient and subsequently, the hydrodynamic pressure. To determine the effect of sediment property on the reflection coefficient, sediment of elastic modulus = 35 000 MPa with varying density is considered. It is observed from Fig. 20 that the mathematical modeling as given in Eqs. (5)–(7) would be appropriate to evaluate the reflection coefficient with varying density. The results show that the increase in reflection coefficient is more in the sediments having low density. Hence, higher damping at the reservoir bottom will occur with the presence of low density sediments. The effect of thickness of sediment layer on the variation of magnitude of reflection coefficient is further studied and shown in Fig. 21. Depths of two sediment layers are considered in this particular case as equal. The moduli of elasticity of sediments for first and second layer are marked as E1 and E2 respectively in the figure. It is interesting to note that the reflection coefficient is decreased and thereby damping is increased if the underlying layer consist of soft sediment (i.e., if E2 < E1). This implies that the property of the bottom

layer is significant in ascertaining the value of reflection coefficient. The frequency dependent reflection coefficient (Eq. (8)) is further studied considering the viscoelastic modeling of the sediment layer lying over a rigid rock. Here, the sediment is considered to be of density 1962 kg/m3 and elastic modulus of 5.12 · 105 MPa. It is seen from Fig. 22 that at lower excitation frequency the variation with sediment depth is almost negligible. In the figure, n is the excitation frequency in Hz. But for higher frequencies, the variation in sediment height is pronounced. Hence, if the frequency content of an earth-

1 0.8 n

0.6

10 20 30 40 50

0.4 0.2 0 0

10

20 Depth in m

30

40

Fig. 22. Variation of frequency dependent reflection coefficient.

I. Gogoi, D. Maity / Advances in Water Resources 29 (2006) 1515–1527

quake excitation is high, it would necessitate incorporating the effect of frequency dependent reflection coefficient in the analysis procedure.

4. Conclusions For effective implementation of a finite element model in the dynamic analysis of a dam–reservoir system, a new truncation boundary is proposed to model the infinite reservoir as a finite one. To get accurate results, a non-reflecting boundary condition is developed and imposed at the truncated boundary of the reservoir. The truncation boundary condition is complex valued and obtained as a derivation of the exact solution that can be used for determining hydrodynamic pressure at the upstream face of a rigid dam. The proposed technique is easy to implement in a finite element solution as the truncation boundary coefficient can easily be incorporated. It is observed from the results that the proposed boundary condition produces accurate results for all ranges of excitation. It is observed from the results that the convergence is attained at a very short distance away from the dam–reservoir interface (i.e., L = 0.02H to 0.2H) for all Tc/H of interest for the case when the upstream face of the dam is vertical. However, for the irregular geometry of the computational domain having inclined upstream face or inclined reservoir bed, the convergence would depend on inclination, reflection coefficient and excitation frequency. The present TBC is found to perform efficiently in all frequency range of interest when used in a time domain analysis procedure. Since the truncation boundary condition is developed from the case when the dam face is vertical and the reservoir bed profile is horizontal, an increase in percentages of error is observed for higher inclinations. For the case of inclined reservoir bed, the error percentage is found to be less for an increase in reservoir bed slope (/) in negative direction with lower reflection coefficient. The effect of reservoir bottom absorption is found to be very small for higher values of Tc/H. But this effect cannot be neglected for values of Tc/H less than and nearer to 4. Further, the effect of gravity waves is found to be considerable for higher values of Tc/H. A study of the frequency dependent reflection coefficient reveals that for lower excitation frequency the depth of sediment layer does not have much effect on

1527

the reflection coefficient. But for an increasing excitation frequency it is observed that an increase in depth of sediment layer would reduce the reflection coefficient considerably. If the frequency content of an earthquake excitation is high, it would necessitate incorporating the effect of frequency dependent reflection coefficient in the analysis procedure.

References [1] Bouaanani N, Paultre P, Proulx J. A closed-form formulation for earthquake-induced hydrodynamic pressure on gravity dams. J Sound Vib 2003;261:573–82. [2] Bouaanani N, Paultre P, Proulx J. Two-dimensional modeling of ice cover effects for the dynamic analysis of concrete gravity dams. Earthquake Eng Struct Dynam 2002;31:2083–102. [3] Bougacha S, Tassoulas JL. Effect of sediment material on the response of concrete gravity dams. Earthquake Eng Struct Dynam 1991;20:849–58. [4] Chwang AT. Hydrodynamic pressures on sloping dams during earthquakes. Part 2. Exact theory. J Fluid Mech 1978;78(2): 343–8. [5] Chuhan Z, Chengda Y, Guanglun W. Numerical simulation of reservoir sediment and effects on hydrodynamic response of arch dams. Earthquake Eng Struct Dynam 2001;30:1817–37. [6] Domı´nguez J, Gallego R, Japo´n BR. Effects of porous sediments on seismic response of concrete gravity dams. J Eng Mech 1997;123(4):301–11. [7] Fenves G, Chopra AK. A computer program for earthquake analysis of concrete gravity dams, Report No. UCB/EERC-4/11, EAGD-84, 1984. [8] Givoli D. High-order local non-reflecting boundary conditions: a review. Wave Motion 2004;39:319–26. [9] Hall JF, Chopra AK. Two-dimensional analysis of concrete gravity and embankment dams including hydrodynamic effects. Earthquake Eng Struct Dynam 1982;10:305–32. [10] Hatami K. Effect of reservoir bottom on earthquake response of concrete dams. Soil Dynam Earthquake Eng 1997;16:407–15. [11] Maity D, Bhattacharyya SK. Time-domain analysis of infinite reservoir by finite element method using a novel far-boundary condition. Finite Elem Anal Des 1999;32:85–96. [12] Pavic A, Zivanovic S, Reynolds P. Errors in numerical solution of equation of motion of lightly damped SDOF system near resonance. Int J Struct Stabil Dynam 2005;5(1):135–42. [13] Sharan SK. Efficient finite element analysis of hydrodynamic pressure on dams. Comput Struct 1992;42(5):713–23. [14] Sharan SK. Finite element modeling of infinite reservoirs. J Eng Mech 1985;111(12):1457–69. [15] Westergaard HM. Water pressure on dams during earthquakes. Trans ASCE 1933;98:418–72. [16] Zienkiewicz OC, Newton RE. Coupled vibration of a structure submerged in a compressible fluid. In: International symposium on finite element techniques, Stuttgart, Germany, 1969.