- Email: [email protected]
A direct boundary-element method for three-D wave diffraction and radiation problems
Ocean Engng, Vol. 15, No 5, pp. 431-455, 1988. Printed in Great Britain.
0029-8018/88 $3.00 + .00 Pergamon Press pie
A DIRECT BOUNDARY-ELEMENT METHOD FOR THREE-D WAVE DIFFRACTION AND RADIATION PROBLEMS CHIH-KANG LEE* *Veritas Technical Services Inc., Houston, TX, U.S.A.
and JACK Y. K. Lout tOcean Engineering Program, Texas A&M University, College Station, TX77843, U.S.A. Abstract--Three-dimensional analytical procedures are developed for determining the waveexciting and motion-induced hydrodynamic forces for fixed and floating ocean structures in a fluid of finite or infinite depth. The fluid potential is expressed by means of a Helmholtz integral, which involves the normal fluid velocities at the fluid-structure interface. Upon discretizing the integral equation by using constant or linear elements and by enforcing the velocity continuity condition at the discrete points of the interface, the fluid potential can be determined. A computer program based on the direct boundary-element method is developed to calculate the wave-exciting and motion-induced hydrodynamic forces. Effects of structural flexibility are not considered in the present paper. Comparisons of the results with those obtained by several previous investigators reveal good agreements.
INTRODUCTION THE TOPIC of fluid loading of ships and large structures has received considerable attention in recent years. For design applications one must assess, with a high degree of certainty, the magnitude of the load induced by the extreme waves associated with the design storm conditions and the load imposed by the motion of the structure. Thus, even though the numerical procedures in free-surface hydrodynamics tend to be rather involved, practical needs exist for dealing with the hydrodynamics of large structures. Of particular concern has been the development of the wave-exciting force and the socalled added mass and damping coefficients. Wave-body interaction at infinitesimal amplitudes is a richly developed subject that, using the so-called indirect boundary-element method (or source distribution method), has served the field of naval architecture with great success. Its applications in offshore technology have stimulated further advances in recent years. Ursell (1949a,b) studied the hydrodynamics of two-dimensional semi-circular or eliptical cylinders in an infinite fluid. He introduced nonorthogonal expansions (multiples) whose terms satisfy the Laplace's equation and the free surface condition, but vanish at infinity. A source or dipole is then added to satisfy the radiation condition. The expansion coefficients and strengths of the singular sources were determined by imposing the body boundary condition at a finite number of collocation points. Frank (1967) used an integral equation approach and the strip method for solving ship wave problems. He assumed that the velocity potential could be represented by a distribution of two-dimensional sources over the contour of a ship segment. The 431
-~32
CIIlH-KANGL~-. and JACKY. K. Lot
unknown source strength is determined by satisfying the boundary condition at the midpoint of the segment. Frank's method can be applied to a two-dimensional section of arbitrary shape. However, it appears difficult to generalize the approach for threedimensional bodies. For the infinite fluid problem in three dimensions, Lewis (1929) computed exactly the added mass of a vibrating spheroid by the method of separation of variables in spheroidal coordinates. Hess and Smith (1962) solved the three-dimensional infinite fluid problem for a body of arbitrary shape in steady motion by solving the Fredholm integral equation of the second kind. For the three-dimensional free surface problem, Havelock (1955) used the multipole expansion method to solve the problem of a heaving sphere. As in the two-dimensional expansion method, he introduced a velocity potential represented by a set of higherorder singularities that satisfy the free surface condition and radiation condition trivially, and a wave source or multipole at the origin within the body that satisfies the radiation condition. The unknown coefficients for the polynomials and the wave sources or multipoles are determined from the kinematic boundary condition on the interface. The linear algebraic equations for the unknown coefficients are solved numerically. An integral equation formulation for calculating wave forces on arbitrarily shaped large bodies was developed by John (1950). The boundary-value problem is posed in terms of velocity potential functions, which are equivalent to a distribution of singularities on suitable boundaries of the fluid region, with the strengths of the singularities determined by the boundary conditions. Solution of the resulting integral equations is approximated by choosing a discrete set of points to obtain a set of matrix equations. The method has received widespread usage in free surface flow problems and various forms of the Green's function for wave potentials can be found in the work of Wehausen and Laitone (1960). Their research has been extended and boundary integral equations have now been frequently used to solve wave-diffraction and radiation problems; for instance, in the work of Garrison (1974, 1978), Hogben and Standing (1974), Faltinsen and Michelsen (1974). In this paper, an approach using a direct boundary-element method is developed which has more physical relation to the wave-body interaction problem. Like the finiteelement method, the direct boundary-element method also divides the external surface of a domain into a series of elements over which the functions under consideration can vary in different ways. This capability is important as the boundary integral equation formulations generally were restricted to constant source strength over the elements and the sources were assumed to be concentrated at a series of points on the external surface of the boundary. FORMULATION OF THE BOUNDARY-VALUE PROBLEMS A right-handed Cartesian coordinate system with the origin on the free surface is used as shown in Fig. 1. The angle o~ is the wave-heading angle with respect to the positive x-axis, the xy-plane coincides with the free surface when the fluid is at rest and the z-axis is taken directly opposite to the force of gravity. The structure is considered as a large rigid body oscillating sinusoidally in a long crested regular wave. The amplitudes of the motions of the structure, as well as of the wave, are assumed to be small. Since the fluid is assumed to be ideal and irrotational, there exists a velocity potential that governs the kinematics of the fluid and the velocity
3-D Wave diffraction
433
Z
v
I
s~
s3
I I I L
~'-~l~#~'~I~'#'~'l ss FIG. 1.
s~ 3-D Coordinate system and various boundaries.
V may be represented by the gradient of the potential function 4. Furthermore, we assume that the motion is sinusoidal in time. The total velocity potential is written as , ( x , y , z ; t ) = qJ~ + Oa = R (dO(x,y,z) e ''t }
(1)
where 4, and ~d are the r~diation potential and the potential for incident and diffracted waves respectively; 1:1 represents the real part of a complex quality and d~(x,y,z) is the spacial component of 4. Since a free floating body can move in a wave with six degreesof-freedom, the velocity potential is conveniently broken into components relating to the incident and diffracted wave fields and the motions of the body 6
+(x,y,z) = +,, + +~ + ~] +~, ~=!
(2)
434
()IIHt-KANG Ll~t~ and JACk. Y. K, Lot,
where +o is the complex amplitude of the incident wave potential, d)v the diffracted wave potential and ~bj the velocity potential due to the forced motion of the body in the jth mode; j is taken to be 1. . . . 6 for surge, sway, heave, roll, pitch and yaw. as shown in Fig. 1. Radiation p r o b l e m
Consider a rigid body oscillating with a frequency to on or beneath the free surface of a given fluid, which is otherwise at rest. The velocities of the body are described by U~ = F~ { - i t o ~ e
i~"}=~{Ui,,e-i"
j = l . . . . 6,
}
(3)
where ~i is the complex amplitude of the motion of the body in the jth mode. Due to the forced motion of the body in the jth mode, the complex velocity potential ¢b~(x,y,z) must satisfy the Laplace equation, 7 2 +~ = 0
in the fluid region, ~ ,
(4)
together with the following boundary conditions hdm
"~J On 0+j
ito~n i
v ¢b~ = 0
j = 1. . . . 6
on z = 0
on the body surface Sl,
at the free surface S~,
Oz
on rigid wall $3 and sea bottom $6,
On
lim
0~
Rr O+j - iv+j \ Or
. . . . . v d~j = 0 Oz
)
= 0
at the distant boundary S4,
(7)
(8)
on the imaginary boundary
S_s in case of infinite water depth,
o) 2
(6)
-
0+j = 0
~r--+a:
(5)
(9)
where v = - - ; r = (x 2 + y2)½; R is the wave number which satisfies the dispersion g relation, K tanh k d = v; S~ . . . . $6 are the categorized boundaries of the domain, as shown in Fig. 1 and the derivative in equation (5) is taken with respect to the outward normal n of the body surface; n~-n6 are the generalized direction cosines on Sl. The set of equations (4-9) constitutes a boundary-value problem. To solve for the radiated wave potential for each mode of motion, a boundary-element method can be
3-D Wave
diffraction
435
employed. The total hydrodynamic pressure can be obtained from the linearized Bernoulli's equation (Newman 1978),
I
p=-p~-=-pFl
- i ~ o e -'°''
~+,
.
(10)
i=l
The oscillating hydrohynamic forces (i=1,2,3) and moments (i=4,5,6) acting on the body surface can be determined by integrating the fluid pressure (10) over the wetted surface S~ and the result can be written as 6
~=-
~(a~+bo~).
(11)
j=l
The coefficients aij and b 0- are known as the added mass and damping coefficient respectively.
Diffraction problem The velocity potential for the diffracted waves is computed as if the structure were restrained from moving. Similar to the radiated wave potential, the diffracted wave potential should satisfy the Laplace's equation (4) and the boundary conditions (5-9). Since the body is assumed to be restrained, the boundary condition on the body surface S~ becomes
0+7 _ O~
a+,, O~
(12)
'
where the incident wave potential is given by
d?o -
i gtoa,, cOShcoshk(Zkh + h) eitK,x+K2y),
(13)
K1 = g COSOt, I(.2 = K sinc~ and a., K and to are the amplitude, wave number and frequency of the incident wave, respectively. Thus we have a set of equations that constitutes a wave diffraction boundary-value problem. Again, a boundary-element method can be employed to solve the diffracted wave potentials and the diffracted pressure can be obtained from the Bernoulli's equation. The resulting wave-exciting forces and moments, Fj, j = 1 . . . . 6, can be written as
14, I
436
CHIH-KANG LEE
and JACK Y.
K. L()L~
BOUNDARY-ELEMENT METHOD Finite difference, finite element (domain methods), and the boundary-element (boundary method) methods are frequently used to solve the above differential equations numerically. The difference between the domain methods and the boundary methods is significant. The domain methods discretize the domain, while the boundary methods discretize the boundary. Thus, the boundary method reduces the dimension by one. This method is well suited to problems in which the limits of the domain are infinite or difficult to define, in that the solution is applied to the boundary rather than the domain. There are two types of boundary element methods, the "indirect" and the "direct". In the "indirect boundary-element" method, which is the "source" method developed by Kellogg for potential problems, the unknowns are not the physical variables of the problem. On the other hand, the "direct boundary-element" method not only allows us to solve the problem in terms of physical variables but also serves as the first step towards a better understanding of the technique and its relationship to other approximate methods, in particular, mixed finite-element methods. The direct boundary-element method involves the transformation of the governing differential equation within the domain under consideration to an integral equation defined on the surface of the domain. The surface then may be discretized into a number of elements over which a polynomial form of the solution is assumed. This enables the evaluation of the relevant integrals, usually by some numerical process, resulting in a system of linear algebraic equations. As only the surface of the domain needs to be discretized, the resulting systems of equations are considerably smaller than those involving domain type solutions, e.g. finite elements, finite differences, and considerable savings also can be achieved in the time required to prepare data. The key to this method, in potential problems, is the adoption of a fundamental solution used to eliminate the domain integral in the formulation. Such solutions can be found in Brebbia (1981).
Formulation of boundary integral equations The direct boundary-element method used to solve boundary-value problems can be formulated rigorously, using either an approach based on Green's theorem or a particular case of the weighted residual method. With the advantages of being more powerful and the ability to relate boundary solutions to other more classical engineering techniques, the weighted residual method is preferred. For the radiation and diffraction problems, which involve more general impedance boundary conditions of the Robin type, we can write the following error functions; for the radiation problem ~ = V2 +j a+j e~ - On
~2 = --
0+~ _ O~F/
in ~,
Ui° ni
on S~,
+~
on $2,
v
on
$3,
(15)
3-D Wave diffraction
0~
¢4 = ff~n - i v ~b~
0,~
~-s = ff~n + v d~
on 84, on S.s,
0+j ~-6 :
~
$6.
on
On
437
We then minimize the errors (15) by distributing a weighting function tb*. The resulting weighted residual statement can be written as
fn(V%3~*dt~=fsl(O~/-U~(a~)'b*dSfs~(O~J-v~J) On
+ ( O*i,'dS+ Js3 On
k On
~ (O*'-iv,O,*dS+
Js 4 k On
fs ( O * Y + v * , ) * * d S ~ k On
+ fs~** dS.
(16)
Integrating by parts the terms on the left-hand side of equation (16) twice will yield,
-
+
+f~4i~+*dS-~.~*~+*dS-~+,O--~ndS.
(17)
.
We now take the simple fundamental solution of V:~b* + ~i = 0 as the weighting function, i.e.
4" -
1 4~rr
(18) '
where r = X/(x - xi) 2 + (y - yi) 2 + (z - zi)2; (x, y, z) and (Xi, Yi, Zi) are the coordinates of the field and source points respectively. Substituting equation (18) into equation (17) and rearranging, we obtain
~ = -
On f,l \l, '.°**
)
Uj,, nj dp* dS
438
(hm,:-K,~.~G LEE and JACK Y, K. Lou
_i,,~(~,~,~:~_,~,~,~,,/dS _
0n
_f,,(,~e,:t~s_f,g,.o~:_ ~an a~ /
.
i v d~ +*) dS
~,
_f~.~(%o+* . an + v +~+*) dS
({+.~:~ On.]dS"
(19)
-L. ~ '
Equation (19) implies an integral relationship between the potential at the field point i inside the domain l~ and its values on the boundary S of the domain. When this point is moved to the boundary S, the integrals involved in equation (19) become singular at i and must be evaluated in the Cauchy Principle Value sense (Brebbia, 1978). Thus, for any point i on the boundary surface, the boundary integral equation can be written as
c'~=-
i,.(, ~-o~~* -u~,,~,~* t d S - i,(~ ¢'~~ - d Z - , ~ * t dS -I~,(~,,~-~*-~dS-f,.~I~,,~*-',,~,,~,*)dS _f, 5 (,,~,, a n + v m , + * >d S - j s ~~ \ t,~.o~:t~s , an/ •
On /
On
,~0,
where the free term c' is due to the singularity of the boundary point. For smooth boundaries, ci = ½. Similar to the radiation problem derived above, the potential for the diffraction problem can be written as
' O~b*+ c' ,l,~ ~ - I'; . [\*~ o,~ 02_~,) d S
_i~,~/.~m~ ~fi o** -. m~ m*~ dS
_ ~,~~ o**~ _,~,,**)~, o~,~,-~,(,~°*" On .
439
3-D Wave diffraction
~(~°'~ + 7 ~~.. -n
-
+ v+7+*
~dS
•
- I~ (,~ °*'/~,.0~ ~
(21)
NUMERICAL METHOD FOR INTEGRAL EQUATIONS
Discretization and modeling of element geometry The basic problem of the present method is the numerical solution of equations (20) and (21). These two equations can be rewritten in a discretized form as
•
c'+i+ ~=~ Z
-~ % ~ - n - #+* dF k=Nl+l
+
E
,k*~
dF+ E
k=N2+l
+
~ k=N4+l
,
k=N3+l
(i(
+
+ Z k=N5+l
,k+~
dr
= 0
(22)
.
where ~t = U ~on~
for radiation problem, j = 1,2 .... 6.
o=-(~:)
for wave-diffraction problem,
I~ is the area of the kth element on the boundary, (Nj - Nj-I) is the number of elements for the boundary Sj and N6 (= N) is the total number of elements. Equation (22) represents, in discrete form, the relationship between the potential at the node i at which the fundamental solution is applied and the potentials of all the elements on the boundary.
440
CHIH-KANG LEE and JACK Y. K. Lou
The discretization requires an approximate evaluation of the relevant integral and an approximate representation of the body surface. Of several possible procedures, the one chosen for this study consists of approximating the boundary surface by a large number of plane quadrilateral elements, as suggested by Hess and Smith (1962), over each of which the potential is assumed to be either constant or varied linearly.
Interpolation functions As in the finite-element method, interpolation functions are used to describe the geometry or the variation of potential ~ and its normal derivative 0_~ for the element. On For geometry, only the linear interpolation function in intrinsic coordinates is used on plane quadrilateral elements through the entire paper, i.e. 4
4
X= E Mk(4,,42)Xk, k=l
Y=
4
~ M k ( 4 1 , 4 2 ) Y k'
Z=
k=l
E
M k ( 4 , , 4 2 ) z~
(23)
k=l
in which (x, y, z) are the coordinates of any point on the element, (xk, yk, zk) are the coordinates of four corners of the element (see Fig. 2) and M~(41, 4e) are the interpolation functions given by 1
M,(41,42) = ~(1 - 4 , ) ( 1 - 42),
M2(4,, 42) = 3 (1 + 41)(1
M3(4., 42) =~(1 + 4,)(1 + 42),
M4(4,, 42) = 14(1 - 40(1 + 42).
-
-
42),
(24)
Higher-order interpolation functions, i.e. quadratic, cubic, etc. can also be used to better represent the geometry of the body, the variables and their normal derivatives, but they are not within the scope of this study.
Constant element
a+
We assume the potential function ~b and its normal derivative ~
to be constant
within each element. Thus, the general expression for the boundary-value problem can be written as
~+'+ k=~ Z . + ~=~+~ E
~kfinn dF 0,. "~ ~
dF
+
" (i
~ ", 1 k=N~
k
"k \
On
(es)
3-D Wave diffraction
/
,,._...7,,.,~.~, ~¥
xJ
441
-,~ (~:1,~'~,~'~)
A plane quadrilateral element and its centroid
~, (-~,~)
(~,~)
_-.~
(~,-~)
(-1, -1
L~trinsic coordinates FI6. 2.
(~
N4
E
k=N3+ 1
A plane quadrilateral element and its centroid.
))
(O+* - iv+* d r
.~ \ On
E
Od~* •~ ~-n d r
k=N5+ 1
~-~(i.~/o°*+~*l~)~o~
~
k=N4+ I
N6
+
4,~+
/
~bk
= ~:,~
t
~* dF q-~', "~
where q-k :
k njk U~o
q-~___(~,/~
for radiation problem, j = 1,2 . . . . 6 for wave-diffraction problem.
442
CmH-KAN¢i I~i~t a n d JAC'K Y. K. Lotr
Notice that for constant elements, the boundary is always "'smooth," hence the coefficient c i has been replaced by ½. Introducing the notations, Hi/~ =
io .
dF and Gik
=
"k
i
(26)
6*dF. "/,
which relates the ith node to the kth element over which the integral is to be performed, equation (25) can be written as N 1
N2
N3
Z H,k +k +
.y_., (/-/,k -- "a,.O +k +
.y_., /4,,,+"
k--1
k-Nl÷
I
k
N4
+
2 k
N5
( H i k - - i P G i k ) ~ bk
N3 ~, 1
~.
+
(H,k + vG,/,) Ok
2 k-N4-,
,'v6
+
N2+l
1
,'v~
H , k ~ = ~. Gi~q~,
k--Ns+l
(27)
k=l
where Hik
=
Hik
for i 4= k, 1
Hi, = H~k + }
f o r / = k.
Finally, equation (27) can be written in a N × N complex matrix form as A X = F,
(28)
where X contains all the unknowns, A is the coefficient matrix (which is dense and nonsymmetric) and the elements of F are known. The solution of equation (28) will give the values of scattered potential for the diffraction problem or radiated potential for the radiation problem. In general, the integrals in equation (26) are computed numerically, except the one corresponding to the node under consideration. For these cases, analytical results for I-I, and G, must be employed (see Lee 1987). Gaussian quadrature formulae suggested by Stroud and Secrest (1966) can easily be applied to the nonsingular integrals with accuracy to the order of ( 2 m - l ) , where m is the number of the integration points. Also, a simple formula exists for approximating the integrals by using source and quadrupole for most cases, as shown in Appendix ,~. If the boundaries of the domain are symmetric, only the nonredundant portion of the boundaries needs to be approximated by quadrilateral elements. The remainder of the boundary is taken into account by suitable reflections, either symmetric or antisymmetric, which results in either saving much computational effort for a given accuracy, or a significant increase in accuracy with the same amount of computational effort.
3-D Wave diffraction
443
Linear element In contrast to the constant element, the potential and its normal derivative will now be assumed to vary linearly within each element. The interpolation function (25), which was previously employed to describe the geometry, is used to represent the linear variation of potential and its normal derivative. This allows us to write 4
4, = ~ Mz(~, ~2) 4,', 1=1 4. (~4,/.t,
(29)
O~ = ~,~ Mt(~l, ~z) ~On] On t=~
(30)
where the superscript I denotes the corner number of the quadrilateral element. Thus,
o4,
the values of 4, and ~nn at any point on the element are defined in terms of their nodal
o4,
values. The integrals, involving the unknowns 4, and ~nn in equations (20) and (21), become
I .k 04,* [hikl 4, ~ - n d r =
hik2 hik3 hik4.
,2
~
(31)
and
Onj
(32)
where
hik~=f.kM~ndF
(33)
/=1,2,3,4
and
(
gikt
=|
Ji "k
M~ 4,* dF
1 = 1,2,3 ,~4
.
(34)
444
('IIIH-KANG L~
a n d JACK Y . K. L o t '
Substituting equations (33) and (34) for all the kth elements into equations (20) and (21), we obtain the following equation for node i,
Y. k~l
h,.+' /
I
+
( h , . - ,,g,..)+' k=,%' I ~ 1
/=1
+ k
N2~" I
I
I
~ k
(ha,,- ivga,,)+'
N 3 ~ I
I~ I
~,
(h,,, + vg,~,)+ z
k = N 4 ~ 1
+
I~ I
"t
~'~ k
N51 1
h,kdb z : /-I
t
~'~ g,kzq z , "~1
/=l
(35)
or in an N × N matrix form, C Y = G,
(36)
where, as (28) in the constant element, Y contains all the unknowns, C is the coefficient matrix which is dense and nonsymmetric and the elements of G are known. The solution of equation (36) will give the values for either the scattered potential or radiation potential for the wave-diffraction or radiation problem, respectively. As for the constant element, the integrals in equations (33) and (34) are computed numerically, except when the collocation node i is inside the element under consideration. The nonsingular integrals hikz and gi~-z can be obtained by the use of Gaussian quadrature integration formulae and can be written as
01~'5:'(~11,,~.~Zq) IJ(~,,, G2q)1%, to'Z'
hik' • Z p
1 q t~t
g,~, = Z p
1
(37)
~'~
t~t
Z +*(~,,,, ~zq)IJ(~,,,, ~e,,)l ~,, ~q, I q
(381
I
where m is the total number of integration points, to, and tog are the weighting factors and (~i,, ~2q) are coordinates of the integration point.
3-D Wave diffraction
445
When the collocation point i belongs to the element under consideration, the kernels in equations (33) and (34) become weakly singular as the integration point approaches point i. Therefore, the integration over the vicinity of point i needs special treatment. Analytical calculation is formidable because of the presence of the Jacobian. However, these singular integrals can be computed by several different ways, which are discussed in Lee (1987).
Discontinuous element When continuous elements are employed, nodes are located at the element edges. When the element edge coincides with a discontinuity in the problem definition, modeling ambiguities will result. Clearly, at an edge of corner of the boundary in a three-dimensional problem, the normal to the boundary is not defined. The problem can be improved by using discontinuous elements, introduced by Patterson and E1-Sebai (1982) and Danson (1982), in which the interior nodes are used to define the boundary elements. Definitely, such elements present difficulties when maintaining inter-element continuity. Inter-element continuity may be advantageous for the boundary-element method but is not necessary. Indeed, the viability of constant elements demonstrates this fact. In the study by Patterson and Sheikh (1984), it was shown that a need for special strategies exists to cope with problems associated with the discontinuity of normal derivatives. The following linear interpolation functions that describe the variation of potential
04,
function 4, and its normal derivative 0n for discontinuous elements are used
M,(~,, ~ ) = ( ~ - ~ , ) ( ~ -
~2), M2(~l,~2)= (~+ ~,)(12-- ~2), (39)
M3(~,,~)= (12+~1)(12 + ~2), .M~(~,,~)= (~-~,)(~ + ~). The nodes are taken at the four interior points (see Lee, 1987). Higher-order interpolation functions can be found in Patterson and Sheikh (1984). This linear discontinuous element will be used to improve the result of the previous linear element and some case studies will be shown in the next section. NUMERICAL RESULTS AND DISCUSSION Simple examples were used to show the convergence rate of the direct boundaryelement method by using different types of elements and integration techniques. The direct method developed was also tested for accuracy and validity for a variety of rigid structures. Results show good agreement with other published data.
Accuracy and convergence of direct boundary-element method Added mass for the heave motion of a cube submerged in an infinite fluid was computed to demonstrate the convergence and accuracy of the direct boundary-element method. It was found that (1) when a four-point Gaussian integration technklue was
446
('ttlII-KANG Le~ and JACK Y. K. Lo;~ 0.75
~ 0.70 m
Const:ant element Linear element Discontinuous linear element
~-~
~ ~ ~
, Stellon and Mavis ~ ' ~ . _
0.85 --
~ 0*50
~
~°°~e
~ • • • • • • •
0.55 ~ iI
I I I I
I 0
~0
20
30
40
50
50
70
80
gO
loo
Number of elements per facet
Fro. 3.
Convergence rates of constant and linear elements.
used, an error ranging from 5 to 15% was introduced if the singularities of the integrals were simply disregarded. However, this error can be minimized for increasing the number of integration points. (2) The approximate method for evaluating the singularities as presented in Appendix A yields very good results even for relatively coarse meshes. (3) The convergence rate of the linear element is much slower than that of constant element (see Fig. 3). The reason for this is not yet known. It is probably due to the discontinuity of the normal derivative of the potential at the edges or corners of the boundary, since using discontinuous linear elements which lead to continuous normal derivatives has shown improvement in convergence rate. Detailed results are presented in Lee (1987) and will not be repeated here.
Radiation problems Added mass of a sphere in an infinite fluid. Again, both constant and linear elements were used to compute to heave added mass coefficients, while plane elements were used to approximate the curved spherical surface. The results are presented in Table 1. It was found that the constant element yields more accurate results. The errors were less than 3 and 1% respectively, when 24 and 40 constant elements were used in the computations. The results also show that considerable saving in computer time can be achieved by taking advantage of the symmetry of the problem.
3-D Wave diffraction TABLE 1.
447
HEAVE ADDED MASS COEFFICIENTS OF A SPHERE IN AN INIZlNITE FLUID BY USING CONSTANT AND LINEAR ELEMENTS WITH SYMMETRICAL PLANES
Number of elements per quarter
Cases
a33/pV
0.448 *
CPU
0.42
t 0.26
a33/pV 0.96
t 0.78
40
0.486
0.497
t 2.24
~ 1.00
*
28.5
0.438 $ 0.67
analytical result Note:
*
7.50
0.291 *
CPU
~ 0.21
24
~ 12.0
t 6.04
t 6.40
~ 2.57 0.441
$ 3.81
* 29.4
t 12.7
$ 8.90
0.5
a33/pV denotes added mass coefficient. Computing CPU based on the case a with 24 elements per quarter and two symmetrical planes as a unit. *denotes the case without plane of symmetry. tdenotes the case with one plane of symmetry. :]:denotes the case with two planes of symmetry. Case a uses constant elements. Case b uses linear elements.
A d d e d mass of a box in semi-infinite water. The box has a length and beam of 90 m with draft of 40 m. The water depth is infinite. The calculated added mass and damping coefficients of surge, heave and pitch motions are presented in Figs 4-6. The numerical results based on indirect boundary-element method and experimental results of Faltinsen and Michelsen (1974), are also presented in these figures for comparison purposes. Agreement between these two data sets is generally good, with the exception of damping coefficients in pitch and surge motions. This could be due to (1) too coarse mesh size on the free surface ($2) and (2) too short a distance between the radiation boundary ($4) and the structure. Diffraction problems Wave forces on a submerged hemisphere on the sea bottom.Horizontal and vertical wave forces on a submerged hemisphere of a radius a in a water depth of 3a were computed. The incident sea wave amplitude is ao. The nondimensionalized forces are shown in Fig. 7, where K is the wave number that satisfies the dispersion relation. Data presented by Chakrabarti and Naftzger (1974) are also included for comparison purposes. Good agreements for both horizontal and vertical forces were obtained by employing quadrilateral constant elements. Wave forces on the submerged storage tank and tower. Both the horizontal and vertical forces were calculated by using constant elements. The results are presented in Fig. 8. Comparison with the data obtained by Hogben and Standing (1974), who used
448
(?HIH-KANGLE~-;and JA(K Y. K. Lo~: 2.0
2.0
~
FalfJn~en
and
Mlchelsen
-~.
Faltinsen
and
Idlchel~en experiment
~-
Pre~ent
method
damping
added
=l > 1.oo q. "--~
mass
.~•
O. ~ I ~
~
~"~--~_
"~ ~
oo
,
•
~
~ ' ~
•• .,~*~
0.5--
-.
1.0
>
~.
•
-- 0.5
..
0.0
8
I 10
N
I 12
I
I 14 Period
FIG. 4.
~
,
16
I 18
0.0
20
(sec.)
Surge added mass and damping coefficientsfor a floating box (L x B×d=90m ×90m × 40m) in semiinfinite water.
the indirect boundary-element method, substantiates the validity and accuracy of the present method. Wave forces on a box in semi-infinite water. By using the same structure in the radiation problem, wave-exciting forces and moments on the floating box, which has length and beam of 90 m, at a draft of 40 m were calculated using constant elements. The results, together with the data obtained by Faltinsen and Michelsen (1974), who used indirect boundary-element method, are shown in Figs 9-11. Again, agreement is satisfactory.
CONCLUSIONS AND RECOMMENDATIONS The main purpose of this research is the development of a numerical method for studying fluid-structure interaction problems and to verify its applicability and accuracy. A variety of problems, including the wave-diffraction and radiation problems, were studied using the method developed. Most results obtained by this method agree with those obtained by previous investigators using other methods, including analytical, numerical or experimental. It would be of further interest to extend this method to study the wave resistance of
3-D Wave diffraction
449
2.0
2.0
~ -l. ~-
F~ltln~en and Idi~.hel~en Falt~;-,~,-; and Ml~hl~en experiment Pr~ent method ~
1.5
-,.0:l add~
8~ . . . .
•- ~ . . . ~ - , ~
mass
. . . . .
~.___
" ' " - . e . . .
e.o.
e
O.~m
-
0.S
damping e-...e...
~.
.-e
... o
0,0
0.0
8
10
12
14,
18
18
20
Period (sec.) FI~. 5.
Heave added mass and damping coefficients for a floating box infinite water.
(L×B×d=9Om×9Omx4Om) in
semi-
a surface ship with forward speed, although it may be difficult to express the radiation condition in moving coordinates. It should be pointed out that, in principle, any boundary geometries, no matter how complicated, can be handled by this method. Modeling the boundary surface by using quadratic elements is recommended for future studies. Not only do quadratic elements represent the curved boundary surface better, but they can also improve the variation of the unknown function and its normal derivative on the boundary. The boundary-element method is suited for solving problems that have a high ratio of domain volume to boundary surface area, especially, for problems having infinite or semi-infinite domains. However, special techniques are needed for solving a large number of full, nonsymmetric systems of equations. Further research is needed to obtain an in-depth undersianding of the present technique. Although Green's function for an infinite domain and the boundary-element method can be used conjunctively to solve problems with complicated boundaries, it would result in a large-sized system of equations. To solve these equations would require an excessive computer time. Thus, it would be much more advantageous to use the special Green's function for the boundaries under consideration in some tractible problems. For instance, for problems of wave diffraction and radiation in semi-infinite water, an
450
CHffI-KAN(; LEE and JACK ¥, K. Lou 0,20
~ -~.
~
0.020
F'a~n~en and MIchel~en Fa~Jn~e~ and M l ~ e l ~ e x
0.15 ~
-- 0.015
~ ~ .~ o.o-
"~
~1~ x'~,N,N, x 0.0~ ~
~ ~
~
~
~
x
~
....
0.~
"
~ ' ~ - ~ . ~
I
1
I
1
I
10
12
14
18
18
Pe~od FI~. 6.
"¢
~ --
o.~
-o.o,o I ~
~
.
o.~o
20
(~ec.)
Pitch added mass and damping coefficients for a floating box (L x B x d=90m x 90m x40m) in semi-infinite water.
1o0 , , |
~'~
~
Chakrabartl and Naf~gm"
0.8---2
~
I~
~
i,
~
~
0.4--
0.,-/]
¢__
~
"\
-_,
-~.~.,
0.0
0
0.0
0.1
0.2
0.3
0.4.
0.5
0.8
0.7
0.8
0.9
1,0
ka
FI6. 7.
The horizontal and vertical forces on a submerged hemisphere on the sea bottom.
3-D Wave diffraction
451
~,~
- - 0.3
l•
o.,,~--t ,',
~
/,'
o
~
",~'~,.
w~, / \ °=1 ~i,,,,, ::=~.:,,,~~,, ..tZo.~ ~ ~. / l'/,' ',',,., "X .~L.J, o=1/ ":,x,....: ~;;~,,,,, 1; ',",//
I , , ~,
-- 0.2
':': ~~-.. ._
o..
--0.1
~-~'a ~,,~ ~,
~.~
[
I
I
1
~
~
~
3
0.0
k~ FI~. 8.
The horizontal and vertical forces on the submerged tank and tower.
- - F o l t l n l e n 0nd Mlchell~n -- - I~m~t rn~tho4
2.0~
1.5--
/
.,.
o
~o+ 1.0-~.
/ / s
0.5~
0.0 10
12
14 Pe~
FI6. 9.
Surge wave force amplitudes for a floating box
18
1~
(se~.)
(LxBxd=9Omx9Omx4Om)in
semi-infinite water.
452
CItItt-KANG LEE and J,,~cK Y. K, L,otr
~
FaltJn~en
and Mlchelsen
- - -- Prewent m e t h o d
2.0--
1.5 o
~ '~ I "~°1"° ~.
--
0.5--
0,0
I
10
I
12
14. Period
FIG. 10.
16
18
20
(S~.)
Heave wave force amplitudes for a floating box
(LxBxd=9Omx9Omx4Om) in semi-infinite water
0.20 - -
0.15 --
/
:°oo,o_
I / / /
/
8
10
- -
Fdtln~
~-
P~t~ent m e t h o d
12
14. Period
FIG. 11.
~nd ~ l ¢ h ~ t ~ n
18
18
20
(S~.)
Pitch wave moment amplitudes for a floating box water.
(L×Bxd=9Omx9Omx4Om) in
semi-infinite
3-D Wave diffraction exact expression for the Green's
f u n c t i o n s o l u t i o n c a n b e f o u n d in W e h a u s e n
453 and
L a i t o n e (1960) as
u* = -1 + P V f f k + ~
e k~+~) J o ( k R ) d k + i 2"rrl~e~ + t ' ) Jo(~xR).
T h e first a u t h o r has e m p l o y e d this s o l u t i o n to s o l v e s o m e o f t h e test c a s e s c i t e d p r e v i o u s l y . P r e l i m i n a r y r e s u l t s i n d i c a t e t h a t significant i m p r o v e m e n t c a n b e a c h i e v e d ; for example, the surge and pitch dampings can be more accurately predicted with a s m a l l e r - s i z e d s y s t e m o f e q u a t i o n s as c o m p a r e d to t h e p r e v i o u s results.
REFERENCES BREB~IA, C. A. 1978. The Boundary Element Method for Engineers. Pretech Press, London. BREBBIA, C. A. 1981. Progress in Boundary Element Methods-1. Wiley, New York. CnAK~,~SAR~, S. K. and NAVrZ~ER, R. A. 1974. Non-linear wave forces on half-cylinder and hemisphere. J. Watways Harbors Coastal Eng.ng. Div., Am. Soc. Chem. Engrs 1~, 189-204. DANSO~, D. J. 1982. A boundary element analysis system. In Boundary Element Methods in Engineering. Edited by B~rss~A, C. A., pp. 557-575. Springer, Berlin. F,~LW~SrN, O. M. and MICnrLS~, F. C. 1974. Motions of large structures in waves of zero Froude number. Proc. Int. Syrup. Dynamics of Marine Vehicles and Structures in Waves, pp. 99-114. Institute of Mechanical Engineering, London. F ~ K , W. 1967. Oscillation of cylinder inor below the free surface of deep fluids. Report No. 2375. Naval Ship Research and Development Center, Bethesda, MD. G~lSON, C. J. 1974. Hydrodynamics of large objects in the sea, part I--Hydrodynamic analysis. J. Hydronautics 8, 5-12. GAte,SOY, C. J. 1978. Chapter 3---Hydrodynamic loading of large offshore structures: three-dimensional source distribution methods. In Numerical Methods in Offshore Engineering. Edited by ZIEr~K~EW~CZ,O. C., LEWIS, R. W. and Sv~,~, K. G., pp. 87-140. Wiley, New York. HAVELOCK,T. H. 1955. Waves due to a floating sphere making periodic heaving oscillations. Proc. R. Soc. Lond. 230A, 1-7. HESS, J. L. and SranrH, A. M. O. 1962. Calculation of nonlifting potential flow about arbitrary threedimensional bodies. Report No. E.S. 40622. Douglas Aircraft Division, Long Branch, CA. Ho~aEr~, N. and S~A~D~G, R. G. 1974. Wave loads and large bodies. In Proc. Int. Symp. Dynamics of Marine Vehicles and Structures in Waves, pp. 258-277. Institute of Mechanical Engineering, London. JoHr~, F. 1950. On the motion of floating bodies II. Commun pure appl. Math. 3, 45-101. KELLO6~, O. D. 1929. Foundations of Potential Theory. Springer, Berlin. LEE, C. K, 1987. A coupled boundary and finite-element method for three dimensional fluid-structure interaction problems. P h . D . dissertation, Texas A&M University, T. LEW~S, F. M. 1929. The inertia of water surrounding a vibrating ship. Trans. SNAME 37, 1-20. NEWMAN, J. N. 1978. Marine Hydrodynamics. The MIT Press, Cambridge, Mass. PAttErSON, C. and EL-SEaA~, N. A. S. 1982. A regular boundary method using non-conforming elements for potential problems in three dimensions. In Boandary Element Methods in Engineering. Edited by BRE~mA, C. A., pp. 137-152. Springer, Berlin. PA~E~SOt~, C. and SHEIKH, M. A. 1984. Interelement continuity in th eboundary element method. In Topics in Boundary Element Research--1. Edited by B~,EBa~A,C. A., pp. 123-141. Springer, Berlin. SrELSOS, T. E. and MAVIS, F. T. 1955. Virtual mass and acceleration in fluids. Proc. ASCE 81, 670. Sxaot~), A, H. and SECaESX,A. H. 1966. Gaussian Quadrature Formulae. Prentice-Hall, New York. URSELL, F. 1949a. On the heaving motion of a circular cylinder in the surface of a fluid. Q. J. Mech. appl. Math. 2 (2), 218--231. URSELL, F. 1949b. The rolling motion of a circular cylinder in the surface of a fluid. Q. J. Mech. appl. Math. 2 (3), 335-353. WEnAUSE~, J. V. and LAIrO~r, E. V. 1960. Surface waves. Encyclop. Phys. 9, 446-778.
454
CHIH-KANG LEE and JACK Y. K. Lou node p(x,y,z)
' [ '
c~, ~ ) l
'Z:__%_ ~
l/0.w ,,,"' --/
--,,,,-,=,,-,---I~ X
,~
(~2,r~2)
Fro. A . I .
DERIVATION
OF
Local element system.
APPENDIX A APPROXIMATE FORMULAE COEFFICIENTS
FOR
INFLUENCE
Consider a potential point P(x, y, z) and a quadrilateral source element in the xy-plane, as shown in Fig. l A. 1. The integrand, i.e. 3' in coefficient Gij in equation (33), can be expanded in a power series in ~ and -q about the origin as 1 r
l
V/~ - ~)2 + (y - Xl)2 + z2
(A.l) l 1 = ~ - w..{ - ~o.~n + ~ to.. ~" + ~..,, { n + ~ ~o.,..,,n-',
where
tO =
1 r( i
r. = ~r 2 + y2 + z 2. It can be show that an approximation for the influence coefficient G,, can be given as G
1
3-D Wave diffraction =
1
A to - (o.,~,, M,~ - o%, Ms) + ~ (to.~ I~,~ + 2 o.~r l~y + tO.ry Iyy)
455
1
(A.2)
,
where
Z:I0 , lyy = f ~q2d~d~.
Ly=fK~dKd~,
In equation (A.2), the first term can be interpreted as the potential of a point source of strength A. The terms in the first parathesis are the dipoles, the axes of which are along the x-axis and y-axis, respectively. The remaining terms are quadrupoles. In the present application, since the origin of the coordinate system is at the centroid of the element, all the first moments, M, and My, vanish; thus, only the source and quadrupole terms exist. 1 From equation (A.1), we can obtain the normal derivative of the fundamental solution, , as r
1
On - - t ° ~ - ° ' ~ z ~ - ° > ' ~ + ~
°~xx~2
+OxYz~+~Oyyz~2"
(A.3)
Thus, we have
H,~ = 4~r
1[
dr
1
l
= 4 ~ A o , - (o.~z Mx - o r, M r ) + } (oxxz lxx + 2 ~Oxyz lxy + o~yz !,,) •
(A.4)
As in equation (A.2), only the first and third terms exist, corresponding to the source and quadrupole. Equations (A.2) and (A.4) can be used in two forms. In the first, both the source and quadrupole are used for intermediate distant points. In the second, only the source is used for a far distant point. According to the study by Hess and Smith (1962), the form of the approximation used for a particular point depends on the ratio of the distance between the point and the centroid of the quadrilateral to the major diagonal of the quadrilateral. If this ratio is less than 2.45, we use the exact method; if this ratio is greater than four, the second form is adequate; otherwise the first form may be employed. In addition, Hess and Smith's results showed that the maximum error owing to the approximate forms is less than 0.1%. Thus, the approximate forms are quite successful in computation without losing accuracy and these forms were used in the present work to save computing time.
0029-8018/88 $3.00 + .00 Pergamon Press pie
A DIRECT BOUNDARY-ELEMENT METHOD FOR THREE-D WAVE DIFFRACTION AND RADIATION PROBLEMS CHIH-KANG LEE* *Veritas Technical Services Inc., Houston, TX, U.S.A.
and JACK Y. K. Lout tOcean Engineering Program, Texas A&M University, College Station, TX77843, U.S.A. Abstract--Three-dimensional analytical procedures are developed for determining the waveexciting and motion-induced hydrodynamic forces for fixed and floating ocean structures in a fluid of finite or infinite depth. The fluid potential is expressed by means of a Helmholtz integral, which involves the normal fluid velocities at the fluid-structure interface. Upon discretizing the integral equation by using constant or linear elements and by enforcing the velocity continuity condition at the discrete points of the interface, the fluid potential can be determined. A computer program based on the direct boundary-element method is developed to calculate the wave-exciting and motion-induced hydrodynamic forces. Effects of structural flexibility are not considered in the present paper. Comparisons of the results with those obtained by several previous investigators reveal good agreements.
INTRODUCTION THE TOPIC of fluid loading of ships and large structures has received considerable attention in recent years. For design applications one must assess, with a high degree of certainty, the magnitude of the load induced by the extreme waves associated with the design storm conditions and the load imposed by the motion of the structure. Thus, even though the numerical procedures in free-surface hydrodynamics tend to be rather involved, practical needs exist for dealing with the hydrodynamics of large structures. Of particular concern has been the development of the wave-exciting force and the socalled added mass and damping coefficients. Wave-body interaction at infinitesimal amplitudes is a richly developed subject that, using the so-called indirect boundary-element method (or source distribution method), has served the field of naval architecture with great success. Its applications in offshore technology have stimulated further advances in recent years. Ursell (1949a,b) studied the hydrodynamics of two-dimensional semi-circular or eliptical cylinders in an infinite fluid. He introduced nonorthogonal expansions (multiples) whose terms satisfy the Laplace's equation and the free surface condition, but vanish at infinity. A source or dipole is then added to satisfy the radiation condition. The expansion coefficients and strengths of the singular sources were determined by imposing the body boundary condition at a finite number of collocation points. Frank (1967) used an integral equation approach and the strip method for solving ship wave problems. He assumed that the velocity potential could be represented by a distribution of two-dimensional sources over the contour of a ship segment. The 431
-~32
CIIlH-KANGL~-. and JACKY. K. Lot
unknown source strength is determined by satisfying the boundary condition at the midpoint of the segment. Frank's method can be applied to a two-dimensional section of arbitrary shape. However, it appears difficult to generalize the approach for threedimensional bodies. For the infinite fluid problem in three dimensions, Lewis (1929) computed exactly the added mass of a vibrating spheroid by the method of separation of variables in spheroidal coordinates. Hess and Smith (1962) solved the three-dimensional infinite fluid problem for a body of arbitrary shape in steady motion by solving the Fredholm integral equation of the second kind. For the three-dimensional free surface problem, Havelock (1955) used the multipole expansion method to solve the problem of a heaving sphere. As in the two-dimensional expansion method, he introduced a velocity potential represented by a set of higherorder singularities that satisfy the free surface condition and radiation condition trivially, and a wave source or multipole at the origin within the body that satisfies the radiation condition. The unknown coefficients for the polynomials and the wave sources or multipoles are determined from the kinematic boundary condition on the interface. The linear algebraic equations for the unknown coefficients are solved numerically. An integral equation formulation for calculating wave forces on arbitrarily shaped large bodies was developed by John (1950). The boundary-value problem is posed in terms of velocity potential functions, which are equivalent to a distribution of singularities on suitable boundaries of the fluid region, with the strengths of the singularities determined by the boundary conditions. Solution of the resulting integral equations is approximated by choosing a discrete set of points to obtain a set of matrix equations. The method has received widespread usage in free surface flow problems and various forms of the Green's function for wave potentials can be found in the work of Wehausen and Laitone (1960). Their research has been extended and boundary integral equations have now been frequently used to solve wave-diffraction and radiation problems; for instance, in the work of Garrison (1974, 1978), Hogben and Standing (1974), Faltinsen and Michelsen (1974). In this paper, an approach using a direct boundary-element method is developed which has more physical relation to the wave-body interaction problem. Like the finiteelement method, the direct boundary-element method also divides the external surface of a domain into a series of elements over which the functions under consideration can vary in different ways. This capability is important as the boundary integral equation formulations generally were restricted to constant source strength over the elements and the sources were assumed to be concentrated at a series of points on the external surface of the boundary. FORMULATION OF THE BOUNDARY-VALUE PROBLEMS A right-handed Cartesian coordinate system with the origin on the free surface is used as shown in Fig. 1. The angle o~ is the wave-heading angle with respect to the positive x-axis, the xy-plane coincides with the free surface when the fluid is at rest and the z-axis is taken directly opposite to the force of gravity. The structure is considered as a large rigid body oscillating sinusoidally in a long crested regular wave. The amplitudes of the motions of the structure, as well as of the wave, are assumed to be small. Since the fluid is assumed to be ideal and irrotational, there exists a velocity potential that governs the kinematics of the fluid and the velocity
3-D Wave diffraction
433
Z
v
I
s~
s3
I I I L
~'-~l~#~'~I~'#'~'l ss FIG. 1.
s~ 3-D Coordinate system and various boundaries.
V may be represented by the gradient of the potential function 4. Furthermore, we assume that the motion is sinusoidal in time. The total velocity potential is written as , ( x , y , z ; t ) = qJ~ + Oa = R (dO(x,y,z) e ''t }
(1)
where 4, and ~d are the r~diation potential and the potential for incident and diffracted waves respectively; 1:1 represents the real part of a complex quality and d~(x,y,z) is the spacial component of 4. Since a free floating body can move in a wave with six degreesof-freedom, the velocity potential is conveniently broken into components relating to the incident and diffracted wave fields and the motions of the body 6
+(x,y,z) = +,, + +~ + ~] +~, ~=!
(2)
434
()IIHt-KANG Ll~t~ and JACk. Y. K, Lot,
where +o is the complex amplitude of the incident wave potential, d)v the diffracted wave potential and ~bj the velocity potential due to the forced motion of the body in the jth mode; j is taken to be 1. . . . 6 for surge, sway, heave, roll, pitch and yaw. as shown in Fig. 1. Radiation p r o b l e m
Consider a rigid body oscillating with a frequency to on or beneath the free surface of a given fluid, which is otherwise at rest. The velocities of the body are described by U~ = F~ { - i t o ~ e
i~"}=~{Ui,,e-i"
j = l . . . . 6,
}
(3)
where ~i is the complex amplitude of the motion of the body in the jth mode. Due to the forced motion of the body in the jth mode, the complex velocity potential ¢b~(x,y,z) must satisfy the Laplace equation, 7 2 +~ = 0
in the fluid region, ~ ,
(4)
together with the following boundary conditions hdm
"~J On 0+j
ito~n i
v ¢b~ = 0
j = 1. . . . 6
on z = 0
on the body surface Sl,
at the free surface S~,
Oz
on rigid wall $3 and sea bottom $6,
On
lim
0~
Rr O+j - iv+j \ Or
. . . . . v d~j = 0 Oz
)
= 0
at the distant boundary S4,
(7)
(8)
on the imaginary boundary
S_s in case of infinite water depth,
o) 2
(6)
-
0+j = 0
~r--+a:
(5)
(9)
where v = - - ; r = (x 2 + y2)½; R is the wave number which satisfies the dispersion g relation, K tanh k d = v; S~ . . . . $6 are the categorized boundaries of the domain, as shown in Fig. 1 and the derivative in equation (5) is taken with respect to the outward normal n of the body surface; n~-n6 are the generalized direction cosines on Sl. The set of equations (4-9) constitutes a boundary-value problem. To solve for the radiated wave potential for each mode of motion, a boundary-element method can be
3-D Wave
diffraction
435
employed. The total hydrodynamic pressure can be obtained from the linearized Bernoulli's equation (Newman 1978),
I
p=-p~-=-pFl
- i ~ o e -'°''
~+,
.
(10)
i=l
The oscillating hydrohynamic forces (i=1,2,3) and moments (i=4,5,6) acting on the body surface can be determined by integrating the fluid pressure (10) over the wetted surface S~ and the result can be written as 6
~=-
~(a~+bo~).
(11)
j=l
The coefficients aij and b 0- are known as the added mass and damping coefficient respectively.
Diffraction problem The velocity potential for the diffracted waves is computed as if the structure were restrained from moving. Similar to the radiated wave potential, the diffracted wave potential should satisfy the Laplace's equation (4) and the boundary conditions (5-9). Since the body is assumed to be restrained, the boundary condition on the body surface S~ becomes
0+7 _ O~
a+,, O~
(12)
'
where the incident wave potential is given by
d?o -
i gtoa,, cOShcoshk(Zkh + h) eitK,x+K2y),
(13)
K1 = g COSOt, I(.2 = K sinc~ and a., K and to are the amplitude, wave number and frequency of the incident wave, respectively. Thus we have a set of equations that constitutes a wave diffraction boundary-value problem. Again, a boundary-element method can be employed to solve the diffracted wave potentials and the diffracted pressure can be obtained from the Bernoulli's equation. The resulting wave-exciting forces and moments, Fj, j = 1 . . . . 6, can be written as
14, I
436
CHIH-KANG LEE
and JACK Y.
K. L()L~
BOUNDARY-ELEMENT METHOD Finite difference, finite element (domain methods), and the boundary-element (boundary method) methods are frequently used to solve the above differential equations numerically. The difference between the domain methods and the boundary methods is significant. The domain methods discretize the domain, while the boundary methods discretize the boundary. Thus, the boundary method reduces the dimension by one. This method is well suited to problems in which the limits of the domain are infinite or difficult to define, in that the solution is applied to the boundary rather than the domain. There are two types of boundary element methods, the "indirect" and the "direct". In the "indirect boundary-element" method, which is the "source" method developed by Kellogg for potential problems, the unknowns are not the physical variables of the problem. On the other hand, the "direct boundary-element" method not only allows us to solve the problem in terms of physical variables but also serves as the first step towards a better understanding of the technique and its relationship to other approximate methods, in particular, mixed finite-element methods. The direct boundary-element method involves the transformation of the governing differential equation within the domain under consideration to an integral equation defined on the surface of the domain. The surface then may be discretized into a number of elements over which a polynomial form of the solution is assumed. This enables the evaluation of the relevant integrals, usually by some numerical process, resulting in a system of linear algebraic equations. As only the surface of the domain needs to be discretized, the resulting systems of equations are considerably smaller than those involving domain type solutions, e.g. finite elements, finite differences, and considerable savings also can be achieved in the time required to prepare data. The key to this method, in potential problems, is the adoption of a fundamental solution used to eliminate the domain integral in the formulation. Such solutions can be found in Brebbia (1981).
Formulation of boundary integral equations The direct boundary-element method used to solve boundary-value problems can be formulated rigorously, using either an approach based on Green's theorem or a particular case of the weighted residual method. With the advantages of being more powerful and the ability to relate boundary solutions to other more classical engineering techniques, the weighted residual method is preferred. For the radiation and diffraction problems, which involve more general impedance boundary conditions of the Robin type, we can write the following error functions; for the radiation problem ~ = V2 +j a+j e~ - On
~2 = --
0+~ _ O~F/
in ~,
Ui° ni
on S~,
+~
on $2,
v
on
$3,
(15)
3-D Wave diffraction
0~
¢4 = ff~n - i v ~b~
0,~
~-s = ff~n + v d~
on 84, on S.s,
0+j ~-6 :
~
$6.
on
On
437
We then minimize the errors (15) by distributing a weighting function tb*. The resulting weighted residual statement can be written as
fn(V%3~*dt~=fsl(O~/-U~(a~)'b*dSfs~(O~J-v~J) On
+ ( O*i,'dS+ Js3 On
k On
~ (O*'-iv,O,*dS+
Js 4 k On
fs ( O * Y + v * , ) * * d S ~ k On
+ fs~** dS.
(16)
Integrating by parts the terms on the left-hand side of equation (16) twice will yield,
-
+
+f~4i~+*dS-~.~*~+*dS-~+,O--~ndS.
(17)
.
We now take the simple fundamental solution of V:~b* + ~i = 0 as the weighting function, i.e.
4" -
1 4~rr
(18) '
where r = X/(x - xi) 2 + (y - yi) 2 + (z - zi)2; (x, y, z) and (Xi, Yi, Zi) are the coordinates of the field and source points respectively. Substituting equation (18) into equation (17) and rearranging, we obtain
~ = -
On f,l \l, '.°**
)
Uj,, nj dp* dS
438
(hm,:-K,~.~G LEE and JACK Y, K. Lou
_i,,~(~,~,~:~_,~,~,~,,/dS _
0n
_f,,(,~e,:t~s_f,g,.o~:_ ~an a~ /
.
i v d~ +*) dS
~,
_f~.~(%o+* . an + v +~+*) dS
({+.~:~ On.]dS"
(19)
-L. ~ '
Equation (19) implies an integral relationship between the potential at the field point i inside the domain l~ and its values on the boundary S of the domain. When this point is moved to the boundary S, the integrals involved in equation (19) become singular at i and must be evaluated in the Cauchy Principle Value sense (Brebbia, 1978). Thus, for any point i on the boundary surface, the boundary integral equation can be written as
c'~=-
i,.(, ~-o~~* -u~,,~,~* t d S - i,(~ ¢'~~ - d Z - , ~ * t dS -I~,(~,,~-~*-~dS-f,.~I~,,~*-',,~,,~,*)dS _f, 5 (,,~,, a n + v m , + * >d S - j s ~~ \ t,~.o~:t~s , an/ •
On /
On
,~0,
where the free term c' is due to the singularity of the boundary point. For smooth boundaries, ci = ½. Similar to the radiation problem derived above, the potential for the diffraction problem can be written as
' O~b*+ c' ,l,~ ~ - I'; . [\*~ o,~ 02_~,) d S
_i~,~/.~m~ ~fi o** -. m~ m*~ dS
_ ~,~~ o**~ _,~,,**)~, o~,~,-~,(,~°*" On .
439
3-D Wave diffraction
~(~°'~ + 7 ~~.. -n
-
+ v+7+*
~dS
•
- I~ (,~ °*'/~,.0~ ~
(21)
NUMERICAL METHOD FOR INTEGRAL EQUATIONS
Discretization and modeling of element geometry The basic problem of the present method is the numerical solution of equations (20) and (21). These two equations can be rewritten in a discretized form as
•
c'+i+ ~=~ Z
-~ % ~ - n - #+* dF k=Nl+l
+
E
,k*~
dF+ E
k=N2+l
+
~ k=N4+l
,
k=N3+l
(i(
+
+ Z k=N5+l
,k+~
dr
= 0
(22)
.
where ~t = U ~on~
for radiation problem, j = 1,2 .... 6.
o=-(~:)
for wave-diffraction problem,
I~ is the area of the kth element on the boundary, (Nj - Nj-I) is the number of elements for the boundary Sj and N6 (= N) is the total number of elements. Equation (22) represents, in discrete form, the relationship between the potential at the node i at which the fundamental solution is applied and the potentials of all the elements on the boundary.
440
CHIH-KANG LEE and JACK Y. K. Lou
The discretization requires an approximate evaluation of the relevant integral and an approximate representation of the body surface. Of several possible procedures, the one chosen for this study consists of approximating the boundary surface by a large number of plane quadrilateral elements, as suggested by Hess and Smith (1962), over each of which the potential is assumed to be either constant or varied linearly.
Interpolation functions As in the finite-element method, interpolation functions are used to describe the geometry or the variation of potential ~ and its normal derivative 0_~ for the element. On For geometry, only the linear interpolation function in intrinsic coordinates is used on plane quadrilateral elements through the entire paper, i.e. 4
4
X= E Mk(4,,42)Xk, k=l
Y=
4
~ M k ( 4 1 , 4 2 ) Y k'
Z=
k=l
E
M k ( 4 , , 4 2 ) z~
(23)
k=l
in which (x, y, z) are the coordinates of any point on the element, (xk, yk, zk) are the coordinates of four corners of the element (see Fig. 2) and M~(41, 4e) are the interpolation functions given by 1
M,(41,42) = ~(1 - 4 , ) ( 1 - 42),
M2(4,, 42) = 3 (1 + 41)(1
M3(4., 42) =~(1 + 4,)(1 + 42),
M4(4,, 42) = 14(1 - 40(1 + 42).
-
-
42),
(24)
Higher-order interpolation functions, i.e. quadratic, cubic, etc. can also be used to better represent the geometry of the body, the variables and their normal derivatives, but they are not within the scope of this study.
Constant element
a+
We assume the potential function ~b and its normal derivative ~
to be constant
within each element. Thus, the general expression for the boundary-value problem can be written as
~+'+ k=~ Z . + ~=~+~ E
~kfinn dF 0,. "~ ~
dF
+
" (i
~ ", 1 k=N~
k
"k \
On
(es)
3-D Wave diffraction
/
,,._...7,,.,~.~, ~¥
xJ
441
-,~ (~:1,~'~,~'~)
A plane quadrilateral element and its centroid
~, (-~,~)
(~,~)
_-.~
(~,-~)
(-1, -1
L~trinsic coordinates FI6. 2.
(~
N4
E
k=N3+ 1
A plane quadrilateral element and its centroid.
))
(O+* - iv+* d r
.~ \ On
E
Od~* •~ ~-n d r
k=N5+ 1
~-~(i.~/o°*+~*l~)~o~
~
k=N4+ I
N6
+
4,~+
/
~bk
= ~:,~
t
~* dF q-~', "~
where q-k :
k njk U~o
q-~___(~,/~
for radiation problem, j = 1,2 . . . . 6 for wave-diffraction problem.
442
CmH-KAN¢i I~i~t a n d JAC'K Y. K. Lotr
Notice that for constant elements, the boundary is always "'smooth," hence the coefficient c i has been replaced by ½. Introducing the notations, Hi/~ =
io .
dF and Gik
=
"k
i
(26)
6*dF. "/,
which relates the ith node to the kth element over which the integral is to be performed, equation (25) can be written as N 1
N2
N3
Z H,k +k +
.y_., (/-/,k -- "a,.O +k +
.y_., /4,,,+"
k--1
k-Nl÷
I
k
N4
+
2 k
N5
( H i k - - i P G i k ) ~ bk
N3 ~, 1
~.
+
(H,k + vG,/,) Ok
2 k-N4-,
,'v6
+
N2+l
1
,'v~
H , k ~ = ~. Gi~q~,
k--Ns+l
(27)
k=l
where Hik
=
Hik
for i 4= k, 1
Hi, = H~k + }
f o r / = k.
Finally, equation (27) can be written in a N × N complex matrix form as A X = F,
(28)
where X contains all the unknowns, A is the coefficient matrix (which is dense and nonsymmetric) and the elements of F are known. The solution of equation (28) will give the values of scattered potential for the diffraction problem or radiated potential for the radiation problem. In general, the integrals in equation (26) are computed numerically, except the one corresponding to the node under consideration. For these cases, analytical results for I-I, and G, must be employed (see Lee 1987). Gaussian quadrature formulae suggested by Stroud and Secrest (1966) can easily be applied to the nonsingular integrals with accuracy to the order of ( 2 m - l ) , where m is the number of the integration points. Also, a simple formula exists for approximating the integrals by using source and quadrupole for most cases, as shown in Appendix ,~. If the boundaries of the domain are symmetric, only the nonredundant portion of the boundaries needs to be approximated by quadrilateral elements. The remainder of the boundary is taken into account by suitable reflections, either symmetric or antisymmetric, which results in either saving much computational effort for a given accuracy, or a significant increase in accuracy with the same amount of computational effort.
3-D Wave diffraction
443
Linear element In contrast to the constant element, the potential and its normal derivative will now be assumed to vary linearly within each element. The interpolation function (25), which was previously employed to describe the geometry, is used to represent the linear variation of potential and its normal derivative. This allows us to write 4
4, = ~ Mz(~, ~2) 4,', 1=1 4. (~4,/.t,
(29)
O~ = ~,~ Mt(~l, ~z) ~On] On t=~
(30)
where the superscript I denotes the corner number of the quadrilateral element. Thus,
o4,
the values of 4, and ~nn at any point on the element are defined in terms of their nodal
o4,
values. The integrals, involving the unknowns 4, and ~nn in equations (20) and (21), become
I .k 04,* [hikl 4, ~ - n d r =
hik2 hik3 hik4.
,2
~
(31)
and
Onj
(32)
where
hik~=f.kM~ndF
(33)
/=1,2,3,4
and
(
gikt
=|
Ji "k
M~ 4,* dF
1 = 1,2,3 ,~4
.
(34)
444
('IIIH-KANG L~
a n d JACK Y . K. L o t '
Substituting equations (33) and (34) for all the kth elements into equations (20) and (21), we obtain the following equation for node i,
Y. k~l
h,.+' /
I
+
( h , . - ,,g,..)+' k=,%' I ~ 1
/=1
+ k
N2~" I
I
I
~ k
(ha,,- ivga,,)+'
N 3 ~ I
I~ I
~,
(h,,, + vg,~,)+ z
k = N 4 ~ 1
+
I~ I
"t
~'~ k
N51 1
h,kdb z : /-I
t
~'~ g,kzq z , "~1
/=l
(35)
or in an N × N matrix form, C Y = G,
(36)
where, as (28) in the constant element, Y contains all the unknowns, C is the coefficient matrix which is dense and nonsymmetric and the elements of G are known. The solution of equation (36) will give the values for either the scattered potential or radiation potential for the wave-diffraction or radiation problem, respectively. As for the constant element, the integrals in equations (33) and (34) are computed numerically, except when the collocation node i is inside the element under consideration. The nonsingular integrals hikz and gi~-z can be obtained by the use of Gaussian quadrature integration formulae and can be written as
01~'5:'(~11,,~.~Zq) IJ(~,,, G2q)1%, to'Z'
hik' • Z p
1 q t~t
g,~, = Z p
1
(37)
~'~
t~t
Z +*(~,,,, ~zq)IJ(~,,,, ~e,,)l ~,, ~q, I q
(381
I
where m is the total number of integration points, to, and tog are the weighting factors and (~i,, ~2q) are coordinates of the integration point.
3-D Wave diffraction
445
When the collocation point i belongs to the element under consideration, the kernels in equations (33) and (34) become weakly singular as the integration point approaches point i. Therefore, the integration over the vicinity of point i needs special treatment. Analytical calculation is formidable because of the presence of the Jacobian. However, these singular integrals can be computed by several different ways, which are discussed in Lee (1987).
Discontinuous element When continuous elements are employed, nodes are located at the element edges. When the element edge coincides with a discontinuity in the problem definition, modeling ambiguities will result. Clearly, at an edge of corner of the boundary in a three-dimensional problem, the normal to the boundary is not defined. The problem can be improved by using discontinuous elements, introduced by Patterson and E1-Sebai (1982) and Danson (1982), in which the interior nodes are used to define the boundary elements. Definitely, such elements present difficulties when maintaining inter-element continuity. Inter-element continuity may be advantageous for the boundary-element method but is not necessary. Indeed, the viability of constant elements demonstrates this fact. In the study by Patterson and Sheikh (1984), it was shown that a need for special strategies exists to cope with problems associated with the discontinuity of normal derivatives. The following linear interpolation functions that describe the variation of potential
04,
function 4, and its normal derivative 0n for discontinuous elements are used
M,(~,, ~ ) = ( ~ - ~ , ) ( ~ -
~2), M2(~l,~2)= (~+ ~,)(12-- ~2), (39)
M3(~,,~)= (12+~1)(12 + ~2), .M~(~,,~)= (~-~,)(~ + ~). The nodes are taken at the four interior points (see Lee, 1987). Higher-order interpolation functions can be found in Patterson and Sheikh (1984). This linear discontinuous element will be used to improve the result of the previous linear element and some case studies will be shown in the next section. NUMERICAL RESULTS AND DISCUSSION Simple examples were used to show the convergence rate of the direct boundaryelement method by using different types of elements and integration techniques. The direct method developed was also tested for accuracy and validity for a variety of rigid structures. Results show good agreement with other published data.
Accuracy and convergence of direct boundary-element method Added mass for the heave motion of a cube submerged in an infinite fluid was computed to demonstrate the convergence and accuracy of the direct boundary-element method. It was found that (1) when a four-point Gaussian integration technklue was
446
('ttlII-KANG Le~ and JACK Y. K. Lo;~ 0.75
~ 0.70 m
Const:ant element Linear element Discontinuous linear element
~-~
~ ~ ~
, Stellon and Mavis ~ ' ~ . _
0.85 --
~ 0*50
~
~°°~e
~ • • • • • • •
0.55 ~ iI
I I I I
I 0
~0
20
30
40
50
50
70
80
gO
loo
Number of elements per facet
Fro. 3.
Convergence rates of constant and linear elements.
used, an error ranging from 5 to 15% was introduced if the singularities of the integrals were simply disregarded. However, this error can be minimized for increasing the number of integration points. (2) The approximate method for evaluating the singularities as presented in Appendix A yields very good results even for relatively coarse meshes. (3) The convergence rate of the linear element is much slower than that of constant element (see Fig. 3). The reason for this is not yet known. It is probably due to the discontinuity of the normal derivative of the potential at the edges or corners of the boundary, since using discontinuous linear elements which lead to continuous normal derivatives has shown improvement in convergence rate. Detailed results are presented in Lee (1987) and will not be repeated here.
Radiation problems Added mass of a sphere in an infinite fluid. Again, both constant and linear elements were used to compute to heave added mass coefficients, while plane elements were used to approximate the curved spherical surface. The results are presented in Table 1. It was found that the constant element yields more accurate results. The errors were less than 3 and 1% respectively, when 24 and 40 constant elements were used in the computations. The results also show that considerable saving in computer time can be achieved by taking advantage of the symmetry of the problem.
3-D Wave diffraction TABLE 1.
447
HEAVE ADDED MASS COEFFICIENTS OF A SPHERE IN AN INIZlNITE FLUID BY USING CONSTANT AND LINEAR ELEMENTS WITH SYMMETRICAL PLANES
Number of elements per quarter
Cases
a33/pV
0.448 *
CPU
0.42
t 0.26
a33/pV 0.96
t 0.78
40
0.486
0.497
t 2.24
~ 1.00
*
28.5
0.438 $ 0.67
analytical result Note:
*
7.50
0.291 *
CPU
~ 0.21
24
~ 12.0
t 6.04
t 6.40
~ 2.57 0.441
$ 3.81
* 29.4
t 12.7
$ 8.90
0.5
a33/pV denotes added mass coefficient. Computing CPU based on the case a with 24 elements per quarter and two symmetrical planes as a unit. *denotes the case without plane of symmetry. tdenotes the case with one plane of symmetry. :]:denotes the case with two planes of symmetry. Case a uses constant elements. Case b uses linear elements.
A d d e d mass of a box in semi-infinite water. The box has a length and beam of 90 m with draft of 40 m. The water depth is infinite. The calculated added mass and damping coefficients of surge, heave and pitch motions are presented in Figs 4-6. The numerical results based on indirect boundary-element method and experimental results of Faltinsen and Michelsen (1974), are also presented in these figures for comparison purposes. Agreement between these two data sets is generally good, with the exception of damping coefficients in pitch and surge motions. This could be due to (1) too coarse mesh size on the free surface ($2) and (2) too short a distance between the radiation boundary ($4) and the structure. Diffraction problems Wave forces on a submerged hemisphere on the sea bottom.Horizontal and vertical wave forces on a submerged hemisphere of a radius a in a water depth of 3a were computed. The incident sea wave amplitude is ao. The nondimensionalized forces are shown in Fig. 7, where K is the wave number that satisfies the dispersion relation. Data presented by Chakrabarti and Naftzger (1974) are also included for comparison purposes. Good agreements for both horizontal and vertical forces were obtained by employing quadrilateral constant elements. Wave forces on the submerged storage tank and tower. Both the horizontal and vertical forces were calculated by using constant elements. The results are presented in Fig. 8. Comparison with the data obtained by Hogben and Standing (1974), who used
448
(?HIH-KANGLE~-;and JA(K Y. K. Lo~: 2.0
2.0
~
FalfJn~en
and
Mlchelsen
-~.
Faltinsen
and
Idlchel~en experiment
~-
Pre~ent
method
damping
added
=l > 1.oo q. "--~
mass
.~•
O. ~ I ~
~
~"~--~_
"~ ~
oo
,
•
~
~ ' ~
•• .,~*~
0.5--
-.
1.0
>
~.
•
-- 0.5
..
0.0
8
I 10
N
I 12
I
I 14 Period
FIG. 4.
~
,
16
I 18
0.0
20
(sec.)
Surge added mass and damping coefficientsfor a floating box (L x B×d=90m ×90m × 40m) in semiinfinite water.
the indirect boundary-element method, substantiates the validity and accuracy of the present method. Wave forces on a box in semi-infinite water. By using the same structure in the radiation problem, wave-exciting forces and moments on the floating box, which has length and beam of 90 m, at a draft of 40 m were calculated using constant elements. The results, together with the data obtained by Faltinsen and Michelsen (1974), who used indirect boundary-element method, are shown in Figs 9-11. Again, agreement is satisfactory.
CONCLUSIONS AND RECOMMENDATIONS The main purpose of this research is the development of a numerical method for studying fluid-structure interaction problems and to verify its applicability and accuracy. A variety of problems, including the wave-diffraction and radiation problems, were studied using the method developed. Most results obtained by this method agree with those obtained by previous investigators using other methods, including analytical, numerical or experimental. It would be of further interest to extend this method to study the wave resistance of
3-D Wave diffraction
449
2.0
2.0
~ -l. ~-
F~ltln~en and Idi~.hel~en Falt~;-,~,-; and Ml~hl~en experiment Pr~ent method ~
1.5
-,.0:l add~
8~ . . . .
•- ~ . . . ~ - , ~
mass
. . . . .
~.___
" ' " - . e . . .
e.o.
e
O.~m
-
0.S
damping e-...e...
~.
.-e
... o
0,0
0.0
8
10
12
14,
18
18
20
Period (sec.) FI~. 5.
Heave added mass and damping coefficients for a floating box infinite water.
(L×B×d=9Om×9Omx4Om) in
semi-
a surface ship with forward speed, although it may be difficult to express the radiation condition in moving coordinates. It should be pointed out that, in principle, any boundary geometries, no matter how complicated, can be handled by this method. Modeling the boundary surface by using quadratic elements is recommended for future studies. Not only do quadratic elements represent the curved boundary surface better, but they can also improve the variation of the unknown function and its normal derivative on the boundary. The boundary-element method is suited for solving problems that have a high ratio of domain volume to boundary surface area, especially, for problems having infinite or semi-infinite domains. However, special techniques are needed for solving a large number of full, nonsymmetric systems of equations. Further research is needed to obtain an in-depth undersianding of the present technique. Although Green's function for an infinite domain and the boundary-element method can be used conjunctively to solve problems with complicated boundaries, it would result in a large-sized system of equations. To solve these equations would require an excessive computer time. Thus, it would be much more advantageous to use the special Green's function for the boundaries under consideration in some tractible problems. For instance, for problems of wave diffraction and radiation in semi-infinite water, an
450
CHffI-KAN(; LEE and JACK ¥, K. Lou 0,20
~ -~.
~
0.020
F'a~n~en and MIchel~en Fa~Jn~e~ and M l ~ e l ~ e x
0.15 ~
-- 0.015
~ ~ .~ o.o-
"~
~1~ x'~,N,N, x 0.0~ ~
~ ~
~
~
~
x
~
....
0.~
"
~ ' ~ - ~ . ~
I
1
I
1
I
10
12
14
18
18
Pe~od FI~. 6.
"¢
~ --
o.~
-o.o,o I ~
~
.
o.~o
20
(~ec.)
Pitch added mass and damping coefficients for a floating box (L x B x d=90m x 90m x40m) in semi-infinite water.
1o0 , , |
~'~
~
Chakrabartl and Naf~gm"
0.8---2
~
I~
~
i,
~
~
0.4--
0.,-/]
¢__
~
"\
-_,
-~.~.,
0.0
0
0.0
0.1
0.2
0.3
0.4.
0.5
0.8
0.7
0.8
0.9
1,0
ka
FI6. 7.
The horizontal and vertical forces on a submerged hemisphere on the sea bottom.
3-D Wave diffraction
451
~,~
- - 0.3
l•
o.,,~--t ,',
~
/,'
o
~
",~'~,.
w~, / \ °=1 ~i,,,,, ::=~.:,,,~~,, ..tZo.~ ~ ~. / l'/,' ',',,., "X .~L.J, o=1/ ":,x,....: ~;;~,,,,, 1; ',",//
I , , ~,
-- 0.2
':': ~~-.. ._
o..
--0.1
~-~'a ~,,~ ~,
~.~
[
I
I
1
~
~
~
3
0.0
k~ FI~. 8.
The horizontal and vertical forces on the submerged tank and tower.
- - F o l t l n l e n 0nd Mlchell~n -- - I~m~t rn~tho4
2.0~
1.5--
/
.,.
o
~o+ 1.0-~.
/ / s
0.5~
0.0 10
12
14 Pe~
FI6. 9.
Surge wave force amplitudes for a floating box
18
1~
(se~.)
(LxBxd=9Omx9Omx4Om)in
semi-infinite water.
452
CItItt-KANG LEE and J,,~cK Y. K, L,otr
~
FaltJn~en
and Mlchelsen
- - -- Prewent m e t h o d
2.0--
1.5 o
~ '~ I "~°1"° ~.
--
0.5--
0,0
I
10
I
12
14. Period
FIG. 10.
16
18
20
(S~.)
Heave wave force amplitudes for a floating box
(LxBxd=9Omx9Omx4Om) in semi-infinite water
0.20 - -
0.15 --
/
:°oo,o_
I / / /
/
8
10
- -
Fdtln~
~-
P~t~ent m e t h o d
12
14. Period
FIG. 11.
~nd ~ l ¢ h ~ t ~ n
18
18
20
(S~.)
Pitch wave moment amplitudes for a floating box water.
(L×Bxd=9Omx9Omx4Om) in
semi-infinite
3-D Wave diffraction exact expression for the Green's
f u n c t i o n s o l u t i o n c a n b e f o u n d in W e h a u s e n
453 and
L a i t o n e (1960) as
u* = -1 + P V f f k + ~
e k~+~) J o ( k R ) d k + i 2"rrl~e~ + t ' ) Jo(~xR).
T h e first a u t h o r has e m p l o y e d this s o l u t i o n to s o l v e s o m e o f t h e test c a s e s c i t e d p r e v i o u s l y . P r e l i m i n a r y r e s u l t s i n d i c a t e t h a t significant i m p r o v e m e n t c a n b e a c h i e v e d ; for example, the surge and pitch dampings can be more accurately predicted with a s m a l l e r - s i z e d s y s t e m o f e q u a t i o n s as c o m p a r e d to t h e p r e v i o u s results.
REFERENCES BREB~IA, C. A. 1978. The Boundary Element Method for Engineers. Pretech Press, London. BREBBIA, C. A. 1981. Progress in Boundary Element Methods-1. Wiley, New York. CnAK~,~SAR~, S. K. and NAVrZ~ER, R. A. 1974. Non-linear wave forces on half-cylinder and hemisphere. J. Watways Harbors Coastal Eng.ng. Div., Am. Soc. Chem. Engrs 1~, 189-204. DANSO~, D. J. 1982. A boundary element analysis system. In Boundary Element Methods in Engineering. Edited by B~rss~A, C. A., pp. 557-575. Springer, Berlin. F,~LW~SrN, O. M. and MICnrLS~, F. C. 1974. Motions of large structures in waves of zero Froude number. Proc. Int. Syrup. Dynamics of Marine Vehicles and Structures in Waves, pp. 99-114. Institute of Mechanical Engineering, London. F ~ K , W. 1967. Oscillation of cylinder inor below the free surface of deep fluids. Report No. 2375. Naval Ship Research and Development Center, Bethesda, MD. G~lSON, C. J. 1974. Hydrodynamics of large objects in the sea, part I--Hydrodynamic analysis. J. Hydronautics 8, 5-12. GAte,SOY, C. J. 1978. Chapter 3---Hydrodynamic loading of large offshore structures: three-dimensional source distribution methods. In Numerical Methods in Offshore Engineering. Edited by ZIEr~K~EW~CZ,O. C., LEWIS, R. W. and Sv~,~, K. G., pp. 87-140. Wiley, New York. HAVELOCK,T. H. 1955. Waves due to a floating sphere making periodic heaving oscillations. Proc. R. Soc. Lond. 230A, 1-7. HESS, J. L. and SranrH, A. M. O. 1962. Calculation of nonlifting potential flow about arbitrary threedimensional bodies. Report No. E.S. 40622. Douglas Aircraft Division, Long Branch, CA. Ho~aEr~, N. and S~A~D~G, R. G. 1974. Wave loads and large bodies. In Proc. Int. Symp. Dynamics of Marine Vehicles and Structures in Waves, pp. 258-277. Institute of Mechanical Engineering, London. JoHr~, F. 1950. On the motion of floating bodies II. Commun pure appl. Math. 3, 45-101. KELLO6~, O. D. 1929. Foundations of Potential Theory. Springer, Berlin. LEE, C. K, 1987. A coupled boundary and finite-element method for three dimensional fluid-structure interaction problems. P h . D . dissertation, Texas A&M University, T. LEW~S, F. M. 1929. The inertia of water surrounding a vibrating ship. Trans. SNAME 37, 1-20. NEWMAN, J. N. 1978. Marine Hydrodynamics. The MIT Press, Cambridge, Mass. PAttErSON, C. and EL-SEaA~, N. A. S. 1982. A regular boundary method using non-conforming elements for potential problems in three dimensions. In Boandary Element Methods in Engineering. Edited by BRE~mA, C. A., pp. 137-152. Springer, Berlin. PA~E~SOt~, C. and SHEIKH, M. A. 1984. Interelement continuity in th eboundary element method. In Topics in Boundary Element Research--1. Edited by B~,EBa~A,C. A., pp. 123-141. Springer, Berlin. SrELSOS, T. E. and MAVIS, F. T. 1955. Virtual mass and acceleration in fluids. Proc. ASCE 81, 670. Sxaot~), A, H. and SECaESX,A. H. 1966. Gaussian Quadrature Formulae. Prentice-Hall, New York. URSELL, F. 1949a. On the heaving motion of a circular cylinder in the surface of a fluid. Q. J. Mech. appl. Math. 2 (2), 218--231. URSELL, F. 1949b. The rolling motion of a circular cylinder in the surface of a fluid. Q. J. Mech. appl. Math. 2 (3), 335-353. WEnAUSE~, J. V. and LAIrO~r, E. V. 1960. Surface waves. Encyclop. Phys. 9, 446-778.
454
CHIH-KANG LEE and JACK Y. K. Lou node p(x,y,z)
' [ '
c~, ~ ) l
'Z:__%_ ~
l/0.w ,,,"' --/
--,,,,-,=,,-,---I~ X
,~
(~2,r~2)
Fro. A . I .
DERIVATION
OF
Local element system.
APPENDIX A APPROXIMATE FORMULAE COEFFICIENTS
FOR
INFLUENCE
Consider a potential point P(x, y, z) and a quadrilateral source element in the xy-plane, as shown in Fig. l A. 1. The integrand, i.e. 3' in coefficient Gij in equation (33), can be expanded in a power series in ~ and -q about the origin as 1 r
l
V/~ - ~)2 + (y - Xl)2 + z2
(A.l) l 1 = ~ - w..{ - ~o.~n + ~ to.. ~" + ~..,, { n + ~ ~o.,..,,n-',
where
tO =
1 r( i
r. = ~r 2 + y2 + z 2. It can be show that an approximation for the influence coefficient G,, can be given as G
1
3-D Wave diffraction =
1
A to - (o.,~,, M,~ - o%, Ms) + ~ (to.~ I~,~ + 2 o.~r l~y + tO.ry Iyy)
455
1
(A.2)
,
where
Z:I0 , lyy = f ~q2d~d~.
Ly=fK~dKd~,
In equation (A.2), the first term can be interpreted as the potential of a point source of strength A. The terms in the first parathesis are the dipoles, the axes of which are along the x-axis and y-axis, respectively. The remaining terms are quadrupoles. In the present application, since the origin of the coordinate system is at the centroid of the element, all the first moments, M, and My, vanish; thus, only the source and quadrupole terms exist. 1 From equation (A.1), we can obtain the normal derivative of the fundamental solution, , as r
1
On - - t ° ~ - ° ' ~ z ~ - ° > ' ~ + ~
°~xx~2
+OxYz~+~Oyyz~2"
(A.3)
Thus, we have
H,~ = 4~r
1[
dr
1
l
= 4 ~ A o , - (o.~z Mx - o r, M r ) + } (oxxz lxx + 2 ~Oxyz lxy + o~yz !,,) •
(A.4)
As in equation (A.2), only the first and third terms exist, corresponding to the source and quadrupole. Equations (A.2) and (A.4) can be used in two forms. In the first, both the source and quadrupole are used for intermediate distant points. In the second, only the source is used for a far distant point. According to the study by Hess and Smith (1962), the form of the approximation used for a particular point depends on the ratio of the distance between the point and the centroid of the quadrilateral to the major diagonal of the quadrilateral. If this ratio is less than 2.45, we use the exact method; if this ratio is greater than four, the second form is adequate; otherwise the first form may be employed. In addition, Hess and Smith's results showed that the maximum error owing to the approximate forms is less than 0.1%. Thus, the approximate forms are quite successful in computation without losing accuracy and these forms were used in the present work to save computing time.