HIGHLY PRECISE APPROXIMATION OF FREE SURFACE GREEN FUNCTION AND ITS HIGH ORDER DERIVATIVES BASED ON REFINED SUBDOMAINS

The infinite depth free surface Green function (GF) and its high order derivatives for diffraction and radiation of water waves are considered. Especially second order derivatives are essential requirements in high-order panel method. In this paper, concerning the classical representation, composed of a semi-infinite integral involving a Bessel function and a Cauchy singularity, not only the GF and its first order derivatives but also second order derivatives are derived from four kinds of analytical series expansion and refined division of whole calculation domain. The approximations of special functions, particularly the hypergeometric function and the algorithmic applicability with different subdomains are implemented. As a result, the computation accuracy can reach 10 in whole domain compared with conventional methods based on direct numerical integration. Furthermore, numerical efficiency is almost equivalent to that with the classical method.


Introduction
The wave loads on fixed structures and the oscillatory motions of vessels free to response to the waves are common problems to be solved in ocean engineering.For ideal flow-field and based on sufficiently small motion assumption, the free surface GF is of remarkable significance in solving those problems with boundary element method.
As to frequency domain, there are generally different categories of Green function according to whether or not the harmonically time dependent unit source beneath a free surface is translating, and whether the water depth is finite or infinite.It should be noted that one concentrates on only the frequency domain infinite depth GF without translating.The evaluation of free surface GF and its derivatives is a complicated mathematical issue, especially the second order derivatives, which are very necessary in high-order panel method.
Free surface GF was known to us because of the work of John [1].There are several different style of mathematical representations [2] for GF.Noblesse [3] advocated the two parts of the GF, which were the so-called near-field and far-field representations.On the basis of that, linear table interpolation fast method was proposed by Ponizy et al. [4], which gave a precision of 10 -5 .In 2017, Wu et al. [5] further proposed simple approximations to the local flow components of GF and its first order derivatives without discussion about the calculation error.
Firstly one considers a pulsing source point ( , , ) q    , an image point ( , , )  q     of the source point relative to the free surface, and a field point ( , , ) p x y z as Fig. 1 showing.The source point and the field point both are lying in the negative half-plane.Z is the vertical distance between field point and image source point.r is the horizontal distance between these two points,

R r Z 
The following expression [1] is the complex GF of infinite depth water in frequency domain Then these two coordinates of X and Y may take on all positive values, one quadrant of a twodimensional plane should be considered.The equation ( 2) may be written in the form Where The elementary singularity 11 can be implemented maturely with numerical or analytic method.So the issue has been translated from the evaluation of ( , ) G p q  and its derivatives to those of ( , ) F X Y .Taking further treatment to infinite integration [6,8] of equation, yields, From equation (5), one can derive the following first and second order partial derivatives of function ( , ) F X Y with respect to independent variable X andY : ( 2 ) Where Here one can find that ( , ) F X Y X  respectively.So if the calculation of GF and its derivatives are based on Equation ( 5), the calculation emphasis is the evaluation of ( , ) Besides, the special function of Struve function, Bessel function of first order and Bessel function of second order have the following identity [15] which will be implemented in the latter manipulation.Furthermore, without loss of generality and ambiguity, here one assumes that the ( , ) F X Y represents the GF in the latter discussion.

Series Expansion Method
The division of the whole domain of physical importance is the important precondition of numerical evaluation of GF, in this paper which is derived from different kinds of Series Expansion Method (SEM).In this part, considering the different location of X and Y in the XY Highly Precise Approximation of Free Surface Green Penghao Shan, Jiemeng Wu Function and Its High Order Derivatives Based on Refined Subdomains 57 two-dimensional plane, one primarily outlines four different kinds of SEM.On the basis of these kinds of SEM, the refined boundary of the whole domain is discussed in the latter part.

SEM1
When X is relatively small beside Y axis,   0 J Xt can be expanded as the following identity in even powers of


