INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL.

21, 1659-1669 (1985)

FINITE ELEMENT ANALYSIS OF UNBOUNDED AND

INCOMPRESSIBLE FLUID DOMAINS

SHAILENDRA K. SHARAN’

School of Engineering, Laurentiun University, Sudbury, Canadu

SUMMARY
A finite element technique of analysing the hydrodynamic pressure, resulting from the vibration of a structure
submerged in an unbounded fluid domain, is presented. A suitable boundary condition is proposed for the
surface where the unbounded fluid domain is truncated. Pressure is assumed to be the nodal unknown and the
fluid is treated as being incompressible and inviscid. The proposed method is computationally very efficient
and its implementation in a finite element program is quite straightforward. The efficiency of the method is
demonstrated by analysing some two-dimensional problems.

INTRODUCTION
For an accurate analysis of fluid loading on submerged structures, such as dams and offshore
structures, having irregular geometries, the fluid domain is usually treated as an assemblage of
finite elements.’-3 In many problems, the fluid may be assumed to be incompressible and the fluid
loading may be considered as an ‘added mass’ in the dynamic response analysis of structures4
In the finite element analysis, difficulties arise as a result of the large extent of the fluid domain
which, for many problems, is practically unbounded. In such a case, the infinite domain must be
truncated at a certain distance away from the structure (Figure 1). Accuracy in results may be
improved by truncating the domain at a larger distance. However, this results in an increased cost
of computation. In addition, if the effect of fluid structure interacton is included in the analysis, the
cost of computation could become prohibitive.
Several techniques have been suggested in the past to circumvent this computational difficulty.
These may be broadly classified as, (a) imposition of a boundary condition along the truncation
~ (b) coupling the finite element discretization with other types of discretization such
~ u r f a c e , ”and
as ‘infinite’ elements,6 boundary elernent~;~ or with continuum solutions.8 The first 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.’ However, for an incompressible fluid, this or other similar boundary
condition5takes the form of that for a rigid stationary boundary. The behaviour of the fluid motion
at the truncation boundary is, therefore, not truly represented. As a result, a very large extent of the
fluid domain is required to be included in the analysis.
The objective of the present paper is to propose a suitable boundary condition along the
truncation surface in the finite element analysis of unbounded and incompressible fluid domains.

‘Assistant Professor

0029-5981/85/091659-11$01.10 Received 3 M a y 1984

0 1985 by John Wiley & Sons, Ltd. Revised 6 November 1984
1660 S. K. SHARAN

STRUCTURE---\
\\ FREE SURFACE
-
1
I I \' 1
V

FLUID I

Figure 1. A structure submerged in an unbounded fluid domain

EQUATIONS OF FLUID MOTION

If the fluid is assumed to be incompressible and inviscid, the hydrodynamic pressure p (in excess of
the hydrostatic pressure) resulting from the vibration of a submerged structure satisfies the Laplace
equation'

