# A novel far-boundary condition for the finite element analysis of infinite reservoir

Applied Mathematics and Computation 170 (2005) 1314–1328

www.elsevier.com/locate/amc

A novel far-boundary condition for the ?nite element analysis of in?nite reservoir

Damodar Maity

Department of Civil Engineering, Indian Institute of Technology, Guwahati 781 039, India

Abstract The present paper deals with the ?nite element analysis of the reservoir of in?nite extent using a novel far-boundary condition. The equations of motion are expressed in terms of the pressure only assuming water as inviscid and incompressible. The truncation boundary condition is developed numerically from the classical wave equation. Comparative studies show that the proposed far-boundary condition is numerically e?cient and accurate over the existing ones, available in the literature. The e?ect of the geometry of the reservoir bed and the adjacent structure on the development hydrodynamic pressure has been studied. The results show that the geometry of the reservoir bed and as well as the adjacent structure has considerable e?ect on the development of hydrodynamic pressure at the dam–reservoir interface. ? 2005 Elsevier Inc. All rights reserved.

Keywords: In?nite reservoir; Far-boundary condition; Finite element method; Hydrodynamic pressure; Radiation boundary

E-mail address: damodar@iitg.ernet.in 0096-3003/$ - see front matter ? 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.01.020

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1315

1. Introduction A large variety of structures subjected to predominantly ?uid loading need to be analysed precisely for external exciting forces. The estimation of precise hydrodynamic forces on dam faces due to earthquakes is an important aspect of the analysis and design of dams. The hydrodynamic pressure on vertical rigid structure subjected to ground motion was ?rst solved analytically by Westergaard [1]. In most of the practical problems, it is di?cult to obtain closed form analytical solutions due to complex geometries of the dam–reservoir systems. The ?nite element method (FEM) is recognised as a powerful numerical tool for solving such practical problems. In the ?nite element analysis of such problems, di?culties arise mainly because of the large extent of the ?uid domain, where the ?uid is unbounded. There are a variety of far-boundary conditions reported in the literature. These may be broadly classi?ed as (i) imposition of a boundary condition on the truncation surface [2–8], and (ii) coupling the ?nite element discretisation with other type of discretisation such as ?in?nite elements? [9], ?boundary elements? [10,11]. There is a common belief that boundary element method (BEM) is superior to the FEM for the modeling of in?nite or semi-in?nite domains. However, in the reported literature [12], the e?ciency of BEM in time domain analysis is not ascertained. This is because of the presence of the convolution integral and singularity of the kernels of the formulation which requires large storage space and computational time for the evaluation of the e?ect of past time history and numerical integration of the kernels. The distance between the structure and truncation surface is not considerably less while adopting the in?nite elements [9]. On the other hand, the ?nite element approach has the distinct advantage of being straightforward in implementation. The most commonly used boundary condition along the truncation surface is the Sommerfeld radiation condition [2,3]. However, for an incompressible ?uid, Sommerfeld or other similar boundary condition takes the form of that for a rigid stationary boundary. As a result, the behavior of the ?uid motion at the truncation boundary is not represented truly. Hence a large extent of ?uid domain is required to be included in the analysis. All the above-mentioned cases, the boundary conditions have been developed considering the structure as rigid [1–7,12–14]. The focus of the present paper is on the development of an e?cient truncation boundary condition to model the in?nite reservoir, which includes the radiation e?ects properly and can be adopted in the ?nite element formulation in a simple form. Most of the far-boundary conditions proposed to date are approximate and require the distance of the truncated boundary su?ciently away from the structure. To obtain the e?ect of the unbounded water at the structure–reservoir interface, an equation along the truncation boundary has been developed. The coe?cients along the truncation surface are determined numerically. A ?nite element technique of analyzing the hydrodynamic

1316

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

pressure developed in the reservoir resulting from the vibration of the structure is presented. In the present work, the equation of motion governing the ?uid is expressed in terms of pressure alone considering the ?uid as incompressible, inviscid and irrotational. Comparative studies show clearly the e?ciency and e?ectiveness of the proposed far-boundary condition.