Substituting the above equation into equation ( 4), and after successive partial integration, yields, Here, all these equations are consist of a double infinite series with positive powers of X and negative powers of Y , one can implement equations ( 12)~ ( 17) to approximate the GF ( , ) F X Y and its derivatives.
Penghao Shan, Jiemeng Wu Highly Precise Approximation of Free Surface Green Function and Its High Order Derivatives Based on Refined Subdomains 58

SEM2
To some extent, this SEM works when X and Y is moderate value, relative to the other 3 SEM.The following identity exists [14]: So the evaluation will be very convenient if one can makes the integral part of equation ( 5)~( 10) to be expressed with   . So ( 6) and ( 8) can be transformed as follows.
As to the integral part of equation ( 6), Making twice trigonometric substitutions ( tan tX   ) and partial integration, after some tedious manipulations, yields,


Substituting the above equation into equation ( 6), thus, The same careful manipulations can be taken with equation ( 8), it becomes So one arrives at the required expression about ( , ) Xt  . In this SEM, kernel of the issue is the evaluation of the equation ( 5), (19) and (20).Furthermore, one introduces the following identity:    , the integral range about X and Y that is divided into the following four parts [14], is illustrated as fig. 2.
One makes substitution t r X  to the integral part of the above equation, and implements the Taylor expansion to 2   1 Then utilizing the definition of Gauss hypergeometric function is the Gamma function.Meanwhile, one introduces the following two identities.

SEM3
When X is not very small and X/Y is less than 2, one implements the Taylor expansion to   Substituting the above equation into the integral part of (5), and exchanging the order of integration and summation, yields Where Successive partial integration of the right integrals of Substituting equation ( 31) into ( 5), one has From equation ( 32 ) Finally, the evaluation can be implemented with the required expressions (32) (33) (34) ( 7) )(9) (10).

SEM4
When X and Y are all not very small, One makes substitutions


Substituting the above equation into (35), and exchanging the order of integration and summation, yields Similar to the subdomain 3, one has ( ) Finally, one obtains the required expressions (37) (38) (39) (7) (9) (10) with which the evaluation can be implemented.
Highly Precise Approximation of Free Surface Green Penghao Shan, Jiemeng Wu Function and Its High Order Derivatives Based on Refined Subdomains 63

Numerical Results and Discussion
From final identities of all the above four kinds of SEM, one concludes that there are two evaluating emphases to be accomplished.One is accurate and efficient approximations of those special functions, which are 0 () YX , 1 () YX , 0 () JX, 1 () JX, 0 () HX, 1 () z , the other is the exact division of the whole domain of physical importance.For convenient comparison, there are four kinds of different methods.The first kind method [14] is called M1, which will be used to the comparison of computational precision.M2 stands for the classical fast method [6].The third method derived from equations (5) (6) (7) (8) (9) (10), called Mpre, is the direct numerical method, which can provide double precision result employing Romberg integral method.The last method from this paper is called Mnew.One benchmarks the corresponding numerical errors of M1 and Mnew against the Mpre results, with the representations E1=|M1-Mpre| and Enew=|Mnew-Mpre|.
The notations T1, T2, Tpre and Tnew represent the computational time of M1, M2 and Mpre, Mnew respectively.Although X and Y may take on all positive values, without loss of generality one assumes that for error analysis, calculation range of X and Y is up to 40, in which the step length of X and Y is 0.2, and that for efficiency analysis, calculation interval of X and Y is (0, 500] with variable step lengths, which are 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 respectively.The constant range of X-Y plane with changeable step length is corresponding to one certain wetted surface with different number of discrete panel elements.The following numerical calculations were performed on a desktop computer (Microsoft Windows 7, 64 bit, Intel Core i5-5200 2.2 GHz/ 4 core, 4 GB implementation memory).

A note on the approximations of special functions
Efficient approximations of 0 () JX are derived from rational functions [16], which are theoretically accurate to at least 18 significant decimal digits.The fast calculation method about 0 ()

