Application of Chebyshev Polynomials to Efficient Approximation of Frequency-Domain Infinite Depth Green Function

-,-,-,-cn

(1.School of Naval Architecture,Ocean&Civil Engineering,Shanghai Jiao Tong University,Shanghai 200240,China;2.Marine Design and Research Institute of China,Shanghai 200011,China;

3.Shanghai Key Laboratory of Ship Engineering,Shanghai 200011,China)

Abstract: New research to fast evaluate the Green function and its derivatives without translating in frequency domain is presented. Concerning the classical representation, composed of a semi-infinite integral involving a Bessel function and a Cauchy singularity, the Green function is divided into two parts： relatively simple terms with analytic expressions and slowly-varying residual functions. Economized Chebyshev polynomials are derived in three modified regions to approximate the residual function,and subsequently the Clenshaw’s method is used to fast evaluate the resultant nested multiplications.To guarantee 6D accuracy,the evaluations of derivatives also are implemented concurrently with almost the same polynomial coefficients, more precisely, with several extra or reduced terms. Numerical comparisons show that agreement of computational results between this paper and the famous HydroStar is fairly well, and that the algorithm in this paper can virtually obtain the same efficiency as HydroStar.In conclusion,this method can yield results with high accuracy and efficiency to the practical application.

Key words:Green function;Chebyshev polynomials;potential theory;frequency domain

0 Introduction

Under the assumption of a‘high frequency-low speed’and rigid or hydroelastic body,the infinite depth Green function without translating has been widely applied to the seakeeping analysis and marine structure design in recent several decades. Because of the enlargement trend of ship’s scale in recent years and the onerous analysis task involving structural spectral analysis, a realistic description of wave-body interactions dictates that computations of the Green function between 5×107and 5×1011are required to analyze an individual ship, which is regarded as the main time-consuming part in performing the hydrodynamic analysis.

The evaluation of Green function and its derivatives have always been a complicated mathematical issue. Numerous practical works have been done on this subject. Noblesse[1]proposed the representation of Green function with the so-called near-field and far-field parts. As to a representation with a semi-infinite integral involving a Cauchy singularity, for the first time Newman[2]advocated the classical fast method combining analytic expansions with two-dimensional Chebyshev polynomial approximations,however,the variety of diverse algorithms and associated branches precludes effective use on every single subdomain.A more systematic numerical approach[3]was derived from approximation of the residual functions in all four rectangular regions with economized polynomials. This approach is the basis of the famous analysis code WAMIT.To some extent,Chen[4]utilized the philosophy of polynomial approximations[3]and the expression of residual functions from Noblesse[1]and developed a new efficient method, which was applied in the notable analysis code HydroStar. In addition, there are also some novel implementations of other representations.Clement[5]outlined the method using classical fourth-order Runge-Kutta method to solve a second-order ordinary differential equation of Green function. Similarly, in 2015, Shen et al[6]investigated the method combinging Numerov type method with Power-series method to solve this equation. New series expansions for different subdomains were proposed by Duan et al[7]. On the basis of that, Shan[8]further proposed highly precise approximation on refined subdomains.All in all,to the best of authors’knowledge,in view of computational efficiency and practical application,the above algorithms[3-4]might be the most effective scheme. However, there are still rooms for research on the selection of the residual function, on the refinement of subdomain division based on economized polynomials methods and on the solution to the resultant Chebyshev polynomial sum.

This paper presents the implementations of the efficient approximations of free-surface Green function without translating in frequency domain. Firstly, computations are divided into evaluations of relatively simple terms and slowly-varying residual functions, whose expressions are introduced in detail. After the approximation and error analysis of Green function, the efficiency and accuracy of the code are verified by comparing the results with those from Hydro-Star. It is concluded that the algorithms in this paper are accurate and efficient for the wavebody interactions.

1 Efficient algorithm of infinite depth free-surface Green function

We consider a field pointp(x,y,z),source pointq(ξ,η,ζ)and an image pointqˉ(ξ,η,-ζ)of the source point relative to the free surface. The following expression is the complex Green function of infinite depth water corresponding to the time exponential item exp(-iωe t),ωeis the encounter frequency.

F(X,Y)and its first derivatives are defined as follows：