2. Theoretical formulation 2.1. Governing equations for ?uid and boundary conditions Assuming the water in the reservoir to be inviscid and incompressible and its motion to be of small amplitude, the hydrodynamic pressure equation will be: r2 p ? 0;

2

?1?

in which, $ is the Laplacian operator and p is the hydrodynamic pressure in the ?uid. The hydrodynamic pressure distributions within the domain may be obtained by solving Eq. (1) with the following boundary conditions (Fig. 1). (a) At the free surface (Cf): Neglecting the e?ects of surface waves of the water, the boundary condition of the free surface is: p ? 0: (b) At the structure–?uid interface (Cs): op ? ?qan ; on

y y=H

f

?2?

?3?

p=0

? p/? n = ρan

Γs Rigid Structure Γb a x=L Γt

Truncation Surface p=0

x=∞

x

? p/? n = 0

Fig. 1. Problem geometry and boundary condition ?uid domain.

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1317

where q is the mass density of ?uid, on is the outwardly directed normal to the elemental surface along the interface, and an is the normal acceleration of the solid surface in the direction of on. (c) At the reservoir bed–?uid interface (Cb): Assuming the reservoir ?oor to be rigid, the condition adopted is op=on ? 0: ?4? (d) At truncation boundary (Ct): The speci?cation of the far-boundary condition is one of the most important features in the development of the reservoir model due to the fact that the hydrodynamic pressure on the structure is highly sensitive to the behavior of the in?nite region of the reservoir. Most of the far-boundary conditions proposed to date are approximate and requires the distance of the truncated boundary su?ciently away from the structure. In order to consider the e?ect of radiation damping, it is assumed that Eq. (2) is satis?ed at an in?nitely large distance away from the structure. If the unbounded ?uid domain is truncated at a su?ciently large distance away from the region of interest, Sommerfeld radiation boundary [2] is normally used at the truncation surface. For an incompressible ?uid, Sommerfeld radiation boundary condition satis?es Eq. (4). However, the Sommerfeld boundary at the truncation surface is the same as that for a rigid stationary boundary. Sharan [4] proposed the truncation boundary as follows which seems better than the Sommerfeld radiation boundary, but is an approximate one. pp op=on ? ? ; ?5? 2H where, H is the depth of the ?uid domain.

2.2. Proposed far-boundary condition: The proposed boundary condition along the truncation surface is derived on the basis of the following assumptions: (i) The bottom of the ?uid domain is horizontal and rigid. (ii) The ?uid–structure interface is vertical. (iii) The ?uid domain extends to in?nity and its motion is two-dimensional. The general solution of Eq. (1) satisfying Eqs. (2)–(4) and the radiation condition, at any point (x, y) is given by Westergaard [1] as: p ? 2an qH

1 y X ??1?n?1 x e??kn H ? cos kn ; 2 H kn n?1

?6?

1318

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

where, kn ? ?2n ? 1?p ; 2 ?7?

an is the acceleration of the ?uid–structure interface in the normal direction. The partial derivative of the hydrodynamic pressure with respect to x considering Eq. (6) will be 1 y X ??1?n?1 op x ? ?2an q e??kn H ? cos kn : ?8? ox H kn n?1 From Eqs. (6) and (8), the following equation is obtained. op op p ? ? ? f; on ox H ?9?

where, p is the hydrodynamic pressure of the ?uid domain and f is given by ? y? P1 ??1?n?1 ??kn x ? e H cos kn H n?1 kn f?P ?10? ? y?: 1 ??1?n?1 ??kn x ? e H cos kn H 2 n?1 k

n