21
, , , F a b c z is a challenging task [19] [20].Here one introduces a fast and accurate computation method of this function [19].
Furthermore, the computational efforts of all these special functions are also acceptable, outlined in Table 1 (the consumed time is for one time, which is the average computational efforts from 10 5 times calculation of the special function).It can be found that the calculation of Gauss hypergeometric function is a little time-consuming, compared with those of other special functions.
4.2 The refined division of four subdomains Generally, exact division of the domain is closely related to the calculation precision and efficiency under appropriate SEM.In this paper, the algorithmic applicability with different (X, Y) and constant truncated number of each SEM is outlined on the X-Y two-dimensional plane.
Penghao Shan, Jiemeng Wu Highly Precise Approximation of Free Surface Green Function and Its High Order Derivatives Based on Refined Subdomains 64 Furthermore, the exact division and different truncated number of SEM is derived from the balance between calculation precision and efficiency.Firstly, one derives the representative numerical E new of ( , ) F X Y with constant truncated numbers (N 1 =15, N 2 =35, N 3 =13, N 4 =15 respectively) of each SEM.
From Fig 3(1), it can be clearly found that SEM1 is available when Y/X is more than 2.0 in one quadrant of the two-dimensional X-Y plane.Fig. 3(2) illustrates that the approximation of the square region (when X is less than 10 and Y is less than 11) can be derived from SEM2.Fig. 3(3) shows that SEM3 is suitable when X/Y is more than 2.0 and X is more than 6.5.From Fig. 3(4) one can conclude that SEM4 can be applied to the whole zone except when X is less than 11.0 and Y is less than 31.0.
(1) Numerical Enew of SEM1 with N1=15 (2) Numerical Enew of SEM2 with N2=35 (3) Numerical Enew of SEM3 with N3=13 (4) Numerical Enew of SEM4 with N4=15 Fig. 3 Algorithmic applicability of each SEM with constant truncated number However, from the computational efficiency of special functions presented earlier, one believes that efficiency sorting order of these four SEM is: SEM1  SEM3> SEM4> SEM2, as to exact division of the whole domain, the above coarse result is far from enough.Consequently based on the general consideration for three aspects, which are the efficiency difference between these SEM, the computational error control of transitional region (especially the exact boundary of the SEM2), the balance between the calculation precision and efficiency of the whole domain, one outlines the refined boundary and truncated number of infinite series for each SEM as follows.In Fig. 4 of X and Y coordinate plane, the D i is shorted for subdomain i (i=1, 2, 3, 4), N i stands for the truncated number of infinite series for subdomain i (i=1, 2, 3, 4).The boundaries of the whole calculating domain are Y/X=2, Y=15, X=9.5, X=6.5, and X/Y=2.With the desired accuracy as the precondition, for the sake of efficiency, D 1 , D 2 , D 4 is re-divided into different zones.It must be pointed out that this introduction of the re-division for D 1 , D 2 and D 4 with different truncated number did not increase computational efforts in practice compared those derived from selection of evaluation algorithm.D 1 is re-divided into two zones with different truncated number of infinite series, its following boundary is Y=14 and Y=17, with truncated number N 1 =15 and N 1 =19 respectively.D 2 is re-divided into four zones with different truncated number of infinite series, its following boundary is Y=6 and x=4.5, Y=8 and X=7.5, and Y=11.The truncated number of these four zones are N 2 =26, N 2 =30, N 2 =36, and N 4 =42 respectively.D 4 is re-divided into two zones with boundary X=14, the truncated number of these two zones are N 4 =20, N 4 =15 respectively.

