- Email: [email protected]
A B-spline-based method for radiation and diffraction problems
ARTICLE IN PRESS
Ocean Engineering 33 (2006) 2240–2259 www.elsevier.com/locate/oceaneng
A B-spline-based method for radiation and diffraction problems Ranadev Datta, Debabrata Sen Department of Ocean Engineering and Naval Architecture, IIT-Kharagpur, Kharagpur 721302, India Received 18 May 2005; accepted 4 December 2005 Available online 9 March 2006
Abstract An open uniform B-spline-based panel method is developed for solution of potential flow problems. In this method, both geometry as well as the field variables are represented by the same open uniform B-spline basis function. The method is initially applied for the radiation problem in unbounded fluid. Computed results for a spheroid of different aspect ratio are found to be in excellent agreement with analytical results. The method is then applied for diffraction problem formulated based on the transient (time-domain) Green’s function. Computed results for a hemisphere and Wigley hull are compared with published results and the comparison shows good agreement. r 2006 Elsevier Ltd. All rights reserved. Keywords: B-spline; Basis function; Panel method; Radiation problem; Diffraction problem
1. Introduction Computation for the motion of a ship when it traverses through a wave field and the corresponding force it experiences is of fundamental interest in marine hydrodynamics. Over the years, several researchers have studied this problem. Most of the earlier methods were based on two-dimensional strip theory. Over the recent past, more advanced three-dimensional (3D) methods are being developed. Corresponding author. Fax: +03222 282284.
E-mail addresses: [email protected] (R. Datta), [email protected] (D. Sen). 0029-8018/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.oceaneng.2005.12.002
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2241
A majority of these solutions are based on panel methods in some form or other. Based on surface description and unknown function representations, there can be different versions of panel methods which can be generally referred to as lower and higher order panel methods. In a lower (0th or 1st order) panel method, planar panels are used and the flow variables (potential, normal velocity and/or singularity strengths) are assumed either constant or linearly varying over individual panels. In the higher order panel methods, some higher order functions like quadratic or cubic Lagrangian interpolation are used to describe the surface geometry as well as the unknown field variables. Higher order methods are in general relatively more efficient and accurate compared to the lower order methods, but often these methods suffer from some undesirable characteristics. For example, in the higher order method of Hess (1979), use of cubic panel shape with quadratic source density seems to result in an extremely complicated algorithm. Also, the higher order method using Lagrangian interpolation function gives discontinuous normal vector and spatial derivatives of the potential across the panel boundaries. Most recently, B-spline based panel methods are being developed for wave– structure interaction and free-surface flow problems. Some of the relevant works in this context are as follows. Hsin et al. (1994) developed a 2D higher order panel method based on B-splines in which the geometry and singularity distribution were defined by B-spline and Green’s theorem was then applied at selected collocation points to solve for the unknown flow parameters. Maniar (1995) used B-splines to develop a higher order 3D panel method, where the potential and the geometry of the body are allowed any degrees of continuity. The developed method is then successfully applied for the zero-speed frequency-domain radiation–diffraction problem. Kim and Nam (1999) proposed a B-spline panel method where the surface is described by NURBS which provides for more accurate geometry modeling, specially for the conic surface. Zhao and Zou (1999) described the surface as well as unknown potential by NURBS to study the potential flow around the body with/ without free surface. Chen et al. (2001) proposed bi-quadratic patch method to predict the wave-induced ship motions and load in frequency domain. Lee et al. (1998) and Danmeier (1999) have presented a geometry independent higher order method in which the geometrical and hydrodynamic representations are decoupled. Here the velocity potential is represented by a set of B-spline basis functions, but the geometry can be modeled by any regular parameterization. In this way the hydrodynamic problem is completely separated from the geometric model. Qiu and Hsiung (2002) introduced a panel free method and applied it to solve the radiation problem in time domain, and subsequently extended the method for the diffraction problem as well (Qiu et al., 2004). In this method, the body boundary was represented by NURBS surface and the integral equations based on source-strength distribution were globally discretized over the body surface by Gaussian quadrature. Due to a de-singularization procedure adopted, there was no need to express the source strength and velocity potential by spline-representation. As can be seen from above, there have been quite some interest in developing Bspline-based solution schemes for a variety of wave–structure and free-surface flow problems. The motivations of these developments are that (i) B-splines provide an
ARTICLE IN PRESS 2242
R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
accurate description of the geometry. Such B-spline-based geometry modeling are presently used extensively in CAD application for ship hull design and thus the geometry description from such hull-development modules can be directly used for the hydrodynamic computations. (ii) A substantial reduction in computational efforts (time as well as required memory space) compared to lower order methods for the same level of accuracy can be achieved. One of the important problems in ship hydrodynamic in which (ii) above becomes very significant is the solution of the 3D forward speed ship motion problem in timedomain using the transient Green function. Lin and Yue (1990) have used first-order panel method to solve this problem, but the method is extremely time consuming due to the convolution nature of the integral and the large number of panels needed for accurate description of the hull. For example, Bingham et al. (1996) studied the motions of a ship in waves using a version of the 3D transient Green function method and applying a low-order panel method for the discretization of the integral relations. For a series 60, C b ¼ 0:7 hull, computation carried out with 256 unknowns and upto 300 time-steps required 22 CPU hours on the entry-level SGI Indigo (MIPS 3000 processor). The memory requirement was very high. In fact, the memory space required in these methods are of the order OðN 2p N t Þ and the computation effort is OðN 2p N 2t Þ, where Np is the number of panel and Nt is the number of time steps. For practical computations, it may be required to run a simulation for substantially higher number of time steps, in particular if behaviors in irregular waves are to be studied. Clearly, too much of computation time and memory is involved in the process: a typical problem of this type for a realistic ship geometry may require a duration of the O(15) hours in todays PC for the solution by lower order panel methods. Motivation for developing a B-spline-based method for this application is therefore substantial. In the present work, an open uniform B-spline-based method is developed with the final goal of applying this method for the 3D forward speed ship motion problem based on transient Green function. A lower order panel method for this problem was developed earlier (Sen, 2002). Despite several developments on B-spline-based methods for the general wave–body interactions, there appear to be a lack of application of B-spline-based methods for the 3D forward speed ship motion problem. In order to check the performance of the basic spline-based integral-equation solver, initially the radiation problem for a fully submerged body is considered. Here the problem is formulated based on a mixed source-sink distribution, and both geometry and field variables (velocity potential and normal velocity) are represented by the same open-uniform B-spline basis functions. Numerical integration is then performed at selected collocation points. Results in terms of added mass for some regular geometric shapes are found to be in excellent agreement with analytical solutions. Next the solution of the zero-speed diffraction problem of a surfacepiercing floating body is considered, where the problem is formulated based on transient Green function, and the integral equation is based on source distribution (as in Lin and Yue, 1990). For evaluation of the Rankine part (the singular part) of the integration, the adaptive subdivision scheme proposed by Telles (1987) is applied
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2243
for near field, while Gaussian integration is performed for far field. The regular (memory) part of the Green’s function is integrated using a 4-point quadrature rule. Trapezoidal rule is applied for the convolution (time) integration. The computed results for diffraction forces for Wigley hull are compared with published results, and are found to be in good agreement. In the following, we first discuss the radiation problem in unbounded fluid, followed by the zero-speed diffraction problem for a ship hull. (1) Radiation problem: For the radiation problem of a body in unbounded fluid, the potential function Fk ðp; tÞ satisfies the Laplace equation and the boundary conditions, as follows: r2 Fk ðp; tÞ ¼ 0 in O; ðFn¯ Þk ¼ n¯ k for k ¼ 1; 2; 3 on Sb , ðFn¯ Þk ¼ ð¯r n¯ Þk3 for k ¼ 4; 5; 6 on
(1.1)
Sb .
ð1:2Þ
In the above, k ¼ 1; 2; 3 refer to linear modes of motions (surge, sway, heave), k ¼ 4; 5; 6 are the rotational modes of motions (roll, pitch, yaw), n¯ ¼ ðnx ; ny ; nz Þ, r¯ ¼ ðrx ; ry ; rz Þ and Sb is the body surface. It is well known that from application of Green’s theorem, the potential function can be represented by a distribution of sources and strengths over the boundary surface as follows: ZZ qGðp; qÞ qFðqÞ aðpÞFðpÞ ¼ FðqÞ Gðp; qÞ dS. (1.3) qnq qnq S The above integral relation represents potential at a field point p by distributing sources G of local source strength qFðqÞ=qn and the dipole qGðqÞ=qn of moment FðqÞ at the point qðx; Z; zÞ on the surface S of the fluid domain. G here is the Green’s function, which for the present problem can be taken as a simple source, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G ¼ 1=R ¼ 1= ðx xÞ2 þ ðy ZÞ2 þ ðz zÞ2 . Also, aðpÞ is the solid angle on the surface at a point pðx; y; zÞ. For an unbounded domain, S reduces to Sb, the body surface, since the contribution over the infinite outer surface vanishes. q=qn represents the normal derivative and the subscript q denotes that this derivative is to be taken at the ‘source point’ q. In order to solve the problem using B-splines, the body surface Sb can be divided into number of large patches. Over each patch, the surface as well as the field variables F, qF=qn are represented by open uniform B-spline basis functions as follows. Let X¯ ¼ ðx; y; zÞ are the Cartesian coordinates. These can be written parametrically as x ¼ X ðu; wÞ, y ¼ Y ðu; wÞ, z ¼ Zðu; wÞ where u,w are the parameters. Any point X¯ ðpÞ on the surface can then be described as X¯ ðpÞ ¼
nþ1 m þ1 X X i¼1 j¼1
Ax¯ ði; jÞN ki ðuÞM lj ðwÞ.
(1.4)
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2244
Here N ki ðuÞ and M lj ðwÞ are the open uniform B-spline basis functions of order k and l, respectively. The same representation is applied over all the patches describing the full surface. In (1.4), Ax¯ ði; jÞ ¼ Ax ði; jÞ; Ay ði; jÞ; Az ði; jÞ are called the vertices of the defining polygon net, also known as control points. ðn þ 1Þ and ðm þ 1Þ are the total number of control points along u and w parametric directions, respectively. The surface is thus described by ðn þ 1Þ ðm þ 1Þ number of control points. Once the basis functions are specified, depending on the order of the spline used, the problem of representing the surface reduces to finding Ax ði; jÞ; Ay ði; jÞ; Az ði; jÞ. ðn þ 1Þ ðm þ 1Þ number of equations are needed to solve for these unknowns, These can be obtained by simply making the surface go through that many collocation points, i.e by fitting a B-spline surface passing through a set of known points. The collocation points in the parametric space have been determined based on chord line approximation (see Rogers and Adams, 1990). Specifically, for r ¼ n þ 1 physical data points in the u parametric direction, the parameter value at the lth data point is u1 ¼ 0;
ul umax
Pl g¼2 jDg;s Dg1;s j ; l ¼ 2; n. ¼ Pr g¼2 jDg;s Dg1;s j
(1.5)
Similarly, the lth data point in the w parametric direction is Pl wl g¼2 jDr;g Dr;g1 j w1 ¼ 0; ; ¼ Ps wmax g¼2 jDr;g Dr;g1 j
l ¼ 2; m
(1.6)
In the above, umax and wmax are the maximum value of the appropriate knot vectors. D is defined by Dp;q ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxp xq Þ2 þ ðyp yq Þ2 þ ðzp zq Þ2 .
(1.7)
The unknown potential F and normal velocity Fn are expressed by the same basis functions: F¼
nþ1 X mþ1 X
Af ði; jÞN ki ðuÞM lj ðwÞ,
(1.8)
i¼1 j¼1
Fn ¼
nþ1 X mþ1 X
Afn ði; jÞN ki ðuÞM lj ðwÞ.
(1.9)
i¼1 j¼1
If the surface is described by Pn number of patches, then the discretized form of Eq. (1.3) for the sth collocation point (i.e, in parametric domain, for
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2245
u ¼ us1 ; w ¼ ws2 ) is given by " ( Pn nþ1 X mþ1 X X aðps Þ Af ðr; i; jÞN ki ðus1 ÞM lj ðws2 Þ r¼1
i¼1 j¼1 nþ1 X mþ1 X
Z
f
u¼0
Pn Z umax X r¼1
Z
wmax
A ðr; i; jÞ
i¼1 j¼1
¼
umax
Z
u¼0
wmax
w¼0
w¼0
nþ1 m þ1 X X
fn
A
)#
N ki ðuÞM lj ðwÞG n ðus1 ; ws2 ; u; wÞjJj du dw !
ðr; i; jÞN ki ðuÞM lj ðwÞ
Gðus1 ; ws2 ; u; wÞjJj du dw.
i¼1 j¼1
ð1:10Þ The normal n¯ can be found from the following relations: X¯ u ¼
nþ1 X mþ1 X
k
Ax¯ ði; jÞN 0 i ðuÞM lj ðwÞ,
i¼1 j¼1
X¯ w ¼
nþ1 m þ1 X X
l
Ax¯ ði; jÞN ki ðuÞM 0 j ðwÞ,
(1.11)
i¼1 j¼1
n¯ ¼
X¯ u X¯ w . X¯ u X¯ w
The Jacobian is given by jJj ¼ X¯ u X¯ w .
(1.12)
(1.13)
fn
In (1.10), A in the right-hand side is determined directly from application of the body boundary condition (1.2), ðFn¯ Þk ¼ n¯ k , k ¼ 1; 2; 3; ðFn¯ Þk ¼ ð¯r n¯ Þk3 , k ¼ 4; 5; 6. To obtain a solution for the unknowns Af , u and w are divided into n m points resulting into ðn þ 1Þ ðm þ 1Þ rectangular segments in the parametric domain. Each segment corresponds to a surface panel in the physical domain. Collocation method is used to get the linear system of equations with ðn þ 1Þ ðm þ 1Þ unknowns. This results into a system of nn ¼ ðn þ 1Þ ðm þ 1Þ linear equation: ½BfAf g ¼ fCg,
(1.14)
where [B] is a nn nn matrix of influence coefficients containing the integrals over individual panels in the parametric domain as given in the left-hand side of (1.10). All integrations are carried out using a 3 3 Gaussian quadrature scheme. However, when the integration is performed over the so-called self-influencing panels, the integrand contains a singularity. There are various schemes available to handle such a singularity. In the present work, we use the adaptive subdivision scheme proposed by Telles (1987). In this, the self-influence panel is subdivided successively to isolate the singular point, and a 3 3 Gaussian quadrature rule is applied to the subdivided parts of the panel. For the small subdivision part of the surface over which the singular point lies, analytical integration is performed, assuming that the surface is
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2246
Table 1 Absolute relative error No of subdivision
Relative error
5 10 12 15 20 30
0.085400 0.002776 0.000692 0.000085 0.000002 0.000000
(0,1)
(1,1)
w
(0,0)
u
(1,0)
Fig. 1. Subdivision of panel.
planar. The subdivision is continued till the integrand converges within a relative absolute accuracy of 0:00005 for a unit square patch in the parametric domain. Relative absolute accuracy is defined as jðas ac Þ=ac j where as is the computed value at the sth level of subdivision and ac is the value at a high level of subdivision where the result is assumed to have converged (e.g. 30 subdivisions as per Table 1 above). It is found that usually about 15–20 subdivisions are sufficient to achieve this accuracy level. The contribution of the small isolated panel to the integral is extremely small, and can even be ignored without any difference for this level of accuracy. This subdivision scheme is shown in Fig. 1 and illustrative results are presented in Table 1. Variety of methods exists to solve the system of Eqs. (1.14). Once the unknowns Af’s are found, F at any desired location can be determined directly from (1.8). The results obtained for this radiation problem are presented in terms of the familiar added masses (mi;j ) defined by ZZ mi;j ¼ r fi nj ds; ðm1;1 for surge; m3;3 for heaveÞ. (1.15) dO
In discretized form, this can be written as " # ! pn Z u¼umax Z w¼wmax X nþ1 m þ1 X X f k l mi;j ¼ r A ði; jÞN i ðuÞM j ðwÞ n¯ :jJj du dw . r¼1
u¼0
w¼0
i¼1 j¼1
(1.16) where the normal n¯ and Jacobian jJj can be determined from (1.12) and (1.13).
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2247
2. Results Numerical computations are performed for a spheroid of different aspect ratio, since for this geometry the analytical results are available. Due to symmetry, only one octant of the spheroid is needed to be modeled as one patch. We use 5 5 ¼ 25 control points to represent this patch. Fig. 2 shows the geometry of one octant of the spheroid and the full body. Tables 2 and 3 present results of non-dimensional heave and surge added masses (m ¯ i;j ¼ mi;j =M, M ¼ mass of the body) for the range of aspect ratios a=b ¼ 0:221. It
Fig. 2. Geometry of spheroid.
Table 2 Added mass for heave b=a
Numerical
Analytical
% of error
Abs error
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.91070 0.82350 0.76080 0.70390 0.65320 0.60750 0.56730 0.53130 0.49990
0.89426 0.82587 0.76189 0.70421 0.65294 0.60757 0.56741 0.53176 0.50000
1.838 0.272 0.143 0.044 0.039 0.011 0.019 0.086 0.020
0.0166 0.0023 0.0011 0.0003 0.0003 0.0001 0.0001 0.0005 0.0001
Table 3 Added mass for surge b=a
Numerical
Analytical
% of error
Abs error
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.06000 0.10520 0.15630 0.21080 0.26570 0.32280 0.38120 0.44040 0.49950
0.05912 0.10542 0.15626 0.21002 0.26576 0.32294 0.38120 0.44028 0.50000
1.488 0.209 0.026 0.371 0.022 0.043 0.000 0.027 0.100
0.0009 0.0002 0.00004 0.0007 0.00006 0.0001 0.0000 0.0001 0.0005
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2248
0.9
3D B-spline Analytical results
m33/ (4/3)πρab2
0.8
0.7
0.6
0.5 0.2
0.4
0.6 b/a
0.8
1.0
Fig. 3. Heave added mass.
may be noted that heave added mass for an aspect ratio a=b is equivalent to surge added mass for the aspect ratio b=a. This means, computational results presented cover linear added masses for a range of aspect ratio a=b ¼ 0:225, i.e. the geometry is varied from a shallow oblate spheroid to a fairly long prolate spheroid. The analytical values as well as the percentage and absolute errors are also shown. The results are plotted in Figs. 3 and 4. As can be seen, the results are in excellent agreement and are graphically indistinguishable, except marginally at the ends. The agreement somewhat deteriorates as the body gets elongated. However, the differences are still less than 0.4% up to 0:3pa=bp3:33. Only for the last case of a=b ¼ 0:2 (or equivalently a=b ¼ 5), the error is somewhat large at around 2%. This may be due to the fact that for elongated shape, a single patch may not be sufficient to represent the geometry with the required level of precision. For such cases, a larger number of patches may be used to improve the results. (2) Diffraction problem: We next consider the diffraction problem of a floating free surface piercing body, formulated based on transient Green function in the time domain. Consider an arbitrary 3D body floating on the free surface of an incompressible ideal fluid. Let Oxyz be the right-handed co-ordinate system with Oxy plane coincident with the mean water level z ¼ 0. The z-axis is positive upwards, and sinusoidal (linear) waves traveling in a direction making an angle b to x-axis is incident on the body. The initial boundary value problem for the diffraction potential Fd ðp; tÞ is specified by the following governing equation, boundary
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2249
0.5 3D B-spline Analytical results
m11/ (4/3)πρab2
0.4
0.3
0.2
0.1
0.2
0.4
0.6 b/a
0.8
1.0
Fig. 4. Surge added mass.
condition, initial conditions and far field conditions: r2 Fðp; tÞ ¼ 0 in O;
(2.1)
q2 F qF ¼ 0 on z ¼ 0, þg qt2 qz
(2.3)
qFd qFi ¼ qn qn
(2.4)
on
S0
Fd ! 0 at t ! 0,
(2.5)
qFd !0 qt
at t ! 0,
(2.6)
rFd ! 0
as r ! 1.
(2.7)
In the above, Fi is the incident wave potential, S0 is the mean wetted body surface, n¯ is the normal on S0 directed outwards to the fluid domain (inward to body), and g is the acceleration due to the gravity. The solution of the above initial boundary-value problem is based on the transient Green function approach of Lin and Yue (1990). The diffracted potential can be expressed by the following convolution integral by distributing transient Green’s
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2250
function of unknown strength over the body surface, 8 2 39 > > ZZ Z t ZZ = < 1 6 7 Fd ðp; tÞ ¼ sG0 dS þ dt4 sG ft dS5 , > 4p > 0 ; : S 0 ðtÞ
(2.8)
S 0 ðtÞ
where q ¼ qðx; Z; zÞ represents a ‘source point’ on the body surface and sðq; tÞ is the associated source strength at time step t. G is the transient Green’s function given by Gðp; t; q; tÞ ¼ G 0 þ Gf
for paq; t4t
(2.9)
where 1 1 ðRankine partÞ, r r0 Z 1 pffiffiffiffiffiffi f ½1 cos gkðt tÞ ekðzþxÞ J 0 ðkRÞ dk ðfree surface memory partÞ, G ¼2 G0 ¼
0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ðx xÞ2 þ ðy ZÞ2 þ ðz zÞ2 , r0 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xÞ2 þ ðy ZÞ2 þ ðz þ zÞ2 ,
R¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xÞ2 þ ðy ZÞ2 ,
k¼
o2 g
is the wave number,
J 0 ¼ Bessel function of order zero.
(2.10)
Eq. (2.9) satisfies all conditions except the body boundary condition. Application of the body boundary condition to (2.8) yields the following integral relation: 9 8 > > ZZ Z t ZZ = < qFd ðp; tÞ 1 f 0 (2.11) ¼ sGnp dS dt sG np t dS . > qnp 4p > 0 ; : S 0 ðtÞ
S 0 ðtÞ
Detail derivation of the above is available in several sources, e.g. Lin and Yue (1990). Once the diffracted potential is obtained, pressure on the body can be determined from Bernoulli equation, and the linear wave forces and moments can be easily obtained by integrating pressure over the body surface: 2 3 ZZ ZZ qFi qFd 6 7 F i ðtÞ ¼ r4 dS þ dS 5, (2.12) ni ni qt qt S 0 ðtÞ
S 0 ðtÞ
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2251
where ni is the generalized normal and ¼ ni for i ¼ 1; 2; 3 and ¼ ð~ r~ nÞi3 for i ¼ 1; 2; 3. For the numerical solution, as in the radiation problem earlier, an open uniform B-spline basis function is used to describe the geometry and field variables. The surface is represented by Eq. (1.4) while, for an instant of time t, the unknown potential Fd and source strength s are defined as Fd ¼
nþ1 m þ1 X X
Af ðt; i; jÞN ki ðuÞM lj ðwÞ,
(2.13)
i¼1 j¼1
s¼
nþ1 m þ1 X X
As ðt; i; jÞN ki ðuÞM lj ðwÞ,
(2.14)
i¼1 j¼1
where NðuÞ and MðwÞ are the open uniform B-spline basis functions. It may be noted that this formulation is based on source distribution unlike the radiation formulation based on a mixed source and dipole distribution, and thus here the unknown source strengths are expressed by B-spline basis functions. Due to the convolution nature of the integrals, the unknown control nets Af and As are also functions of time. The convolution integral in (2.11) is treated by discretizing time into equal time-steps. Thus, (2.11) in discretized form for the Mtth time step and for sth collocation point can be written as " Pn X nþ1 Mþ1 X X As ðr; M t ; i; jÞ r¼1
Z
i¼1 J¼1 umax u¼0
Z
wmax
w¼0
# N ki ðuÞM lj ðwÞG 0np ðus1 ; ws2 ; u; wÞjJj du dw
qF ¼ 4p qnp " n (Z t M 1 P X X er 1 Dt r1 ¼1
r¼1
umax
u¼0
)#
G fnp ðu; wÞjJ j du dw
,
Z
wmax
w¼0
nþ1 m þ1 X X
! A
s
ðr; r1 ; i; jÞN ki ðuÞM lj ðwÞ
i¼1 j¼1
ð2:15Þ
where u ¼ us1 ; w ¼ ws2 are the parameters representing the sth collocation point, and t ¼ M t Dt, Dt is the time step size. Applying Eq. (2.15) to selected collocation points as in previous case, a linear system of equations for the unknowns As for time step M t can be formed: ½BfAs ðM t Þg ¼ fDg.
(2.16)
Matrix [B] consists of the integration of Rankine part G 0 of the Green’s function, which has a singularity when integrated over self-influence panels, as in the previous case. This is treated in the similar manner as in the radiation problem, by applying
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2252
a 3 3 quadrature rule for the non-singular part and using the adaptive subdivision for the integral containing singularity. For the integrals containing the regular memory part G f of the Green function, we use a 2 2 quadrature rule since this is found adequate. The convolution time integration, as shown in (2.15), is based on a simple trapezoidal rule. Efficient and accurate computation of G f and its derivatives follow the method described in Sen (2002). Once the unknown As for a time step is determined from a solution of (2.16), the unknown Af for that time step can be found from (2.8). The discretized form of (2.8) for the Mtth time step and sth collocation point (i.e, in parametric domain at u ¼ us1 ; w ¼ ws2 ), is " # pn X nþ1 m þ1 X X f t k l A ðr; M ; i; jÞN i ðus1 ÞM j ðws2 Þ r¼1
i¼1 j¼1
" n Z ( ) # Z p nþ1 X mþ1 1 X umax wmax X As ðr; M t ; i; jÞN ki ðuÞM lj ðwÞ G 0 ðus1; ws2; ; u; wÞjJj du dw 4p r¼1 u¼0 w¼0 i¼1 j¼1 " n (Z ! )# t M 1 P nþ1 X mþ1 umax Z wmax X X X s f k l e r1 A ðr; r1 ; i; jÞN i ðuÞM j ðwÞ G np ðu; wÞ du dw . Dt ¼
r1 ¼1
r¼1
u¼0
w¼0
i¼1 j¼1
ð2:17Þ Once again, applying (2.17) successively at the collocation points, one gets the following system of equations: fAf ðM t Þg ¼ ½CfAs ðM t Þg.
(2.18)
The velocity potential F can be now determined in a straightforward from (2.13). Time derivative of this potential needed for evaluating pressure is found from numerical differentiation, and forces, moments are evaluated from (2.12). 3. Results Computations are first performed for a floating hemisphere due to availability of results for this geometry. Considering flow symmetry, only one quadrant of the surface needs to be modeled. Computations are carried out using both (5 5) and (7 7) control points to define the quadrant. Fig. 5 shows the surface as modeled by (5 5) control points. Fig. 6 shows the total surface of the hemisphere (using symmetry). Computations cover a range of frequencies and the results are presented in terms of force amplitudes in the frequency domain. In all cases, the time history is extremely regular (a typical time history is shown later for the next computation). The force amplitudes are determined by averaging the peak-to-peak values over the steady part of the time-history. Fig. 7 shows the non-dimensional heave diffraction force amplitude jX 3 j=rgpR2 against non-dimensional frequency kR, where k is the wave number and R is the radius of the hemisphere. Fig. 8 shows the nondimensional sway diffraction force jX 2 j=rgpR2 against kR (sway here is defined along the direction in which the incident wave travels). Results are compared with
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2253
Fig. 5. One octant of a hemisphere.
Fig. 6. Hemisphere.
the results of Beck and King (1989), results obtained by Qiu et al. (2004) using their panel-free-method (PFM), and results obtained using Haskind relations and by Cohen. These last two results are as reported in Qui et al. (2004). As can be seen, computations with (5 5) and (7 7) points produce almost identical results. This can also be observed in Table 4. Convergence appears to be best at the mid-range of the computed frequency-range with some minor variations at the high and low end of the frequencies. With regard to the present computations in comparison to the other available results, the matching can be stated to be very good. It is interesting to
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2254
0.30 (5x5) B-spline (7x7) B-spline Beck & King PFM Haskind Relation
Nondimensional heave force
0.24
0.18
0.12
0.06
0.00 0
1
2 kR
3
4
Fig. 7. Heave diffraction force.
0.25
Nondimensional sway force
0.20
0.15
0.10 (5x5) B-spline (7x7) B-spline Beck & King PFM Haskind Relation Cohen
0.05
0.00 0
1
2 kR
Fig. 8. Sway diffraction force.
3
4
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2255
Table 4 Convergent characteristics kR
Hemisphere Sway
0.5 0.75 1.0 1.6 2.0 2.5 3.0 3.5 4.0
Heave
55
77
% Er.
Abs er.
55
77
% Er.
Abs er.
0.15669 0.16359 0.14230 0.01280 0.11160 0.20957 0.24165 0.23633 0.20260
0.15672 0.16361 0.14235 0.01280 0.11160 0.20957 0.24160 0.23600 0.19500
0.019 0.012 0.035 0.000 0.000 0.000 0.020 0.139 3.751
0.00003 0.00002 0.00005 0.00000 0.00000 0.00000 0.00005 0.00033 0.00760
0.2239 0.2468 0.2670 0.2673 0.2260 0.1354 0.0600 0.0280 0.0229
0.2367 0.2538 0.2693 0.2673 0.2260 0.1354 0.0600 0.0280 0.0229
5.40 2.57 0.85 0.00 0.00 0.00 0.00 0.00 0.00
0.0128 0.0070 0.0023 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
note that the areas where some minor differences are observed, e.g. heave at kR around 1.5–2.0, present results produce highly convergent results from the two level of computations. Next, computations are carried out for a Wigley hull, due to availability of this result in King (1987), where the results from strip theory computations were also presented. The hull is defined by the formula: y=b ¼ ð1 ðz=TÞ2 Þð1 ð2x=LÞ2 Þð1 þ 0:2ð2x=LÞ2 Þ þ ðz=TÞ2 ð1 ðz=TÞ8 Þð1 ð2x=LÞ2 Þ4 .
One half of the surface is modeled as a single patch, and computations are performed using both 5 5 and 7 7 control points to define this patch. The hull surface is shown in Fig. 9 when modeled using 5 5 control points. For zero wave heading, the non-dimensional heave force and pitch moment are shown in Fig. 10 and 11 (surge result is not shown since this mode is relatively less important, and no other results could be found for comparison). Results from computation using a linear code developed earlier (Sen, 2002) is also shown for comparison. Table 5 shows the convergence characteristics from the two levels of computation. Although the percentile differences appear to be somewhat high at the high-frequency limit, the absolute differences are fairly small (less then 0.05). This is evident from the graphs where the difference between the two computations are very small, much less than their differences with other computation. As can be seen, the overall agreement between the various results is fairly good. The present results are in closest agreement with strip theory, and lies between the linear-code results of Sen (2002) and King (1987). It may be noted that Beck and King’s code is also based on a linear panel method. The linear-code computations of Sen(2002) were performed using 90 panels over one half of the surface. In Fig. 12, the time histories of heave exciting force using the linear code with 90 and 120 panels over one half, and the results from the present B-spline code, are plotted. It is seen that increase in the number of panels causes the linear-code results to converge towards the B-spline code. Also to be
ARTICLE IN PRESS 2256
R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
Fig. 9. Model of Wigley Hull.
observed is the steadiness of the peak-to-peak values. While Beck and King’s results are somewhat lower, particularly for pitch moment, considering the present level of computational maturity in 3D ship motion computations, the agreement can be termed good.
4. Conclusion An open uniform B-spline-based panel method is developed with the aim of applying the method for the 3D sea keeping problem in time-domain based on transient Green’s function. The main motivation behind this development is to reduce the computational time and memory that is associated with solution of this problem using lower order panel methods making routine application of lower order panel method for 3D ship motion computations in a PC environment difficult. In the present method, the body surface as well as the field variables (potential, normal velocity, source strength) are described by the same basis function, and then
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2257
Nondimensional heave exciting force
16 (5x5) B-Spline (7x7) B-Spline 1st order(Sen) King Strip Theory
12
8
4
0 0
4
12
8
16
kL Fig. 10. Heave exciting force for Wigley Hull.
Nondimensional pitch exciting moment
3 (5x5) B-Spline (7x7) B-Spline 1st order (Sen) King Strip Theory 2
1
0 0
2
4
6
8
10
kL Fig. 11. Pitch exciting force for Wigley Hull.
12
14
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2258
Table 5 Convergent characteristics kL
Wigley Hull Heave
2 4 6 8 10 11 12 14 16
Pitch
55
77
% Er.
Abs er.
55
77
% Er.
Abs er.
14.500 8.7940 4.7800 1.3822 0.9257 1.1595 1.1622 0.5656 0.3200
14.650 8.8013 4.7900 1.3822 0.9258 1.1597 1.1624 0.5757 0.3211
1.024 0.083 0.208 0.000 0.011 0.017 0.017 1.754 0.342
0.1500 0.0073 0.0100 0.0000 0.0001 0.0002 0.0002 0.0101 0.0011
1.660 2.360 2.210 1.400 0.551 0.272 0.100 0.030 0.080
1.640 2.346 2.120 1.363 0.550 0.270 0.089 0.018 0.075
1.204 0.593 4.072 2.642 0.181 0.735 11.00 40.00 6.250
0.020 0.014 0.090 0.037 0.001 0.002 0.011 0.012 0.005
90 panel on one half 120 panel on one half B-spline
1000
Heave exciting force
750 500 250 0 -250 -500 -750 -1000
0
10
20
30
40
50 t
60
70
80
90
100
Fig. 12. Time history of heave exciting force of a Wigley hull at L=l ¼ 1:2. The no. of panels mentioned refer to the linear-code computations.
collocation method is used for solution of the unknown control points. The method is first applied for the radiation problem in unbounded fluid, where the panel scheme is formulated using a mixed singularity (source and dipole) distribution. Results obtained, in term of added masses, show excellent agreement with analytical results, confirming the accuracy and robustness of the B-spline formulation. Next this methodology is developed for the zero-speed diffraction problem of a floating body, where the formulation is based on transient Green function. The basic solution scheme is in time domain following the approach of Lin and Yue (1990).
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2259
Computational results for a floating hemisphere and Wigley hull show very good agreement with available numerical results based on both 2D strip theory and 3D theory of other workers. Further work to extend this method for the forward speed diffraction problem and finally the forward speed ship motion problem is in progress. References Beck, F.R., King, B., 1989. Time domain analysis of wave exciting forces on floating bodies at zero forward speed. Applied Ocean Research 11, 19–25. Bingham, H., Korsmeyer, F., Newman, J., 1996, Prediction of the sea keeping characteristics of ships. In: Proceedings of 20th ONR Symposium on Naval Hydrodynamics, National Academy Press, Washington, DC, pp. 27–47. Chen, X.B, Diebeld, L., Doutreleau, Y., 2001. New green function method to predict wave induced ship motion and Loads. In: 23rd ONR Symposium on Naval Hydrodynamics, National Academy Press, Washington, DC, pp. 66–81. Danmeier, D.G., 1999. A higher order panel method for large-amplitude simulation of bodies in waves. Ph. D. Thesis, Massachusetts Institute of Technology, Massachusetts. Hess, J.L., 1979. A higher order panel method for three-dimensional potential flow, Report No. NADC77166-30. Hsin, C.Y., Kerwin, J.E., Newman, J.N., 1994, A higher order panel method based on B-spline. In: Proceedings of the 6th International Conference on Numerical Ship Hydrodynamics, National Academy Press, Washington D.C., pp. 133–151. Kim, B.K., Nam, J.H., 1999. A B-spline method using NURBS surface. In: 14th International Workshop on Water Waves and Floating Bodies, Port Huron, USA, pp. 76–79. King, B., 1987. Time domain analysis of wave exciting forces on the ships and bodies, Ph. D. Thesis, Dept. of Naval Arch. and Marine Engineering, The University of Michigan, USA. Lee, C.H., Farina, L., Newman, J.N., 1998, A geometry independent higher order panel method and its application to wave–body interactions. In: Proceedings of Engineering Mathematics and Applications Conference, Adelaide. Lin, W.M., Yue, D., 1990. Numerical Solution for large amplitude ship motions in the time domain. In: Proceedings of 18th Symposium on Naval Hydrodynamics, National Academy Press, pp. 41–66. Maniar, H., 1995. Three-dimensional higher order panel method based on B-splines, Ph. D. Thesis, Massachusetts Institute of Technology, Massachusetts. Qiu, W., Hsiung, C.C., 2002. A panel free method for time domain analysis of the radiation problem. Ocean Engineering 29, 1555–1567. Qiu, W., Chuang, J.M., Hsiung, C.C., 2004. Numerical solution of wave diffraction problem in the time domain with panel free method. Journal of Offshore Mechanics and Artic Engineering 126, 1–8. Rogers, D.F., Adams, J.A., 1990. Mathematical Elements for Computer Graphics. McGraw-Hill, New York. Sen, D., 2002. Time domain computation of large amplitude 3D ship motion with forward speed. Ocean Engineering 29, 973–1002. Telles, J.C.F., 1987. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals. International Journal for Numerical Methods in Engineering 24, 371–381. Zhao, C.B., Zou, Z., 1999. A three-dimensional potential flow computing method based on NURBS. In: 14th International Conference on Water Waves and Floating Bodies, Port Huron, USA, pp. 171–174.
Ocean Engineering 33 (2006) 2240–2259 www.elsevier.com/locate/oceaneng
A B-spline-based method for radiation and diffraction problems Ranadev Datta, Debabrata Sen Department of Ocean Engineering and Naval Architecture, IIT-Kharagpur, Kharagpur 721302, India Received 18 May 2005; accepted 4 December 2005 Available online 9 March 2006
Abstract An open uniform B-spline-based panel method is developed for solution of potential flow problems. In this method, both geometry as well as the field variables are represented by the same open uniform B-spline basis function. The method is initially applied for the radiation problem in unbounded fluid. Computed results for a spheroid of different aspect ratio are found to be in excellent agreement with analytical results. The method is then applied for diffraction problem formulated based on the transient (time-domain) Green’s function. Computed results for a hemisphere and Wigley hull are compared with published results and the comparison shows good agreement. r 2006 Elsevier Ltd. All rights reserved. Keywords: B-spline; Basis function; Panel method; Radiation problem; Diffraction problem
1. Introduction Computation for the motion of a ship when it traverses through a wave field and the corresponding force it experiences is of fundamental interest in marine hydrodynamics. Over the years, several researchers have studied this problem. Most of the earlier methods were based on two-dimensional strip theory. Over the recent past, more advanced three-dimensional (3D) methods are being developed. Corresponding author. Fax: +03222 282284.
E-mail addresses: [email protected] (R. Datta), [email protected] (D. Sen). 0029-8018/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.oceaneng.2005.12.002
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2241
A majority of these solutions are based on panel methods in some form or other. Based on surface description and unknown function representations, there can be different versions of panel methods which can be generally referred to as lower and higher order panel methods. In a lower (0th or 1st order) panel method, planar panels are used and the flow variables (potential, normal velocity and/or singularity strengths) are assumed either constant or linearly varying over individual panels. In the higher order panel methods, some higher order functions like quadratic or cubic Lagrangian interpolation are used to describe the surface geometry as well as the unknown field variables. Higher order methods are in general relatively more efficient and accurate compared to the lower order methods, but often these methods suffer from some undesirable characteristics. For example, in the higher order method of Hess (1979), use of cubic panel shape with quadratic source density seems to result in an extremely complicated algorithm. Also, the higher order method using Lagrangian interpolation function gives discontinuous normal vector and spatial derivatives of the potential across the panel boundaries. Most recently, B-spline based panel methods are being developed for wave– structure interaction and free-surface flow problems. Some of the relevant works in this context are as follows. Hsin et al. (1994) developed a 2D higher order panel method based on B-splines in which the geometry and singularity distribution were defined by B-spline and Green’s theorem was then applied at selected collocation points to solve for the unknown flow parameters. Maniar (1995) used B-splines to develop a higher order 3D panel method, where the potential and the geometry of the body are allowed any degrees of continuity. The developed method is then successfully applied for the zero-speed frequency-domain radiation–diffraction problem. Kim and Nam (1999) proposed a B-spline panel method where the surface is described by NURBS which provides for more accurate geometry modeling, specially for the conic surface. Zhao and Zou (1999) described the surface as well as unknown potential by NURBS to study the potential flow around the body with/ without free surface. Chen et al. (2001) proposed bi-quadratic patch method to predict the wave-induced ship motions and load in frequency domain. Lee et al. (1998) and Danmeier (1999) have presented a geometry independent higher order method in which the geometrical and hydrodynamic representations are decoupled. Here the velocity potential is represented by a set of B-spline basis functions, but the geometry can be modeled by any regular parameterization. In this way the hydrodynamic problem is completely separated from the geometric model. Qiu and Hsiung (2002) introduced a panel free method and applied it to solve the radiation problem in time domain, and subsequently extended the method for the diffraction problem as well (Qiu et al., 2004). In this method, the body boundary was represented by NURBS surface and the integral equations based on source-strength distribution were globally discretized over the body surface by Gaussian quadrature. Due to a de-singularization procedure adopted, there was no need to express the source strength and velocity potential by spline-representation. As can be seen from above, there have been quite some interest in developing Bspline-based solution schemes for a variety of wave–structure and free-surface flow problems. The motivations of these developments are that (i) B-splines provide an
ARTICLE IN PRESS 2242
R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
accurate description of the geometry. Such B-spline-based geometry modeling are presently used extensively in CAD application for ship hull design and thus the geometry description from such hull-development modules can be directly used for the hydrodynamic computations. (ii) A substantial reduction in computational efforts (time as well as required memory space) compared to lower order methods for the same level of accuracy can be achieved. One of the important problems in ship hydrodynamic in which (ii) above becomes very significant is the solution of the 3D forward speed ship motion problem in timedomain using the transient Green function. Lin and Yue (1990) have used first-order panel method to solve this problem, but the method is extremely time consuming due to the convolution nature of the integral and the large number of panels needed for accurate description of the hull. For example, Bingham et al. (1996) studied the motions of a ship in waves using a version of the 3D transient Green function method and applying a low-order panel method for the discretization of the integral relations. For a series 60, C b ¼ 0:7 hull, computation carried out with 256 unknowns and upto 300 time-steps required 22 CPU hours on the entry-level SGI Indigo (MIPS 3000 processor). The memory requirement was very high. In fact, the memory space required in these methods are of the order OðN 2p N t Þ and the computation effort is OðN 2p N 2t Þ, where Np is the number of panel and Nt is the number of time steps. For practical computations, it may be required to run a simulation for substantially higher number of time steps, in particular if behaviors in irregular waves are to be studied. Clearly, too much of computation time and memory is involved in the process: a typical problem of this type for a realistic ship geometry may require a duration of the O(15) hours in todays PC for the solution by lower order panel methods. Motivation for developing a B-spline-based method for this application is therefore substantial. In the present work, an open uniform B-spline-based method is developed with the final goal of applying this method for the 3D forward speed ship motion problem based on transient Green function. A lower order panel method for this problem was developed earlier (Sen, 2002). Despite several developments on B-spline-based methods for the general wave–body interactions, there appear to be a lack of application of B-spline-based methods for the 3D forward speed ship motion problem. In order to check the performance of the basic spline-based integral-equation solver, initially the radiation problem for a fully submerged body is considered. Here the problem is formulated based on a mixed source-sink distribution, and both geometry and field variables (velocity potential and normal velocity) are represented by the same open-uniform B-spline basis functions. Numerical integration is then performed at selected collocation points. Results in terms of added mass for some regular geometric shapes are found to be in excellent agreement with analytical solutions. Next the solution of the zero-speed diffraction problem of a surfacepiercing floating body is considered, where the problem is formulated based on transient Green function, and the integral equation is based on source distribution (as in Lin and Yue, 1990). For evaluation of the Rankine part (the singular part) of the integration, the adaptive subdivision scheme proposed by Telles (1987) is applied
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2243
for near field, while Gaussian integration is performed for far field. The regular (memory) part of the Green’s function is integrated using a 4-point quadrature rule. Trapezoidal rule is applied for the convolution (time) integration. The computed results for diffraction forces for Wigley hull are compared with published results, and are found to be in good agreement. In the following, we first discuss the radiation problem in unbounded fluid, followed by the zero-speed diffraction problem for a ship hull. (1) Radiation problem: For the radiation problem of a body in unbounded fluid, the potential function Fk ðp; tÞ satisfies the Laplace equation and the boundary conditions, as follows: r2 Fk ðp; tÞ ¼ 0 in O; ðFn¯ Þk ¼ n¯ k for k ¼ 1; 2; 3 on Sb , ðFn¯ Þk ¼ ð¯r n¯ Þk3 for k ¼ 4; 5; 6 on
(1.1)
Sb .
ð1:2Þ
In the above, k ¼ 1; 2; 3 refer to linear modes of motions (surge, sway, heave), k ¼ 4; 5; 6 are the rotational modes of motions (roll, pitch, yaw), n¯ ¼ ðnx ; ny ; nz Þ, r¯ ¼ ðrx ; ry ; rz Þ and Sb is the body surface. It is well known that from application of Green’s theorem, the potential function can be represented by a distribution of sources and strengths over the boundary surface as follows: ZZ qGðp; qÞ qFðqÞ aðpÞFðpÞ ¼ FðqÞ Gðp; qÞ dS. (1.3) qnq qnq S The above integral relation represents potential at a field point p by distributing sources G of local source strength qFðqÞ=qn and the dipole qGðqÞ=qn of moment FðqÞ at the point qðx; Z; zÞ on the surface S of the fluid domain. G here is the Green’s function, which for the present problem can be taken as a simple source, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G ¼ 1=R ¼ 1= ðx xÞ2 þ ðy ZÞ2 þ ðz zÞ2 . Also, aðpÞ is the solid angle on the surface at a point pðx; y; zÞ. For an unbounded domain, S reduces to Sb, the body surface, since the contribution over the infinite outer surface vanishes. q=qn represents the normal derivative and the subscript q denotes that this derivative is to be taken at the ‘source point’ q. In order to solve the problem using B-splines, the body surface Sb can be divided into number of large patches. Over each patch, the surface as well as the field variables F, qF=qn are represented by open uniform B-spline basis functions as follows. Let X¯ ¼ ðx; y; zÞ are the Cartesian coordinates. These can be written parametrically as x ¼ X ðu; wÞ, y ¼ Y ðu; wÞ, z ¼ Zðu; wÞ where u,w are the parameters. Any point X¯ ðpÞ on the surface can then be described as X¯ ðpÞ ¼
nþ1 m þ1 X X i¼1 j¼1
Ax¯ ði; jÞN ki ðuÞM lj ðwÞ.
(1.4)
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2244
Here N ki ðuÞ and M lj ðwÞ are the open uniform B-spline basis functions of order k and l, respectively. The same representation is applied over all the patches describing the full surface. In (1.4), Ax¯ ði; jÞ ¼ Ax ði; jÞ; Ay ði; jÞ; Az ði; jÞ are called the vertices of the defining polygon net, also known as control points. ðn þ 1Þ and ðm þ 1Þ are the total number of control points along u and w parametric directions, respectively. The surface is thus described by ðn þ 1Þ ðm þ 1Þ number of control points. Once the basis functions are specified, depending on the order of the spline used, the problem of representing the surface reduces to finding Ax ði; jÞ; Ay ði; jÞ; Az ði; jÞ. ðn þ 1Þ ðm þ 1Þ number of equations are needed to solve for these unknowns, These can be obtained by simply making the surface go through that many collocation points, i.e by fitting a B-spline surface passing through a set of known points. The collocation points in the parametric space have been determined based on chord line approximation (see Rogers and Adams, 1990). Specifically, for r ¼ n þ 1 physical data points in the u parametric direction, the parameter value at the lth data point is u1 ¼ 0;
ul umax
Pl g¼2 jDg;s Dg1;s j ; l ¼ 2; n. ¼ Pr g¼2 jDg;s Dg1;s j
(1.5)
Similarly, the lth data point in the w parametric direction is Pl wl g¼2 jDr;g Dr;g1 j w1 ¼ 0; ; ¼ Ps wmax g¼2 jDr;g Dr;g1 j
l ¼ 2; m
(1.6)
In the above, umax and wmax are the maximum value of the appropriate knot vectors. D is defined by Dp;q ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxp xq Þ2 þ ðyp yq Þ2 þ ðzp zq Þ2 .
(1.7)
The unknown potential F and normal velocity Fn are expressed by the same basis functions: F¼
nþ1 X mþ1 X
Af ði; jÞN ki ðuÞM lj ðwÞ,
(1.8)
i¼1 j¼1
Fn ¼
nþ1 X mþ1 X
Afn ði; jÞN ki ðuÞM lj ðwÞ.
(1.9)
i¼1 j¼1
If the surface is described by Pn number of patches, then the discretized form of Eq. (1.3) for the sth collocation point (i.e, in parametric domain, for
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2245
u ¼ us1 ; w ¼ ws2 ) is given by " ( Pn nþ1 X mþ1 X X aðps Þ Af ðr; i; jÞN ki ðus1 ÞM lj ðws2 Þ r¼1
i¼1 j¼1 nþ1 X mþ1 X
Z
f
u¼0
Pn Z umax X r¼1
Z
wmax
A ðr; i; jÞ
i¼1 j¼1
¼
umax
Z
u¼0
wmax
w¼0
w¼0
nþ1 m þ1 X X
fn
A
)#
N ki ðuÞM lj ðwÞG n ðus1 ; ws2 ; u; wÞjJj du dw !
ðr; i; jÞN ki ðuÞM lj ðwÞ
Gðus1 ; ws2 ; u; wÞjJj du dw.
i¼1 j¼1
ð1:10Þ The normal n¯ can be found from the following relations: X¯ u ¼
nþ1 X mþ1 X
k
Ax¯ ði; jÞN 0 i ðuÞM lj ðwÞ,
i¼1 j¼1
X¯ w ¼
nþ1 m þ1 X X
l
Ax¯ ði; jÞN ki ðuÞM 0 j ðwÞ,
(1.11)
i¼1 j¼1
n¯ ¼
X¯ u X¯ w . X¯ u X¯ w
The Jacobian is given by jJj ¼ X¯ u X¯ w .
(1.12)
(1.13)
fn
In (1.10), A in the right-hand side is determined directly from application of the body boundary condition (1.2), ðFn¯ Þk ¼ n¯ k , k ¼ 1; 2; 3; ðFn¯ Þk ¼ ð¯r n¯ Þk3 , k ¼ 4; 5; 6. To obtain a solution for the unknowns Af , u and w are divided into n m points resulting into ðn þ 1Þ ðm þ 1Þ rectangular segments in the parametric domain. Each segment corresponds to a surface panel in the physical domain. Collocation method is used to get the linear system of equations with ðn þ 1Þ ðm þ 1Þ unknowns. This results into a system of nn ¼ ðn þ 1Þ ðm þ 1Þ linear equation: ½BfAf g ¼ fCg,
(1.14)
where [B] is a nn nn matrix of influence coefficients containing the integrals over individual panels in the parametric domain as given in the left-hand side of (1.10). All integrations are carried out using a 3 3 Gaussian quadrature scheme. However, when the integration is performed over the so-called self-influencing panels, the integrand contains a singularity. There are various schemes available to handle such a singularity. In the present work, we use the adaptive subdivision scheme proposed by Telles (1987). In this, the self-influence panel is subdivided successively to isolate the singular point, and a 3 3 Gaussian quadrature rule is applied to the subdivided parts of the panel. For the small subdivision part of the surface over which the singular point lies, analytical integration is performed, assuming that the surface is
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2246
Table 1 Absolute relative error No of subdivision
Relative error
5 10 12 15 20 30
0.085400 0.002776 0.000692 0.000085 0.000002 0.000000
(0,1)
(1,1)
w
(0,0)
u
(1,0)
Fig. 1. Subdivision of panel.
planar. The subdivision is continued till the integrand converges within a relative absolute accuracy of 0:00005 for a unit square patch in the parametric domain. Relative absolute accuracy is defined as jðas ac Þ=ac j where as is the computed value at the sth level of subdivision and ac is the value at a high level of subdivision where the result is assumed to have converged (e.g. 30 subdivisions as per Table 1 above). It is found that usually about 15–20 subdivisions are sufficient to achieve this accuracy level. The contribution of the small isolated panel to the integral is extremely small, and can even be ignored without any difference for this level of accuracy. This subdivision scheme is shown in Fig. 1 and illustrative results are presented in Table 1. Variety of methods exists to solve the system of Eqs. (1.14). Once the unknowns Af’s are found, F at any desired location can be determined directly from (1.8). The results obtained for this radiation problem are presented in terms of the familiar added masses (mi;j ) defined by ZZ mi;j ¼ r fi nj ds; ðm1;1 for surge; m3;3 for heaveÞ. (1.15) dO
In discretized form, this can be written as " # ! pn Z u¼umax Z w¼wmax X nþ1 m þ1 X X f k l mi;j ¼ r A ði; jÞN i ðuÞM j ðwÞ n¯ :jJj du dw . r¼1
u¼0
w¼0
i¼1 j¼1
(1.16) where the normal n¯ and Jacobian jJj can be determined from (1.12) and (1.13).
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2247
2. Results Numerical computations are performed for a spheroid of different aspect ratio, since for this geometry the analytical results are available. Due to symmetry, only one octant of the spheroid is needed to be modeled as one patch. We use 5 5 ¼ 25 control points to represent this patch. Fig. 2 shows the geometry of one octant of the spheroid and the full body. Tables 2 and 3 present results of non-dimensional heave and surge added masses (m ¯ i;j ¼ mi;j =M, M ¼ mass of the body) for the range of aspect ratios a=b ¼ 0:221. It
Fig. 2. Geometry of spheroid.
Table 2 Added mass for heave b=a
Numerical
Analytical
% of error
Abs error
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.91070 0.82350 0.76080 0.70390 0.65320 0.60750 0.56730 0.53130 0.49990
0.89426 0.82587 0.76189 0.70421 0.65294 0.60757 0.56741 0.53176 0.50000
1.838 0.272 0.143 0.044 0.039 0.011 0.019 0.086 0.020
0.0166 0.0023 0.0011 0.0003 0.0003 0.0001 0.0001 0.0005 0.0001
Table 3 Added mass for surge b=a
Numerical
Analytical
% of error
Abs error
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.06000 0.10520 0.15630 0.21080 0.26570 0.32280 0.38120 0.44040 0.49950
0.05912 0.10542 0.15626 0.21002 0.26576 0.32294 0.38120 0.44028 0.50000
1.488 0.209 0.026 0.371 0.022 0.043 0.000 0.027 0.100
0.0009 0.0002 0.00004 0.0007 0.00006 0.0001 0.0000 0.0001 0.0005
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2248
0.9
3D B-spline Analytical results
m33/ (4/3)πρab2
0.8
0.7
0.6
0.5 0.2
0.4
0.6 b/a
0.8
1.0
Fig. 3. Heave added mass.
may be noted that heave added mass for an aspect ratio a=b is equivalent to surge added mass for the aspect ratio b=a. This means, computational results presented cover linear added masses for a range of aspect ratio a=b ¼ 0:225, i.e. the geometry is varied from a shallow oblate spheroid to a fairly long prolate spheroid. The analytical values as well as the percentage and absolute errors are also shown. The results are plotted in Figs. 3 and 4. As can be seen, the results are in excellent agreement and are graphically indistinguishable, except marginally at the ends. The agreement somewhat deteriorates as the body gets elongated. However, the differences are still less than 0.4% up to 0:3pa=bp3:33. Only for the last case of a=b ¼ 0:2 (or equivalently a=b ¼ 5), the error is somewhat large at around 2%. This may be due to the fact that for elongated shape, a single patch may not be sufficient to represent the geometry with the required level of precision. For such cases, a larger number of patches may be used to improve the results. (2) Diffraction problem: We next consider the diffraction problem of a floating free surface piercing body, formulated based on transient Green function in the time domain. Consider an arbitrary 3D body floating on the free surface of an incompressible ideal fluid. Let Oxyz be the right-handed co-ordinate system with Oxy plane coincident with the mean water level z ¼ 0. The z-axis is positive upwards, and sinusoidal (linear) waves traveling in a direction making an angle b to x-axis is incident on the body. The initial boundary value problem for the diffraction potential Fd ðp; tÞ is specified by the following governing equation, boundary
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2249
0.5 3D B-spline Analytical results
m11/ (4/3)πρab2
0.4
0.3
0.2
0.1
0.2
0.4
0.6 b/a
0.8
1.0
Fig. 4. Surge added mass.
condition, initial conditions and far field conditions: r2 Fðp; tÞ ¼ 0 in O;
(2.1)
q2 F qF ¼ 0 on z ¼ 0, þg qt2 qz
(2.3)
qFd qFi ¼ qn qn
(2.4)
on
S0
Fd ! 0 at t ! 0,
(2.5)
qFd !0 qt
at t ! 0,
(2.6)
rFd ! 0
as r ! 1.
(2.7)
In the above, Fi is the incident wave potential, S0 is the mean wetted body surface, n¯ is the normal on S0 directed outwards to the fluid domain (inward to body), and g is the acceleration due to the gravity. The solution of the above initial boundary-value problem is based on the transient Green function approach of Lin and Yue (1990). The diffracted potential can be expressed by the following convolution integral by distributing transient Green’s
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2250
function of unknown strength over the body surface, 8 2 39 > > ZZ Z t ZZ = < 1 6 7 Fd ðp; tÞ ¼ sG0 dS þ dt4 sG ft dS5 , > 4p > 0 ; : S 0 ðtÞ
(2.8)
S 0 ðtÞ
where q ¼ qðx; Z; zÞ represents a ‘source point’ on the body surface and sðq; tÞ is the associated source strength at time step t. G is the transient Green’s function given by Gðp; t; q; tÞ ¼ G 0 þ Gf
for paq; t4t
(2.9)
where 1 1 ðRankine partÞ, r r0 Z 1 pffiffiffiffiffiffi f ½1 cos gkðt tÞ ekðzþxÞ J 0 ðkRÞ dk ðfree surface memory partÞ, G ¼2 G0 ¼
0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ðx xÞ2 þ ðy ZÞ2 þ ðz zÞ2 , r0 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xÞ2 þ ðy ZÞ2 þ ðz þ zÞ2 ,
R¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xÞ2 þ ðy ZÞ2 ,
k¼
o2 g
is the wave number,
J 0 ¼ Bessel function of order zero.
(2.10)
Eq. (2.9) satisfies all conditions except the body boundary condition. Application of the body boundary condition to (2.8) yields the following integral relation: 9 8 > > ZZ Z t ZZ = < qFd ðp; tÞ 1 f 0 (2.11) ¼ sGnp dS dt sG np t dS . > qnp 4p > 0 ; : S 0 ðtÞ
S 0 ðtÞ
Detail derivation of the above is available in several sources, e.g. Lin and Yue (1990). Once the diffracted potential is obtained, pressure on the body can be determined from Bernoulli equation, and the linear wave forces and moments can be easily obtained by integrating pressure over the body surface: 2 3 ZZ ZZ qFi qFd 6 7 F i ðtÞ ¼ r4 dS þ dS 5, (2.12) ni ni qt qt S 0 ðtÞ
S 0 ðtÞ
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2251
where ni is the generalized normal and ¼ ni for i ¼ 1; 2; 3 and ¼ ð~ r~ nÞi3 for i ¼ 1; 2; 3. For the numerical solution, as in the radiation problem earlier, an open uniform B-spline basis function is used to describe the geometry and field variables. The surface is represented by Eq. (1.4) while, for an instant of time t, the unknown potential Fd and source strength s are defined as Fd ¼
nþ1 m þ1 X X
Af ðt; i; jÞN ki ðuÞM lj ðwÞ,
(2.13)
i¼1 j¼1
s¼
nþ1 m þ1 X X
As ðt; i; jÞN ki ðuÞM lj ðwÞ,
(2.14)
i¼1 j¼1
where NðuÞ and MðwÞ are the open uniform B-spline basis functions. It may be noted that this formulation is based on source distribution unlike the radiation formulation based on a mixed source and dipole distribution, and thus here the unknown source strengths are expressed by B-spline basis functions. Due to the convolution nature of the integrals, the unknown control nets Af and As are also functions of time. The convolution integral in (2.11) is treated by discretizing time into equal time-steps. Thus, (2.11) in discretized form for the Mtth time step and for sth collocation point can be written as " Pn X nþ1 Mþ1 X X As ðr; M t ; i; jÞ r¼1
Z
i¼1 J¼1 umax u¼0
Z
wmax
w¼0
# N ki ðuÞM lj ðwÞG 0np ðus1 ; ws2 ; u; wÞjJj du dw
qF ¼ 4p qnp " n (Z t M 1 P X X er 1 Dt r1 ¼1
r¼1
umax
u¼0
)#
G fnp ðu; wÞjJ j du dw
,
Z
wmax
w¼0
nþ1 m þ1 X X
! A
s
ðr; r1 ; i; jÞN ki ðuÞM lj ðwÞ
i¼1 j¼1
ð2:15Þ
where u ¼ us1 ; w ¼ ws2 are the parameters representing the sth collocation point, and t ¼ M t Dt, Dt is the time step size. Applying Eq. (2.15) to selected collocation points as in previous case, a linear system of equations for the unknowns As for time step M t can be formed: ½BfAs ðM t Þg ¼ fDg.
(2.16)
Matrix [B] consists of the integration of Rankine part G 0 of the Green’s function, which has a singularity when integrated over self-influence panels, as in the previous case. This is treated in the similar manner as in the radiation problem, by applying
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2252
a 3 3 quadrature rule for the non-singular part and using the adaptive subdivision for the integral containing singularity. For the integrals containing the regular memory part G f of the Green function, we use a 2 2 quadrature rule since this is found adequate. The convolution time integration, as shown in (2.15), is based on a simple trapezoidal rule. Efficient and accurate computation of G f and its derivatives follow the method described in Sen (2002). Once the unknown As for a time step is determined from a solution of (2.16), the unknown Af for that time step can be found from (2.8). The discretized form of (2.8) for the Mtth time step and sth collocation point (i.e, in parametric domain at u ¼ us1 ; w ¼ ws2 ), is " # pn X nþ1 m þ1 X X f t k l A ðr; M ; i; jÞN i ðus1 ÞM j ðws2 Þ r¼1
i¼1 j¼1
" n Z ( ) # Z p nþ1 X mþ1 1 X umax wmax X As ðr; M t ; i; jÞN ki ðuÞM lj ðwÞ G 0 ðus1; ws2; ; u; wÞjJj du dw 4p r¼1 u¼0 w¼0 i¼1 j¼1 " n (Z ! )# t M 1 P nþ1 X mþ1 umax Z wmax X X X s f k l e r1 A ðr; r1 ; i; jÞN i ðuÞM j ðwÞ G np ðu; wÞ du dw . Dt ¼
r1 ¼1
r¼1
u¼0
w¼0
i¼1 j¼1
ð2:17Þ Once again, applying (2.17) successively at the collocation points, one gets the following system of equations: fAf ðM t Þg ¼ ½CfAs ðM t Þg.
(2.18)
The velocity potential F can be now determined in a straightforward from (2.13). Time derivative of this potential needed for evaluating pressure is found from numerical differentiation, and forces, moments are evaluated from (2.12). 3. Results Computations are first performed for a floating hemisphere due to availability of results for this geometry. Considering flow symmetry, only one quadrant of the surface needs to be modeled. Computations are carried out using both (5 5) and (7 7) control points to define the quadrant. Fig. 5 shows the surface as modeled by (5 5) control points. Fig. 6 shows the total surface of the hemisphere (using symmetry). Computations cover a range of frequencies and the results are presented in terms of force amplitudes in the frequency domain. In all cases, the time history is extremely regular (a typical time history is shown later for the next computation). The force amplitudes are determined by averaging the peak-to-peak values over the steady part of the time-history. Fig. 7 shows the non-dimensional heave diffraction force amplitude jX 3 j=rgpR2 against non-dimensional frequency kR, where k is the wave number and R is the radius of the hemisphere. Fig. 8 shows the nondimensional sway diffraction force jX 2 j=rgpR2 against kR (sway here is defined along the direction in which the incident wave travels). Results are compared with
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2253
Fig. 5. One octant of a hemisphere.
Fig. 6. Hemisphere.
the results of Beck and King (1989), results obtained by Qiu et al. (2004) using their panel-free-method (PFM), and results obtained using Haskind relations and by Cohen. These last two results are as reported in Qui et al. (2004). As can be seen, computations with (5 5) and (7 7) points produce almost identical results. This can also be observed in Table 4. Convergence appears to be best at the mid-range of the computed frequency-range with some minor variations at the high and low end of the frequencies. With regard to the present computations in comparison to the other available results, the matching can be stated to be very good. It is interesting to
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2254
0.30 (5x5) B-spline (7x7) B-spline Beck & King PFM Haskind Relation
Nondimensional heave force
0.24
0.18
0.12
0.06
0.00 0
1
2 kR
3
4
Fig. 7. Heave diffraction force.
0.25
Nondimensional sway force
0.20
0.15
0.10 (5x5) B-spline (7x7) B-spline Beck & King PFM Haskind Relation Cohen
0.05
0.00 0
1
2 kR
Fig. 8. Sway diffraction force.
3
4
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2255
Table 4 Convergent characteristics kR
Hemisphere Sway
0.5 0.75 1.0 1.6 2.0 2.5 3.0 3.5 4.0
Heave
55
77
% Er.
Abs er.
55
77
% Er.
Abs er.
0.15669 0.16359 0.14230 0.01280 0.11160 0.20957 0.24165 0.23633 0.20260
0.15672 0.16361 0.14235 0.01280 0.11160 0.20957 0.24160 0.23600 0.19500
0.019 0.012 0.035 0.000 0.000 0.000 0.020 0.139 3.751
0.00003 0.00002 0.00005 0.00000 0.00000 0.00000 0.00005 0.00033 0.00760
0.2239 0.2468 0.2670 0.2673 0.2260 0.1354 0.0600 0.0280 0.0229
0.2367 0.2538 0.2693 0.2673 0.2260 0.1354 0.0600 0.0280 0.0229
5.40 2.57 0.85 0.00 0.00 0.00 0.00 0.00 0.00
0.0128 0.0070 0.0023 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
note that the areas where some minor differences are observed, e.g. heave at kR around 1.5–2.0, present results produce highly convergent results from the two level of computations. Next, computations are carried out for a Wigley hull, due to availability of this result in King (1987), where the results from strip theory computations were also presented. The hull is defined by the formula: y=b ¼ ð1 ðz=TÞ2 Þð1 ð2x=LÞ2 Þð1 þ 0:2ð2x=LÞ2 Þ þ ðz=TÞ2 ð1 ðz=TÞ8 Þð1 ð2x=LÞ2 Þ4 .
One half of the surface is modeled as a single patch, and computations are performed using both 5 5 and 7 7 control points to define this patch. The hull surface is shown in Fig. 9 when modeled using 5 5 control points. For zero wave heading, the non-dimensional heave force and pitch moment are shown in Fig. 10 and 11 (surge result is not shown since this mode is relatively less important, and no other results could be found for comparison). Results from computation using a linear code developed earlier (Sen, 2002) is also shown for comparison. Table 5 shows the convergence characteristics from the two levels of computation. Although the percentile differences appear to be somewhat high at the high-frequency limit, the absolute differences are fairly small (less then 0.05). This is evident from the graphs where the difference between the two computations are very small, much less than their differences with other computation. As can be seen, the overall agreement between the various results is fairly good. The present results are in closest agreement with strip theory, and lies between the linear-code results of Sen (2002) and King (1987). It may be noted that Beck and King’s code is also based on a linear panel method. The linear-code computations of Sen(2002) were performed using 90 panels over one half of the surface. In Fig. 12, the time histories of heave exciting force using the linear code with 90 and 120 panels over one half, and the results from the present B-spline code, are plotted. It is seen that increase in the number of panels causes the linear-code results to converge towards the B-spline code. Also to be
ARTICLE IN PRESS 2256
R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
Fig. 9. Model of Wigley Hull.
observed is the steadiness of the peak-to-peak values. While Beck and King’s results are somewhat lower, particularly for pitch moment, considering the present level of computational maturity in 3D ship motion computations, the agreement can be termed good.
4. Conclusion An open uniform B-spline-based panel method is developed with the aim of applying the method for the 3D sea keeping problem in time-domain based on transient Green’s function. The main motivation behind this development is to reduce the computational time and memory that is associated with solution of this problem using lower order panel methods making routine application of lower order panel method for 3D ship motion computations in a PC environment difficult. In the present method, the body surface as well as the field variables (potential, normal velocity, source strength) are described by the same basis function, and then
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2257
Nondimensional heave exciting force
16 (5x5) B-Spline (7x7) B-Spline 1st order(Sen) King Strip Theory
12
8
4
0 0
4
12
8
16
kL Fig. 10. Heave exciting force for Wigley Hull.
Nondimensional pitch exciting moment
3 (5x5) B-Spline (7x7) B-Spline 1st order (Sen) King Strip Theory 2
1
0 0
2
4
6
8
10
kL Fig. 11. Pitch exciting force for Wigley Hull.
12
14
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2258
Table 5 Convergent characteristics kL
Wigley Hull Heave
2 4 6 8 10 11 12 14 16
Pitch
55
77
% Er.
Abs er.
55
77
% Er.
Abs er.
14.500 8.7940 4.7800 1.3822 0.9257 1.1595 1.1622 0.5656 0.3200
14.650 8.8013 4.7900 1.3822 0.9258 1.1597 1.1624 0.5757 0.3211
1.024 0.083 0.208 0.000 0.011 0.017 0.017 1.754 0.342
0.1500 0.0073 0.0100 0.0000 0.0001 0.0002 0.0002 0.0101 0.0011
1.660 2.360 2.210 1.400 0.551 0.272 0.100 0.030 0.080
1.640 2.346 2.120 1.363 0.550 0.270 0.089 0.018 0.075
1.204 0.593 4.072 2.642 0.181 0.735 11.00 40.00 6.250
0.020 0.014 0.090 0.037 0.001 0.002 0.011 0.012 0.005
90 panel on one half 120 panel on one half B-spline
1000
Heave exciting force
750 500 250 0 -250 -500 -750 -1000
0
10
20
30
40
50 t
60
70
80
90
100
Fig. 12. Time history of heave exciting force of a Wigley hull at L=l ¼ 1:2. The no. of panels mentioned refer to the linear-code computations.
collocation method is used for solution of the unknown control points. The method is first applied for the radiation problem in unbounded fluid, where the panel scheme is formulated using a mixed singularity (source and dipole) distribution. Results obtained, in term of added masses, show excellent agreement with analytical results, confirming the accuracy and robustness of the B-spline formulation. Next this methodology is developed for the zero-speed diffraction problem of a floating body, where the formulation is based on transient Green function. The basic solution scheme is in time domain following the approach of Lin and Yue (1990).
ARTICLE IN PRESS R. Datta, D. Sen / Ocean Engineering 33 (2006) 2240–2259
2259
Computational results for a floating hemisphere and Wigley hull show very good agreement with available numerical results based on both 2D strip theory and 3D theory of other workers. Further work to extend this method for the forward speed diffraction problem and finally the forward speed ship motion problem is in progress. References Beck, F.R., King, B., 1989. Time domain analysis of wave exciting forces on floating bodies at zero forward speed. Applied Ocean Research 11, 19–25. Bingham, H., Korsmeyer, F., Newman, J., 1996, Prediction of the sea keeping characteristics of ships. In: Proceedings of 20th ONR Symposium on Naval Hydrodynamics, National Academy Press, Washington, DC, pp. 27–47. Chen, X.B, Diebeld, L., Doutreleau, Y., 2001. New green function method to predict wave induced ship motion and Loads. In: 23rd ONR Symposium on Naval Hydrodynamics, National Academy Press, Washington, DC, pp. 66–81. Danmeier, D.G., 1999. A higher order panel method for large-amplitude simulation of bodies in waves. Ph. D. Thesis, Massachusetts Institute of Technology, Massachusetts. Hess, J.L., 1979. A higher order panel method for three-dimensional potential flow, Report No. NADC77166-30. Hsin, C.Y., Kerwin, J.E., Newman, J.N., 1994, A higher order panel method based on B-spline. In: Proceedings of the 6th International Conference on Numerical Ship Hydrodynamics, National Academy Press, Washington D.C., pp. 133–151. Kim, B.K., Nam, J.H., 1999. A B-spline method using NURBS surface. In: 14th International Workshop on Water Waves and Floating Bodies, Port Huron, USA, pp. 76–79. King, B., 1987. Time domain analysis of wave exciting forces on the ships and bodies, Ph. D. Thesis, Dept. of Naval Arch. and Marine Engineering, The University of Michigan, USA. Lee, C.H., Farina, L., Newman, J.N., 1998, A geometry independent higher order panel method and its application to wave–body interactions. In: Proceedings of Engineering Mathematics and Applications Conference, Adelaide. Lin, W.M., Yue, D., 1990. Numerical Solution for large amplitude ship motions in the time domain. In: Proceedings of 18th Symposium on Naval Hydrodynamics, National Academy Press, pp. 41–66. Maniar, H., 1995. Three-dimensional higher order panel method based on B-splines, Ph. D. Thesis, Massachusetts Institute of Technology, Massachusetts. Qiu, W., Hsiung, C.C., 2002. A panel free method for time domain analysis of the radiation problem. Ocean Engineering 29, 1555–1567. Qiu, W., Chuang, J.M., Hsiung, C.C., 2004. Numerical solution of wave diffraction problem in the time domain with panel free method. Journal of Offshore Mechanics and Artic Engineering 126, 1–8. Rogers, D.F., Adams, J.A., 1990. Mathematical Elements for Computer Graphics. McGraw-Hill, New York. Sen, D., 2002. Time domain computation of large amplitude 3D ship motion with forward speed. Ocean Engineering 29, 973–1002. Telles, J.C.F., 1987. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals. International Journal for Numerical Methods in Engineering 24, 371–381. Zhao, C.B., Zou, Z., 1999. A three-dimensional potential flow computing method based on NURBS. In: 14th International Conference on Water Waves and Floating Bodies, Port Huron, USA, pp. 171–174.