v2p = 0 (1)
If the effects of surface waves and viscosity of the fluid are neglected, the following boundary
conditions apply:
1. At the free surface
p=0 (24
(2) At the fluid-solid interface
JP
-= - P*,
an
where p is the density of the fluid, an is the outwardly-directed normal to the elemental surface
along the interface, and a, is the acceleration of the solid surface in the direction of an.
In order to consider the effect of radiation damping, it is assumed that at an infinitely large
distance away from the structure, equation (2a) is satisfied. If the unbounded fluid domain is
truncated at a sufficiently large distance from the region of interest, Sommerfeld radiation
boundary3 is normally used at the truncation surface. For an incompressible fluid this takes the
following form:

This boundary condition at the truncation surface is the same as that for a rigid stationary
boundary.

PROPOSED TRUNCATION BOUNDARY

The proposed boundary condition along the truncation surface is derived on the basis of the
following additional assumptions:
 UNBOUNDED AND INCOMPRESSIBLE FLUID DOMAINS 1661

't TRUNCATION SURFACE

P=O
'
Y=h

aP/ax V2P=0 p= 0
"-Pan

I1
II
11

*
I! *
Figure 2. Assumptions made for the derivation of a boundary condition along the truncation surface

1. The fluid domain extends to infinity and its motion is two dimensional.
2. The fluid-structure interface is vertical.
3. The submerged structure is rigid; its height is not less than the depth of the fluid; and it vibrates
in the direction normal to the fluid-structure interface.
4. The bottom of the fluid domain is rigid and horizontal.
Use of symmetry leads to the boundary value problem shown in Figure 2. The solution for the
hydrodynamic pressure at any point ( x , y ) is obtained as

(3)
where
(2n - 1)7t
An =
2
and h is the depth of the fluid domain.
For sufficiently large value of x/h, the second and subsequent terms in the summation series of
equation (3) may be neglected and p may be approximated as

(4)

Differentiation of the above equation partially, with respect to x, gives

Use of equation (4) in equation (5) leads to the following relationship between the pressure and its
gradient along the direction which is normal to the truncation surface (Figure 2):
1662 S. K. SHARAN

It should be noted that if the truncation boundary is located at a relatively large distance away
from the structure, the value of p is approximately equal to zero and equation (6) is reduced to
equation (2c). However, this is not true if the fluid domain is truncated at a relatively short distance
away from the structure.

FINITE ELEMENT IMPLEMENTATION

The fluid domain is discretized as an assemblage of finite elements, assuming pressure to be the
nodal unknown. The pressure at any point inside an element is given by
P = “1 {P>“ (7)
where { p } “ is the vector of pressures at the element nodes and [N] is matrix of shape functions.
By the use of the Galerkin process, the discretized form of equation (1) is obtained as9
[GI {PI = { B ) (8)
where {p} represents the vector of nodal pressures for the fluid domain. The elements of [GI and
{ B } are obtained by summing the contributions from each element as follows:

[GI =q( (&)“l’,x“l

a a
+-“IT-“
aY
(3
aY
a
aZ
+-“1’-“1 aZ” dv (9)

In the above equations v, and s, represent, respectively, the volume and the boundary surface of an
element.
Through equation (10) all the boundary conditions for the normal pressure gradients are
incorporated. { B } may be separated into the following parts for different boundaries:
{B}= P, } + {Bs}+ {BJ (1 1)
where the subscripts f , s and E stand for the free surface, the fluid-solid interfaces and the
truncation surface, respectively. For the present paper, fluid-solid interfaces represent both fluid-
structure and fluid-foundation interfaces.
For the nodal points at the free surface, the elements of { B , } are unknown but the pressure is
assumed to be zero. Hence, for the purpose of obtaining the solutions for the remaining unknown
nodal pressures, the nodal points at the free surface may be deleted from equation (8).In such a case
(4)= lo> (1 2)
For the fluid-solid interface, if { u } is the vector of nodal accelerations along the directions of
generalized co-ordinates, equations (2b) and (10) are used to express { B , } as
(13)
where

(13a)

In the above equation, [R] is the matrix which transforms the generalized accelerations of a point
on the fluid-solid interface to the acceleration in the direction which is normal to the interface.
Similarly, [N,] is the matrix of shape functions used to interpolate the generalized acceleration of
 UNBOUNDED AND INCOMPRESSIBLE FLUID DOMAINS 1663

any point on the fluid-solid interface in terms of the generalized nodal accelerations of an element.
If the motion of the fluid domain is to be coupled with the finite element model of a structure, [ N , ]
represents the matrix of shape functions of the structural element at the fluid-structure interface.
The integrations and summation in equation (13a) are performed only for the element surfaces on
the fluid-solid interface.
The truncation surface is assumed to be a plane which is normal to the free surface. Substitution
of the proposed equation (6) in equation (10) leads to

where

[TI = C
Lt [NITIN]ds

The integrations and summation in equation (14a) are performed only for the element surfaces on
(1 4 4

the truncation boundary. [TI is symmetrical and has very few non-zero terms, depending upon the
number of nodal points on the truncation surface.
Substitution of equations (ll), (12), (13) and (14) in equation (8) leads to the following equation:
[GI {PI = - Pcsl{a> (1 5 )
where
71
[GI = CCl + 5 c7-I (1 5 4

and h is the depth of fluid at the truncation boundary.

If the fluid domain is finite and there is no truncation boundary, { B , ) = ( 0 ) and hence
[GI = [GI (15b)
For any prescribed acceleration of the fluid-solid interface, equation (15) may be directly used to
solve for the hydrodynamic pressure. Alternatively, if the effect of fluid-structure interaction is to
be included in the analysis, incorporation of equation (1 5) in the discretized equation of motion of
the structure is straightf~rward.~
It should be noted that [GI is symmetrical and highly banded and an identical form is
maintained for [C,] as well. Furthermore, the additional computational effort required for the
inclusion of the proposed boundary condition at the truncation surface is negligible.

TEST EXAMPLES
The proposed truncation boundary condition was incorporated in a computer program for the
two-dimensional finite element analysis of fluid domains using linear triangular elements. In order
to test the accuracy, convergence and efficiency of the proposed method, some example problems
were solved. These examples were related to the analysis of hydrodynamic pressure distribution on
two-dimensional submerged structures, such as dams, subjected to a uniform acceleration a,. The
following assumptions were made for all the cases analysed: the structure is rigid; water in the
reservoir is incompressible; the reservoir bed is rigid, horizontal, and extends to infinity; the
direction of vibration of the structure is horizontal and normal to its longitudinal direction. These
assumptions were made in order to compare the present results with those obtained by
classical'0*" or other" methods. The geometry of the infinite reservoir was varied by varying the
1664 S. K. SHARAN

shape of the fluid-structure interface. These different example problems were analysed by
assuming the interface to be vertical, fully inclined, or partially inclined.

Example 1
In the first example, the hydrodynamic pressure on a dam, having a vertical upstream face, was
obtained. The infinite reservoir was analysed for four different locations of the truncation
boundary, each resulting in a different size 1 of an equivalent finite reservoir. Each size of the finite
reservoir was then discretized in two different sizes of the finite element mesh. Results for pressure
distribution were obtained for two types of boundary condition along the truncation surface:
(a) a p / d n = 0, i.e. a rigid stationary boundary, and (b) the proposed boundary. Figure 3 shows a
typical finite element model of the infinite reservoir.
All the different cases analysed are summarized in Table I along with the results for the pressure
coefficient, C, = po/pa,h, at the bottom of the dam face. For this example, C, is the maximum

Y
4

-a

RESERVOIR

cc X
a0

Figure 3. A typical finite element model of an infinite reservoir impounded by a dam having a vertical upstream face
(l/h = 0.2, N , = 10 and N , = 2)

Table I. Convergence test for the analysis of hydrodynamic pressure on a dam having a vertical upstream face

d p p n = 0 at x = 1 Proposed boundary at x = 1
Mesh size Percentage Percentage
Ilh N" x N h c, error+ in C, c, error+in C,

lox 1 5.0250 577 07079 4.66

0.1
20 x 2 5.03 12 578 0.7072 4.75
10 x 2 2.5625 245 0.7298 1.71
0.2
20 x 4 2.5656 246 0.7294 1.76
10x5 1.1645 57 0.7410 0.20
0.5
2 0 x 10 1.1658 57 0.74 15 0.13
10 x 10 0.8 142 10 0.7410 0.20
1.o
20 x 20 0.8153 10 0.742 1 0.05

Compared with the value of C , = 0.7425, computed by classical method for l/h = co.
 UNBOUNDED AND INCOMPRESSIBLE FLUID DOMAINS 166.5

pressure coefficient. Results for C,, thus obtained, were compared with Westergaard's'" classical
solution and are presented in Table I. It should be noted that results obtained by both types of
truncation boundary approach the exact solution as the distance of the truncation boundary from
the dam is increased and as the finite element mesh is further subdivided. However, in the casc of the
proposed boundary, a relatively very short length of the reservoir is required to be considered. For
example, the proposed boundary produced reasonably accurate results (i.e. an error of about 5 per
cent) for 1 = 0.1 h, whereas for the other boundary, such a level of accuracy was not achieved
(resulting in an error of about 10 per cent) even for I = h.
Similar conclusions were drawn from a comparison of the total hydrodynamic pressure
coefficient, C, = total hydrodynamic pressure/total hydrostatic pressure on the face of the dam.
Figure 4 shows the results for C, obtained for different locations of the truncation boundary and
are compared with the exact result. It should be noted that the result for 1 = 0.1 h with the proposed
boundary has the same order of accuracy as that for 1 = 2h with the rigid truncation boundary.

Example 2
This example was similar to the previous one, with the only difference that the upstream face of
the dam was considered to be inclined. Hydrodynamic pressure distributions were obtained for
different angles of inclination 8 of the dam face. Five different slopes were considered by varying 0
from 15" to 7.5". For each slope, different finite element models, similar to the previous example,
were analysed. Figure 5 shows a typical finite element model. Results of C, computed for a typical
value of 8= 30" are summarized in Table I1 and are compared with the exact solution obtained by
Chwang.' Again, results obtained by using either truncation boundary condition converged to the

7 -
? - CLASSICAL (L.00)
I
6 - 1 -+- F.E.M. (ap/an = 0 AT x =L)
I
I x EE.M. (PROPOSED BOUNDARY AT x = L )
I
5 - I
I
u- I
I
4- \

3 -
a \
\
\
\
2-
\
\
&---
- \

---+------ *
1-
I I

0 ; I 1 I 1
0.0 0.5 1.0 1.5 2 .o
L/h
Figure 4. Total hydrodynamic pressure coeficients for a dam having a vertical upstream face
1666 S. K. SHARAN

--m

RIGID DAM

Figure 5. A typical finite element model of an infinite reservoir impounded by a dam having a fully inclined upstream
face (0 = 30",l/h = 0.2, N , = 10 and N , = 2).

Table 11. Convergence test for the analysis of hydrodynamic pressure on a dam having an inclined upstream
face (angle of inclination = 30")

a p p n = 0 at x = 1 Proposed boundary at x = 1
Mesh size Percentage Percentage
llh N" x N h Co error+in C, c, error+in C ,

lox 1 1.1431 145 0.482 1 3.15

0.1
20 x 2 1.1331 142 0.4765 1.95
10 x 2 0.8993 92 0.4728 1.16
0.2
20 x 4 0.8953 92 0-4705 0.66
10x5 0.6114 31 0.4683 0.19
0.5
2 0 x 10 0.6 104 31 0.4611 0.06
10 x 10 0.4957 6 0.4677 0.06
1.o
20 x 20 0.4954 6 0.4675 0.02

'Compared with the value of C , = 0.4674, computed by classical method for l/h = m.

exact solution as the ratio l / h was increased and as the finite element mesh was further subdivided.
However, as in the previous example, results for 1 = 0.1 h with the proposed boundary were more
accurate (having an error of about 2 per cent) than those for 1 = h with the rigid truncation
boundary (having an error of about 6 per cent).
It is interesting to note that although the proposed boundary condition is derived assuming the
face of the structure to be vertical, the results for an inclined face are relatively more accurate than
those for a vertical face. This demonstrates that the proposed boundary condition is equally valid
for inclined faces.
Figure 6 shows the distribution of pressure, obtained by the proposed method, for different
values of 8. Results are presented for two different finite element models: (a) 1 = 0 2 h, number of
elements in the vertical direction N , = 20, and number of elements in the horizontal direction
N , = 4;(b) 1 = h, N , = 20 and N , = 20. These results are compared with the exact results obtained
by Chwang. If the proposed truncation boundary is located at 1 = h, the finite element results
are almost identical to the exact results. The agreement between the results is found to be excellent
even for 1 = 0-2h.
 UNBOUNDED A N D INCOMPRESSIBLE FLUID DOMAINS 1667

--?---

h
+@I

-f-- RESERVOIR

“i‘
Figure 7. A typical finite element model of an infinite reservoir impounded by a dam having a partially inclined upstream
face (l/h = 0.2, N , = 20 and N , = 4)

Example 3
In both the previous examples, it was possible to compare the present results with exact
solutions. The third example was for a partially inclined surface of the dam-reservoir interface for
which, to the author’s knowledge, no analytical solution has yet been reported. This example was
selected in order to compare the present results with the experimental results obtained by
1668 S. K. SHARAN

1.o

0.8

0.6

r
a
\

0.4

0.2

0
0 02 0.2 0.3 0.4 0.5 0.6 0.7
p/Pa,h

Figure 8. Distribution of hydrodynamic pressure on fully inclined and partially inclined upstream faces of dams

Zangar,” who used the method of electrical analogy. The surface selected is shown in Figure 7
along with a typical finite element model for the reservoir.
Figure 8 shows the distribution of pressure obtained for two different finite element models: (a)
1 = 0.2 h, N , = 20 and N , = 4; (b) 1 = h, N , = 20 and N , = 20. Both of these models were analysed by
using the proposed truncation boundary. It should be noted, as in the previous examples, that with
the use of the proposed method, the hydrodynamic pressure on dams could be computed very
accurately by considering the length of the reservoir to be as small as equal to 0.2h. Figure 8 also
shows a comparison of the present results with those obtained by Zangar.” Although the
agreement is very good, it was found that Zangar’s results could not be reproduced, even by
increasing the length of the reservoir and by further subdivision of the finite element mesh. In order
to explain the cause of the minor discrepancy, Zangar’s result for a fully inclined surface having
H = 60” was compared with the exact solution,” which was almost identical to the present result
(Figures 6 and 8). From such a comparison, it was concluded that the insignificant error could be in
Zangar’s experimental modelling.

CONCLUSIONS
A finite element technique has been proposed for the analysis of hydrodynamic pressure generated
due to the vibration of a submerged structure in an infinite and incompressible fluid domain.
In the proposed method, a boundary condition has been derived to represent the surface where
the infinite fluid domain is truncated to obtain an equivalent finite domain model.
 UNBOUNDED A N D INCOMPRESSIBLE F L U I D DOMAINS 1669

The effectiveness of the proposed method has been demonstrated by analysing hydrodynamic
pressures on rigid dams of varying surfaces. It was found that the results converged to the exact
solutions as the lengths of the equivalent finite reservoirs were increased and as the finite
element meshes were further subdivided.
There are two main merits of the proposed technique. First, the truncation boundary may be
located very close to the structure, resulting in great computational advantages. For example,
for the dam problem analysed, the length of the equivalent finite reservoir to be considered was
found to be relatively very small as compared to that required in the conventional method in
which the truncation boundary is assumed to have a zero normal pressure gradient. Secondly,
incorporation of the proposed technique in a finite element program is quite straightforward
and the additional computational effort required to include the proposed truncation boundary
condition is insignificant.

ACKNOWLEDGEMENT

The financial support provided by the Natural Sciences and Engineering Research Council of
Canada is gratefully acknowledged.

REFERENCES
1. 0.C. Zienkiewicz and R. E. Newton, ‘Coupled vibrations of a structure submerged in a compressible fluid’, Int. Symp.
on Finite Element Techniques, Stuttgart, 1969.
2. S. K. Sharan and G. M. L. Gladwell, ‘A general method for the dynamic response analysis of fluid -structure systems’,
Comp. Struct. (in press).
3. 0. C. Zienkiewicz, P. Bettess and D. W. Kelly, ‘The finite element method for determining fluid loadings on rigid
structures: two- and three-dimensional formulations’, Chapter 4 of Numerical Methods in Offshore Engineering ( 0 . C .
Zienkiewicz, R. W. Lewis and K. G. Stagg, Eds.), Wiley, Chichester, U.K., 1978.
4. 0.C. Zienkiewicz and P. Bettess, ‘Dynamic fluid-structure interaction. Numerical modelling of the coupled problem’,
ChaDter 5 of Numerical Methods in Offshore Enqineerina (0.C. Zienkiewicz, R. W. Lewis and K. G. Stam, -_ Eds.), Wiley,
ChiEhester, U.K., 1978.
5. J. Humer and M. Roufaiel. ‘Finite element analvsis of reservoir vibration’. Proc. A.S.C.E.. 109 (EM1).
, I _
215-230 (1983).
, ,
6. S. S. Saini, P. Bettess and 0.C. Zienkiewicz,‘Co;pled hydrodynamic response ofconcrete gravity dams using finite and
infinite elements’, Earthquake Engng Struct. Dyn., 6, 363-374 (1978).
7. C. A. Felippa, ‘Interfacing finite element and boundary element discretizations’, Appl. Math. Modell., 5,383-386 (I98 1).
8. J. F. Hall and A. K. Chopra, ‘Dynamic analysis of arch dams including hydrodynamic effects’, Proc. A.S.C.E.,
109(EM1), 149-167 (1983).
9. 0. C. Zienkiewicz, The Finite Element Method, 3rd edn, McGraw-Hill, London, 1977.
10. H. M. Westergaard, ‘Water pressure on dams during earthquakes’, Trans. A.S.C.E., 98, 418-472 (1933).
11. A. T. Chwang, ‘Hydrodynamic pressures on sloping dams during earthquakes. Part 2. Exact theory’, J . Fluid Mech.,
87(2), 343-348 (1978).
12. C. N. Zangar, ‘Hydrodynamic pressures on dams due to horizontal earthquakes’, Proc. Soc. E x . Stress Anal. 10(2),
93-102 (1953).