Computations of Wave Loads Using a BSpline Panel Method
C.H.Lee, H.Maniar, J.Newman, X.Zhu (Massachusetts Institute of Technology, USA)
Abstract
This paper describes the development of a higherorder panel method based on Bsplines, and its application to problems involving linear and nonlinear wave interactions with floating or submerged bodies. The use of Bsplines permits a very accurate representation of the body geometry, and a continuous representation of the velocity potential on the body surface. These features result in more efficient and robust solutions relative to the conventional loworder panel method.
Several applications are described to illustrate the computational advantages of this method. These include firstorder solutions for a floating hemisphere, a tensionleg platform, and a long periodic array which is representative of very large platforms proposed for floating airports and cities. Also considered are second and thirdorder wave loads, where fundamental difficulties exist in evaluating the freesurface forcing functions using the loworder method. This is particularly important at third order, where we present preliminary results from an investigation aimed at predicting the occurrence of ringing.
1.
Introduction
Computations based on panel methods have achieved widespread use among practicing engineers and hydrodynamicists for externalflow potential problems. In applications concerning the wave loads and motions of offshore platforms, and other stationary vessels, the firstorder analysis in the frequency domain is now routine (Korsmeyer et al. [10], Herfjord & Nielsen [6]). Panel methods can be used not only for single bodies but also for multiple interacting vessels, with important applications to marine operations. Success also has been achieved in extending this method to the analysis of secondorder wave effects (Molin & Chen [18]; Eatock Taylor & Chau [3]; Lee & Newman [12]). However, notwithstanding the popularity and success of this technique, there are various special problems where the computational burden is severe, or where quantities are required such as gradients of the velocity field, which cannot be computed robustly from the conventional ‘loworder' panel method.
In loworder methods the submerged surface of the body is represented by a large number N of flat quadrilateral panels, and the velocity potential (or the source strength) is assumed piecewise constant on each panel. Only sources were used in the pioneering work of Hess & Smith [7], with the unknown source strength determined from an integral equation derived from the boundary condition on the body. In the complementary procedure where sources and normal dipoles are used, the potential is determined more directly from an integral equation based on Green 's theorem. In either of these approaches, the integral equation is replaced by a system of linear equations, with the number of unknowns equal to N.
The use of flat panels simplifies the evaluation of the integrated singularity distribution on each panel, but a relatively large number of panels are required to achieve accurate representations of the geometry and the potential. A more fundamental restriction is that the potential derived from Green's theorem cannot be differentiated analytically to derive the fluid velocity on or near the body surface; this difficulty can be circumvented in a limited manner by using a distribution of sources only, but that method also fails if gradients of the velocity field are required.
Higherorder panel methods have been devel
oped to overcome these difficulties. Most of the higherorder methods are based on piecewise polynomial approximations of the geometry and potential on each panel, usually restricted to linear or quadratic representations using local polynomials of first or seconddegree, respectively. In this paper we describe a different higherorder approach, where Bspline basis functions are used to represent the geometry and potential. This offers the possibility of a more continuous representation, with greater geometrical flexibility and numerical efficiency. Following a preliminary investigation of this technique in two dimensions (Hsin et al., [8]), we have developed a threedimensional panel program which will be referred to as ‘HIPAN'. A more detailed description of the method is contained in the thesis of Maniar [16].
HIPAN was developed initially to solve linearized radiationdiffraction problems in the frequency domain. Thus the body is fixed, or performing small oscillatory motions about a fixed mean position, and plane progressive waves are incident from a prescribed direction. The fluid is either of constant finite depth, or infinitely deep. The body may be floating on the surface or submerged. For this case the freesurface Green function can be used effectively. The analogous loworder panel method WAMIT, which is described by Lee [11], can be used to test the accuracy and relative computational efficiency of HIPAN. The efficiency of HIPAN is most apparent for relatively complicated body shapes. This will be illustrated by considering a large array of floating cylinders similar to a floating bridge. In that particular problem we find not only that substantially larger arrays can be analyzed, but also that the results display very interesting features closely related to the occurrence of trapped waves in a channel.
In many applications important nonlinear effects must be analyzed, particularly the secondorder sum and differencefrequency loads which occur at relatively high and low frequencies compared to the firstorder wave spectrum. Highfrequency loads are important in cases where resonant structural response is encountered. Examples include hull deflections of long slender ships, bending of vertical monotowers, and vertical motions of tensionleg platforms. Conversely, lowfrequency loads are important for rigidbody motions where the restoring forces are relatively weak. Examples include the horizontal oscillations of moored and towed vessels, and the vertical response of vessels with small waterplane areas. The analysis of these secondorder loads can be performed using a loworder panel method, but the computational burden is quite large and much care is necessary to ensure robust results. The extension of HIPAN to include secondorder loads is currently in progress, and we are not able to show complete computations, but preliminary results which have been obtained demonstrate the improvements that may be achieved using the Bspline methodology in the secondorder analysis.
It should be noted that we refer to ‘order' in two completely different contexts here, one specifying the numerical approximation of the geometry and velocity potential on the body surface, and the other referring to the order of the perturbation expansion in terms of the amplitude of the waves and body motions. Loworder and higherorder panel methods are distinguished by the type of numerical approximations used, whereas firstand higherorder wave loads are defined with respect to the corresponding powers of the wave amplitude.
Recent attention in the offshore community has been directed toward higherorder nonlinear loads which are thought to cause ‘ringing', a hydrodynamic/structural resonance which has been observed for large platforms in extreme wave conditions. The relevant frequencies suggest that ringing may be caused by thirdharmonic wave loads, and this has motivated two recent studies which are restricted to the simplest possible geometrical configuration, a vertical circular cylinder. In the work of Faltinsen et al. [4] a perturbation scheme is employed appropriate to the regime where the cylinder radius and wave amplitude are of comparable magnitude, and both are small compared to the wavelength. In the complementary work of Malenica & Molin [15], a conventional Stokes expansion is used with the wave amplitude assumed small compared to the cylinder radius and wavelength and without restricting the radius/wavelength ratio. Partly to clarify the differences between these two works, and also to develop a more general computational tool for analyzing practical bodies, we have extended HIPAN to include thirdorder loads. Here too our present results are incomplete, but they do indicate that the Stokes expansion may be more appropriate for contemporary platforms.
In Sections 2–3 we review the boundaryvalue problems to be addressed, and the formulation of the panel method for solving these problems based on Bsplines. Solutions of several illus
trative problems are described in Section 4. In the concluding Section 5 other possible applications and extensions of HIPAN are discussed. For simplicity we restrict our attention to diffraction problems where the body is fixed. Except where noted the fluid depth is assumed infinite.
2.
The BoundaryValue Problem
Cartesian coordinates x=(x, y, z) are defined with z=0 the plane of the undisturbed free surface and z<0 the fluid domain. In accordance with the assumptions of potential theory, and a perturbation solution in terms of a small parameter (e.g. the incidentwave slope), the velocity potential is expanded in the form
(1)
where t denotes time.
Assuming a discrete spectrum with frequency components ω_{n}, the firstorder potential can be expressed as
(2)
where ω_{n}>0. The second and thirdorder potentials, which result from quadratic and cubic interactions among firstorder solutions, can be expressed as
(3)
(4)
where Re denotes the real part of the complex quantity.
The complete higherorder solutions contain a variety of frequency components. For example, the secondorder solution contains the secondharmonic, mean and sum and differencefrequency components. The thirdorder solution contains more diverse combinations of frequency components. To simplify the discussion here we consider the case of monochromatic firstorder incident waves. Specifically we consider the secondharmonic and mean components of the secondorder solution and the thirdharmonic component of the thirdorder solution.
The diffraction potentials can be decomposed into incident and scattering components. In general, it is convenient to solve for these potentials separately, due to their different asymptotic behavior at the farfield. However in the case of monochromatic incident waves in infinite water depth, the second and thirdorder incident wave potentials vanish and it is not necessary to distinguish between the diffraction and scattering components at these orders.
The diffraction potentials are subject to the freesurface condition
(5)
where j=1, 2, 3. The forcing functions at each order j are evaluated from the following equations, with z=0:
(6)
Each of the diffraction potentials also satisfies the condition of zero normal velocity on the body surface, and vanishes at large depths. In addition, they are subject to appropriate radiation conditions in the farfield. The scattering part of the firstorder potential satisfies a Sommerfeld radiation condition. The farfield behavior of the second and thirdorder potentials include outgoing ‘free' waves, which satisfy a Sommerfeld radiation condition, and ‘locked' waves which are in phase with the [15].
Integral equations suitable for evaluating the diffraction potentials are derived by applying Green's theorem to the volume of fluid in the domain which is bounded by S_{b}, S_{f}, and a closure surface which is in the far field and at large depths below the free surface. For this purpose we use
the wave source potential as the Green function (cf. [19], equation 7.9). The resulting equations are Fredholm integral equations for the unknown ^{(j)} over the domain S_{b}. The firstorder diffraction potential is obtained from
(7)
and the higherorder potentials from
(8)
In the forcing function on the righthandside of the latter equation the integrals over S_{b} and the farfield closure vanish due to the boundary conditions of the diffraction potentials as stated above.
The principal difficulty in the computation of the higherorder potentials lies on the evaluation of the integral over the free surface. The integrands are oscillatory, and slowly decaying. Truncation at a finite distance from the body is generally not adequate. Instead, following Lee & Newman [12], S_{f} is divided into two regions separated by a circle which encloses the body, and is of sufficiently large radius so that the effect of the evanescent modes of the scattering waves is negligible outside the circle, compared to the propagating modes. In this infinite outer region, the Green function and the velocity potentials can be expanded in FourierBessel series and the integrals can be reduced to the sums of onedimensional integrals of the products of the corresponding Bessel functions. The latter integrals are nontrivial, but they can be evaluated numerically with appropriate algorithms. Direct numerical quadrature is used for the evaluation of the surface integral over the finite inner region. The details of this procedure, and the utilization of Bspline basis functions are discussed in Section 3.3.
The integrated wave loads or exciting forces are evaluated by integrating the pressure due to the diffraction potentials over the wetted surface of the body. The final expressions for the first, second and thirdorder forces are
(9)
(10)
and
(11)
Here ρ is the water density, C_{w} denotes the intersection of the body with the mean free surface, and n is the normal vector to the body surface. For the corresponding moments n is replaced by (x×n).
3.
The Bspline Methodology
We consider a threedimensional body (or ensemble of separate bodies) where the submerged body surface consists of a finite set (p=1,…, ) of continuous surfaces, referred to as patches. The Cartesian points x=(x, y, z) on each patch are mapped to a rectangular parametric domain (u, v) by a continuous mapping function, which is assumed here to be described by a Bspline tensor product expansion of the form
(12)
Here and are Bsplines of order (degree −9) in the parametric variables u and v respectively, and are prescribed vertices or coefficients. The total number of ^{p} and splines used for the approximation are _{p} and _{p} respectively.
For a detailed description of Bsplines we refer to DeBoor [2]. It suffices here to note that Bsplines are polynomial functions which are nonzero over a finite span (support), and that their shape and continuity are dependent on the order and on a monotonic increasing sequence of real numbers
(13)
which are called knots. The rectangular parametric space which is mapped by (12) onto the physical surface is
(14)
This is referred to as the usable space.
As an example, Figure 1 shows one quadrant of a tensionleg platform (TLP), represented by 35 patches. On the pontoons a total of 32 patches are used, and outlined separately in the figure. Three larger patches are used to represent the column, including one on the bottom circular disk, one on the side to represent the complete circular cylinder above the pontoons, and one on the lower outside surface of the cylinder between the pontoons. (While not relevant to the representation of the geometry, each of these larger patches is subdivided for hydrodynamic purposes into smaller panels which we define later. These panels are outlined in Figure 1.)
Utilizing the parametric space provided by each patch, we also approximate the potential by a Bspline tensor product expansion. Thus, the potential on patch p is represented in the form
(15)
Here and are Bsplines of order k in the same parametric variables u and v, dependent on the knot vectors
(16)
and are unknown complex coefficients to be solved for.
Since different Bsplines may be used in (12) and (15) the method is not restricted to be isoparametric. To ensure that the usable parametric space implicit in (15) matches that of the geometry (12), we require that Multiple knots may be used in (12–13) to represent special features of the patch geometry. Multiple knots are not allowed in the usable space for the potential splines. In the absence of multiple knots for potential splines, we define M_{p}=I_{p}–k+1 and N_{p}=J_{p}–k+1, which correspond to the total number of nonzero intervals between consecutive knots of the usable space in the u and v directions respectively.
3.1
Features of the method
By allowing piecewise continuous bodies to be described in a patchwise fashion, one can represent a wide variety of complex structures with a very accurate geometric description of each patch.
While the potential is approximated in the same parametric variables as the surface, these two representations are otherwise independent. Consequently, increasingly accurate potential approximations can be sought without reconstructing the geometry.
The representations (12) and (15) are local approximants due to the finite support of the splines. Further, implicit in these equations is the possibility of a nonuniform distribution of the splines over the parametric space. Both of these features are desirable for approximating the solution in the vicinity of flow singularities (e.g., in the vicinity of body edges and the free surface).
An approximant based on Bsplines of order k has (k–2) degrees of continuity everywhere on a patch, unlike many higherorder methods where continuity is ensured only over a panel or element. Continuous descriptions of parametric derivatives of the potential can be obtained simply by analytic differentiation. Cartesian derivatives of the potential can also be derived analytically, using the parametric derivatives of the patch geometric description (see Bingham & Maniar [1]). Not only first but also higherorder derivatives of the
potential can be derived, provided the order of the Bsplines is appropriately selected.
No explicit conditions are imposed on the potential at the common boundaries of the patches. This absence of ‘connectivity' between patches greatly simplifies the implementation of the method. No special considerations are required if a patch boundary coincides with a surface discontinuity where the potential is singular, as in the case of the corner at the lower edge of the TLP column in Figure 1.
3.2
The integral equation
The unknown coefficients which appear in (15) may be determined by imposing Green's theorem (7–8). For the sake of generality we define the forcing function which appears on the right side of this equation by , and also the residual r, by
(17)
i.e., the difference between the left and righthand sides of the integral equation. To obtain a linear system of equations for the complex coefficients, a Galerkin procedure is adopted. Specifically, we minimize the residual with respect to each of the (potential) Bspline tensor products directly over their support in the usable parametric space. Thus,
(18)
for i=1, 2, …, I_{p}; j=1, 2, …, J_{p}; p= 1, 2, …, . Each patch provides the same number of equations as the number of unknowns, resulting in a square complex dense system of equations of dimension U×U where
(19)
As a result of the Galerkin procedure, the influence coefficients in the linear system are double spatial integrals of ∂G/∂n. These coefficients are evaluated by a semidiscrete Galerkin approach where the outer integral in (18) is replaced by a fixed order (N_{g}) product Gauss rule, applied to each of the intervals between the knots of the support of the (potential) splines in usable space. For this purpose GaussLegendre quadratures are efficient, and numerical experiments indicate that N_{g}=2 and 3 suffice for linear and cubic Bsplines, respectively. With this algorithm, the field points of the Green function in the inner integration are specified by the Gauss nodes in parametric space.
Special attention is required for the inner integrals, to account for the Rankine singularity of the Green function. When the field point lies within the domain of integration, the singularity is analytically removed and the resulting integrand expanded in a series with algebraic terms which are integrable in closed form. For other cases, and for the integrals of the nonsingular part of the Green function, GaussLegendre quadratures are used. If the field point is near the integration domain, adaptive subdivision is used for the nearlysingular Rankine component to preserve the accuracy of the quadrature scheme. The contribution to the Green function from the free surface is evaluated using the algorithms described by Newman [19]. For further details see Maniar [16].
Similar procedures are used for the integrals of the forcing function, which are of the form
(20)
For the firstorder diffraction potential, = 4π_{I}, and (20) is straightforward to evaluate. In other cases the forcing function is of the form
(21)
where f=∂/∂n for firstorder radiation problems and the integral is over the body surface, and for second or thirdorder diffraction problems with the integral over the free surface. In these cases f is approximated on each panel by a series expansion in the parametric variables, and it follows that
(22)
where (u_{e},v_{e}) is a central expansion point on the parametric space of the panel. The truncation index N_{t} is determined such that the first neglected term is of magnitude consistent with the other local errors.
3.3
Nonlinear methodology
When the integral equation (8) is used for the second and thirdorder potentials, the principal modification is to evaluate the forcing function on the righthandside which involves integration over the free surface. As noted in §2, an appropriate procedure is to use numerical integration inside a circle of finite radius, and a semianalytic analysis in the far field outside the same circle.
Since the inner integration is the most difficult computational task in the loworder panel method, we concentrate on it here and truncate the free surface at a specified partition radius r, and neglect the far field contribution in the forcing function. The resulting solutions are incomplete, but serve to demonstrate the advantages of HIPAN in solving the higherorder boundaryvalue problems. Another simplification which
is made here is to neglect the contribution to the thirdorder forcing function from the secondorder potential. This is motivated by the longwavelength approximation of Faltinsen et al. [4]. Both of these simplifications can be overcome in a relatively straightforward manner, using the same procedures as in the secondorder version of WAMIT.
After solving for the firstorder diffraction potential using (7), the same equation is used (with the factor 4π in the first term) to evaluate ^{(1)} at specified points on the free surface. A leastsquares procedure is then used on an appropriate set of freesurface patches to approximate ^{(1)} with Bspline expansions throughout the domain of the truncated free surface. The required firstand second derivatives in the horizontal plane are then evaluated analytically. The required vertical derivatives are derived using the freesurface condition and Laplace's equation.
Since the above procedure results (locally) in a polynomial form for the functions the integrals on the righthandside of equation (8) can be evaluated using the same algorithms described in Section 3.2. The subsequent procedure to solve (8) is essentially the same as for the first order potential.
A fundamental advantage of this procedure is that it is relatively robust in the vicinity of the body waterline, by comparison to the loworder
panel method. To demonstrate this we first show results from WAMIT where the diffraction potential and its radial derivative are evaluated on the free surface, and compared with the analytical solution for a vertical circular cylinder. (These correspond more physically to the freesurface elevation and slope.)
The relative errors from the loworder computations are shown in Figures 2 & 3. In most cases the error is small, and can be reduced by increasing the number of panels, but close to the waterline intersection the derivative cannot be evaluated in a robust manner, as is true more generally near the body surface (Zhao & Faltinsen [22]). (To achieve maximum accuracy from WAMIT the radial derivative is evaluated using the source formulation, and along the normal to the body boundary at the midpoint of waterline elements.) Clearly this procedure cannot be used to evaluate the second derivatives near the body. In the secondorder solution this problem can be minimized by using Gauss' theorem [11], but this approach cannot be extended to third order. In Figures 4 & 5 we show the corresponding results from HIPAN. Using a similar number of unknowns as in WAMIT, the relative error of the potential is about one order smaller near the body and decays more rapidly in the radial direction. The relative error of the radial derivative is
substantially smaller as well. Near the waterline intersection the error is larger than away from the waterline intersection, but the problem is far less serious than in the low order panel method. By increasing the number of panels or the order of the potential on the body, the error of the potential decreases. While a loss of accuracy is inherent in obtaining derivatives by analytic differentiation of the polynomial representation of ^{(1)}, it can be reduced with a suitable refinement of the Bspline characteristics on the body and on the free surface. The results shown in Figure 5 for use a Bspline representation where the potential is fit at 55×18 points on one quadrant of the free surface in the radial interval 1.0<r/a<1.5. This sector is subdivided into 12×4 panels and Bsplines of order 4 are used. This representation is sufficiently accurate so that the errors shown in Figure 5 are associated only with the firstorder potential representation on the body surface. Similar results have been confirmed for the second derivative of the potential.
The angular variation of the error using HIPAN is shown in Figure 6. It is interesting to note that the error is oscillatory, with maximum errors at the Gauss nodes and at the midpoints between these nodes. The error diminishes with increasing distance from the waterline intersection of the body, both radially and vertically. (Only the real part of the error is shown here since the imaginary component is relatively small.)
The freesurface Green function has a logarithmic singularity when the source and field points coincide on the free surface. This is presumed to be the dominant contribution to the errors discussed above. It is possible to treat this singularity in a semianalytic manner, analogous to the technique described by Newman & Sclavounos [21]. This extension of HIPAN should further reduce the local error at the intersection.
4.
Results
A variety of examples are presented to demonstrate the robustness and efficiency of HIPAN. For the results presented, Bspline representation of patches are constructed by fitting points on the surface by a linear leastsquares procedure. The geometric representations can be considered exact in relation to the approximation of the solutions for the velocity potentials. Generally we start with ‘coarse' representations of the potential using the same knot vectors (13) as in the geometry representation. More accurate potential representations are obtained using knot vectors which subdivide the intervals between consecutive geometry knots. We define panels as the physical surface corresponding to the parametric
Table 1: Definition of commonly used symbols.
Notation 
Quantity 
A 
wave amplitude 
a 
cylinder radius 
d 
halfdistance between bodies 
g 
acceleration due to gravity 
T 
cylinder draft 
λ 
wavelength 
v 
wavenumber 
ρ 
fluid density 
ω 
radian frequency 
space between consecutive knots of the potential Bsplines. Therefore, the number of panels refers to the product, M_{1}×N_{1} for a body described by one patch, or the sum of products ∑_{p}M_{p}×N_{p} for a body described by several patches.
The most significant errors in this scheme are associated with the order (N_{g}) of the product GaussLegendre rule for the outer integration (Galerkin), and the limitations of the potential Bspline approximation (15). The latter is controlled by the total number of splines used, their degree, and their distribution over the region. Since these factors are independent, a simple but effective indicator of the approximating ability of (15) and the ensuing computational effort is the number of unknowns (19) involved. Secondary sources of error include the truncated normal derivative expansion, the amplification of errors in the solution due to an illconditioned linear system of equations, inaccuracies in the Bspline surface patches, and lack of continuity across patches.
Table 1 defines the notation used in presenting the results below.
4.1
Firstorder results
As the first example we evaluate the surge addedmass and damping coefficients for a floating hemisphere, where the benchmark results of Hulme [9] are available for comparison. One quadrant of the hemisphere is modelled as a patch (=1). Symmetry is used to restrict the unknown solution to this patch, where the potential is represented by a cubic (k=4) Bspline expansion. The expansion (22) for the normal derivative is truncated at degree N_{t}=4. The outer integration is an N_{g}×N_{g}=3×3 Gauss rule. Figure 7 compares the HIPAN results with those
of Hulme [9]. Two sets of results are plotted, corresponding to discretizations with M_{1}×N_{1}=2×2 and 3×3 panels, corresponding to 5×5 and 6×6 unknowns. The corresponding curves are practically indistinguishable, except in the vicinity of the irregular frequencies. The bandwidth of the irregular frequencies is very narrow, particularly for the finer discretization.
To illustrate the application to a more complicated problem we next consider the surge exciting force on the Snorre TLP shown in Figure 1. Each quadrant of the TLP is represented by 35 patches. These computations are performed in a finite water depth (For the dimensions and depth see [12]). Two discretizations are used for the potential, labelled ‘coarse' and ‘fine'. The coarse discretization is shown in Figure 1, with panels as large as the entire pontoon side. The fine discretization corresponds to subdividing each panel of the coarse set into four smaller panels. Figure 8 compares the firstorder surge exciting force on the TLP as computed by the present method and by the loworder panel method WAMIT. The graphical agreement with WAMIT is clear. The use of linear potential Bsplines for the coarse case and cubic splines for the fine case, involve 164 and 989 unknowns per quadrant respectively.
We observe from Figure 8 that the maximum value of the surge exciting force occurs near vd/π=1, where the wavelength λ is equal to the column spacing 2d, and minimum values occur near the points where λ/2d is equal to 0.5 and 1.5. These correspond respectively to the condition where the local phase of the incident waves is the same at both pairs of columns, and where the waves are out of phase. This simple explanation, which is well known, essentially follows from considering only the FroudeKrilov exciting force and neglecting the scattered pressure field. The importance of the latter, and the small shift of the maximum force to vd/π0.86 corresponding to λ/2d1.16, is illuminated by considering a long array of N circular cylinders in a single row, with equal spacing 2d between adjacent cylinders.
The first application of HIPAN to a long array was performed for the case of bottommounted circular cylinders in a finite water depth (Maniar & Newman [17]), to permit comparison with the semianalytic interaction theory of Linton & Evans [13]. Computations were performed with up to 100 cylinders, revealing a singularly large surge exciting force (parallel to the array) acting locally on each cylinder, at the critical wavenumber where trapped waves occur for a single cylinder in a channel (Linton & Evans [14]). For N=100 and a/d=0.5 the local surge force on cylinders near the center of the array is about 35 times the magnitude of the force on a single isolated cylinder!
To further illustrate this phenomenon we consider the headsea diffraction by a periodic array of N truncated circular cylinders, of radius a and draft T=a, with the centers separated by five diameters (2d=10a). Figure 9 shows the magnitude of the exciting force on the first cylinder, for N=1, 3, 5, 9. Sharp peaks are evident when vd/π is equal to an integer or an integer plus one half. The peak values increase with increasing N, whereas the bandwidth decreases. Similar behavior has been observed with arrays of floating hemispheres and spheroids.
Figure 10 shows the magnitude of the exciting force on each cylinder in a long array with N=100 elements, at the critical wavenumber near vd/π=1. In this case the force increases along the first 20% of the array, and diminishes thereafter, with substantial sheltering only present near the leeward end. The maximum force is amplified by a factor of four times the force on a single cylinder, and total integrated force on the entire array is about two times the simple estimate based on one cylinder.
Computations of the freesurface elevation show that local ‘sloshing modes' of large amplitude exist at the same wavenumbers where the exciting force on each cylinder is large. The first mode, near vd/π=1/2, corresponds to the trapped mode for a single cylinder in a channel
with rigid walls; in this case the phase of the sloshing mode and local force differs by 180° between adjacent cylinders, and the total integrated force on the array is relatively small. The second mode, near vd/π=1, appears to correspond to a trapped mode for a single cylinder in a channel with ‘Dirichlet walls' where the potential is zero; here the motion at adjacent cylinders is in phase and the total integrated force is large. These trapped modes correspond to eigensolutions of the beamsea diffraction problem in the limit N → ∞. It appears that similar ‘trapped' modes occur for the higher wavenumbers, corresponding to the successive peaks in Figure 9.
To address the issue of computational efficiency, Table 2 compares the number of unknowns required and execution time, for equivalent accuracy, using HIPAN and using the conventional loworder panel method. Two applications are included in this comparison, the heave exciting force on a single truncated circular cylinder and the vertical mean drift force on a submerged sphere. The times reported in the two cases are for a spectrum of 8 and 6 wavenumbers respectively. For a 1.0%–0.1% relative error, the Bspline higherorder method is 10–200 times faster than the constant panel method. If greater accuracy is required, or for more complicated bodies,
Table 2: Comparison of the number of unknowns and run times required to obtain equivalent accuracy, for HIPAN (H) and the loworder panel program WAMIT (W). The results in the upper table are for the surge exciting force on a truncated circular cylinder, and in the lower table for the vertical secondorder drift force on a submerged sphere in finite depth. The figures in parenthesis are for cubic Bsplines; the other HIPAN results are for linear Bsplines. The last column shows the typical execution time, per wavenumber, on a DEC Alpha 3000/700 workstation.
Relative Error 
# of unknowns per quadrant 
Time [seconds] 