To get the e?ect of unbounded ?uid domain in the truncation surface, the coe?cient f is determined numerically. The signi?cance of the proposed truncation boundary condition of the in?nite ?uid domain is that the value of normal derivative of the pressure is not constant at all, though previous researchers [2,4] had assumed it as constant. 2.3. Finite element implementation: Here the ?uid domain is discretised as an assemblage of ?nite elements, assuming pressure to be the nodal unknown. The pressure at any point inside an element is given by p ? ?N ?fg; p ?11? where fg are the pressures at the element nodes and [N] is the interpolation p function. By the use of Galerkin process, the discretised form of Eq. (1) becomes ?G?fpg ? fBg; ?12? where, {p} represents the vector of nodal pressures for the ?uid domain. The elements of [G] and {B} are obtained by summing the contributions from each element as follows: XZ o o T o T o ?G? ? ?N ? ?N ? ? ?N ? ?N ? dX; ?13? ox oy oy X ox

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1319

fBg ?

XZ

C

?N ?T

op dC; on

?14?

X and C represent the ?uid domain and ?uid boundary respectively. Through Eq. (14) all the boundary conditions for the normal pressure gradients are incorporated. {B} may be separated into the following parts for di?erent boundaries: fBg ? fBf g ? fBs g ? fBb g ? fBt g: ?15? Here the subscripts f, s, b and t stand for the free surface, the structure–?uid interfaces, reservoir bed–?uid interface and the truncation surface respectively. According to the boundary conditions for the ?uid domain, pressure at free surface is zero and may be written in the following form: fBf g ? 0: ?16? At the structure–?uid interface, if {a} is the vector of nodal accelerations of generalised coordinates, Eqs. (3) and (14) are used to express {Bs} as fBs g ? ?q?S?fag; where, ?S? ? XZ

Cs

?17?

?N ? ?T ??N s ? dC:

T

?18?

In Eq. (18), [T] is the matrix which transforms the generalised accelerations of a point on the structure–?uid interface to the acceleration in the direction which is normal to the interface. Similarly, [Ns] is the matrix of shape functions used to interpolate the generalised acceleration of any point on the ?uid–structure interface in terms of the generalised nodal accelerations of an element. As the reservoir bed–?uid interface is considered as rigid, Eqs. (4) and (14) will become as follows: fBf g ? 0; f ?R?fpg; H ?19? Substituting Eq. (10) in Eq. (14) the following relation may be obtained: fBt g ? ? where, ?R? ? ?20?

XZ

Ct

?N ?T ?N ?dC:

?21?

Substitution of Eqs. (15)–(17), (19) and (20) in Eq. (12) leads to the following equation: ?Gt ?fpg ? ?q?S?fag; ?22?

1320

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

where, f ?R?: ?23? H For any prescribed acceleration of the ?uid–structure interface, Eq. (23) may be used to solve the hydrodynamic pressure in the reservoir. ?Gt ? ? ?G? ?

3. Numerical results Few examples have been solved in order to examine the feasibility and the accuracy of the proposed truncation boundary condition. The examples are related to the analysis of the hydrodynamic pressure distribution on two-dimensional structures exposed to in?nite ?uid media, such as shore protection walls, dams, etc. The structure in this particular case is considered as rigid and its direction of vibration at its base is horizontal and normal to its longitudinal direction. 3.1. Structure with vertical upstream face The hydrodynamic pressures at the vertical upstream face of the structure are obtained and compared with those obtained by classical [1] and other [2,4] methods. The in?nite reservoir is truncated at four di?erent locations (L/H) and di?erent truncation boundary conditions have been adopted for comparison. The ?nite reservoir obtained after truncation in each case is discretised with two di?erent sizes of the ?nite element mesh. The numbers of elements in vertical and horizontal directions are represented by Nv and Nh respectively. A typical ?nite element discretisation is shown in Fig. 2, considering four noded isoparametric quadrilateral elements. Here, L and H represent the distance of truncation boundary from the structure and the reservoir height respectively as shown in Fig. 2. All the di?erent cases analysed, are summarised in a tabular form along with the results for the maximum pressure coe?cient, c0 = p0/qaH at the bottom of the ?uid–structure interface. Here, p0 is the developed hydrodynamic pressure at the bottom of the ?uid–structure interface. Results are compared with Westergaard?s [1] classical solution. Table 1 shows the e?ciency and accuracy of the proposed truncation boundary. The proposed boundary is capable of producing reasonably accurate results even at a relatively very short length of the reservoir. It is interesting to note from Table 1 that the results converge to the exact solution as with the ?ner meshes and is practically independent of the ratio L/H, when proposed truncation boundary is adopted. Fig. 3 shows the distribution of hydrodynamic pressure coe?cient, c (the ratio of hydrodynamic pressure to hydrostatic pressure = p/qaH) along the upstream face of the structure for di?erent truncation boundary conditions.

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1321