wherer=[R2+(z-ζ)]1/2,r1=[R2+(z+ζ)]1/2,ν=ωe2g-1,X=νR,R=[(x-ξ)2+(y-η)]1/2is the horizontal distance between p andq,Y=ν|z+ζ|,R1=(X2+Y2)2,gis the gravitational acceleration. Jn(·) is then-th order Bessel function of the first type,Yn(·)is then-th order Bessel function of the second type,and Hn(·)is then-th order Struve function.

The evaluations of the Rankine sources 1/r+1/r1and the analytic part 2πνexp(-Y)J0(X)can be easily solved. ThusF(X,Y) can temporarily represent the Green function. The main task is reduced to the solution toFandFXin Eq.(2)throughout the quadrant.

Here we utilize the‘heavy-handed’application of residual function[3], which means that computations of Green function are divided into evaluation of relatively simple terms and slowly-varying residual functions. Unlike the division of Newman’s work with four rectangular regions,we re-derive the modified regions in view of the cyclical variation trend ofFwith continued amplitude attenuation, with one newly-combined region from the original Rregions 2 and 4,and the other two are identical to the original.Furthermore,the slowly-varying residual functionsf(X,Y)andfX(X,Y)ofFandFXrespectively and the rest analytic expression for three regions are introduced explicitly in Eqs.(3～5).

Region 1：0＜X≤3,0＜Y≤4

Region 2：0＜X≤3,Y≥4

Region 3：X≥3,Y＞0

In view ofXandY,it is necessary to define the maximum ofXandYto be 1 800 and 100 respectively if a ship of 400 meters in length is sailing in the heading angle,with the draft of 20 m and a velocity of 20 knots. Numerous tests of numerical computations indicate that special modified partitions atX=1,2,3,4,5,7,10,15,24,40,60,100,200,400,700,1 000 andY=1,2, 3, 4, 5, 7, 10, 15, 22, 35, 60 are explicitly applied. In the final 123 subdomains (non-uniform with respect to the partitions), economized doublepolynomials are adopted to evaluate the residual functionsfwith 6D accuracy.FandFXin Eqs.(3-5) can be evaluated by the Romberg integral. However, considering the efficiency of computing Chebyshev coefficients (to be introduced later),means developed by Shan[8]with 9D accuracy is adopted in this paper.It should be stressed that the inefficient series expansion[3]is also adopted in case the calculation pair(X,Y) is beyond the range of [0, 1 800]×[0, 100], even though it may not happen.

Fig.1 The division of the whole(X,Y)quadrant

Besides, considering the origin of (X,Y) quadrant is the singularity of this issue,so whenX→0+,Y→0+,FandFXcan be defined to be zero without influences on results. Furthermore, whenX→0+andY≠0,meaning the field point coincides with the source point only in horizontal direction, it is clear that the above Eqs.(3-5)cannot be adopted directly because of the singularity ofYn(X)andX-1.But there exist the following supplementary relations through simple manipulation[12].

Finally, a comprehensive algorithm covering the whole (X, Y) quadrant is derived for the evaluation of infinite depth case (see Fig.1), in other words, with 5 regions. Moreover, the approximations of all involved special functions J0(·),J1(·),Y0(·),Y1(·),Ei(·)are also implemented efficiently[9-11].The errorεin this paper to determine whetherXorYequals to 0 is 10-5in the final Fortran subroutine.

The method to expand a function into a series of Chebychev polynomials needs to be outlined briefly. Based on Chebyshev polynomialsTm(x) of the first kind which lie in the interval[-1,1],and the givenf(X,Y)defined on[a,b]×[c,d],we compute Chebyshev coefficientsai j

wherezk,zlare Chebyshev interpolation nodes on [-1,1],wklis the evaluation off(X,Y) at approximation nodes on[a,b]and[c,d]intervals,ai j,i=0,…,n1;j=0,…,n2(m≥n1+1,m≥n2+1).

For bothfandfX, the final achieved coefficient orders (n1,n2) vary from (2, 1) associated with 6 coefficients up to(6,4)with 35 coefficients.Moreover,the coefficients should have a precision of several more decimals than needed in the final approximation and be neglected if the magnitudes of them are smaller than 10-9.

At last,arriving at approximation off(X,Y)on[a,b]×[c,d]：

A first possibility for the computation of this nested-multiplication in Eq.(8) is to rewrite the Chebyshev polynomialsTn(x) in terms of powers ofxand then use the Hörner scheme[4].However,this has to be done carefully because for some expansions there is a considerable loss of accuracy due to cancellation effects[13]. In this paper, an alternative and efficient Clenshaw’s method for evaluating this nested-multiplication is adopted. Meanwhile, the derivativefXcan be evaluated as efficiently asf,only with one or two extra or reduced terms(depending on different subdomains)in order to guarantee the same level of accuracy.