W 
1% 
512 
44 
H 
1% 
16 (48) 
1 (3) 
W 
0.1% 
2048 
1140 
H 
0.1% 
32 (72) 
4 (17) 
Relative Error 
# of unknowns 
Time [seconds] 

W 
0.6 ~ 3% 
800/quadrant 
637 
H 
1% 
–(49)/halfbody 
–(10) 
W 
0.05 ~ 1% 
3081/quadrant 
11428 
H 
0.1% 
–(121)/halfbody 
–(108) 
the ratio of the computational effort between the two methods is even larger due to the faster rate of convergence of the higherorder method.
A comparison of the accuracy and convergence of the two programs is shown in Figure 11, for the surge exciting force on the Snorre TLP plotted in Figure 8. This example is intended to provide a more practical illustration of the computational efficiency of HIPAN, but it is somewhat difficult to study the very small errors involved in a systematic manner due to the complexity of the geometry and the difficulty of deriving benchmark values. For WAMIT the benchmarks are established by extrapolation of the computations shown. In the case of HIPAN the benchmark is a separate run with cubic Bsplines and a greater number of unknowns. Uncertainties in the benchmark may be responsible for the inconsistent convergence rates shown. Nevertheless it is clear that the average errors experienced with HIPAN are 10 to 100 times smaller than with WAMIT, with similar numbers of unknowns.
The rate of convergence depends primarily on the order of the (potential) splines, their distri
bution on the body (governed by the use of uniform or nonuniform knot vectors) and the specific physical problem. Our experience to date indicates that for smooth bodies with uniform knots, the error in integrated quantities is proportional to u^{–α} where u is the total number of unknowns (19) and the exponent α 2k–1 (k=order of the splines). In the presence of flow singularities for bodies with corners, the use of higherorder splines based on uniformly distributed knots may result in linear convergence only. In such cases, the convergence rate is expected to improve with the use of nonuniform splines.
4.2
The secondorder solution
Preliminary computations are presented here to verify the effectiveness of HIPAN in solving for the secondorder potential. The problem considered is the secondorder diffraction by a truncated circular cylinder with radius a and draft T=6a. Cosine spacing is used to discretize the cylinder and a uniform discretization is used on the free surface. The freesurface integration in (8) is truncated at a finite radius r=6a to simplify the computations. The resulting potential
Table 3: The secondorder horizontal force as computed by HIPAN and WAMIT. For the results in column 2, 96+86 panels are used on the body and the free surface for HIPAN and 2688+7680 panels for WAMIT. Finer discretizations are used in column 3, with 160+134 panels for HIPAN and 9984+30720 panels for WAMIT.
va 
{F^{(2)}} from HIPAN 