y

H Reservoir ∞ Rigid Structure x

a

L

Fig. 2. A typical ?nite element mesh of an in?nite reservoir impounded by a structure having a vertical upstream face (L/H = 0.2, Nh = 4 and Nv = 20). Table 1 Comparison of hydrodynamic pressure coe?cient for di?erent truncation boundary conditions L/H Mesh size (Nv · Nh) 20 · 1 40 · 2 20 · 2 40 · 4 20 · 4 40 · 8 20 · 10 40 · 20 Sommerfeld boundary condition c0 0.05 0.1 0.2 0.5

a

Sharan (1985) c0 0.6812 0.6811 0.7071 0.7069 0.7294 0.7292 0.7417 0.7416 % error in c0a ?8.26 ?8.27 ?4.77 ?4.79 ?1.76 ?1.79 ?0.11 ?0.12

Proposed boundary condition c0 0.74272 0.74252 0.74272 0.74252 0.74270 0.74252 0.74255 0.74251 % error in c0a 0.030 0.003 0.030 0.003 0.027 0.003 0.007 0.001

% error in c0a 1248.5 1248.9 577.6 577.8 245.5 245.6 57.0 57.1

10.0125 10.0156 5.0313 5.0328 2.5656 2.5664 1.1659 1.1662

Compared with the value of c0 = 0.7425, computed by classical method for L/H = 1.

Similar conclusions may be drawn from Fig. 4, which shows the comparison of total hydrodynamic pressure coe?cient, ct (ratio of total hydrodynamic pressure to total hydrostatic pressure) for di?erent truncation boundary condition at the ?uid–structure interface. Fig. 5 shows the variation of the coe?cient f, derived from Eq. (10). It is interesting to note that this value has di?erent magnitude for di?erent position of the truncation boundary. Previous researchers had considered it as either zero or some constant. For example, the coe?cient of f is considered as p/2 by Sharan [4]. Hence signi?cant amount of error

1322

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

Fig. 3. Comparison of hydrodynamic pressure for di?erent truncation boundary along the upstream face of the structure.

develops while truncation surface is considered at a closer distance from the structural interface. 3.2. Structure with inclined upstream face The hydrodynamic pressure distribution along the reservoir–structure interface having di?erent inclination of the upstream face (Fig. 6) is presented in graphical form (Fig. 7). Results are presented for two di?erent ?nite element meshes: (i) L = 0.1H, Nv = 20, Nh = 2; (ii) L = 0.2H, Nv = 20, Nh = 4. The results are compared with the classical results obtained by Chwang [13] for different slopes of the ?uid–structure interface (h = 15°, 30°, 45°, 60° and 75°). It is interesting to note that as the proposed far-boundary condition is derived assuming the face of the structure to be vertical, the results for a vertical face are relatively more accurate than those for an inclined one. However, if proposed boundary is located at L = 0.2H, the results become almost identical to the exact solution. 3.3. Structure with partially inclined upstream face The example considered herein is for a partially inclined surface of the ?uid– structure interface for which, to the author?s knowledge, no analytical solution

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1323

Fig. 4. Convergence of total hydrodynamic pressure for di?erent truncation boundary conditions.