Numerical results of Green function and its derivatives
Fig. 5 Numerical E1 of ( , )  F X Y with M1 Fig. 6 Numerical Enew of ( , ) F X Y with Mnew Fig. 5 shows numerical E 1 for ( , ) F X Y with M 1 .One finds that M 1 can obtain the accuracy of 10 -9 in some area, not covering the whole domain.The precision of some sizeable areas is about 10 -6 ~10 -8 .The numerical results also shows that there are actually some areas with the precision of 10 -2 ~10 -5 , to some extent, which is not sufficiently accurate for computational Penghao Shan, Jiemeng Wu Highly Precise Approximation of Free Surface Green Function and Its High Order Derivatives Based on Refined Subdomains 66 application, such as for the wave-induced motions, forces, and resistance [21].However, Fig. 6 gives numerical E new for ( , ) F X Y .It can be obviously concluded that precision of all the evaluating points with different X and Y from the whole domain can reach at least 10 -9 , which is due to the proper SEM, the highly precise approximation of special function and refined subdomain division.The same results can be found in Table 2, which is the comparison of numerical results between M 1 and M new .As a result, Fig. 7 and 8 show that evaluations of the first and second order derivatives of GF can also obtain a precision of 10 -9 at least under M new presented earlier.Furthermore, one can rationally expect satisfactory calculation results of GF other derivatives due to the relationship of equations (6) (7) (8) (9) (10).
In addition, the computation time is also acceptable.Efficiency of the new method is compared with that of the rest three methods.From Table 3, one finds that T 1 , T 2 , T new have the same order of magnitude, which are much less than T pre derived from the direct integration.Except from T pre , T new is more approximate to T 2 .9 Change trends of GF and its derivatives when X and Y is up to 40 Penghao Shan, Jiemeng Wu Highly Precise Approximation of Free Surface Green Function and Its High Order Derivatives Based on Refined Subdomains 68 Fig. 9 shows that along the X direction, which is also the propagation direction of the surface wave, the variation trend generally meets the cyclical variation law with continued amplitude attenuation, and that along the Y direction, which is the increasing direction of water depth, wave amplitude sharply attenuates to zero.The above law accords with the actual physical phenomenon.The above-described numerical method of evaluating GF is tested here through a floating hemisphere with radius of 25 meters (see Fig. 10).The wetted surface of this hemisphere is represented by an ensemble of connected four-sided panels.A panel degenerates to a triangle when the coordinates of two vertices coincide.Total number of panel is 1008.The centroid of each panel is used as field point or source point.Here one investigates GF between a fixed field point located at (-19.71,-0.98,-2.93)and all the points regarded as unit sources.So 1007 times calculations need to be solved except for the situation when the field point coincides with the source point.In addition, considering the range of the subdomains, it is reasonable to assume wave number equals to 0.8.Then the resulting (X, Y) is outlined in Fig. 11.
Highly Precise Approximation of Free Surface Green Penghao Shan, Jiemeng Wu Function and Its High Order Derivatives Based on Refined Subdomains 69 One chooses four points from XY plane in each subdomains to make the verification.Table 4 depicts these results obtained by M new and M pre .As to the hemisphere, it can be concluded that evaluations of the first and second order derivatives of GF can obtain a precision of 10 -9 at least from all the listed data with accuracy of 11 decimal digits.

Conclusion
This paper explored in detail the classical representation of free surface GF, composed of a semi-infinite integral involving a Bessel function and a Cauchy singularity.For the evaluation of the function itself and its high order derivatives, only four kinds of analytical series expansion was developed, and the involved special functions, especially the hypergeometric function were evaluated with high precision and efficiency.In addition, one put forward the refined subdomains of the whole domain.The numerical results showed that the new method could acquire a high degree of precision and an acceptable and practical efficiency.

pqR
is the distance between them, pq R is the distance between field point and image source point.

0 JFig. 1
Fig. 1 Layout and notations of field and source point

YX
is the m-th order derivatives of Bessel function of the second type.() () m n HX is the m-th order derivatives of Struve function.() () m n JX is the m-th order derivatives of Bessel function of the first type.

Fig. 2
Fig. 2 Four integral parts of u-t two-dimensional plane For the convenient evaluation of the above equation, One introduces the coordinate substitution cos ur   , sin tr   27)Substituting equation (26) into(19), (27) into (5), (26) (27) into )(20) respectively and using the relation of special functions presented earlier, one gets

X
), one can achieves the derivatives of function ( andY , and simplifies those equations with the former special function identity, thus

HighlyFig. 4
Fig. 4 Refined boundary and truncated number of infinite series for each SEM.

4. 4
Application to a floating hemisphere

Fig. 10
Fig. 10 The wetted surface of the hemisphere Fig. 11 Different location of X and Y of the assumed situation

Table 1
Computational efforts of special functions Special function Consumed time/ s/once Special function Consumed time/ s/once 2F1

Table 3
Comparison of the computational efforts between four methods

Table 4
Comparison of the hemisphere results