0.10 
–0.07969 
–0.07969 
0.15 
–0.32143 
–0.32145 
0.20 
–0.47781 
–0.47778 
va 
{F^{(2)}} from WAMIT 

0.10 
–0.0760 
–0.0783 
0.15 
–0.3107 
–0.3177 
0.20 
–0.4568 
–0.4699 
va 
{F^{(2)}} from HIPAN 

0.10 
–0.05968 
–0.05967 
0.15 
0.00932 
0.00931 
0.20 
0.36915 
0.36921 
va 
{F^{(2)}} from WAMIT 

0.10 
–0.0460 
–0.0532 
0.15 
0.0321 
0.0208 
0.20 
0.3946 
0.3811 
is incomplete, but since the integration over the truncated freesurface is the most computationally expensive task, we are able to evaluate the effectiveness of HIPAN without the additional effort associated with the farfield integration.
Table 3 shows the resulting secondorder horizontal force, for three wavenumbers and two discretizations. Also shown are the corresponding results evaluated from the loworder panel program WAMIT. For the two Bspline discretizations the HIPAN results are converged to 4–5 decimals. The WAMIT results converge more slowly, although the number of panels is much larger.
4.3
The thirdorder solution
As noted in the Introduction, the problem of ‘ringing' has focused attention on thirdorder wave loads, particularly the highfrequency components in a spectrum or the thirdharmonic in a regular wave system. Two complementary studies of the thirdharmonic force on a circular cylinder have been made by Faltinsen et al. [4] (‘FNV') and by Malenica & Molin [15] (‘M&M'), using different assumptions and boundary conditions. While the magnitude of the thirdharmonic force is similar in these two works, the phase is opposite
Table 4: Real (top) and imaginary (bottom) parts of the thirdharmonic force based on the longwavelength approximation. The numerical solution is obtained by truncating the free surface at the radius r/a.
va 
FNV 
r/a=6 
r/a=12 
0.025 
–9.817E04 
–8.721E04 
–8.661E04 
0.05 
–3.927E03 
–3.610E03 
–3.616E03 
0.10 
–1.571E02 
–1.582E02 
–1.606E02 
0.15 
–3.534E02 
–3.894E02 
–3.931E02 
0.20 
–6.283E02 
–7.321E02 
–7.204E02 
0.025 
0.00 
4.017E07 
3.480E07 
0.05 
0.00 
7.600E06 
4.447E06 
0.10 
0.00 
1.101E04 
–6.292E05 
0.15 
0.00 
9.586E05 
–1.426E03 
0.20 
0.00 
–2.178E03 
–8.090E03 
over the most important regime (0.1<va<0.4).
We are using HIPAN to develop a more general solution of the thirdorder problem which is applicable to practical body shapes including the TLP. This computational approach can be developed following the boundary conditions and assumptions of either FNV or M&M, so that it is possible to develop independent results for comparison with their more analytical solutions. As a first step we simplify the solution of the thirdorder potential by neglecting both the farfield forcing on the free surface (as in the secondorder results described above) and also the contribution to the thirdorder forcing function due to the secondorder potential. Both of these simplifying assumptions are justified asymptotically in the longwavelength regime, as shown by FNV.
The results presented here are for the truncated cylinder with draft T=6a. This is considered to be sufficiently deep so that comparison can be made with the infinitely deep cylinder of FNV.
First, following the longwavelength approximations of FNV, only the second term in the forcing function (6) is retained, and the Green function G=1/r+1/r′ is used to solve for the thirdorder potential. The resulting thirdharmonic force is tabulated in Table 4 for two values of the truncation radius on the free surface, and compared with the corresponding force from the FNV analysis. The computations are relatively insensitive to the truncation radius, since the forcing is confined to the near field in the longwavelength approximation. The agreement
Table 5: The real (top) and the imaginary (bottom) part of thirdharmonic force based on the diffraction analysis.
va 
r/a=6 
r/a=9 
r/a =12 
0.025 
–1.182E03 
–1.161E03 
–1.163E03 
0.05 
–6.469E03 
–6.317E03 
–6.113E03 
0.10 
1.013E03 
–3.141E03 
–3.354E03 
0.15 
3.607E02 
4.447E02 
3.031E02 
0.20 
7.451E02 
7.217E02 
9.383E02 
0.025 
1.158E04 
7.391E05 
6.085E05 
0.05 
3.701E03 
3.805E03 
3.668E03 
0.10 
3.748E02 
3.410E02 
3.907E02 
0.15 
5.654E02 
7.300E02 
6.284E02 
0.20 
9.317E02 
6.980E02 
6.644E02 
with the FNV results is reasonable, considering the finite draft of the cylinder. The imaginary part of the force is zero in the FNV analysis, and relatively small in the computed results.
To correspond more closely with the diffraction analysis of M&M, the procedure described above has been extended by including all of the contributions from the firstorder potential to the thirdorder forcing function (6), and by using the complete freesurface Green function. The results obtained in this manner are listed in Table 5. In this case the convergence with truncation radius is somewhat less satisfactory, except for the longest wavelength (va=0.025). Referring to the results for the largest radius, in the last column, and comparing these with the corresponding results in 4, it appears that the real component changes sign in the same manner as in the results of M&M, and the imaginary component is significant for all but the longest wavelength. These preliminary results appear to support the conclusion of M&M that diffraction effects are significant in the regime va>0.1.
5.
Conclusions
The higherorder panel method described here makes use of Bsplines to describe both the body geometry and the velocity potential on the body surface. The accuracy and efficiency of Bsplines means that the description of the geometry can be practically exact, and the representations of the velocity potential and its derivatives are continuous. Thus wave radiation and diffraction problems can be analyzed more efficiently than with the loworder panel method, and it is possible to develop second and thirdorder solutions with a more robust evaluation of the inhomogeneous forcing function on the free surface.
A variety of examples have been used to illustrate the results of this method, based on the program HIPAN. For firstorder wave loads comparisons of run time and accuracy are made with the loworder program WAMIT. One indication of accuracy is the very narrow bandwidth of the irregular frequencies in the results for a floating hemisphere. The computational efficiency of HIPAN is evident especially for applications where a relatively large number of panels are required. The long periodic array with up to 100 separate cylinders is used to illustrate the feasibility of analyzing a large complex structure which would be difficult or impractical with loworder programs.
The tests for long periodic arrays were originally intended to test the computational capabilities of HIPAN, but the results have given us a more complete understanding of the role of trapped wave modes. In particular, finite arrays of bodies are shown to experience relatively large wave loads at certain critical wavenumbers. These correspond to the existence of trapped modes, which have been associated in the past with wave diffraction by a single body in a channel of finite width.
The higherorder method is particularly advantageous in the analysis of second and thirdorder wave loads. In these problems the computational cost is relatively high, and thus the efficiency of the Bspline technique is more significant. More fundamentally, the continuity of the resulting solution makes it possible to evaluate the freesurface forcing effects in a robust manner, particularly in the vicinity of the body/freesurface intersection. Another possible application where this feature may be useful is in the analysis of wavedrift damping [5].
One more possible advantage of the continuous Bspline representation is in applying the resulting wave loads to a structuralanalysis code. In finiteelement structural analysis it is necessary to use a refined discretization, for example with smaller elements in regions of stress concentration. This mismatch between appropriate discretizations for the hydrodynamic and structural analyses can be overcome easily if the wave loads are described continuously on the body surface.
The problems described here are based on perturbation expansions in the frequency domain,
and on the use of the special Green functions which satisfy the linearized freesurface boundary condition with harmonic time dependence. Similar higherorder techniques can be applied in the context of other applications of panel methods. One example is the solution of the doublebody flow past a ship hull, where HIPAN has been used by Bingham & Maniar [1] to evaluate the socalled mterms which involve second derivatives of the velocity potential. Another example is in the timedomain analysis of wavebody interactions, where parallel work is now underway using the transient freesurface Green function. Similar techniques may be applied also to linear and nonlinear wave problems where Rankine Green functions are distributed on both the body and free surface.
An important practical consideration is the difficulty of using HIPAN. Unlike loworder panel methods where the geometry is described by a set of vertices, or points in Cartesian space situated on the body surface, HIPAN requires a set of Bspline control points which are relatively small in number, but more abstract to interpret. In practice, specialized preprocessor programs must be used in both cases, and we anticipate that suitable software will be developed so that Bsplines can be used by practicing engineers with no more effort than is presently required to input the coordinates of panel vertices.
One feasible extension which offers substantial benefit in terms of ‘userfriendliness' is to automate the process of testing for numerical convergence. In conventional panel methods there is no method for determining a priori if the panel representation of the body is sufficiently accurate in relation to the hydrodynamic parameters of interest; this question can only be answered by repeating the analysis with smaller panels and testing the results for convergence. Several examples of this procedure are described by Newman & Lee [20]. Assuming the Bspline description of a body surface is practically exact, it is relatively straightforward to automatically subdivide patches, or use Bsplines of increasing order to represent the potential, and to test these results for convergence within the program. With such a procedure implemented it would be feasible for a user to specify the required precision of the hydrodynamic parameters of interest, and the program would then ensure that this precision is achieved without the need for the user to perform laborious convergence tests.
In the initial development of HIPAN we were attracted by the flexibility and elegance of Bsplines. However they are but one of a large number of possible basis functions which might be used. An obvious extension is to consider rational Bsplines with nonuniform knotvectors (NURBS). In HIPAN the equivalence of Bsplines to polynomials is exploited in the evaluation of the influence functions, but it appears that this restriction can be removed without substantial complication or loss of computational efficiency. It now seems feasible to develop a more general panel method where the description of the body surface is completely general, provided only that it is defined by a set of algorithms which correspond respectively to subregions or ‘patches’. The only restriction is that the Cartesian coordinates of points on each patch can be transformed by a continuous mapping function onto a rectangular parametric coordinate space. The complete submerged surface of the body would then be described by the union of these patches. In this manner we envisage an extension of HIPAN which is applicable to any continuouslydefined body shape, retaining the Bspline representation only for the solution on each patch and leaving it to the user to define the body geometry in whatever manner is most convenient and exact.
Acknowledgment
This work was conducted under a Joint Industry Project. The sponsors have included the Chevron Petroleum Technology Company, Conoco, David Taylor Research Center, Exxon Production Research, Mobil Oil Company, National Research Council of Canada, Norsk Hydro, Offshore Technology Research Center, Petrobrás, Saga Petroleum, Shell Development Company, Statoil, and Det Norske Veritas. Additional support was provided by the National Science Foundation, Grant 9416096CTS.
References
[1] Bingham, H.B. & Maniar, H., “Calculating the double body mterms with a higher order Bspline based panel method”, Eleventh International Workshop on Water Waves and Floating Bodies . Hamburg, 1996.
[2] DeBoor, C., A Practical Guide to Splines, Applied Mathematical Sciences, 27, SpringerVerlag, 1978.
[3] Eatock Taylor R. & Chau F.P. “Wave diffraction theory—some developments in linear and nonlinear theory,” J. of Offshore Mechanics and Arctic Engineering, 114, 1992, pp. 185–194.
[4] Faltinsen, O.M., Newman, J.N. & Vinje, T., “Nonlinear wave loads on a slender vertical cylinder.” J. Fluid Mech., 289, 1995, pp. 179–198.
[5] Grue, J., & Palm, E., “The mean drift force and yaw moment on marine structures in waves and current.” J. Fluid Mech., 250, 1993, pp. 121–142.
[6] Herfjord, K. & Nielsen, F.G., “A comparative study on computed motion response for floating production platforms: Discussion of practical procedures.”, Proc. 6th. International Conf. Behaviour of Offshore Structures (BOSS '92), Vol. 1, London, 1992.
[7] Hess, J.L. & Smith, A.M.O., “Calculation of nonlifting potential flow about arbitrary threedimensional bodies.”, J. Ship Research, 8, 1964, pp. 22–44.
[8] Hsin, C.Y., Kerwin, J.E. & Newman, J.N., “A HigherOrder Panel Method Based on Bsplines”, Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics, Iowa City, 1993.
[9] Hulme, A., “The wave forces acting on a floating hemisphere undergoing forced periodic oscillations.”, J. Fluid Mech., 121, 1982, pp. 443–463.
[10] Korsmeyer, F.T., Lee, C.H., Newman, J.N. &; Sclavounos, P.D., “The analysis of wave interactions with TensionLeg Platforms”, OMAE Conference, Houston, 1988.
[11] Lee, C.H., WAMIT Theory Manual, Report No. 95–2, Dept. of Ocean Engg., Massachusetts Institute of Technology, 1995.
[12] Lee, C.H. & Newman, J.N., “SecondOrder wave effects on offshore structures.”, Proc. 7th. International Conf. Behaviour of Offshore Structures (BOSS '94) , Vol. 2, Pergamon, 1994.
[13] Linton, C.M. & Evans, D.V., “The interaction of waves with arrays of vertical circular cylinders ”, J. Fluid Mech., 215, 1990, pp. 549–569.
[14] Linton, C.M. & Evans, D.V., “The radiation and scattering of surface waves by a vertical circular cylinder in a channel”, Phil. Trans. R. Soc. Lond. A, 338, 1992, pp. 325– 357.
[15] Malenica, Š. & Molin, B., “Thirdharmonic wave diffraction by a vertical cylinder.” J. Fluid Mech., 302, 1995, pp. 203–229.
[16] Maniar, H., A three dimensional higher order panel method based on Bsplines, Ph.D. thesis, Massachusetts Institute of Technology, 1995.
[17] Maniar, H. & Newman, J.N., “Wave diffraction by a long array of circular cylinders”, Eleventh International Workshop on Water Waves and Floating Bodies . Hamburg, 1996.
[18] Molin B. & Chen X.B. 1990 Vertical resonant motions of tension leg platforms, Unpublished report, Institut Français du Pétrole, Paris.
[19] Newman, J.N., “The approximation of freesurface Green functions”, in Wave Asymptotics, P.A.Martin & G.R.Wickham, editors, Cambridge University Press, 1992, pp. 107–135.
[20] Newman, J.N., & Lee, C.H., “Sensitivity of wave loads to the discretization of bodies.”, Proc. 6th. Conf. on the Behaviour of Offshore Structures (BOSS '92), Vol. 1, London, 1992.
[21] Newman, J.N., & Sclavounos, P.D., “The computation of wave loads on large offshore structures,” Proc. 5th. Conf. on the Behaviour of Offshore Structures (BOSS ‘88), Trondheim, 1988.
[22] Zhao, R. & Faltinsen, O., “A discussion of the mterms in the wavebody interaction problem”, Fourth International Workshop on Water Waves and Floating Bodies. Øystese, 1989.
DISCUSSION
L.J.Doctors
University of New South Wales, Australia
I would like first to say how much I enjoyed the presentation of this work. I believe the careful analysis of the errors, such as that shown in Figures 2 through 6, is a most important and useful aspect of such research because it provides confidence in the method and the computer program.
Figure 7 shows the addedmass and damping coefficients for a floating hemisphere. The wellknown problem of the misbehavior of the panel method in the neighborhood of the irregular frequencies is displayed there. Could the authors comment on how they eliminate this difficulty in practical applications of their work? In my own twodimensional panelmethod analysis of ship sections for a striptheory shipmotion computer program, I have found the use of a “lid” on the internal free surface of a surfacepiercing section to be both easy to program and very effective. Either a stationary or a moving lid seems to be equally good. This is detailed in the publication: Doctors, L.J., “Application of the BoundaryElement Method to Bodies Oscillating near a Free Surface,” Computational Fluid Dynamics, Proc. Intl Symposium on Computational Fluid Dynamics, Sydney, Elsevier Science Publishers B.V., Amsterdam, pp. 377–386 (1988).
AUTHORS' REPLY
We thank Dr. Doctors kindly for his comments.
While the irregular frequencies are a set of isolated frequencies in the continuous problem, their effects in the discrete problem appear as erroneous solutions in the vicinity of those frequencies. The frequency bandwidth of the erroneous solution decreases as the discrete formulation of the problem is made more exact, either by increasing the number of panels in a loworder method or by using higherorder algorithms. The latter is illustrated by Figure 7, where the bandwidth is very narrow compared to typical results from loworder panel methods. Thus, the results in this Figure demonstrate the more accurate approximation to the continuous problem which can be achieved using the higherorder panel method.
There are several known techniques to eliminate the effects of irregular frequencies. Our own experience with the loworder panel method also suggested that putting a “lid” on the interior free surface is effective from the computational point of view. We have refined this technique and used it in the loworder panel code WAMIT to remove the irregular frequency effects from the first and secondorder nonlinear solution. This work was reported in Lee, Newman, and Zhu (1996). The same technique can be easily applied to the higherorder panel method but we have not done so, partly to illustrate the much narrower bandwidth that is present in the latter method.
ADDITIONAL REFERENCE
[24] Lee, C.H., Newman, J.N., and Zhu, X. “An extended boundary integral equation for the removal of the irregular frequency effect,” Int. J. for Numerical Methods in Fluids (in print).
DISCUSSION
W.W.Shultz
University of Michigan, USA

Higherorder panels are known to be more susceptible to singularities. Have you noticed such problems at the corners? What constraints do you put on the patches?

The lines midway between cylinders look like appropriate locations to apply periodic boundary conditions if the wavenumber in the inline cylinder direction is properly chosen. Then, couldn't waves come in many directions?
AUTHORS' REPLY

We do not place any constraints on the solution at the boundaries between contiguous patches. We find that the potential is practically continuous at these boundaries, even in cases where there is an external corner flow. Figure 12, reproduced from [16], illustrates this in the case of the streaming flow past a cube (without a free surface). The equipotential lines in this figure appear to be smooth and continuous within graphical accuracy, except for a small discontinuity which is evident near the corner where the three edges meet.

For a long array, the solution is nearly periodic along the array, with a constant phase shift between adjacent cylinders. Away from the resonant peaks, this phase shift is governed by the longitudinal component of the incidentwave wavenumber and by the cylinder spacing, as assumed by Linton & Evans [23]. However, at the resonant peaks, where the “sloshing modes” are dominant, the phase shift is either zero (Dirchlet modes), or 180° (Neumann modes), in accordance with the requirement that the sloshing modes be continuous between adjacent cylinders. In the latter case, the direction of the incident waves is irrelevant, as suggested by Professor Schultz.
ADDITIONAL REFERENCE
[23] Linton, C.M. & Evans, D.V., “The interaction of waves with a row of circular cylinders,” J. Fluid Mech., 251, 1993, pp. 687–708.