has yet been reported. The geometry of the ?uid–structure interface is shown in Fig. 8. The example is selected in order to compare the present results with the experimental results obtained by Zangar [15] who used electrical analogy method. The hydrodynamic pressure distribution at the ?uid–structure interface is shown in Fig. 9 for two di?erent ?nite element meshes: (i) L = 0.1H, Nv = 40, Nh = 4; (ii) L = 0.5H, Nv = 40, Nh = 20. Fig. 9 shows a comparison of the present results with those obtained by Zangar [15], in which minor discrepancy is observed. The di?erence in results may be attributed to di?erent reasons such as assumption in experiments, imperfections, etc. The results obtained by Zangar [15] for a fully inclined surface having = 60° is compared with the analytical solution [15] which is almost identical to the present results as observed from Figs. 7 and 9 and a similar discrepancy is observed. It may be mentioned that, similar discrepancy was also observed by Sharan [4], who pointed out that this small error may be due to Zangar?s [15] experimental modeling. 3.4. Structure with inclined reservoir bed The hydrodynamic pressure distribution along the upstream face of the structure with di?erent inclination of the reservoir bed is also investigated.

1324

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

Fig. 5. Coe?cient f for di?erent position of the truncation boundary along the height of the structure.

Fig. 6. A typical ?nite element discretisation of an in?nite reservoir impounded by a structure having an inclined upstream face. (L/H = 0.2, Nh = 4 and Nv = 20).

A typical ?nite element discretisation of the reservoir is shown in Fig. 10 in which Li represents the length of inclined bed. L represents the length between the structure and the truncation boundary and Lh is the length beyond the inclined part and up to the truncation boundary. a is the angle of inclination of the reservoir bed. Number of vertical elementsNv, assumed for the problem

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1325

Fig. 7. Comparison of pressure distribution along the upstream sloping faces of the structure.

Fig. 8. Finite element mesh of an in?nite reservoir having partially inclined upstream face (L/H = 0.1, Nh = 4 and Nv = 40).

are 10; number of horizontal elements Nhi, to represent the inclined bed are 10 and number of horizontal elements Nhh, to represent the horizontal bed of

1326

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

Fig. 9. Hydrodynamic pressure distribution on fully and partially inclined upstream faces of the structure.

length Lh are 3. Results are presented for di?erent inclination angle of the reservoir bed. The classical solutions of such problems for the above-mentioned con?guration is not available in the open literature. Therefore, the results are compared with those obtained adopting the Sommerfeld radiation boundary condition [2] and the boundary condition proposed by Sharan [4] for L/ H = 4.0. The in?nite reservoir is truncated at three di?erent locations and di?erent truncation boundaries have been adopted for comparison. All the di?erent cases analysed, are summarised in a tabular form (Table 2) along with the results for the maximum pressure coe?cient c0 at the bottom of the ?uid–structure interface. Table 2 shows the e?ciency and accuracy of the present truncation boundary condition. The results con?rm that the present far-boundary condition is also e?cient to account for the inclination of the reservoir bed,

Fig. 10. Finite element mesh having inclined reservoir bed.

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

1327

Table 2 Comparison of hydrodynamic pressure coe?cient (c0) for di?erent truncation boundary conditions, having inclined reservoir bed a Lh/H Nv · Nh · Nhh Sommerfeld boundary condition c0 +5 0.05 0.10 0.20 0.05 0.10 0.20 0.05 0.10 0.20 0.05 0.10 0.20 20 · 2 · 1 20 · 2 · 2 20 · 2 · 4 20 · 2 · 1 20 · 2 · 2 20 · 2 · 4 20 · 2 · 1 20 · 2 · 2 20 · 2 · 4 20 · 2 · 1 20 · 2 · 2 20 · 2 · 4 3.3857 2.5705 1.7727 3.3787 2.5614 1.7601 3.3902 2.5763 1.7804 3.3759 2.5577 1.7548 % error in c0a 347.55 239.79 134.33 363.22 251.16 141.30 339.20 233.76 130.65 370.59 256.62 144.68 Sharan (1985) Proposed boundary condition c0 0.7590 0.7579 0.7571 0.7275 0.7285 0.7292 0.7765 0.7745 0.7729 0.7131 0.7150 0.7165 % error in c0a 0.33 0.19 0.08 0.26 0.13 0.02 0.60 0.34 0.13 0.57 0.30 0.09

