'.e
Author.,
A three dimensional higher order panel method
based ion B-splines,.
by
Hiren Dayalal Maniar
BTech, Marine Engineering, D.M.E.T, 'Calcutta, 1984 MS, Civil Engineering, University of Wisconsin-Madison, 1989
Submitted to the Department of Ocean Engineering
in partial fulfillment of the requirements for the degree of
Doctor of Philosophy in Hydrodynamics
at the
MASSACHUSETTS INSTITUTE OF TECHNOLOGY
September 1995
(:) Massachusetts Institute of Technology 1995. All rights reserved.
a woo woe . wo ow ow. ro ow, no
Department of Ocean Engineering
July 25th 1995
Certified by..
..
J. Nicholas Newman
Professor of Naval Architecture
Thesis Supervisor
Accepted by
A (.
c.
n
A,. Douglas Carmichael
Chairman, Departmental Committee
on Graduate Students
...ASSACHUSETTS INS ITNTE.
TCMI2ISCH2 11176115AMEIT Laboratorium voor
Scheepshydromechanical Archie/
Mekelweg 2. 2e28 CD Delft
Acknowledgments
With much gratitude, I wish to thank Professor J. N. Newman for his suggestions, guidance and support over the past few years. I have benefitted, in no small measure, by his expertise and experience. It is a privilege and an honour to have had Professor Newman as a supervisor.
I am also grateful to the other thesis committee members, Professor J. Kerwin and Dr. F. T. Korsmeyer, for their comments and suggestions. Their constructive criticisms aimed at improving this dissertation was wellappreciated. I also wish
to thank Professor P. Sclavounos for being a part of the thesis committee, before departing on a sabbatical leave.
Many thanks are due to Dr. C.H. Lee, Professor C.Y. Hsin and Ming Xue, for sharing with me their experiences with other panel methods.
This work was supported by the Joint Industry Project 'Wave effects on offshore structures'. This liaison with the industry sponsors was enriching, and their contribu-tions and support are gratefully acknowledged. The sponsors have included: Chevron Petroleum Technology Company, Conoco, David Taylor Model Basin, Exxon Produc-tion Research, Mobil Oil Company, the NaProduc-tional Research Council of Canada, Norsk Hydro, Offshore Technology Research Center, Petrobras, Saga Petroleum, the Shell Development Company, Statoil and Det Norske Verita.sResearch.
For their sacrifices and devotion to our cause, I dedicate this dissertation to my parents, Dayalal and Chandreeka Maniar.
Title: Professor of Naval Architecture
A three dimensional higher order panel method based on
Bsplines.
by
Hiren Dayalal Maniar
Submitted to the Department of Ocean Engineering on July 25th 1995, in partial fulfillment of the
requirements for the degree of Doctor of Philosophy in hydrodynamics
Abstract
A higher order method for solving integral equations over the surface of three dimen-sional bodies commonly encountered in potential flows is presented.
This work is motivated by the drawbacks of the constant panel method, some of which are: the need for large number of panels in situations where the geometric configuration or the hydrodynamics is complex (for e.g., large volume structures, multiple body interactions or bodies with sharp corners), difficulties in obtaining accurate pointwise values like run-up and derivatives of the potential on the body.
The method accepts bodies composed of 'large' patches, each a B-spline tensor product representation, and the unknown potential is approximated by a B-spline tensor product expansion on the parametric space of each patch. The B-splines used to model the potential and the geometry are not necessarily thesame. A local series expansion in the parametric variables i used to approximate the normal derivative
of the potential. Using a semi-discrete GaJerkin procedure, a square linear system of equations for the unknown coefficients of the potential expansions is obtained.
Several examples for flow past bodies in an infinite fluid and for floating bodies are considered for validating and numerically characterizing the method. These include, in particular, the severe test case of a translating cube is used to demonstrate the method's ability to approximate by: localizing and selectively distributing the splines, and using higher polynomial representation, to achieve rapid convergence. Based on the examples studied, the method can be said to provide accurate results with rapid convergence. Derived quantities involvi-g derivatives of the potential are also characterized similarly. It is believed that this is a consequence of the higher degree of continuity implicit in the B-spline potential representation and its derivatives.
In relation to constant panel methods, it
is found that the present method is
less susceptible to the influences of irregular frequencies which affect the solution for surface-piercing bodies.Contents
1
2
3
Introduction
The Numerical Method
2.1 The method 2.2 Discussion
The details of the numerical method
3.1 Notations
3.2 The surface description
3.3 The solution as a B-sdline expansion
3.4 Residual niinimization and the system set-up
3.4.1 Some characteristics of the matrices encountered
3.5 The forcing function
3.6 The influence coefficients
3.7 Other points to note
3.7.1 Definition of a panel and its significance
3.7.2 The choice of the solution B-splines
3.7.3 Integrating the 27r0 term exactly
3.8 The computational effort
3.8.1 A theoretical estimate for the Galerkin assembly effort
3.8.2 Relative estimates of the influence coefficient evaluation effort
3.8.3 Analysis of the computational effort through numerical
experi-17 23 93 26 28 98 30 39 34 36 49 43 45 45 45 46 48 49 . .
3.8.4 Comparing run times with a constant panel method. 60
4 The Boundary Value Problem
624.1 The linearized mathematical formulation
64
4.2 The integral equation reformulation 68
5 The Influence Coefficients
715.1 Self influence-coefficients 74
5.1.1 Notations 75
5.1.2 Analytic desingularization of the integrals 76
5.1.3 On the feasibility of a series expansion Si
5.1.4 Rankine approximations 83
5.1.5 The final series approximations 88
5.1.6 Comments on the overall computations 91
5.1.7 Some related analysis 99
5.2 Near field influence coefficient evaluations via adaptive subdivision . 100
5.3 Far field influence-coefficients 104
5.4 The wave part of the influence coefficien's 105
5.5 Approximating the normal and the Jacobian 106
5.5.1 Series expansions 107
5.5.2 On the convergence, truncation and choice of the expansion point 114
6 Results
1156.1 Motion of bodies in an infinite fluid 118
6.1.1 Translational motion of a sphere 119
6.1.2 Translational motion of a circular cylinder 129
6.1.3 Translational motion of a cube 137
6.1.4 Summary of the infinite fluid results 152
6.9 Motion of bodies in the presence of a free surface 154
6.2.1 Definitions of the derived quantities 156
6.2.2 The motion of a submerged spheroid 158
. ... . -.
.
.. . . . . ... .. ., . . ... .... . .. .B Near field influence coefficient evaluation via an algebraic series
ex-E On converting the surface Bspline representation to a polynomial
representation.
2346.2.3 A hemisphere in surge and heave
, ,
, 1,616.2.4 Heave and surge motion of a rigid configuration of 4 cylinders, 167
6.2.5 Heave and surge motions of a tension leg platform 174:
6.2.6 Freesurface elevation on a bottom mounted circular cylinder.. 176
6.2.7 Horizontal mean drift forces bn a bottom, mounted circular
cylinder.
, 6 , .4! _ 1816.2..8 Arrays of Circular cylinders , 187
6.2.9 Summary of the results for floating/submerged, 'bodies
,
,
1907
Discussion and Conclusions,
191A Creating Bspline patches, via surface interpolation.
.194pansion
B.G.10 Notations
B.0.11 Reparametrization of the surface
44-B.0.12 Series expansions for the integrands 4.4 ,
B.0,13 Expressions for the final .sum
,
B.0.14 Recursion relations for Sit, .
B.0.15 Evaluation of the line integrals, ANB,Nri
B.0.16 On evaluating GmNi
B.0.17 Some related analysis ,41
C Multiplication of univariate polynomials.
D Multiplication of bivariate polynomials.
D.1 Multiplication of two homogeneous polynomials of different degrees D.2 Multiplication of two series in a homogeneous polynomial grouping,
;41 198 200 901 204 206, 210 913 216 991 227 230 230 939 . .
List of Figures
2-1 Linking the physical space, parametric space and the solution surface. 95
3-1 A tension leg platform as described by 35 B-spline patches. The lines
on the surface indicate panel boundaries. 31
3-9 The typical condition number, K., for splines of various orders. The
dependence on the specific problem is very weak. 38
3-3 A shaded plot of the real and imaginary parts of the matrix showing the structure of the dominant elements. The potential spline is cubic
and the matrix size is 81 x 81 39
3-4 A shaded plot of the real part of the matrices showing the band struc-ture. The potential splines are quadratic (top) and linear (bottom).
The matrix sizes are 64 x 64 and 49 x 49. 40
3-5 A shaded plot of the real part of the matrices showing the structure of the dominant elements. The potential spline is cubic and the matrix
size is 346 x 346. 41
3-6 The distribution of the computational effort for a single frequency run The evaluation of the Rankine and wave part of the influence coef-ficients are labelled 'Rankine' and 'Wave' repectively. The Galerkin assembly and solution procedures are labelled `Galerkin' and 'Solver' respectively. The potential spline is cubic, the outer integration is a
3 x 3 Gauss rule, and the wave part of the influence coefficients are
evaluated by a 4 x 4 Gauss rule 54
3-7 The quadratic dependence of the computational effort on the number
of panels. The plots from top to bottom are the: Rankine, Galerkin and Wave efforts. The three curves on each plot correspond to the order of the potential splines.
3-8 The dependence of the computational effort on the order of the
poten-tial spline. The plots from top to bottom are the: Rankine, Galerkin and Wave efforts. The three curves on each plot correspond to fixed
discretizations of 24,96 and 384 panels 56
3-9 The distribution of the computational effort for a six frequency run
The top figure shows the effort with the use of cubic splines and a 3 x 3 outer Gauss rule. The bottom figure shows the effort with the use of linear splines with 2 x 2 outer Gauss rule. The wave part of the influence coefficients are evaluated by a 2 x 2 Gauss rule 59
4-1 The physical domain of the boundary value problem 67
5-1 Triangulation of the parametric space and the quadratic transformation
for each triangle. 77
5-9 Definition of the length scales and the subdivision of the surface. . . 102
5-3 Equipartitioning of the parametric space 103
6-1 Discretizations of an octant of a sphere. The lines on the surface
in-dicate panel boundaries. The set of discretizations are: 1 x 1, 2 x 2,
4 x 4. 121
6-9 Error in the added mass of a translating sphere. The unknowns are
per octant of the sphere. Additional discretizations of 3 x 3 for the cubic, and 16 x 16 for the constant cases are also plotted to confirm
their trend 199
6-3 Error in the added mass of a translating sphere. The change in
con-vergence rates when the normal derivative expansion retains: 3-terms
(top) and 4-terms (bottom) 124
55
6-4 Error in the added mass of a translating sphere. Effects of reducing
the order of the Gauss rule to 9 X 9 126
6-5 Error in the added mass of a translating sphere in an infinite fluid
A limiting case where every panel is a patch: Effects of not imposing continuity across patches. The discretizationsequence is 1 x 1, .1(1 x 1),
16(1 x 1) 198
6-6 Uniform discretization of an octant of a circular cylinder
(Diame-ter/Length = 1). The lines on the surface indicate panel boundaries. The discretizations are: 2(3 x 21 and 2(3 x repectively. 130
6-7 Cosine-spaced discretization of an octant of a circular cylinder
(Diam-eter/Length = 1). The lines on the surface indicate panel boundaries. The discretizations are: 2[3 x 2] and 2[3 x 4] repectively. 1'31
6-8 Error in the added mass for a cylinder translating parallel to its axis.
The plots correspond to the uniform (top) and cosine-spaced (bottom)
discretizations. 134
6-9 The distribution of the potential from the center of the disc to the
midpoint of the side. The plots correspond to the use of linear (top) and cubic (bottom) potential splines. The solid line is tlw reference curve corresponding to a 3 x 16 cosine spaced discretization. The remaining curves on each plot, correspond to the use of a 3 x 2 uniform (solid) and cosine spaced (hollow) discretizations 136
6-10 Uniform discretization of an octant of a cube. The lines on the surface indicate panel boundaries. The discretizations shown are: 3(2 x 2] and
3[4 x 4] 138
6-11 Cosine-spaced discretization of an octant of
a. cube. The lines on
the surface indicate panel boundaries. The discretizations shown are:3[2 x 2) and 3[4 x 4] I 39
6-12 Equipotential contours on an octant of a translating cube. The dis-cretization is 3[8 x 8] panels/octant. Below a close--up of the corner
6-13 Error in the added mass for a translating cube. The plots correspond to the uniform (top) and cosinespaced (bottom) discretizations. . . . 144
6-14 Equipotential contours on a translating cube. Effects of the order of the splines: linear (top), quadratic (center) and cubic (bottom). The cliscretization (uniform) is 3[2 x 2] panels/octant. 146
6-1:3 Equipotential contours on a translating cube. Effects of reduced Gauss order, 2 x 2. Discretizations (uniform) are: 3[2 x 2] (top), 3[4 x 4]
(center) and 3[8 x 8] (bottom) 148
6-16 Equipotential contours on a translating cube. Effects of reduced Gauss
order (2 x 2) on potential splines of order:
linear (top), quadratic (center) and cubic (bottom). The discretization (uniform) is 3[2 x 2]panels/octant 149
6-17 The distribution of the potential from the stagnation point to the mid-point of a side. The plots correspond to the use of linear (top) and cubic (bottom) potential splines. The solid line is the reference curve corresponding to a 8 x 8 cosine spaced discretization. The other curve on each plot corresponds to the use of a 2 x 2 cosine spaced
discretiza-tion. 151
6-18 Fitch added mass for a submerged spheroid 159
6-19 The surge added mass/damping for a floating hemisphere. 163
6-20 The heave added mass/damping for a floating hemisphere 164
6-21 Effects of finer disc:retization on the bandwidth reduction at the first irregular frequency in heave. The potential spline is cubic and the
order of the outer Gauss rule is 3 x 3 165
6-22 Irregular frequency bandwidth reduction. Effects of order of spline (top) and order of outer Gauss integration (bottom). The abbrevia-tions 4c and 3c imply 4 x 4 and 3 x 3 outer Gauss rule. The solid curve
is provided for reference 166
. . .
. . ..
6-23 The configuration of four cylinders placed at the corners of a square.
The lines on the surface indicate panel boundaries. Each cylinder is composed of two patches: the cylindrical side and the bottom disc. In the above example, cylindrical side and the bottom disc patch consist
of 4 x 2 and 4 x 1
panels repectively. The upper boundaries of thecylinders coincide with the free surface 169
6-24 Surge exciting forces on the rigid four cylinder arrangement as shown
in figure 6-23. For validation, results from a constant panel method.
WA MIT, are also included. The side and bottom patches are dis-cretizeci equally and indicated in the legend as 2[M x N] 170
6-25 Heave exciting forces on the rigid four cylinder arrangement as shown
in figure 6-23. For validation, results from a constant panel method,
WAMIT, are also included. The side and bottom patches are dis-cretized equally and indicated in the legend as 2[M x 171
6-26 Exciting forces on a tension leg platform. Comparison with WAMIT. 175 6-27 Discretizations of half of a bottom mounted circular cylinder. The
lines on the surface indicate panel boundaries. The set of geometric
discretizations is: 6 x 1, 12 x 9, 94 x 4 177
6-28 Real and imaginary parts of the run-up on a bottom mounted circular
cylinder from the analytic solution. 178
6-29 Convergence rates for the line integral contribution to the total
hori-zontal mean drift force on a bottom mounted circular cylinder. The top and bottom plots show the use 3 x 3 and 4 x 4 Gauss rule respectively. 183 6-30 Convergence rates for the surface integral contribution to the total
horizontal mean drift force on a bottom mounted circular cylinder. The top and bottom plots show the use :3 x 3 and 4 x 4 Gauss rule
respect i \rely 181
. .
.
6-31 Convergence rates for the total horizontal mean drift force on a bottom mounted circular cylinder. The solid curves are results from WAMIT
The top and bottom plots show the use 3 x 3 and 4 x 4 Gauss rule
respectively 185
6-32 A square array of 4 x 4 floating truncated circular cylinders. Each cylinder is composed of two patches: the circular side and flat bottom disc. The cylinders are spaced uniformly with their centers three radii
List of Tables
3.1 Comparing run times to obtain equivalent accuracy between a con-stant pahel method (WAMIT! and the present method (HIPAN3D). For HIPAN3D, the use of linear and cubic splines is shown (cubic case in parenthesis). The tables, from top to bottom, correspond to the evaluation of the: added mass for a floating hemisphere in surge, the vertical diffraction force on a floating truncated circular cylinder, the heave exciting force on a rigid configuration of four cylinders arranged in a square and the vertical mean drift force on a submerged sphere in
finite depth 61
6.1 Added mass for a cylinder translating parallel to its axis, for a uniform
discretization 135
6.9 Added mass for a cylinder translating parallel to its axis, for a
cosine-spaced discretization 135
6.3 Added mass for a translating cube based on uniform discretizations. 143
6.4 Added mass for a translating cube based on cosine-spaced discretizations.143
6.5 The absolute errors in computing the pitch added mass for a submerged
spheroid. The potential splines used are cubic 159
6.6 The absolute errors in computing the vertical exciting forces on a sub-merged spheroid. The top and bottom tables correspond to the real and imaginary parts of the exciting forces. The potential splines used
are cubic. 160
. . .
6.7 Tabular comparison of the surge added mass (top) and damping
(bot-tom) for a floating hemisphere. 163
6.8 Tabular comparison of the heave added mass (top) and damping
(bot-tom) for a floating hemisphere. 164
6.9 Surge exciting forces: Comparison with a constant panel method,
WAMIT. The potential splines used are cubic 170
6.10 Heave exciting forces: Comparison with a -_-.onstant panel method,
WAMIT. The potential splines used are cubic 171
6.11 Surge exciting forces: Comparing the use of potential splines of various
orders 172
6.12 Heave exciting forces: Comparing the use of potential splines of various
orders 173
6.13 Runup on a bottom mounted
circular cylinder. The top and bottomhalves correspond to va = 1,2 respectively. The real and imaginary
parts are tabulated separately 179
6.14 Runup on a bottom mounted circular cylinder. The top and bottom halves correspond to va = 4, 8 respectively. The real and imaginary
parts are tabulated separately 180
6.15 The error in the horizontal mean drift force on a
bottom mountedcircular cylinder. From top to bottom: the tables correspond to the line integral contribution, the surface integral contribution
and the
total force. The panels and unknowns are for half the body. The potential spline is cubic, and a 3 x 3 outer Gauss rule is used. . . . 1866.16 The convergence of the total heave and surge exciting force on a rigid square array of floating truncated circular cylinders (see Figure 6-32).
From top to bottom, the tables correspond
to an array of 1,4 and
16 cylinders/quadrant. The cylindersare identical with a draft = 3x radius and their centers are three radii apart. The unknowns listed in the second and third table are 4 and 16 times those in the topmost table. The unknowns (topmost table) correspond to a sequence of discretizations: (4 x 2 + 4 x1), (8 x 4 + 8 x 2) and (12 x 6+12 x 3), and the use of linear splines for the potential. The outer integration isChapter 1
Introduction
In a range of situations, fluid flows past bodies can be ma.'hematically modelled by Laplace's equation. In engineering practice, this idealization of real flows finds wide application in the analysis of flow past bodies. The motion of a floating body excited by waves (Newman [26]), flows past propellers (Kerwin [19], Kinnas [201) and various
applications in the field of subsonic aerodynamics (Hess [121) are some examples of the application of potential flow theory. Analytic or closed form solutions for three dimensional potential flows exist for a few simple geometries, but for more general shapes numerical solution of Laplace's equation is imperative. Some of the numerical methods commonly used for analyzing potential flows are: finite differences, finite elements and boundary integral equations. For exterior flow problems,
due to the
infinite extent of the fluid domain, methods like finite differences and finite elements can be computationally expensive. In contrast, boundary integral equation methods only require the discretization of some of the bounding surfaces of the domain, and compared to volume discretization methods the computational effort may be smaller. In this chapter we review some representative approaches for the numerical so-lution of integral equations commonly used for potential flows. In particular weexamine: the constaut panel method, a higher order method based on
Lagragian interpolating polynomials, a higher order method based on Taylor expansious, and a two dimensional higher order method based on the use of B-splines. We do nottegral equations. For a comprehensive review of panel methods we refer the reader to the works of Atkinson [3], Hunt [16] and Romate [38].
Constant panel methods stern from the pioneering work of Hess and Smith [10] (see also Smith [39] and Hess [12]). The continuous body surface is approximated by flat triangular or quadrilateral panels, and the unknown singularity strength and normal derivatives are assumed constant on each panel. Collocating at the centroid of each panel, the discrete integral equation gives as many equations as unknowns, and a square linear system of equations is obtained. The coefficients of this system of equations are integrals of the derivatives of the Green function over flat panels. The choice of the Green function can vary from the simple inverse distance function to more complex ones (see Equations (4.17)-(4.18)). The evaluation of these integrals require substantial computational effort, and efficient algorithms must be devised. On solving this linear system of equations, the singularity strength on each panel is obtained. Most quantities of physical interest can be obtained from the singularity distribution on the body. To obtain the first derivatives of the potential on the surface the source formulation is used.
For constant panel methods, increasingly accurate solutions are obtained through a finer and selective discretization of the body. By reducing the size of the panels, the geometric description of the body is more accurate and the error in the assumption that the singularity strength is constant on each panel, is also reduced.
If N is
the total number of panels, the computational effort in setting-up the system of equations and for its (iterative) solution are both 0(N2). Typically, a sequence of runs with increasing N are used to determine whether the results are converged. Derived integrated quantities which are linear functionals of the potential, usually converge at a rate inversely proportional to N.
The advantage of the constant panel method is its simplicity and the relative ease of its implementation. The main drawbacks are: large number of panels are required in situations where the geometric configuration and/or the hydrodynamics is complex, for wave problems the method is susceptible to the influences of irregular frequencies, and obtaining accurate derivatives of the potential is difficult.
The use of a large number of panels leads to an increased computational effort for the set-up of the system of equations and its solution. Large volume structures', multiple bodies and bodies with regions of high curvature are examples where many panels are required to describe the configuration. In many instances, the rapid vari-ation of the solution requires a finer and selective discretizvari-ation of the body. Some examples are: bodies with sharp corners, the pointwise evaluation of the potential (run-up) at the intersection of the body with the free surface, and multiple body interactions.
In some instances, for example, second order problems in the frequency domain or for the evaluation of m-terms for the time domain problems, second derivatives of the potential are required. The procedure of differentiating the integral equation for the potential cannot be used to obtain these.
For wave problems, the use of the Green function (see Chapter 4) results in the in-tegral equations being non-unique at a discrete set of frequencies (called the irregular frequencies) which coincide with the eigenfrequencies of the interior Dirichlet
prob-lem (Kleinman [21]). For a finite band of frequencies centered about this discrete set, constant panel methods give inaccurate results. For problems where derived quanti-ties are rapidly varying functions of frequency, to distinguish the true from spurious behaviour can be problematic. To reduce their corrupting influence requires a filler discretization of the body, and to eliminate their influence requiresthe use of modified integral equations (Lee and Sclavounos [22]).
A promising and fundamentally different method is the '0( N)'approach of Nabors, Korsmeyer, Leighton and White (251. The method is efficient with the usc of very large number of panels (much larger than is possible with conventional panel meth-ods), but it presently is restricted to problems solvable with the useof simple Rankine
singularities.
In an effort to achieve equivalent accuracy for a smaller number of panels, higher
order panel methods have been proposed in the past. The usage 'higher order'
call intended to imply better surface and function representations than the cons t ant panel method. Based on the surface and unknown descriptions, we outline two ap-proaches.
In the first, the surface is discretized by curvilinear triangular (Atkinson [3)) and/or rectangular [46]) panels or elements. Each panel is described by fitting the nodes (points on the panel: panel vertices, boundary midpoints and perhaps, an interior point) by an expansion using interpolating polynomials as their basis set. The number of terms in this expansion is the same as the number of nodes on the panel. These interpolating polynomials are quadratic in the local triangular (Atkinson [31) or rectangular (Xii [46]) parametric space of the panel. At any node all but one of this basis set vanish, therefore the Cartesian coordinates of the panel nodes are the coefficients of the expansion. While there are no 'gaps' at the common nodes of the neighbouring panels, the normal or the Jacobian must be discontinuous across panel boundaries since their evaluation reauires differentiating the surface representation.
The unknown potential over each panel is approximated using the same basis set as used for the panel description. Therefore, the coefficients of these local expansions are the unknown functional values at the nodes of the panels. To determine these unknown coefficients, a collocation scheme is employed where the nodes of the panels are chosen as the collocation points. The resulting square linear system of equations are solved for these coefficients. A similar approach is utilized by Chau [6] for solving the second order diffraction problem.
The use of quadratic interpolating polynomials ensures functional continuity at the common nodes of the neighbouring panels. To obtain the derivatives of the potential, Xii [46] utilizes a finite difference scheme, while Chau [6] obtains the first derivatives by directly differentiating the polynomial approximant. Presumably, on
smooth surfaces, the Cartesian2 derivatives are discontinuous across the elements. Overall, these schemes provide better accuracy than the constant panel method for a fixed number of unknowns or panels.
A different higher order approach is developed in the works of Johnson el al.
[17] and Hess [Ht. In their scheme, paraboloidal' panels are constructed from a mesh of points lying on the true surface by a local weighted least squares fit to the mesh points in the neighbourhood of the panel. On each panel, a linearly varying source distribution and a quadratically varying dipole distribution are proposed in the tangent plane variables. The choice of the degree of these polynomials, for the panel geometry and the singularity distribution, is determined by a local truncation error analysis. The unknown coefficients of the higher order terms in these distributions are analytically related by a weighted least squares scheme to the leading order coefficients
(the singularity strength at the panel `centroid') on that panel and its immediate
neighbours. The discrete intregal equation results in a linear system of equations forthe unknown singularity strength at the centroid of each panel. Hess [11] notes the superiority of the higher order scheme over the constant panel method for interior flows with an example of flow through a duct of varying cross-section. The results of Johnson et al. [17] include a case of random discretization of the surface to illustrate
the insensitivity of the higher order scheme to panel distributions.
Higher order formulations using B-splines are not common in engineering. A more recent investigation of their use in two dimensions is the work of Hsin et al. [14]. An isoparametric scheme is devised where the two dimensional body profile and the unknown potential are all modelled by B-spline expansions in the global parametric space of the body. The use of B-splines in such a fashion gives the potential expansion and its derivatives several degrees of continuity everywhere on the body. Periodicity is imposed for the non-lifting flows by equating some kif the unknown coefficients. The performance of the scheme using B-splines of various orders are presented. In general, the use of cubic B-splines show significant improvement in the results when
compared to those for constant splines.
In this thesis, we present a higher order method for solving integral equations over bodies in three dimensions using B-splines. In doing so we are primarily motivated by (a) the need for obtaining derivatives of the potential on the body and (b) the
need to overcome the several drawbacks of constant panel methods.
In two dimensions, simply connected bodies can numerically always be 'globally' parametrized. In three dimensions, for arbitrarily shaped bodies it may he difficult, if not impossible, to create a global parametric space. Once again, Figure 3-1 (pp. 31) illustrates the futility of such a task. However, complicated bodies can be subdivided into large patches, each of which is parametrically described. All approximations which follow are conducted on the parametric space of each patch.
In particular, the method accepts bodies constructed by patches, each a B-spline tensor product expansion (as constructed by surface interpolation or from some geo-metric modelling program). On the parageo-metric space of each patch the potential is approximated by a B-spline tensor product representation, not necessarily the same as the B-splines used to describe the geometry, and the normal derivative is approx-imated by a series expansion in the parametric variables. To obtain a square system of equations a semi-discrete Galerkin scheme is employed, where the residual is min-imised with respect to the (potential) B-splines directly in the parametric space. Overall, the method can be characterized as being very accurate, efficient and the linear systems which must be solved are an order of magnitude smaller than constant panel methods.
In the chapter to follow, we first outline the method with a broad perspective and discuss the several advantages of the method. In chapter 3, we provide the details of the method, and also discuss the computational effort involved in the method. We apply the method to the frequency domain analysis of the linear wave-body interactions, and the corresponding boundary value problem is described in Chapter
4. In Chapter 5, we describe the evaluation of the Rankine and wave part of the
influence coefficients. Following which, the results for flow past bodies in an infinite fluid and for floating bodies are presented in chapter 6. Based on the results presented, the conclusions drawn are discussed in Chapter 7.
Chapter 2
The Numerical Method
2.1
The method
We consider the numerical solution of the Fredholm second kind integral equation
2ircb
if
an
dS,s
aG
(2.1)
commonly encountered for 'exterior potential flow problems (Atkinson PI, Newman [30]). Above, .cf), is an unknown function (real or complex) on some surface S in three dimensional space, G the Green function which satisfies Laplace's equation and
a
boundary conditions,, except those on. the surface, S., Often it is of the form1
G H
where 1/R is the inverse distance function and H some smooth function (see Chapter 5). And, .F some prescribed forcing, also encountered in the form
=
f
dS (2.2)where f is some function defined on the. surface S..
The method allows the surface S to be composed of several 'large' patches, each
surface S. Each patch is a parameterized surface with its Cartesian components given by a B-spline tensor product expansion (as constructed bysome geometric modelling package or by surface interpolation outlined in appendix A). In practice, these expansions can describe a surface accurately and can be considered 'exact'.
On the two dimensional parametricspace of each patch, the solution, 95, is viewed as another surface, the solution surface, which is modelled by proposing q5
as a
B-spline tensor product expansion in the parametric variables. The unknown coefficients of this expansion are to be solved for via the above integral equation.The 'patch-wise' mapping of the physical surface and the solution surface to a rectangular parametric space is illustrated in Figure 2.1. The entire surface is composed of a circular cylindrical side and a flat bottom disc, but only a quadrant of each of these surfaces shown.
To obtain a system of equations for the unknown coefficients of the solution B-spline expansion, we employ the Galerkin procedure. Denoting the approximate
so-lution by cb. and the (possibly) approximate forcing by we define the residual as aG
r = 27r0
ff
q5dS
-and minimize it with respect to the (solution) B-splines tensor products. Over every patch, the solution is a sum of the product of an unknown coefficient and a B-spline tensor product, therefore we have as many B-splines to minimize the residual with as unknowns. Thus, every patch provides exactly the same number of equations as
unknowns over that patch, resulting in a square system of equations. In practice, the Galerkin approach is a semi-discrete one where the outer integration is replaced by a product Gauss-Legendre rule.
PATCH 2
X = X (u, v ) 2 2 (P (p1= q) (u, v ) 1 V ---/-1 _--- ,,__---X = __---X (u, v)
1 1 (I) q) = 2 cp (u, v)2 U2.2
Discussion
There are several features of the method worth noting explicitly:
Describing piecewise smooth' surfaces in a patchwise fashion not only facilitates their construction, but also allows more accurate descriptions.
The presence of a rectangular parametric space offers several advantages: (a) Besides B-splines it allows for other tensor product expansions for the solution. (b) It allows the choice of solution approximant to be independent from that used for the surface description. I.e. the method does not have to be isoparametric (see also Section 3.3). As a consequence, increasingly accurate approximations for the solution may he sought without re-constructing the surface. (c) The uniform rectangular space could be used for parametric finite difference approximations for derivatives of functions defined on the patch.
For surfaces composed of several patches, no explicit conditions are imposed on the solution at or across the common boundaries of the patches. To impose such
condition? in two dimensions is possibly trivial, but in three dimensions it is not. Es-pecially, if the surface consists of patches topologically rectangular and triangular (see Figure 3-1, P. 31). This absence of 'connectivity' between patches greatly facilitates the implementation of the method.
Implicit in the solution representation is the fact that an expansion based on kth.-order B-splines has (A: 2) degrees of continuity2 ensured everywhereon the patch. Therefore, with an appropriately chosen order, even higher parametric derivatives of the solution retain some degree of continuity everywhere on the patch. Consequently, depending on the degree of the patch descriptions', higher Cartesian derivatives of the solution also retain some degree of continuity on the patch.
An advantage in using the Galerkin procedure is that derived quantities of the form, fjcbgd.S (g is some smooth function), based on Galerkin approxirnants, 0, can
'Surfaces with corners or with abrupt changes in their curvature. 2Assuming no multiple knots. See Section 3.3.
3The Cartesian derivatives of a function defined on a surface in three dimensions, require para-metric derivatives of the surface representation.
exhibit convergence rates higher than theoretically predicted rates (Sloan [36]). The dominant sources of error in the method are: the limitations of the function representation by B--splines, those introduced in replacing the outer Galerkin inte-gration by product Gauss rules, and the right hand side representation. We note that the geometric descriptions can considered exact to single precision accuracy . The
integrals of the Green function are accurate to single precision, with no simplfying assumptions for the surface or its boundaries.
Chapter 3
The details of the numerical
method
3.1
Notations
Some variables we commonly refer to are:
S: The entire surface.
p: Index indicating patch under consideration. Sp: A patch on the surface, S.
u, v: The parametric variables on some patch. u,, ye: Some point on the parametric space. k: The order of the B-splines.
The ith. B-spline of order k in the parametric variable u.
The jth. H-spline of order k in the parametric variable v.
Total number of segments the u-usable parametric space is divided into. Total number of segments the v-usable parametric space is divided into.
M = M k
1= N + k 1
uP ,v P The mth. and nth. knot in the completel u and v parametric space, on patch,
rn
p.
vP: The complete knot vector in the u and v parametric space, on patch, p. At.): The parametric space, [u,+k, x [112+k, Vi+k+1)
The usable' parametric space spanned by the Bsplines tensor product, U,V,.
The above declarations are mainly related to Bsplines in a general way. To draw a distinction between those used for the geometry and the solution, we 'top' those associated with the geometry with a tilde ( ). The parametric variables u and v are exceptions since the usable parametric space is the same for both.
Subscripts on u and v are used to indicate knots. An exception being (ue,t)e) which is used to refer to some point on the parametric space.
Any variable appended with the index p implies that it is specific to patch, S. Su-perscripts used on the Bsplines U and V do not refer to the order of the splines, hut, the patch under consideration.
Boldface type variables imply a vector or matrix, the distinction will be evident from the context it is used in.
3.2
The surface description
The surface cver which the integral equation is to be solved, is described compositely by P patches
S = S, + S, +
. .. +S,
where each patch, S, is described parametrically by a B-spline tensor product ex-pansion (Rogers [37j)
' P Sip
XP(U, V) = E Ex:',0,P00-/:(v)
(3.1)1.1
where .A-47, =A1 + k 1, = k - 1, M,, and /Cip art the non-zero intervals'
between knots in the u and v parametric directions respectively, UP and 17" are B-splines of order k in the parametric variables u and v respectively and XP,, are the known vertices or coefficients.
Figure 3-1 is an example of a piecewise smooth three dimensional surface,
com-posed of 35 patches each described by a B-spline expansion in the form of Equation (3.1). The lines on the surface are not necessarily patch boundaries, but panels which we define later.
The rectangular parametric space which is mapped by (3.1) onto the physical surface is
-P r -P
11,-"k) X [yr,, )
where UP and I are the mth. and nth. knots of the knot vectors, tiP and . We will
refer to this as the usable parametric space. Note that the entire parametric support of the B-splines in Equation (3.1) is
)r
X 7 F+2i_i ).This is also the complete parametric space for the geometry B-splines.
'To keep the discussion simple, we preclude multiple knots in the interior usable space. In our implementation multiple knots for the geometry Bsplines is permitted.
Figure 3-1: A tension leg platform as described by 35 Bspline patches. The lines on the surface indicate panel boundaries.
3.3
The solution
as a Bspline expansion
On the rectangular parametric domain of each patch, S, we approximate the solution by a B-spline tensor product expansion
A 4, Alp
0 ( 1 I)) = E E OP,,11,* (u)1/ (v) (3.3)
)=1
where Alp = M, k 1 and Al, = k
I, M, and
IV, the non--zero intervalsbetween knots in the usable parametricspace of u and V respectively, OP,, the unknown coefficients (real or complex), LI ,11 and 1/,* the B-splines of order k in the parametric variables u and v respectively.
The B-splines used above and those in Equation (3.1) are not necessarily the same, i.e. not necessarily isoparametric. They can differ not only in theirorders, hut also the knot vectors there are based on. The only restriction is that the usable rectangular parametric domain implicit in Equation (3.3) must match that. of the geometry, Equation (3.2), i.e.
P [UPk UPNi k X [V k VNp+k
" ) = [U
k hip+- ) XCertain discontinuities in the derivatives of the B-splines or the B-splines themselves can be introduced by using multiple knots. There is no reason to introduce discon-tinuities in the solution a priori, unless dictated by physical considerations. For the results presented in this thesis, the B-splines used for the solution have no multiple knots in the interior of the usable space.
From a surface interpolation standpoint, how well the proposition, Equation (3.3), approximates the solution surface over the patch depends on the smoothness of the solution, the location of the knots' and the parameters: M,,N, and k. On a patch, with the order of the B-splines fixed, by increasing M, the span over which the
(3.4)
'In regions of rapid solutioi. variation, non-uniform B-splines can be used improve the
approx-imation. Non-uniform B-splines are obtained by using non-uniform knots in the usable space. To
keep the discussion simple, we consider only uniform knot vectors in the usable space of the patch.
In our implementation we allow B-splines based on non-uniform knot vectors.
B-splines are non-zero are effectively reduced (localized). Further, additional B-splines are introduced in the expansion, Equation (3.3).
Both of these lead to a
better approximation of the solution, the former is particularly effective for non-smooth solutions (for example corner flows). An alternative node of approximation
is to keep M, N, fixed and raise the order of the B-splines, the advantage being
higher continuity in the solution and the disadvantage, for low M and Np, is that the expansion (Equation (3.3)) tends to be a global approximant. Often, a combination of both these modes of approximation is the hest approach.3.4
Residual minimization and the
system set-up
For the integral Equation, (2.1), defining the residual as
r(u, v) = 270 +
dS
-dii
we considt. setting up a linear system of equations by minimizing the residual (a) pointwise (collocation) (b) with respect to the solution B-splines (Galerkin).
A commonly used approach to set-up a linear system of Equations is the colloca-tion method. For constant panel methods, an unknown is associated with every panel and cent roidal collocation on each panel provides as many equations as unknowns. In the current method, collocation at the mid-points of the intervals between knots in the usable space, will not provide sufficient number of equations. Specifically, midpoint collocation oil the intervals provide M, x N, equations for (A19+ k- 1) x k - 1) unknowns on each patch, S9. Therefore, unless additional equations are generated by choosing some additional collocation points or some conditions imposed at or across the boundaries of the patches to eliminate some unknowns, a square system cannot be obtained. It is well known (Sloan [361) that the location of collocation points has a bearing on the accuracy of the solution. Numerical evidence suggests that choosing the zeroes of orthogonal polynomials as collocation points provides better accuracy than other locations. Imposing constraints on the solution at or across patch edges requires anticipating the behaviour of the solution qualitatively. If the patch edges coincide with surface discontinuities (e.g. corners), besides functional continuity, con-ceiving additional conditions a priori may be difficult. In addition, if the method must accept arbitrarily shaped bodies, for e.g. Figure 3-1 (p. 31), such an approach is not practical.
In a robust alternative, an overdetermined linear system of equations is obtained by setting the residual to zero at X; x H, collocation points on each of the Ali, x
intervals between knots. A least squares scheme is then employed to solve this system of equations. Numerical experimentation with the location and number (.A/9)
of collocation points indicates that their choice can affect the accuracy of the solution.
Significantly better results are obtained by collocating at the zeroes of Legendre polynomials in each parametric direction. For a cubic solution B-spline expansion, 3 x 3 collocation points on each interval is a practical choice. An advantage of the approach is that the assembly of the rectangular matrix is simple compared to the Galerkin approach. The drawbacks which outweigh it are: the large computational effort required in solving the system of equations and the large storage requirements for the rectangular matrices encountered.
In the Galerkin approach, we minimize the residual with respect to the (solution) B-splines directly in the usable parametric space of its support
Jfr(u, v) U: Vi° du dv = 0 (3.6)
and since Equation (3.3) has as many unknowns as B-splines, each patch provides as many equations as unknowns. The resulting linear system of equations is of size, C x C, where
C = E(N.c+ k - 1) x (N, + k - 1).
(3.7)p=1
In practice, we take a semi-discrete approach where the outer integral is replaced by a product Gauss rule. Using an Ng th. order product Gauss rule on each interval
(space between consecutive knots) of the B-spline support
141,, V,
E EE
r(uns, v,$) U:(urn) V(v) = 0.
(3.8)oP m.1 r.=1
Above, Efe implies a sum over every interval on the usable domain of çl, and
wi,ui, vi are the weights and nodes of the Gauss rule.The error introduced due to the semi-discrete step, Equation (3.8), depends on the smoothness of the residual on the surface, the choice and order of the Gausssian quadrature. By definition, the residual is the difference between the left- and right-hand sides of the integral equation. Therefore, its "behaviour" is prescribed when
(Equation 3.10) are fixed. In theory, it is always possible to impose Equation (3.6),
f f r
V dudv = 0, accurately by a combination of sub-division of the domain of integration and/or raising the order of the Gauss rule. In practice, irrespective of the smoothness of the residual, we apply a Gauss rule of fixed order over every patiel, and thereby possibly introduce errors in imposing f f r I/ V dudv = 0 inaccurately. We refer to these inaccuracies in the system set-up as errors introduced due to the outer Gauss integration.Conclusions drawn from numerical experimentation with Legendre, Gauss-Chehysliev and other open rules is that Gauss-Legendre performs significantly better than other rules. For (solution) cubic B-splines or splines of a. lesser degree, a 3 x 3 product Gauss-Legendre is sufficient (see Section 3.7.3 for further comments).
3.4.1
Some characteristics of the matrices encountered
For a 3 x 3 outer Gauss rule, Figure 3-2 shows the typical condition numbers for the matrices encountered in the present method. The curves shown correspond to the use of various orders of the potential splines. As the order of the potential spline is raised, the condition numbers get larger. The plot also indicates
that that the
condition number is practically independent of the discretization. We also note that the estimates remain unchanged if a higher order Gauss rule is used. However, for a lower order Gauss rule, the estimates will change to even larger condition numbers.Since the matrices encountered are very ill-conditioned, the use of iterative meth-ods for solution of the linear system of equations may seem doubtful. In Section 3.8.3, we note that the linear system of equations are solved using a direct, method (LU de-composition), and that the solution effort constitutes only 1% of the total execution time. Therefore, at the present time, using LU decomposition and backsubstitution for the solution of these ill-conditioned system of equations pose no problems.
To help identify the location of the dominant elements of the matrices encountered, we plot the magnitude of the elements in Figures 3-3-3-5. The intensity of the shades reflect the magnitude of the normalized elements (normalized by max la,j1, where a, are the real or imaginary elements of the matrix). Since the magnitudes of the
elements vary widely, the scales shown at the sides are exponential. The dominant elements show up as dark patches, and the elements which are nearly zero show up as white patches.
In Figure 3-3 we plot the complex matrix encountered for a floating hemisphere case (see Chapter 6). The real (top) and imaginary (bottom) parts are shown sep-arately. A certain band structure is evident for the real matrix, and the number of
bands are presumably related to the order of the potential spline (cubic). For the same physical problem and the choice of linear and quadratic potential splines, Figure 3-4 show the number and spread of the dominant bands (only the real part of the matrices are shown). In general, the matrices encountered are not symmetric, and this lack of symmetry is seen in Figure 3-5. This matrix is encountered for a geometry similar to the tension leg platform, but without the rounded edges at the pontoon corners (see Figure 3-1).
10'
to"
104 10'
1/System Size
Figure 3-2: The typical condition number, n, for splines of various orders. The dependence on the specific problem is very weak.
- - - 0 Unmet 1/r 10-2 10-4 10 104
----va.0 -5 Cubic va.0.5 Cub:1.40c - vs.4.0 Ouadniek va.0-5 Lineal 104REAL 1.2E0 1AE-1 52E2 1.9E1 72E-3 1.1363 33E-4 1.4E-4 51E5 14E-5 72E-6 1.0E-6 IMAGINARY LCE0 33E1 L4E.1 52E1 1.9E-2 7.2E.3 2JE3 17E-4 1.4E4 52E-5 I.965 2.2E4 2JE-6 11E4
Figure 3-3: A shaded plot of the real and imaginary parts of the matrix showing the structure of the dominant elements. The potential spline is cubic and the matrix size
2.7E3
2.7E4
: -10,..-. -...: n5 4... "1.:,..., :.:.:?..:;, 194 40 4! II Ir. F.... 4." REAL 1.0E0' 3.7E-1 1.4E1 5.2E-2 1.9E2 7.2E 27E-3 1.0E-3 3.7E-4 1.4E-4, 5.2E-5 1.9E.5 7.2E-6 2.7E-6 1.0E'43 REAL .0E0 5..1 3.7E-1 1.4E-I 5 5.2E-2 1.9E-2 7.2E-3 IV! 2.7E-3 I.0E-3 3.7E-4 1.4E-4, 1 9E5 72E-6" 2.7E-6 1.0E -6
Figure 3-47 A shaded plot of the real' part of the matrices Showing -Ole 'band structure. The potential splines are quadratic (top) and linear (bottom). The matrix sizes are,
64 x 64 and 49 x 49.
I
'4 .
t Et:
...
mommovisiditarnsfremr Vrof MO tiff.
Figure 3-5: A shaded plot of the real part of the matrices showing the structure of the dominant elements. The potential spline is cubic and the matrix size is 346 x 346.
r-i REAL 1.0E0 3.7E-1 1.4E-1 52E-2 1.9E-2 7.2E-3 2.7E-3 1.0E-3 3.7E-4 1.4E-4 5.2E-5 1.9E-5 7.2E-2.7E-6 1.0E-6 ....;.:,
1 1
=
,f G dS At, Jug3.5
The forcing function
Due to the Galerkin procedure, the right hand side of the systeth Of equations, requires the evaluation of
If .1 is computationally inexpensive to evaluate, product Gauss rules and/or subdi-ViSion Of the integration domain are easily applied for its accurate evaluation,. Corn= pared to the computational effort required to set-up the global matrix, this requires negligible effort.
Instead if is of the form,. Equation (2.2)
PMp 1 Np
fOds
E E Eff f G dS
AP
p= 1 i=0 j=0
then not only is the effort required in evaluating .F non-trivial, but greater care iS required, due to the presen.oe of the singularity in. G for some of the integrals on the right hand side above. In softie instances, it may be possible to replace f locally with a series expansion in, the parametric. variables (see Section, 5.5)
co in
E EfJ tuu,ey (v va)
G (3.10)m=0
which is convenient for the analysis and evaluation of the singular and nearly-singular integrals of G.
For Neumann problems,, f DOlan, the normal, derivative of the potential. In general, for curved surfacesa series expansion for the normal derivative has an infinite number of terms, arid in practice must be truncated (see Section 5.5, p. 112). For a fixed truncation error, the number of terms required in this series approximation is dependent on the size (i.e. the domain of integration in Equation '(3.10)), and the
intrinsic geometry of the surface..
if
.13,7 .7" du dv. (3.9), nP -= = dS3.6
The influence coefficients
The Galerkin procedure requires the evaluation of double surface integrals of the Green function. With the semi-discrete Galerkin approach, Equation (3.8), the Carte-sian coordinates of the Gauss nodes of the outer integrals provide the field points for the inner integrals
u v
-aG dSflUV G dS
If
where the domain of integration is the support of the solution B-splines, U and V, on the curved surfaces. With the definition of panels (see Section 3.7.1, p. 45).
we view the support. of the B-splines on these surfaces as a collection of panels. Further, instead of evaluating the above integrals over the support of the B-splines, we
consider their evaluation over each panel individually. On each panel substituting the appropriate polynomial form of the B-splines, we encounter integrals (the influence coefficients) of the form
f .1(u
- u.)"' (v - vc)"
dSan
ff(u - Tier (v - vc)" G dS
m,n = 0,1,...,(k - 1)
where (u v,) is some point on the panel. The reason for evaluating the integrals
over panels as opposed to the support of the splines is that the B-splines can be
converted to polynomials. For the near field algorithm we adopt, the computational effort in evaluating the integrals with polynomial moments is smaller than the oneswith B-splines.
To maintain the continuity of this chapter, the details of the influence coefficient evaluations are provided in Chapter 5. Here we note that algorithms have been devel-oped for evaluating the dominant Rankine part (1/R)over curved surfaces. The near field cases are evaluated by an adaptive subdivision scheme. For the self-influence
(3.11)
-panded in a series with algebraic terms which are integrable. The integrals of the wave part of the Green function are evaluated using a product Gauss-Legendre rule. For the functional evaluations of the wave part, algorithms as described by Newman
3.7
Other points to note
3.7.1
Definition of a panel and its significance
We define panels as the space between consecutive knots of the solution B-splines (recall that the solution B-splines have no multiple knots in the interior of the usable space, p. 32).
Since the interior of the usable space of the patch is partitioned by (M, 1) and (N, 1) knots of the U and V splines respectively, there are exactly M x N, panels on patch, S. For convenience, we refer to this panelization of the patch or body as its discretization. Thus, when we refer to a finer discretization of the body or patch, we imply raising M, N in Equation (3.3).
The indirect significance of this definition is: (a) For solution splines built on uniform knots in the usable parametric space, a graphical display of thepatch/body with its discretization, will convey the span of the spline supports, and therefore the extent to which they are local. (b) For solution splines built on non-uniform knots in the usable parametric space, a graphical display of the patch/body with
its discretization, will indicate which regions of the patch or body the splines are 'concentrated' and also the span of their support.
Note that the support of a tensor product spline is a panel only if the spline is a constant (k = 1). For splines of order, k > 1, the support of a tensor product spline is not limited to a panel, but k2 panels. This observation only applies to splines which lie entirely in the usable parametric space of the patch, otherwise, the support of the tensor product spline is less than k2 panels. Therefore for a uniformly discretized body, the size of the panel relative to the patch or body reflects the extent to which the splines are local.
3.7.2
The choice of the solution Bsplines
In Chapter 2, we mentioned that the method is not necessarily isoparametric and
-with the restriction (3.4). The geometry and solution B-splines are different if their orders and/or knot vectors differ. In practice it is convenient to choose the solution B-spline knot locations in the usable space the same as the geometry B-spline knot locations in the usable space. In assigning the geometry B-spline knot locations to the solution B-spline knot locations, the multiple geometry knots are considered a single knot.
The reason for this choice is that the self-influence coefficient algorithm(described in Chapter 5, p. 71) requires a polynomial form for the surface representation, Equa-tion (3.1). A B-spline has different polynomial representaEqua-tions on the space between consecutive non-coincident knots of its support. Therefore Equation (3.1) can be converted to a polynomial form between the finite intervals between knots. If the solution B-spline knot locations are arbitrary, on imposing Equation (3.8), separate domains over which the surface has different polynomial forms must be identified. From an implementational standpoint, the bookkeeping that ensues is non-trivial.
A sequence of increasingly accurate approximations can be obtained in the fol-lowing fashion. For the first approximation, choose the solution B-spline knots as described above (we call these knots the basic set). For the next approximation, new
knots are inserted by equipartitioning the space between knots of the basic set. For subsequent approximations, the above procedure is repeated and additional knots are inserted. By inserting new knots in the old knot vector, this procedure raises the in-tegers, N, of the B-spline expansion for the solution (Equation 3.3) and thereby provides a better approximation (see Section 3.3).
Finally, we note that the solution B-spline knots not on the usable space can be assigned arbitrarily, and the order of the solution B-splines can be selected at will.
3.7.3
Integrating the
271-cpterm exactly
We note that on expanding Equation (3.6), the product of the B-spline and the 27.0 term is a polynomial of order 2k in U. and v each. Exact integration of this product, at least requires a Gauss-Legendre rule of order, k.
an-alytically. Due to the illconditioning of the global matrix for large k (see Section 3.4.1, p. 36), it is important to minimize possible sources of error in the evaluation of the matrix elements. If the order of the Gauss rule is to be kept low (to keep the computational effort small), the exact integration alternative is essential for large k.
3.8
The computational effort
3.8.1 A theoretical estimate
for the Galerkin assembly
effort
To estimate the effort involved in the Galerkin assembly it is convenient to use as a starting point the overdetermined rectangular system of equations as obtained by a collocation procedure (see Section 3.4, p. 34). To keep the analysis simple, we consider a body described by
one patch and discretized by M x N panels.
For this discretization and a solution spline of order k, the number of unknowns are (Al + k1)(N + k 1). If
x .Vg collocation points are selected on each panel, theresulting rectangular system has MNA192 rows and (M+k 1)(N+ k 1) columns. We view the Galerkin assembly as a procedure where the M N.AI7 rows of the rectangular matrix are approprately weighted and "converted" into the (final) square matrix. In practice, the assembly procedure differs since the entire rectangular matrix is never actually constructed. But the effort (total number of multiplications and additions) involved is the same.
Consider the rows of the rectangular matrix to be composed of MN blocks each of length .Ar:. The rows of each block are a result of the .AIg2 equations from the collocation points from a panel. The "conversion" of the rectangular matrix into a square one is obtained by the following sequence:
(a) Consider the first block of .11/7 rows corresponding to the .Ar: collocation pointson
a panel. Associated with every collocation point are k` non-zero B-splines. The outer integration requires multiplying these .A/7 rows with one of the ig non-zero B-spline, a pre-computed weight (of the product Gauss rule) and adding them columnwise to obtain a single row. This row is then added on to a row of the (final) square matrix. Since each row is of length (M + k 1)(N + k 1), the effort involved is
(N92) (M + k 1)(N + k 1) additions-multiplications
The procedure is repeated for each of the k2 non-zero B-splines, and thus the effort
-is
k2
(.V9 + 1) (M + k -1)(N + k - 1)
additions-multiplications(b) The procedure in (a) is repeated for each of the MN blocks to complete the entire Galerkin assembly, and the total assembly effort is
k2 N (M + k
-
1)(N + k - 1) additions-multiplications. (3.12) Therefore, the Galerkin assembly effort is directly proportional to the square of the order of the solution spline (k) and the square of the order of the outer Gauss rule(Ac). For M,N >> k, the effort is also directly proportional to the square of the
number of panels (.Ar,), where Ar, = M x N.Note that it
is inappropriate to estimate the effort of the "conversion" of the rectangular matrix into a square matrix by pre-multiplication with a (M+k-1)(N k -1) x MN.A1,2, matrix. For M,N >> k, such a pre-multiplier is largely sparse due to the local support of the B-splines and the effort would be over-estimated.3.8.2
Relative estimates of the influence coefficient
evalua-tion effort
In the following, we compare the total number of influence coefficient evaluations by constant panel methods and the present method. To keep the analysis simple, it is convenient to assume that the surface consists of only one patch. For the more general case, order of magnitude estimates outlined below will be similar. We also note that the Rankine influence coefficient evaluations due to the reflection of the field point above the free surface, below the horizontal bottom and across the planes of symmetry (see Chapter 4) only change the constant of proportionality of the estimates below. This constant of proportionality for constant panel methods and the present method is the same, and for relative estimates no distinction with respect to the evaluation of the reflected field points is required.
is jV, where Arc is the total number of panels on the surface. In the present method, the total number of panels are Arh = M x N (see Section 3.7.1 for definition of a panel). If Al9 x .1\1; is the order of the product Gauss rule (for the outer integration), then the total number of field points are .1V7..Vh ( see Equation (3.3) and Section 3.6). Therefore, the total number of influence coefficient evaluations in the present method is .A/7../1/7. The ratio (R) of the total influence coefficient evaluations by the two methods is
R Arc2
Ar2d1Ph
To obtain the same relative accuracies, 1% 0.1%, we find that the ratio of number of panels required in the two methods' is, Ard.Arh :sz 10 60. Therefore we estimate, R (100 3600)/A/7. If the outer Gauss rule is of order. .A19 = 3, we have, R
10
400. For N
2, the ratio is R
25 900. Therefore, the total number ofinfluence coefficient evaluations are at least an order of magnitude less in the present method.
To estimate the ratio of the total time spent by the two methods in evaluating the influence coefficients is difficult, since the two methods utilize different algorithms for the Rankine influence coefficient evaluations.
For the evaluation of wave part of the influence coefficients, most of the total execution time is spent in the functional evaluation of the wave part, H. Therefore, the total time spent in the evaluation of the wave influence coefficients is directly proportional to the total number of functional evaluations of H. Typically, for
con-stant panel methods a 1 x 1
product Gauss rule is used for evaluating the wave influence coefficient. Therefore, the total number of functional evaluations are Arz.In the present method, aQxQ product Gauss rule is used and the total number
of functional evaluations are Q2.A/92ArZ. The ratio (1?,) of total number of functional evaluations of H isR R
Q 2
'See Table 3.1, p. 61. There we list the number of unknowns, (M + k 1)(N +k 1), hut since
M x N < (M + k 1)(N + k 1) these estimates hold.
=
For .Al = 3 and Q
= 4, It,
0.6 25. If= 2 and Q = 2, then R,
1.7 56.
If the present method and the constant panel method use the same implementation of the algorithms for the functional evaluation of H, the ratio of the total time spent in evaluating the wave influence coefficients can be estimated by the ratio of total number of functional evaluations of H.