2 Numerical results and discussion

Based on the above mentioned algorithms and the boundary value problem with the source formulation,codes are developed in Intel visual Fortran 2013 version.Furthermore,the symmetry of Green function[14]is always used in the code. However, the symmetry of body geometry and elimination of the irregular frequencies are not included. All of the numerical tests are run on a personal computer with a CPU clock of 3.6 GHz.

2.1 Approximations of residual function

Before we evaluate the residual function with economized polynomials, its variation trend need to be confirmed with straightforward method according to Eqs.(3-5). Fig.2 shows that no matter forf(X,Y) orfX(X,Y), Region 1 has more severe variation than the other two regions.And the newly combined Region 3 is rational because of its slow-varying characteristic. It can also be concluded that approximation for fast-varying Region 1 will need more coefficients than the other two.

Fig.2 The slowly-varying parts f and fX for ε＜X≤15 and ε＜Y≤8

When (X,Y) is located in [ε, 3]×[ε, 4], Fig.3 depicts the approximation error of residual functions with economized Chebyshev polynomials in Region 1 with 3×4=12 subdomains. It can be obviously concluded that accuracy of all the evaluating points with differentXandYcan reach at least 6D accuracy. Now taking the subdomain [ε, 1]×[ε, 1] of Region 1 as an example to explain how the polynomials coefficients are used in the approximation. The used coefficient is illustrated in Tab.1,with three kinds of numbers,which are the normal ones,three bold forn2=5,and one italic forn1=6(the null cell stands for the zero value of coefficients when the magnitudes of them are smaller than 10-9). The normal plus the three bold is adopted for the approximation off,whereas the normal plus the one italic is forfX.Although the precision loss is inherent in numerical differentiation, we observe that the required number of coefficients forfis actually two more than that forfXin order to guarantee the same precision of 6D. Certainly opposite case will happen in other subdomains. All in all, coefficients for approximations offandfXis derived from the same coefficient matrix like Tab.1, which is completely different from Wang[12].

Fig.3 Evaluation error of residual functions on[ε,3]×[ε,4]based on Chebyshev polynomials

Tab.1 Coefficients for the approximation of f and fX in subdomain ε≤X≤1 and ε≤Y≤1

2.2 Numerical comparison of a real ship with HydroStar

To demonstrate the above developments, numerical computations have been performed for a very large container vessel with length/width/draft=380/52.4/16.0 in meter. Three different discretization of M1, M2, M3 composed of 1 507, 3 015, 8 910 individual panels on the wetted hull (see Fig.4) are used for comparison. The sailing speed is assumed to be 5 knots. The famous subroutine rdf.exe of HydroStar Version 7.3 (2016) is adopted for comparison,which is also for radiation and diffraction computation.

Fig.4 Different discretization on the same wetted hull

In this section, the accuracy and efficiency of the codes are validated.The nondimensional added-mass of surgeA11, heaveA33, and pitchA55and damping coefficients of surgeB11, heaveB33, and pitchB55are illustrated in the Fig.5. It is observed that for both small and large values of wave frequency,agreement between these two methods is fairly well.

Fig.5 Accuracy comparisons of hydrodynamic coefficients with 29 wave frequencies from 0.1 to 1.5 rad/s and one heading angle 180°(ρ is density of sea water,V is displacement,L is ship length)

From consumed times of 30 regular waves as presented in Tab.2,we can find that the efficiency in this paper is almost equal to that of HydroStar, with 2.5 percent higher efficiency for M1 and about 3 percent lower efficiency for both M2 and M3. For instance, the consumed time of one single regular wave for M2 is only 121/30≈4.03 seconds.

Tab.2 Efficiency comparison within 30 regular waves composed of 3 heading angles(0°, 90°, 180°) and 10 wave frequencies (0.1～1.9 rad/s with step 0.2)

3 Summary

In this paper, new research of efficient algorithm for free-surface Green function in frequency domain is presented in detail. This algorithm can yield results with adequate accuracy to the practical application,and can achieve almost the same efficiency as the famous hydrodynamic code HydroStar. Furthermore, considering the increasing number of cores in CPU as the relevant technology advances, another effective solution to reduce computational time is to adopt high-performance computers alongside code parallelization, which will be the objective for our next phase of work.