c0 0.7372 0.7445 0.7519 0.7063 0.7153 0.7240 0.7546 0.7611 0.7677 0.6921 0.7020 0.7113

% error in c0a 2.55 1.58 0.61 3.17 1.93 0.74 2.25 1.41 0.55 3.49 2.12 0.82

?5

+10

?10

Compared with the values of c0 = 0.7565, 0.7294, 0.7719 and 0.7172 having = +5°, ?5°, +10° and ?10° respectively, computed using Sommerfeld boundary condition for Lh/H = 4.0.

a

considering the truncation boundary at a relatively closer distance from the structure.

4. Conclusion A simple but e?ective numerical truncation boundary condition has been derived to model the in?nite reservoir into an equivalent ?nite reservoir model. The proposed truncation boundary may be located even at a very small distance away from the structure, resulting in great computational advantages. Incorporation of the present truncation boundary in the ?nite element program is quite simple and the additional computational e?ort required to include the proposed truncation boundary condition is insigni?cant.

References

[1] H.M. Westergaard, Water pressure on dams during earthquakes, Trans. ASCE 98 (1933) 418–472. [2] O.C. Zienkiewicz, R.E. Newton, Coupled vibration of a structure submerged in a compressible ?uid, in: International Symposium on Finite Element Techniques, Stuttgart, Germany, 1969.

1328

D. Maity / Appl. Math. Comput. 170 (2005) 1314–1328

[3] O.C. Zienkiewicz, P. Bettes, D.W. Kelly, The ?nite element method for determining ?uid loadings on rigid structures: two- and three-dimensional formulations, in: O.C. Zienkiewicz, R.W. Lewis, K.G. Stagg (Eds.), Numerical Methods in O?shore Engineering, Wiley, Chichester, UK, 1978 (Chapter 4). [4] S.K. Sharan, Finite element analysis of unbounded and incompressible ?uid domains, Int. J. Numer. Meth. Eng. 21 (1985) 1659–1669. [5] S.K. Sharan, Finite element modelling of in?nite reservoirs, J. Eng. Mech. ASCE 111 (1985) 1457–1469. [6] C.S. Tsai, Semi-analytical solution for hydrodynamic pressures on dams with arbitrary upstream face considering water compressibility, Comput. Struct. 42 (1992) 497–502. ? [7] X. Li, M.P. Romo, L.J. Aviles, Finite element analysis of dam–reservoir systems using an exact far-boundary condition, Comput. Struct. 60 (1996) 751–762. [8] D. Maity, S.K. Bhattacharyya, Finite element analysis of ?uid structure system for small ?uid displacement, Int. J. Struct. 17 (1) (1997) 1–18. [9] S.S. Saini, P. Bettess, O.C. Zienkiewicz, Coupled hydrodynamic response of concrete gravity dams using ?nite and in?nite elements, Earthquake Eng. Struct. Dyn. 6 (1978) 363–374. [10] C.A. Felippa, Interfacing ?nite element and boundary discretisations, Appl. Math. Model 5 (1981) 383–386. [11] C.S. Tsai, G.C. Lee, Arch dam-?uid interactions: by FEM–BEM and substructure concept, Int. J. Numer. Meth. Eng. 24 (1987) 2367–2388. [12] R. Yang, C.S. Tsai, G.C. Lee, Explicit time-domain transmitting boundary for dam–reservoir interaction analysis, Int. J. Numer. Meth. Eng. 36 (1993) 1789–1804. [13] A.T. Chwang, Hydrodynamic pressure on sloping dams during earthquakes. Part-2. Exact theory, J. Fluid Mech. 87 (1978) 343–348. [14] D. Maity, S.K. Bhattacharyya, Time domain analysis of in?nite reservoir by ?nite element method using a novel far-boundary condition, Int. J. Finite Element Anal. Des. 32 (1999) 85–96. [15] C.N. Zangar, Hydrodynamic pressures on dams due to horizontal earthquakes, Proc. Soc. Ex. Stress Anal. 10 (1953) 93–102.