Journal of Ship Research, Vol. 32, No. 1, March 1988, pp. 4 4 - 5 3
J. J.
¥L
Baar^ and W. G. Price^
This paper discusses the numerical evaluation of the characteristic Kelvin wavelike disturbance trailing downstream from a translating submerged source. Mathematically the function describing ttie wavelike disturbance is expressed as a single integral with Infinite integration limits and a rapidly oscillatory integrand. Numerical Integration of such integrals is both cumbersome and time-consuming. Attention is therefore focused on two complementary Neumann-series expansions which were originally derived by Bessho [ 1 ] . ^ Numerically stable algorithms are presented for the accurate and efficient evaluation of the two series representations. When used In combination with the Chebyshev expansions for the nonoscillatory near-field component which were recently obtained by Newman [ 2 ] , the present algorithms provide an effective solution to the numerical difficulties associated with the evaluation of the Kelvin wave source potential.
Introduction
I N T H E L I N E A R I Z E D theory of ship waves and wave resistance as discussed by Andrew et al [3], Baar [4] and Baar and Price [5], a solution is sought to the steady disturbance potential of a ship moving at constant speed i n calm water. The Green's function associated w i t h the three-dimensional boundary-value problem for the disturbance potential (that is, the " N e u m a n n - K e l v i n " problem) is known as the Kelvin wave source potential and physically i t represents the potential of a translating submerged source. Solutions to the Neumann-Kelvin problem may be ob-tained by means of a source distribution technique involving the distribution of Kelvin wave sources over the ship's wetted hull surface and waterline contour. Such techniques require fre-quent evaluation of the Kelvin wave source potential, and the development of accurate and e f f i c i e n t algorithms for the evalua-tion of this Green's funcevalua-tion is necessary and of great importance i n any practical apphcations.
As discussed by Noblesse [6] and Baar [4] the Kelvin wave source potential can be decomposed into three components: (i) the potential of a fundamental Rankine source; (ii) the potential of a nonoscillatory near-field disturbance, symmetric upstream and downstream f r o m the source; and (iii) the potential of a Kelvin wavelike disturbance trailing downstream f r o m the mov-ing source (and zero upstream). E f f e c t i v e algorithms f o r the near-field disturbance have recently been derived by Newman [2], and i n the present paper attention is therefore focused on the calculation of the wave ike disturbance. Mathematically the wavelike disturbance is expressed as a single integral w i t h a rapidly oscillatory integrand. The numerical integration of such integrals is cumbersome, time-consuming, and prone to round-ing errors (see Davis and Rabinowitz [7]), and whenever possible should be avoided.
Following a summary of the most important features of the K e l v i n wave source potential, the wavelike disturbance is de-f i n e d i n equation (5) ode-f the paper and subsequendy de-f u r t h e r investigated. The proposed algorithms f o r the evaluation of the wavelike component are based on two complementary Neu-mann series representations (that is, expansions i n terms of Bessel
1 Brunei University of West L o n d o n , Department of Mechanical En-gineering, Uxbridge, Middlesex; presendy w i t h Shell International Pe-t r o l e u m Co., The Hague, Pe-the NePe-therlands.
2 Numbers i n brackets designate References at end of paper. Manuscript received at S N A M E headquarters August 1, 1986; revised manuscript received A p r i l 6, 1987.
functions), w h i c h were originally derived i n a remarkable paper by Bessho [1] (see also Appendix 2). The termwise summation of the two series expansions, given by equations (11) and (12), is prone to numerical instability due to the presence of Bessel functions of large integer order. This d i f f i c u l t y can be over-come by computing ratios of successive Bessel functions rather than b y such functions themselves, i n a manner originally sug-gested by Gautschi [8] and explained i n Appendix 3. Simple recursion relations are derived for the accurate evaluation of Bessho's series representations. The proposed algorithms are well suited to vectorization and enable the wavelike disturbance to be computed i n a very efficient manner.
It appears that the present algorithm can be used almost throughout the entire physical domain except i n the immediate vicinity of the free surface, where the asymptotic series given by equation (12) converges to the wrong answer, as was pointed out by N e w m a n [9], This problem remains yet to be resolved and is related to the singular behavior of the wavelike disturbance along the source track i n the free surface. I n Appendix 1 i t is discussed how this singular behavior can be related to Dawson's integral f u n c t i o n as defined b y Abramowitz and Stegun [10]. The convergence of Bessho's series representations can probably be somewhat iinproved by subtracting similar series expansions of the singular component defined i n equation (22) of Appendix 1 (see also N e w m a n [11]), but this approach is not pursued i n the present paper.
The Kelvin wave source
I t is convenient to introduce nondimensional variables i n terms of the acceleration of gravity g and the constant source speed V. For example, the nondimensional length £ is defined as = gL/V^ = F~ , where L is dimensional and the Froude number Fn = 'V/\lgL. These nondimensional variables are used exclusively i n this paper and a discussion of their advantages over other possible sets of nondimensional variables is given b y Baar and Price [5].
Figure 1 illustrates the Cartesian reference f r a m e Oxyz i n steady rectilinear motion at unit speed i n the positive 0.x-direc-tion. The ideal f l u i d occupies the lower half-space z < 0 and irrotational f l o w is assumed. The Kelvin wave source potential is defined as the Green's function G(x,a) associated w i t h the so-called Neumann-Kelvin problem formulated i n the linearized theory of ship waves and wave resistance [3-5], Physically, the f u n c t i o n G(x,a) represents the potential at the f i e l d point x{x,tj,z
< 0) where the f l o w is observed by a moving observer, due to a translating singularity of unit strength at the source point a{a,b,c < 0). As illustrated i n Fig. 1 the position vectors x and a relate to the moving axis system Oxyz. Specifically the unit outflow stems f r o m either a submerged source, i f c < 0, or a flux across the mean free surface z = 0, i f c = 0, as explained by Ursell [12] and Noblesse [6].
Also shown i n Fig. 1 is the vector X defined by
X = (X,Y >0,Z>0) = {a-x,\b- y\,\c + z\) ( I ) joining the field point x to the free surface mirror image of the source point a (notice that X- < 0 and X > 0 correspond to the upstream and downstream f l o w regimes, respectively). Fur-thermore, the quantities r and R are defined b y r = |x - a| and R = | X | , respectively, and physically r represents the distance be-tween x and a, whereas R is the distance bebe-tween x and the free surface mirror image of a.
A comprehensive survey of alternative expressions f o r the Green's function G(x,a) is beyond the scope of the present paper and the reader may consult Wehausen and Laitone [13], Eggers et al [14], Noblesse [6], Euvrard [15], and Baar [4]. Nowadays i t is generally recognized that an expression originally due to Peters [16] and later m o d i f i e d by Noblesse [6,17] is the most convenient formulation f r o m both physical, mathematical, and numerical points of view (see also Baar [4]). Thus the Kelvin wave source potential is expressed i n the f o r m
47rG(x,a) = - 1 / r + A'(X) + W ( X ) (2) w i t h components:
(i) the potential —1/r of a fundamental Rankine source i n unbounded f l u i d (as i f there were no free surface); (ii) the potential N ( X ) of a nonoscillatory near-field (local)
disturbance, symmetric upstream and downstream f r o m the source; and
(iii) the potential W ( X ) of a Kelvin wavelike (far-field) dis-turbance trailing downstream f r o m the source (and zero upstream).
As can be seen, the near-field and wavelike components i n equa-tion (2) are funcequa-tions only of the dimensionless vector quantity X as defined i n equation (1), and these terms describe the effects associated w i t h the presence of the free surface.
The near-field disturbance N{X) i n equation (2) is given by N ( X ) = 1/R + (2/7r) I m { e x p ( A ) £ i ( A ) } dt (3)
where the complex argument
A = \-Z^Jl - t^ + Yt + i j x I l V l - t^
I m denotes the imaginary part, and £ i represents the complex-valued exponential integral f u n c t i o n as defined b y Abramowitz and Stegun [10]. A simple picture of the behavior of this f u n c -tion can be obtained by expressing it i n the f o r m N ( X ) = Q(X)/R, where according to Noblesse [17] the f u n c t i o n (?(X) behaves like:
g ( x ) r-Hl + 0 ( R ) fR
1-1 + 0 ( 1 / R ) [R 0
(4) and i t is seen that the near-field disturbance can be interpreted as the potential of a Rankine image sink of strength Q{X).
Numerical integration of the integral i n equation (3) or alter-native integral representations [4] is cumbersome and time-con-suming. Recently, however, N e w m a n [2] has derived four com-plementary Chebyshev expansions trivariate i n X , Y, and Z. Using these expansions the near-field disturbance can be evaluat-ed i n an accurate and efficient manner. Figures 2(a) and 2(b) show the f u n c t i o n Q{X) = RJV(X) for Y = 0 and Z = 0, respective-ly, and clearly illustrate the typical l i m i t i n g behavior of this f u n c t i o n f o r small and large values of R as indicated by equation (4). Because of the effectiveness of Newman's algorithms for the near-field disturbance, the remainder of this paper is devoted exclusively to the investigation of the wavelike disturbance.
The wavelike disturbance
The wavelike disturbance W ( X ) i n equation (2) takes the f o r m W ( X ) = - u ( X ) 4 f sin{(X + Yt)^JÏ+7] e x p { - Z ( l -f-1^)] dt
where the Heaviside unit step f u n c t i o n
f x < 0 (i.e., upstream f r o m source) «(X)
X > 0 (i.e., downstream f r o m source) (5)
(6) This formulation clearly illustrates that the waves are f o l l o w i n g the source. The fact that W represents the potential of a Kelvin wavelike disturbance trailing downstream f r o m the source may be verified by investigating the asymptotic behavior of W f o r large values of R. This classic investigation is reported i n detail by Ursell [12], Wehausen and Laitone [13], and E u v r a r d [15],
Fig. 3 Sl<etcti of tfie Kelvin shiip wave pattem
who show that the f u n c t i o n W behaves as
W ( X )
'X < 2 / 2 Y
X 2 / 2 Y (7)
[X > 2 / 2 Y
W ( X ) = u{X)8 e x p ( - % Z ) 2 ' (-1)7;„(X)K„(V2D) cos(n^)
This expression impUes that the dominant Kelvin waves generat-ed by the moving source are confingenerat-ed to a sector making an angle of arccot(2y'2) c^ 1 9 ° 2 8 ' w i t h the downstream sailing line of the source, as illustrated i n F i g . 3. According to equation (7) the energy f l u x is least near the borderlines of this sector. A more detailed analysis shows that w i t h i n the sector there occur t w o distinct wave systems of transverse and diverging wave patterns w i t h a common angle of arctan(-v/2) ^ 5 4 ° 44' near the cusp lines.
I n addition i t can be shown that the f u n c t i o n W satisfies the following properties: W(X,0,0) = M(X)47rYi(X) W(0,Y,Z) = 0 IWI < 4 e x p ( - Z ) V V Z (8) (9) (10) where Yi denotes the usual Bessel f u n c t i o n as defined by Abramowitz and Stegun [10].
A point of some academic interest is the singular behavior of the wavelike disturbance f o r Z = 0. This relates to the "source" situated i n the free surface and, as indicated b y equation (8), the wavelike disturbance becomes singular along the source track Y = 0 = Z; that is, W(X,0,0) =^ - i i ( X ) 8 / X + O ( X l n X ) as X ^ 0. The singularity can be isolated by decomposing the f u n c t i o n W into singular and regular components. This procedure has been applied by Noblesse [18] and N e w m a n [11] to the centerplane potential W(X,0,Z) and is generalized i n Appendix 1. I t appears that the singularity can be related to Dawson's integral [see Appendix 1 equations (23) and (24)]. However, this approach is not f u r t h e r exploited i n the present paper.
The numerical integration of equation (5) or alternative inte-gral representations [4] is again troublesome due to the presence of the rapidly oscillatory integrand as illustrated i n Fig. 4. H o w -ever, Bessho [1] has derived t w o complementary series represen-tations by expanding terms i n the integrand of equation (5) into Neumann series and evaluating the remaining integrals analyti-cally (see Watson [19] for a comprehensive treatise of Neumann series). I t thus follows that (see also Appendix 2)
n=0
(11) M(X)87r exp(-V2Z) ^ ' Y ; „ ( X ) 7 „ ( V 2 D ) cos(n(3) (12)
n=0
I n these expressions the p r i m e d summation sign indicates that the first t e r m i n the series must be halved; / „ , Yn, / „ , and K„ are the usual Bessel functions of integer order n as defined by Abramowitz and Stegun [10]; and a p r i m e on the Bessel functions J2n and Yin denotes derivatives w i t h respect to X . Furthermore the polar coordinates {D>0,0<P^ ^kit) are defined by Z - j - lY = D exp(i^); that is, D = V Y ^ T Z ^ and /3 = arctan(Y/Z).
The two Neumann expansions given by equations (11) and (12) are complementary and w i l l henceforth be referred to as the /K-series and the Y/-series, respectively. The /K-series given by equation (11) is convergent and well suited for small and moder-ate values of X^/4D, as has indeed been c o n f i r m e d numerically by N e w m a n [11] for /3 = 0 (Y = 0), whereas the Yl-series given by equation (12) is i n fact asymptotic and therefore useful f o r large values of X^iD. However, N e w m a n [9] has pointed out that the Y/-series, as well as a related integral representation derived by Bessho [1], appears to be too regular to be u n i f o r m l y valid i n the immediate vicinity of the Z = 0 plane, where the Y/-series converges to the wrong answer (see also Appendix 2). I t is noted that the Y/-series is correct on the source track Y = 0 = Z [compare w i t h equation (8)] and i t appears also vaHd i f the source track is approached along an inclined radius beneath the free surface. I t is therefore recommended that the algorithms dis-cussed below be used only when 0 < jS < YgTr — e, where e is a small positive quantity.
Computational aspects
Due care is required when evaluating Neumann series of the type expressed i n equations (11) and (12). The Bessel functions satisfy three-term recurrence relationships w i t h respect to their degree n as discussed i n Appendix 3, equations (31) and (32). Only the functions Y„ and Kn are numerically increasing f u n c -tions of n and these func-tions can be evaluated without d i f f i c u l t y by means of f o r w a r d recursion. The functions J „ and / „ are numerically decreasing functions of n and are so-called m i n i m a l solutions to their respective three-term recurrence relationships.
- 0 . 1
Fig. 4 integrand of equation (5) for (X.V.Z) = (21,12,1)
as explained for example by Gautschi [8]; recursion is therefore only stable i f applied i n the backward direction. Moreover, it appears that the functions Y„ and K„ increase exponentially i n magnitude w i t h increasing n , whereas the functions ] „ and /„ vanish i n an exponential manner. This feature causes serious cancellation errors and under/overflow problems during the running summation of the Neumann series. The difficulties associated w i t h the calculation of the Bessel functions can be avoided by computing ratios of successive Bessel functions rather than the functions themselves, as suggested b y Gautschi [8] and explained i n Appendix 3.
Baar [4] has shown that effective and numerically stable algo-rithms are obtained by expressing the Neumann series represen-tations given by equations (11) and (12) i n the f o r m :
W ( X ) = M(X)Cj,(/3) + e^iX) (13) where eN(X) denotes the truncation error and Cf^(/3) is the f i n i t e
Fourier cosine series given by
(14)
11=0
Series of this type may be evaluated rapidly by means of the Goertzel-Clenshaw algorithm (see Nonweiler [20]), w h i c h can be expressed as CN(|S) = yziho — bz), where bo and ^2 are obtained recursively f r o m
= fln + 2 C0s(^)fcn+i - i„+2 (15) for n = N,N — 1 , . . . ,0. Here cos(;ö) =Z/D and the initial values are given by btq+2 = 0 = bi^+i.
The Fourier coefficients a„ (n = 0,1, . . . ,IV) corresponding to the two series representations given by equations (11) and (12) can be calculated i n an efficient manner by means of the follow-ing two-step procedure;
1. Using the algorithms outlined i n Appendix 3 the following sequences of Bessel function ratios are computed:
(16)
( 1 7 )
m = / 2 „ + i ( ^ ) / / 2 „ - i ( ^ )
U\D) = K„+,{%D)/K„{%D)
for n = 0 , 1 , . , . , N i n the /K-series given by equation (11); and y„(X) = Y2„+i(X)/Y2„_i(X)
i„{%D)=I„^,i%D)/lM
for n = 0 , 1 , . . . ,1V i n the YTseries given by equation (12). N o -tice thatyo(X) = - 1 = yoiX).
2. By making use of equations (11) and (12), as well as the derivative relations f o r the Bessel functions [see Appendix 3, equations (38)-(40)], it can be shown that the Fourier coefficients
a„ i n equation (14) may be computed recursively by means of
(18) for n = 0 , 1 , . . . ,IV — 1 i n the /K-series, the initial value being ÜQ = - 8 exp(-y2Z)7i(X)Ko(y2D); and
fln+i = a„i„y„{l
-
?/„-H)/(1 - y^) (19) f o r n = 0,1, . . . , N — 1 i n the Y/-series, the initial value being öo = 87rexp(-y2Z)Yi(X)Io(y2D).Extensive numerical experiments have been p e r f o r m e d to assess the range and capability of the present algorithms, as well as to establish the convergence properties of the two series expan-sions. I n these experiments no attempt has been made to evalu-ate the wavelike disturbance at the free surface (that is, only stricdy positive values of Z have been considered). For selected values of (X,D,^) i n the ranges 0 < X < 50, 0 < D < 20, and 0 < jS < 0.487r 86.4 deg a table of accurate benchmark values of W
has been prepared by numerical integration of equation (5). I t appears that the two series expansions are indeed complemen-tary and absolute accuracy e = 10"^ (giving between five and six significant digits) can be achieved throughout the previously indicated domain.
The JK-series given b y equation (11) is best suited for small and moderate values of X^/AD {X^/AD less than about 20) w i t h no more than about 70 terms required, whereas the Y/-series given by equation (12) may be used for large values of X^/AD, w i t h no more than about 20 terms required. The precise transi-tion value of X^/AD between the two regimes of convergence depends i n a rather complicated fashion on both the coordinate X, the distance D f r o m the source track, the angular orientation /3, and the specified tolerance e (as expected because of the asymptotic nature of the Y/-series). The required number of terms i n the JK- and YZ-series appears to be roughly proportional to X V D and D / X , respectively.
The following device to establish which series should be used has proved effective i n computational practice at very litde additional cost. Initially the ratios ;„(X) and kni^iD) as defined by equation (16) are computed f o r n = 0 , 1 , . . . ,A', where IV is chosen sufficiently large ( N = 75, say). The Fourier coefficient a„ corresponding to the /K-series as defined by equation (18) may now be calculated and convergence of the /K-series may easily be monitored by observing the rate of decay i n the magni-tude of the Fourier coefficients. Specifically, the summation of the /K-series may be terminated when \at} < tr, where er represents the desired accuracy. When the /K-series fails to converge, the Y/-series is used and evaluated by making use of equation (17) (with N = 25, say) and equation (19). Conver-gence of the Y/-series can be established by comparing successive convergents of equation (14) u n t i l the desired accuracy has been achieved.
Similar recursion schemes can be derived f o r the calculation of the derivatives of the wavelike disturbance by differentiation of the series expansions given by equations (11) and (12). Using the present algorithms the computing time is reduced by at least a factor 100 i n comparison w i t h direct numerical integration of equation (5). Moreover the presented recursion schemes are suitable for vectorization on a parallel computer and the average CPU-time is approximately 50 microseconds on a Cray-lS m a - , chine f o r one evaluation of W and its gradient. Figures 5a-5d show the computed wavelike disturbance and gradient compo-nents f o r a given value of Z = 1.
Conclusion
The evaluation of the K e l v i n wave source potential is of great importance i n source distribution programs for the analysis of steady ship generated wave patterns and related phenomena. The near-field and wavelike disturbances expressed as compo-nents w i t h i n this Kelvin potential can be computed by means of Chebyshev and Neumann series representations, respectively. Accurate and efficient algorithms f o r the evaluation of the latter expansions have been proposed and it has been shown that an important obstacle i n the linearized theory of ship waves and wave resistance can thus be removed. However, certain ques-tions regarding the behavior of the wavelike disturbance i n the immediate vicinity of the free surface require f u r t h e r investiga-tion.
Acknowledgment
We gratefully acknowledge the encouragement of D r . R. K. Burcher [chief superintendent. A d m i r a l t y Research Establish-ment (ARE), Haslar (U.K.)] and Dr. R. N . Andrew (ARE) during the development of this work, w h i c h has been carried out w i t h
the support of the Procurement Executive Ministry of Defence We are also grateful to Professor J. N . N e w m a n (Department of Ocean Engineering, Massachusetts Institute of Technology) for providing the Chebyshev expansions of the near-field distur-bance component.
References
1 Bessho, M . , " O n the Fundamental Function i n the Theory of the Wavemaking Resistance of Ships," Memoirs. Defense Academy (Ja-pan), V o l . 4, 1964, pp. 99-119.
2 N e w m a n , J. N . , "Evaluation of the Wave-Resistance Green Function: Part 1—The Double I n t e g r a l , " JOURNAL O F SHIP RE-SEARCH, Vol. 3 1 , No. 2, June 1987, pp. 79-90.
3 Andrew, R. N . , Baar, J. J. M . , and Price, W . G., "Prediction of Ship Wavemaking Resistance and Other Steady F l o w Parameters Using N e u m a n n - K e l v i n Theory," Trans. R I N A , V o l . 130, 1988.
4 Baar, J. J. M . , " A Three-Dimensional Linear Analysis of Steady Ship M o t i o n i n Deep Water," Ph.D. Thesis, Brunei University of West London, Uxbridge, U.K., 1986.
5 Baar, J. J. M . and Price, W . G., "Developments i n the Calculation of the W a v e m a k i n g Resistance of Ships," Proceedings of the Royal
Society of London, Series A, Vol. 416, 1988, pp. 115-147.
6 Noblesse, F., " A l t e r n a t i v e Integral Representations f o r the Green F u n c t i o n of the Theory of Ship Wave Resistance," Journal of Engineer¬ i n g Mat/iemaiics, V o l . 15, 1981, pp. 241-265.
7 Davis, P. J. and Rabinowitz, P., Methods of Numerical Integra-tion, 2 n d Ed., Academic Press, N e w York, 1975.
8 Gautschi, W . , "Gomputational Aspects of T h r e e - T e r m Recur-rence Relations," SIAM (Society for Industrial and Applied
Mathemat-ics) Review, Vol. 9, 1967, pp. 24-82.
9 N e w m a n , J. N . , Private communication, 1986.
10 Handbook of Mathematical Functions, M . Abramowitz and I . A. Stegun, Eds., Dover Publications, N e w York, 1972.
11 N e w m a n , J. N . , "Evaluation of the Wave-Resistance Green Function: Part 2—The Single Integral on the Centerplane," JOURNAL OF SHIP RESEARCH, Vol. 3 1 , No. 3, Sept. 1987, pp. 145-150.
12 Ursell, F., " O n Kelvin's Ship Wave Pattern," Journal of Fluid
Mechanics, Vol. 8, 1960, pp. 418-431.
13 Wehausen, J. V . and Laitone, E. V . , "Surface Waves" i n
Hand-iouch der Physik, Springer-Verlag, Berlin, Band 9, 1960, pp. 446-778.
14 Eggers, K. W . H . , Sharma, S. D . , and W a r d , L . W . , " A n Assess-ment of Some ExperiAssess-mental Methods f o r D e t e r m i n i n g the W a v e m a k i n g Characteristics of a Ship F o r m , " Trans. S N A M E , Vol. 75,1967, pp. 112¬ 144.
15 E u v r a r d , D . , "Les M i l l e et Une Faceties de la Fonction de Green d u Probleme de la Resistance de Vagues," Report No. 144, Ecole N a t i o n -ale Superieure de Techniques Avancees, Paris, 1983.
16 Peters, A. S., " A N e w Treatment of the Ship-Wave P r o b l e m , "
Communications on Pure and AppUed Mathematics, Vol. 2, 1949, pp.
123-148.
17 Noblesse, F., " T h e Fundamental Solution i n the Theory of Steady M o t i o n of a Ship," JOURNAL O F SHIP RESEARCH, V o l . 2 1 , No. 2, June 1977, pp. 82-88.
18 Noblesse, F., "The Steady Wave Potential of a U n i t Source, at the Centerplane," JOURNAL O F SHIP RESEARCH, V o l . 22, No. 2, June 1978, pp. 80-88.
19 Watson, G. N . , A Treatise on the Theory of Bessel Functions, 2nd ed., Cambridge University Press, U . K . , 1944.
20 Nonweiler, T. N . R., Computational Mathematics, Ellis Hor-wood Publications, Chichester, U . K . , 1984.
21 ' Gautschi, W . , " E f f i c i e n t Computation of the Complex E r r o r F u n c t i o n , " S J A M / o u r n a / o n Numerical Analysis, V o l . 7,1970, pp. 187¬ 198.
22 Tal)le of Integrals, Series and Products, I . S. Gradshteyn and I . M . Ryzhik, Eds., 4th ed.. Academic Press, N e w York, 1980.
23 Maruo, H . , " N e w Approach to the Theory of Slender Ships w i t h F o r w a r d Velocity," Bulletin of the Faculty of Engineering, Yokohama National University, Japan, V o l . 3 1 , 1982, pp. 85-100.
24 Yeung, R. W . and K i m , S. H . , " A N e w Development i n the Theory of Oscillating and Translating Slender Ships" i n Proceedings, 15th Symposium on Naval Hydrodynamics, H a m b u r g , 1984, pp. 195¬ 212.
25 ^ N e w m a n , J. N . , "Approximations f o r the Bessel and Struve
Fvmc-Wons," Mathematics of Computation, Vol. 43, 1984, pp. 551-556.
26 Allen, E. E., " P o l y n o m i a l Approximations to Some M o d i f i e d Bessel Functions," Mathematical Tables and Other Aids to
Computa-tion, Vol. 10, 1956, pp. 162-164.
Fig. 5 ( a ) The function IVCX.V.I)
Fig. 5 ( b ) The function WJ,X, V,1)
W ( X ) = - u ( X ) 8 e x p ( - Z ) f W s ( X ) + W „ ( X ) 1 (21) where the functions W s and W u are defined b y
W s ( X ) = [ sin(Xi) cos(Y<2) exgl^-Zt^) dt Jo
WB(X) = . f " {sin(xVrT?) cos(YiVr+7)
.0
- sin(X<) cos(yi')) e x p ( - Z < ' ) dt (22) For large values of the integradon variable t the integrand of this .r. j , , i r ,1/ • i " /ti. f expression behaves as sin(Xi)cos(yi2)exp(-Z<2). This behavior suggests respectively I t may be verified that the f u n c t i o n W R is regular (that decomposing W ( X ) into components is, bounded but not analytic).
Appendix 1
Investigation of the singular behavior
By expanding the trigonometric t e r m i n equation (5) the expression f o r the wavelike disturbance potential can be r e w r i t t e n asW ( X ) = - u ( X) 8 [ °
sin(xVrT7) cos(Yi>/rT?)
e x p { - Z ( l + t ' ) ) dt (20)Fig. 5 ( c ) The function 1%(X,V,1)
Fig. 5{d) The function WJ,X,Y^)
T h e singular f u n c t i o n W s ( X ) can be related to the complex-valued Dawson's integral f u n c t i o n F as d e f i n e d by A b r a m o w i t z and Stegun [10], equation 7.4.7. Specifically
WgOQ = (^M/X (23)
where a)((T) = 2 RelVo-F(v'c)] and the complex similarity variable
<r = XV4(Z - l Y ) = (XV4D) exp(i/3) (24) Here polar coordinates (D,/3) d e f i n e d b y Z -I- j Y = D exp(i/3) have been introduced. A n e f f i c i e n t a l g o r i t h m f o r Dawson's integral has been proposed b y Gautschi [21] and f o r the case where Z = 0 the f u n c t i o n co(o-) can be expressed i n terms of real-valued Fresnel integral functions (see Gradshteyn and Ryzhik [22], equation 3.691.6).
W h e n Y = 0 = Z, that is, on the source track i n the f r e e surface, i t may
be v e r i f i e d that the singular and regular components are given b y
W s ( X , 0 , 0 ) = 1 / X
W R ( X , 0 , 0 ) = - ( 7 r / 2 ) Y i ( X ) - 1 / X
(25)
respectively, so that the singularity is now entirely accounted f o r b y the singular f u n c t i o n W j . Physically the component Ws can be interpreted as the potential of a diverging wave system w h i c h is illustrated i n F i g . 6 for Z = 1.
A similar decomposition of the wavelike disturbance has been studied b y E u v r a r d [15] and, f o r Y = 0, b y Noblesse [18] and N e w m a n [11]. I t i s also noted that the f u n c t i o n Ws arises i n t w o recent new approaches to the theory of steady m o t i o n of slender ships as proposed b y M a r u o [23] and Yeung and K i m [24].
Fig. 6 The function W^^X, YA)
Appendix 2
Derivation of Bessho's Neumann series
By p e r f o r m i n g tlie change of integration variable t = sinh(V2s) i n equation (20) and i n t r o d u c i n g the polar coordinates (D,^) d e f i n e d b y Z+ iY = D exp(j/3), the wavelike disturbance can be w r i t t e n as
W ( X ) = u{X)8 e x p ( - V 2 Z ) Ö P ( X ) / d X (26) where the f u n c t i o n P(X) is d e f i n e d b y
P(X) = % Re f cosjX cosh(V2s)l expl-VaD cosh(s - i(3)l ds (27)
F r o m the generating f u n c t i o n and associated series of the Bessel f u n c -tions (see A b r a m o w i t z and Stegun [10], equation (9.1.41) the f o l l o w i n g N e u m a n n series may be obtained:
cos|X cosh(V2s)) = 2 ^ ' ( - l ) " ; 2 n ( ^ ) ™sh(ns) (28) n=0
where the p r i m e indicates that the first t e r m i n the series must be halved. Substitution of equation (28) into the integrand of equation (27) gives after some m a n i p u l a t i o n
P(X) = ^ ' {-irj,„{X)A„{%D) cos(ri^)
vhere
cosh(ns) expj-VaD cosh(s)j ds
(29)
(30)
Equation (30) is recognized as an integral representation of the m o d i f i e d Bessel f u n c t i o n K„{%D) (see A b r a m o w i t z and Stegun [10], equation 9.6.24) and the /K-series given b y equation (11) follows i m m e d i a t e l y upon substituting equations (29) and (30) into equation (26).
According to Bessho [1] the Y/-series given b y equation (12) can be derived i n a similar manner b y expanding the exponential t e r m i n the integrand of equation (27) into a N e u m a n n series i n v o l v i n g the m o d i f i e d Bessel functions / „ ( y 2 D ) (see A b r a m o w i t z and Stegun [10], equation
9.6.33). The resulting integrals can be related to the Bessel functions
Y2„(X). However, the present authors have not been able to produce a
correct derivation of Bessho's Y/-series representation. I t is noted that the f u n c t i o n W is singular i n the v i c i n i t y of the Z = 0 plane as Y approaches zero (see Appendix 1 and also Ursell [12] and E u v r a r d [15]) while the N e u m a n n series given b y equation (12) is regular. This fact seems to explain the convergence of equation ( l 2 ) to the w r o n g answer, as pointed out i n the discussion below equations (11) and (12).
Appendix 3
Evaluation of the Bessel functions
L e t J„{x), Ynix), I„(x), and K„(x) denote the usual Bessel functions of x and integer order n as d e f i n e d by A b r a m o w i t z and Stegun [10]. For n = 0,1 the Bessel functions ]„(x) and Y„(j;) can be computed b y means of the p o l y n o m i a l approximations developed b y N e w m a n [25]. Similar ap-proximations f o r the zeroth and first-order m o d i f i e d Bessel functions J„(x) and K„(x) have been derived by A l l e n [26] and these are quoted b y A b r a m o w i t z and Stegun [10]. ( F r o m a n u m e r i c a l point of view i t is better to compute exp(-x-)J'n(*) and exp{x)K„{x) rather than I„{x) and
K„{x), respectively, since this removes most of the variation i n these
functions.)
W h e n n > 2, the f o l l o w i n g procedures have been f o u n d to be very convenient. F o l l o w i n g Gautschi [8], consider the general three-term recurrence relationship
/ n + l + flJn + V „ - l = 0 (31) where ƒ„ = f„{x) and n > I. For the Bessel functions the coefficients
a„{n,x) and b„(n,x) are given by A b r a m o w i t z and Stegun [10] as
{-2n/x,\) ( 2 n A , - l ) l ( - 2 n A , - l ) 7„(*)=Y„(«) fn = • U^) KnM (32)
F r o m equation (31) i t follows that the ratio r„{x) = f„+i/f„ satisfies the relationship
0 (33) f o r n > 1, provided t l i a t / „ ( x ) 0 (this restriction is not serious f r o m the practical point of view, as explained by Gautschi [8]).
The functions Y„(x) and l<:„(x) are numerically increasing functions of n and the ratios r„ can simply be computed by means o f the f o r w a r d recurrence
(34) f o r n > 1, the i n i t i a l value being ro = / i / / o .
T h e functions / „ ( x ) and /„(x) are numerically decreasing functions of
n [that is, " m i n i m a l " solutions of equation (31)] and therefore the
recur-rence relationship given b y equation (33) is stable only when applied i n the backward direction, That is
^ - 1 (35)
f o r n = A ' , N - 1, . . . , 2 , 1 , where i t is assumed that the ratio )•„ is k n o w n f o r some value n = TV. Gautschi [8] shows that rw is given by the continued f r a c t i o n
= l i m (36)
where a„ and h„ are d e f i n e d i n equation (32). T h e successive conver-gents W]^ m a y be generated recursively as follows:
";t+i = - (fcw+i+i/aw+jiON+A-n)"itl
«i+i = ^kWw - 1) (37)
f o r k = 1,2, . . . , the i n i t i a l values being «1 = 1 and Vi = -bN+i/an+i =
Wl.
The procedures described i n this appendix are very w e l l suited f o r the computation of the Bessel f u n c t i o n ratios d e f i n e d i n equations (16) and (17) (a small but straightforward m o d i f i c a t i o n is required f o r the calcu-lation of the jn's and y„'s). I n order to v e r i f y the recalcu-lations expressed by equations (18) and (19), as well as to derive similar relationships f o r the evaluation of the gradient of the wavelike disturbance, the f o l l o w i n g derivative relations are useful:
c/„ = ƒ „ - ! + r f / „ + i where ƒ„ = df„{x)/dx and
' ( 2 - 1 ) (JnMXM) ( c , d ) = (2,1) i f / „ = ] / „ ( x )
[ ( - 2 . 1 ) [ K „ ( X )
see A b r a m o w i t z and Stegun [10], I n particular (/o,Y'o,/o.K'o) = ( - 7 i , - Y i , / i , - K i )
(38)
(39)
(40)
