Prof Harry Bingham

Dr Richard Porter

Dr Marion Darbas

Dr Gérard Delhommeau

Hui Liang

Youngmyung Choi

Yunxing Zhang

Youngjun Kim

Lei Wang

Kun Zheng

Wenjin Cao

Zhaobin Li

Yang Zhang

Chunmei Xie
email: chunmei.xie@ec-nantes.fr

François Rongère

Alain H Clément

Aurélien Babarit

Xiaobo Chen

Comparison of existing methods for the calculation of the infinite water depth free-surface Green function for the wave-structure interaction problem

Keywords:

First and foremost, I would like to express my sincere gratitude to my supervisor Aurélien Babarit for his continuous support of my Ph.D study and related

List of Figures

6-2

The nondimensional added mass and damping coefficients of a heaving and surging hemisphere with (dashed lines with circle markers) and without (solid lines) irregular frequencies . . . . . . . . . . . . .

6-3

The nondimensional added mass and damping coefficients of a surging hemisphere of inhouse code (solid lines), Hydrostar (green red and blue circle markers), analytical solution (black circle markers) . . digital accuracy between Direct Integral method and ODE-based method for Boxbarge case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017 . . . . . 6.4 Main dimension and hydrostatic characteristics of the KCS . . . . . . 6. 5 The maximum absolute difference of added mass coefficient with 8 digital accuracy between Direct Integral method and ODE-based method for KCS case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017 . . . . . . . . 6. 6 The maximum absolute difference of damping coefficient with 8 digital accuracy between Direct Integral method and ODE-based method for KCS case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017 . . . . . . . .

List of Tables

Chapter 1 Introduction Résumé

Dans le premier chapitre, le contexte et le sujet de la thèse sont présentés. Nous étudions les efforts et les mouvements induits par les vagues dans le cadre de la théorie des écoulements potentiels avec la méthode des éléments frontières. Le contexte de la méthode des éléments de frontière et de la fonction de Green, qui est la principale préoccupation de cette thèse, est donné. Ensuite, les contributions principales et les grandes lignes de la thèse sont présentées.

Wave-induced loads and motions

The oceans cover over 70 percents of the Earth's surface. They have been traveled and explored since ancient times. The development of marine science and technology in the last century has made the extensive exploitation of the oceans' resources. Nowadays, they provide important living and nonliving resources, such as fish; oil and natural gas; minerals and renewable energy. The need to further harvest the ocean resources is growing up with increasing world population and the lack of land resources. Ships for transportation of passengers and goods, fish farms in open sea, oil and gas platforms, platforms for the recovery of shallow-sea and deep-sea minerals, wind energy farms and wave energy converters are the main structures deployed at sea in order to exploit its resources. For the design of all those structures, the understanding of wave-structure interaction effects is critical in order to provide reliable estimates of wave loads.

A lot of researchers have investigated wave-structure interaction effects. The methods to study the wave-induced loads and motions can be split into : theoretical analysis; experiments; numerical calculations. Theoretical analysis is powerful and highly reliable but it is limited by strict assumptions. Only few cases with simple geometries and simple boundary conditions can be solved theoretically. They are usually considered as validation tests for model tests or numerical calculations.

For experiments, full-scale trials are desirable but expensive and difficult to perform under controlled conditions. Therefore, lots of countries have built towing tank and water tunnels in order to perform model tests. In model tests, dimensional analysis is first applied to find the main characteristics and define the parameters. An issue with model test is that it can be challenging to scale tests' results to full scale when viscous hydrodynamic forces matters. Moreover, available equipment at the test facilities may also limit the experimental possibilities. Last but not least, model tests are usually expensive.

Due to the fast advance of computer technology, numerical modelling play an increasingly important role for the calculation of wave-induced loads and motions. Numerical results can usually be achieved at much cheaper cost and faster than experimental results. Nevertheless, model tests are still indispensable as they can provide insight in the physical phenomena for the numerical calculations and also for the validation of the numerical calculations.

In summary, the three methods complement each other in hydrodynamic study. In this thesis, we focus on numerical modelling of wave-structure interaction.

A variety of hydrodynamics problems concerning ships and offshore platforms can be addressed within the context of potential flow theory (in which viscous effects are neglected) [START_REF] Newman | Panel methods in marine hydrodynamics[END_REF]. Those problems include the wave resistance and motion of ships in presence of waves [25], [START_REF] Noblesse | On the theory of potential flow about a ship advancing in waves[END_REF], [START_REF] Liu | Prediction of added resistance of ships in waves[END_REF], the dynamic response of offshore and wave/wind energy converters to wave action [8], [30], [32], [START_REF] Sánchez | Control performance assessment and design of optimal control to harvest ocean energy[END_REF], propeller performance and interactions between adjacent ships maneuvering in close proximity [START_REF] Kagemoto | Interactions among multiple threedimensional bodies in water waves: an exact algebraic method[END_REF], [13], [START_REF] Kashiwagi | Wave drift forces and moments on two ships arranged side by side in waves[END_REF], [START_REF] Xiang | Time domain simulation of two interacting ships advancing parallel in waves[END_REF], ... In all those cases, the Reynolds number is large. Provided that the geometry is streamlined or if the Keulegan-Carpenter number is small, the separation of the boundary layer is avoided. Thus, viscous effects can be neglected and potential flow theory can be used. Within this framework, the most efficient and widely used numerical method is the boundary element method (BEM), also named the panel method.

The background of boundary element method and

Green function

Having been developed for over four decades, the BEM firmly stands in the arena of the numerical methods for solving the partial differential equations, such as Helmholtz equation.

The concept of the BEM method is to convert the volume problem into a problem on the surface of the domain boundaries. It is achieved by converting the problem defined by partial differential equations in the domain into an integral equation over the boundaries of the domain. Compared to other numerical methods such as finite element method, finite difference method, finite volume method which are classified as domain methods, the boundary element method reduces the spatial dimensions. It leads to easier generation of the mesh, usually smaller linear systems, less computer memory requirements and more efficient computation. However, note that the matrices of the linear problem involved in BEM are usually dense and non-symmetrical.

In BEM, the velocity potential is obtained through the distribution of singularities (source and/or normal dipoles) over on panels representing the geometry. Thus, the BEM is also considered as a semi-analytical method due to the integral nature of singularities used in the formulations. This semi-analytical nature makes the BEM highly accurate. Exterior problems defined by unbounded domains but bounded boundaries can be handled as easily as interior problem. Another advantage of the BEM is that the physical values at the boundaries can be obtained directly which can be important for some applications.

The early history of the boundary element method from the 18th century up to the late 1970s is reviewed in [18]. It covers potential theory, Laplace equation, the existence and uniqueness of the solution of boundary value problems, the Gauss and Stokes theorems that allow the reduction of spatial dimensions, the Green's identities, Green's function, the Fredholm integral equations and the extension of Green's formula to acoustics, elasticity, hydrodynamics and other physical problems. The authors also celebrated the pioneers behind these mathematical developments with short biographies. Another comprehensive review up to 1980s was written by Tanaka [START_REF] Tanaka | Some recent advances in boundary element methods[END_REF]. One can refer to those two review papers for the background, the theory and the applications of the BEM.

Numerical implementations of boundary integral equations in hydrodynamics can be classified according to: which singular element is distributed on the boundaries (sources , normal dipoles, mixed of sources and normal dipoles • • • ); which approximation is used for the distribution (constant, linear, or higher order distribution); and whether the used Green function satisfies the free-surface condition.

The use of boundary element method in hydrodynamic problems did not blossom out until the seminal work for three-dimensional bodies with arbitrary geometries by Hess and Smith [START_REF] Hess | Calculation of non-lifting potential flow about arbitrary three-dimensional bodies[END_REF], [START_REF] Hess | Calculation of nonlifting potential flow about arbitrary threedimensional bodies[END_REF], [START_REF] Hess | Calculation of potential flow about arbitrary bodies[END_REF], which was based on their previous works for two-dimensional and axis-symmetric bodies [START_REF] Smith | Exact Solution of the Neumann Problem: Calculation of Non-circulatory Plane and Axially Symmetric Flows about Or Within Arbitrary Boundaries[END_REF], [39]. In their work originally dedicate to aerodynamics, an unbounded domain is considered (no free-surface). The body surface is discretized by quadrilateral panels with a constant source density distribution on each panel. The pulsating sources and/or sinks are placed at the center of the panels. The integral equation for the body surface source density is replaced by a set of linear algebraic equations for the values of the source density on each panel. Once the linear system built, the Seidel iterative procedure can be used to solve it. Then, it is straightforward to obtain the velocity potential, the velocity, the pressure and the forces.

Numerous studied followed Hess and Smith's work. Some uses the source formulation with source density distribution (indirect method) while others apply the potential formulation with a mixed source and normal dipole distribution (direct method). The computational effort involved in building and solving the linear systems of these two methods is similar, albeit there can be differences in the results depending on the applications. According to [START_REF] Newman | Panel methods in marine hydrodynamics[END_REF], for thin bodies, the normal dipoles approach is more stable than sources for representing the cross-flow components. To get the fluid velocity, the gradient of the Green function is needed for the source formulation while for the potential formulation, the second derivatives are required. Since the local velocity induced by the dipole distribution varies rapidly over distances comparable to the panel dimensions, the potential formulation is less accurate for cases in which curved surfaces are approximated by flat panels.

The body surfaces are usually discretized by flat triangles or quadrangles. The density of singularities is constant on each panel in the constant panel method. It is a mature and very efficient method which is widely used in academia [28], [33], [START_REF] Inglis | Comparison of calculated responses for arbitrary shaped bodies using two-and three dimensional theories[END_REF], [START_REF] Korsmeyer | The analysis of wave effects on tension-leg platforms[END_REF] and industry. Almost all commercial seakeeping softwares are based on the constant panel method. They include Hydrostar [2], WAMIT [5], NEMOH [3], AQWA [1], SESAM [4], • • • . Some of them include options to use higher order singularities distributions.

However, for curved body surfaces, the constant distribution is discontinuous between panels. Thus, it may require a large number of panels to achieve sufficient accuracy. To overcome this limitation, Webster [START_REF] Webster | The flow about arbitrary, three-dimentional smooth bodies[END_REF] proposed the linear sources distribution with triangle panels. Higher order panel methods were then proposed in [START_REF] Liu | Comparison of higher-order boundary element and constant panel methods for hydrodynamic loadings[END_REF], [START_REF] Liu | The computation of mean drift forces and wave run-up by higher-order boundary element method[END_REF], [START_REF] Taylor | Wave diffraction theory-some developments in linear and nonlinear theory[END_REF], [START_REF] Taylor | The effect of corners on diffraction/radiation forces and wave drift damping[END_REF], [START_REF] Teng | New higher-order boundary element methods for wave diffraction/radiation[END_REF], [START_REF] Maniar | A three dimensional higher order panel method based on B-splines[END_REF], [START_REF] Maniar | Wave diffraction by a long array of cylinders[END_REF], [START_REF] Lee | Computations of wave loads using a b-spline panel method[END_REF] which essentially represents the state of the art in BEM solutions today. Linear and higher order BEM are expected to be more accurate than the constant panel method for the same discretization level when curved surfaces are involved. However, they are more challenging to implement because of the more complicated mathematics.

Another major contribution of Hess and Smith [START_REF] Hess | Calculation of non-lifting potential flow about arbitrary three-dimensional bodies[END_REF], [START_REF] Hess | Calculation of nonlifting potential flow about arbitrary threedimensional bodies[END_REF], [START_REF] Hess | Calculation of potential flow about arbitrary bodies[END_REF] is their derivation of the analytical expression for the potential and velocity induced by a source distribution for a quadrilateral. Indeed, the direct numerical integration is slow and inaccurate when the field point is near or on the panels. This source only singularity method is also called as the Rankine panel method. The Rankine panel method with sources only was then extended to constant distributions of sources and normal dipoles on panels of arbitrary polygonal shapes [37] and linear distributions on triangular panels [START_REF] Yeung | A singularity-distribution method for free-surface flow problems with an oscillating body[END_REF], [START_REF] Webster | The flow about arbitrary, three-dimentional smooth bodies[END_REF]. Furthermore, the potential and velocity for a distribution of normal dipoles over a quadrilateral panel where the density of singularities is constant, linear, bilinear or of arbitrary polynomial form was given by Newman [START_REF] Newman | Distributions of sources and normal dipoles over a quadrilateral panel[END_REF].

For wave-structure interaction problems, the free-surface condition must be fulfilled. Rankine sources do not satisfy the free-surface condition by themselves. Thus, additional fundamental solutions of Laplace equation subject to the freesurface condition and radiation condition must be added. The combination of the Rankine source and the additional solutions satisfying the free-surface condition is called the Green function. It is named after the British mathematician George Green [35]. This additional fundamental term includes the image source term and the free-surface term in addition to the Rankine source. One important part of numerical burden in the BEM is the evaluation of free-surface Green function. Fortunately, in some cases, exact formulations of the Green function are available such as the Green functions with zero forward speed in infinite or finite water depth and in time domain or frequency domain. In this thesis, only the infinite water depth case (deep water) is considered.

The mathematical formulations of Green function in that case are reviewed in Chapter 3. Despite the exact formulations, the numerical evaluation of the Green function is still very hard because of its highly oscillatory behavior. Lots of researchers have worked on this issue [START_REF] Newman | An expansion of the oscillatory source potential[END_REF], [?], [START_REF] Noblesse | The Green function in the theory of radiation and diffraction of regular water waves by a body[END_REF], [START_REF] Telste | Numerical evaluation of the Green function of water-wave radiation and diffraction[END_REF], [26], [START_REF] Wu | A global approximation to the Green function for diffraction radiation of water waves[END_REF], [14], • • • . Almost all these methods are based on the integral formulations. A different perspective was proposed by Clément [24], [23] in the time domain and in the frequency domain. He found that the time domain Green function is the solution of a fourth order differential equation (ODE) and the corresponding frequency domain Green function is the solution of a second order differential equation. This provides a new perspective to understand the properties of Green function and also to evaluate the Green function.

The used of the ODE to evaluate the Green function has been implemented in a seakeeping code in time domain [22]. However, it was not done in the frequency domain, which is the main objective of this thesis.

The similar works maybe also be made for finite depth or non-zero forward speed Green functions.

Present contributions

• The exact mathematical expressions of the Green function in infinite depth of water and in the frequency domain with zero speed are collected, which could be useful to better understand the numerical behavior of the Green function.

• Several existing algorithms for the evaluation of the Green function in deep water and in the frequency domain with zero speed are implemented and compared with respect to their accuracy and computational time which could be useful for engineers to choose the most appropriate method for their applications.

• A series of ODEs of the Green function and its gradient with respect to time and frequency are given. The ODEs with respect to the spatial variables are first introduced in this thesis which may give another perspective for the calculation of surface integrals or the volume integrals of the Green function.

• It is challenging to solve the second order ODEs of the Green function and its gradient with usual initial conditions at the original point of variable because of the degenerated properties and the singularities. A new way to remove the difficulty by modifying the ODEs associated with new functions free of the singularity is introduced. The new ODEs are then transformed into their canonic form by using a novel definition of vector functions. The canonic form can then be solved with the initial conditions at the original point since all involved terms are finite.

• The expansion method is applied to find analytical solution of Green function for small frequencies. It is based on series of logarithmic function together with ordinary polynomials

• A new ODE-based boundary element method is developed to calculate the hydrodynamic coefficients (added mass, damping coefficient, excitation force) of marine structures.

Thesis outline

This thesis is organized as follows. In Chapter

Basic equations of the linear potential flow theory

The basic equations of the linear potential flow theory are first recalled. The case study is that of floating or submerged bodies without forward speed with free surface in infinite depth of water and in presence of regular incident propagative waves.

Let us define the three-dimensional Cartesian coordinates system Oxyz. The (x, y) plan coincides with the mean free surface level and the Oz axis is oriented positively upwards. The free surface elevation E is defined by:

z = E(x, y, t) (2.1)
The fluid is assumed to be inviscid and incompressible and the flow is irrotational. Therefore the curl of the velocity -→ V is zero:

∇ × - → V = 0 (2.2)
The curl of the gradient of any continuously twice-differentiable scalar field is always the zero vector. Thus, the flow vector velocity -→ V can be expressed as the gradient of a scalar potential Φ(M, t) at the field point M = (x M , y M , z M ) and time t:

- → V = ∇Φ (2.3)
In potential flow theory, the potential Φ is called the velocity potential.

The mass conservation law for an incompressible fluid is:

∇ • - → V = 0 (2.4)
By combining equations (2.3) and (2.4), one can show that the velocity potential satisfies the Laplace equation:

∇ 2 Φ = 0 (2.5)
The momentum equation for an inviscid flow can be written as:

( ∂ ∂t + - → V • ∇) - → V = -∇( P r ρ + gz) (2.6)
in which P r stands for the fluid pressure, ρ is the density of fluid and g is the gravity. The Bernoulli equation can be obtained from the momentum equation expressed. It reads:

P r ρ + gz + ∂Φ ∂t + ∇Φ • ∇Φ 2 = C(t) (2.7)
where C(t) is an arbitrary function of t. It can usually be absorbed into Φ.

The dynamic behaviour of the free-surface is governed by the kinematic and dynamic free-surface boundary conditions. The kinematic condition is established by considering that a particle of the free surface always stays on the surface. This kinematic condition is obtained by taking the material derivative of (2.1):

∂E ∂t + ∂Φ ∂x ∂E ∂x + ∂Φ ∂y ∂E ∂y - ∂Φ ∂z = 0 (2.8)
Furthermore, the free surface dynamic condition is obtained by neglecting the effect of surface tension and by considering that the pressure given by (2.7) is equal to the atmospheric pressure at the free surface:

gE + ∂Φ(x, y, E, t) ∂t + ∇Φ(x, y, E, t) • ∇Φ(x, y, E, t) 2 = 0 (2.9)
Combining the dynamic and kinematic condition, the free-surface condition is:

∂ 2 Φ ∂t 2 + g ∂Φ ∂z + 1 2 ∇Φ • ∇(∇Φ • ∇Φ) + 2∇Φ • ∇ ∂Φ ∂t = 0 (2.10)
The last assumptions are that the steepness of the incident waves is small and that the body motion is small in comparison to its dimensions. Thus, the velocity potential can be expanded according to the perturbation theory:

Φ(M, t) = Φ (1) (M, t) + Φ (2) (M, t) + • • • (2.11)
The superscript (•) denotes the order of the velocity potential and Φ (2) Φ (1) . Other physical quantities, i.e. free surface elevations, dynamic pressures, wave loads and body's motions, can be expanded in similar forms.

Finally, by introducing the expansion (2.11) into (2.5) and (2.10) and by developing Φ(z = E) on the mean free surface (z = 0), one can obtain the boundary value problems for the different orders.

Boundary value problem at first order

At the first order, the linear boundary value problem for the interaction of a body with waves is: [l]

∇ 2 Φ (1) = 0 M ∈ (τ e ) (2.12a)
gΦ (1) z + Φ

(1)

tt = 0 M ∈ S F (2.12b) Φ (1) n = - → X (1) t • n M ∈ S B (2.12c) ∇Φ (1) = 0 M ∈ (C) (2.12d)
where S F is the mean free-surface, S B is the body surface, -→ X t is the velocity of the body, and (C) This boundary value problem (BVP) includes the Laplace equation (2.12a), the linearised free-surface condition (2.12b), the body surface condition (2.12c) and the radiation condition (2.12d). The the existence and the uniqueness of the solution of the BVP (2.12) for a given incident wave and body motion was prescriptive for non-bulbous body by John [START_REF] John | On the motion of floating bodies II. simple harmonic motions[END_REF].

(x 2 + y 2 + z 2 → ∞) is
The first order free surface elevation is given according to the dynamic condition (2.9) by:

E (1) = - 1 g
∂Φ (1) ∂t (2.13) evaluated at z = 0 . The total dynamic pressure on the body surface is given according to the Bernoulli equation by:

P (1) r = -gρX (1)
3 -ρΦ

where X

(1)

3 the vertical displacement. In this thesis, only the first order problem is considered. Therefore, in the following, for sake of simplicity, the first-order quantities such as Φ (1) , E (1) , P (1) r are written as Φ, E, P r .

Boundary element method

Green function

Let us define the source point, field point and other notations as depicted in Figure 2-2. The source point P (x P , y P , z P ) and the field point M (x M , y M , z M ) are both lying on or under the free surface (z P ≤ 0, z M ≤ 0). The image source point P (x P , y P , -z P ) is the mirror of the source point P with respect to the mean free surface. The horizontal distance between the source point P and the field point M is denoted by r. The vertical distance between the image source point P and the field point M is -Z. The distance between the source and field points is denoted by R and the distance between the image source point and the field point is

R 1 = √ r 2 + Z 2 .
The angle θ is defined by cos θ = -Z/R 1 and sin θ = r/R 1 . The relations between the coordinates are given by: To solve the first-order boundary value problem defined by equation (2.12), a fundamental solution which satisfies the following equations can be used:

r = (x M -x P ) 2 + (y M -y P ) 2 Z = z M + z P R = (x M -x P ) 2 + (y M -y P ) 2 + (z M -z P ) 2
∇ 2 G ∞ (M, P, t) = δ(M, P ) M ∈ τ e (2.16a) gG z + G tt = 0 M ∈ S F (2.16b) G = O(1/r) r → ∞ (2.16c)
in which δ is the Dirac function and δ(M, P ) = δ(x M -x P )δ(y M -y P )δ(z M -z P ).

The solution of equation (2.16) is called the Green function. It gives the field of velocity potential at M (x M , y M , z M ) created by a source located at P (x P , y P , z P ).

Boundary integral equation

Let us apply the Green's second formula to the couple of harmonic functions (Φ, G) [36]:

4πΩΦ(M ) = S(P ) ∂Φ(P ) ∂n(P )
G(M, P ) -Φ(P ) ∂G(M, P ) ∂n(P ) dS(P ) (2.17)

4πΩ is the 'solid angle' over which the fluid is viewed from M and the normal vector. The left hand side is the result of the domain integral while the terms on the right hand side come from the transformation of the volume integral into the surface integral on the boundaries according to the formula of Ostrogradsky. The boundary surface S(P ) includes the body surface S B , the mean free surface S F and the surface at infinity C. Both the integrals on the surface at infinity and on the mean free surface are zero according to equations (2.16) and equation (2.12a).

4πΩΦ(M ) = S B ∂Φ(P ) ∂n(P ) G(M, P ) -Φ(P ) ∂G(M, P ) ∂n(P ) dS(P ) (2.18)
with Ω = 1, 1/2, 0 for M ∈ (τ e ), M ∈ S B , M ∈ (τ i ) with flat panels.

Let us now consider a complementary domain (τ i ) inside the body. It is limited by the body boundary S B and the interior free-surface S F i (see Figure ).

∇ 2 Φ = 0 M ∈ (τ i ) (2.19a) gΦ z + Φ tt = 0 M ∈ S F i (2.19b) Φ n = - - → X (1) t • n M ∈ S B (2.19c)
where Φ is the velocity potential of interior domain.

By applying the Green's second formula:

4πΩΦ (M ) = - S B ∂Φ (P ) ∂n(P ) G(M, P ) -Φ (P ) ∂G(M, P ) ∂n(P ) dS(P ) (2.20)
with Ω = 1, 1/2, 0 for M ∈ (τ i ), M ∈ S B , M ∈ (τ e ) in the constant method. It should be noted that thein the right side hand is because of opposite direction of the normal vector for the interior and exterior domain.

For better use the properties of singularities like sources and dipole and to get the gradient of velocity potential around the body easily. We do addition of Eqns.

Wave radiation and diffraction

Let us further assume that the velocity potential is harmonic. Thus, the velocity potential and Green function can be written as:

Φ(M, t) = (ϕ(M )e -iωt ) (2.22a) G(M, P, t) = (G(M, P )e -iωt ) (2.22b)
in which (•) means to take the real part and ω is the wave frequency. . Linearity allows the decomposition of the body-generated wave perturbation potential into the radiation potential ϕ R and the diffraction potential ϕ D . The details of the decomposition can be seen in Wehausen and Laitone's book [START_REF] Wehausen | Surface waves[END_REF]. The velocity potential can thus be expressed as:

ϕ = ϕ R + ϕ D + ϕ I (2.23)
in which ϕ I is potential of the incoming waves.

For the radiation potential, we can further decompose it into six components:

ϕ R = -iω 6 j=1 ξ j ϕ j (2.24)
where ξ j , j = 1, ..., 6 denote the amplitudes of the body displacement in its six body degrees of freedom.

For the diffraction potential and the potential of the incident waves, one can write:

ϕ D = ξ 0 ϕ 7 (2.25a) ϕ I = ξ 0 ϕ 0 (2.25b)
where ξ 0 is the (complex) amplitude of the incident wave. The complex radiation velocity potential components ϕ j , j = 1, ..., 6 and the diffraction velocity potential ϕ 7 satisfy the Laplace equation in the fluid domain.

They also satisfy the free surface condition on the z = 0 plane:

ϕ jz -k 0 ϕ j = 0, j = 1, ...7 (2.26)
with k 0 = w 2 /g the wavenumber.

The radiation and diffraction velocity potential satisfy the Sommerfeld radiation condition at infinity for infinite depth water. Finally, the body boundary conditions are:

ϕ jn = n j , j = 1, ..., 6 (2.27
)

ϕ 7n = -ϕ 0n (2.28)
where n i , i = 1, 2, 3 are the Cartesian components of a unit vector normal to the body surface and pointing out of the fluid domain, and

(n 4 , n 5 , n 6 ) = (x M - x 0 , y M -y 0 , z M -z 0 ) × (n 1 , n 2 , n 3 ) with O(x 0 , y 0 , z 0 ) the reference point of rotation.
The incident wave potential is defined by: ϕ 0 = -ig ω e k 0 z e ik 0 x cos β+ik 0 y sin β (2.29)

Where β is the propagation angle of the incident waves relatively to the positive x-axis.

Source distribution

By introducing (2.22) into (2.21), by choosing ϕ = ϕ on S B and by defining σ = ϕ n -ϕ n as the source density, one obtains:

ϕ = S B -iω 6 j=1 ξ j ϕ j + ξ 0 ϕ 7 GdS + ξ 0 ϕ 0 (2.30a) ϕ j = S B σ j GdS (2.30b)
The source distribution is determined by satisfying the boundary condition on the body surface S B :

4πΩ S B σ j (M ) + S B σ j (P ) ∂G(M, P ) ∂n(M ) dS(P ) = n j (M ) j = 1, 2, • • • , 6 -∂ϕ 0 /∂n(M ) j = 7 (2.31)

Mixed distribution

An alternative way is to consider a mixed sources and dipoles distribution. The dipole density is defined by µ = ϕ -ϕ . Let us define ϕ = 0. Thus, the dipole distribution is the velocity potential itself and the boundary integral equation reads:

4πΩ S B µ j (M ) + S B µ j (M ) ∂G(M, P ) ∂n(P ) dS(P ) =        S B n j (P )G(M, P )dS j = 1, 2, • • • , 6 S B
(-∂ϕ 0 /∂n) (P )G(M, P )dS j = 7

(2.32)

Numerical solution

The boundary element method is a popular numerical computational method for solving linear partial differential equations (2.12) which have been formulated as integral equations (2.31) and (2.32) [6], [7], [START_REF] Partridge | Dual reciprocity boundary element method[END_REF], [12], [37], [26], [START_REF] Lee | Computation of Wave Effects Using the Panel Method[END_REF], [START_REF] Newman | Boundary-element methods in offshore structure analysis[END_REF], [15].

In this thesis, the mixed source and dipole distribution is considered and the body wet surface is discretized in flat panels (quadrangles). The density of source and dipole distributions on each panel is assumed to be constant. Therefore, equation (2.32) can be transformed into two linear systems with influence coefficients depending on the evaluation of Green function and its gradient.

Equation of motion of a floating body

Once the integral equation (2.32) is solved, the first-order pressure can be obtained by equation (2.14). Moreover, the integration of the first terms in equation (2.14) plus the variation of the gravity loads leads to the hydrostatic stiffness coefficients (C kj ) in (2.34). The second term in (2.14) leads to the definition of the added-mass coefficients (A kj ), the damping coefficient (B kj ), and the wave exciting loads F k . They are defined respectively by integration of the pressure due to radiated waves and the pressure due to incoming and diffracted waves:

ω 2 A kj + iωB kj = ω 2 ρ S B ϕ j n k dS and F k = -iωρξ 0 S B (ϕ 0 + ϕ 7 )n k dS (2.33) with k, j = 1, • • • 6
Finally, the motion equation can be formulated as:

6 j=1 -ω 2 (M kj + A kj ) -iω(B kj + B kj ) + C kj + C kj ξ j = F k (2.34)
where M kj is the inertia matrix. (B kj ), (C kj ) are artificial damping coefficient and stiffness coefficients.

Removal of irregular frequencies 2.3.1 The irregular frequencies issue

When the boundary integral formulations (2.31) and (2.32) are used, irregular frequencies arise in the solution of the exterior problem. It was first pointed out by John [START_REF] John | On the motion of floating bodies II. simple harmonic motions[END_REF] that when the body surface boundary intersects the free surface, the integral boundary equation has eigenfunctions for certain values of the frequencies (eigenvalues called irregular frequencies).

It should be noted that the irregular frequencies effect is due to the formulation of the boundary integral but that does not reflect an irregularity in the solution of the BVP. As a matter of fact, the eigenfunctions of the integral boundary equation for exterior problem are induced by these of the interior problem. Indeed, when the solution of the exterior problem is expressed in terms of sources or/and dipoles distributions over the body surface, it assumes an interior fluid motion such that the potential is continuous across the body surface. If the hypothetical interior problem has eigenfunctions, the corresponding source or/and dipole distribution become undetermined. As a consequence, the solution of the exterior problem cannot be expressed by the source or/and dipole density on the body surface. In summary, the irregular frequencies are certain discrete frequencies for which the solution of the integral equation is not unique.

Numerically, it translates into the determinant of the matrix of the discretized equation to be equal to 0. Near the irregular frequencies, the condition number of the matrix becomes large, which translates into large errors in the solution over a substantial frequency band around the irregular frequencies. Note that the width of the 'polluted' band is reduced when the number of panels is increased.

The removal of the irregular frequencies

The detrimental effect of irregular frequencies in applications was not evident until the late sixties. In 1967, Frank [31] developed a two-dimensional panel method with source distributions for the prediction of the hydrodynamic pressure, force and moment on cylinders. His numerical results showed substantial errors in the prediction of the hydrodynamic coefficients near the irregular frequencies. After Frank's [31] work, researchers have investigated the removal of the irregular frequencies from modifying the system of the internal system or the external system to decrease the condition number of the linear system.

In 1975, Ohmatsu [START_REF] Ohmatsu | On the irregular frequencies in the theory of oscillating bodies in a free surface[END_REF] applied an additional rigid 'lid' on the body waterplane area to close the interior problem in two dimension. This approach is called the extended boundary condition method. The presence of the lid eliminates the resonance associated with the interior Dirichlet eigen problem and therefore eliminates the irregular frequencies for the exterior problem. Ohmatsu [START_REF] Ohmatsu | On the irregular frequencies in the theory of oscillating bodies in a free surface[END_REF] also gave the theoretical demonstration for a circular cylinder that the eigenfrequencies of the interior problem is the irregular frequencies in exterior problems. The extended boundary condition method is effective for two-and three-dimensional problems with a source distribution and the Green method at the cost of using additional panels on the body waterplane areas. A variant of this method was proposed in 1998 by Malenica and Chen [START_REF] Melenica | On the irregular frequencies appearing in wave diffraction-radiation solutions[END_REF]. They moved down the interior surface below the free-surface plane in order to avoid the problem of logarithmic singularity of Green function. Unfortunately, this approach cannot remove all the irregular frequencies.

Ogilvie and Shin [START_REF] Ogilvie | Integral-equation solutions for time-dependent free-surface problems[END_REF] removed the first irregular frequency in two dimensions by placing a point wave source on the interior free surface acting as an absorber of the energy of the corresponding resonant Dirichlet eigensolution. This approach follows Ursell's suggestion [START_REF] Ursell | Short surface waves due to an oscillating immersed body[END_REF] which was to add a source at the origin with the strength selected to absorb the energy associated with a possible sloshing mode. This approach is called the modified Green function method. Ogilvie and Shin discussed the symmetric or asymmetric problems for which we just need to modify the Green function. This approach is simpler than Ohmatsu's method [START_REF] Ohmatsu | On the irregular frequencies in the theory of oscillating bodies in a free surface[END_REF] as it involves only a minor change in the Green function to avoid the eigenfrequencies near the frequencies under consideration. In 1981, Ursell [START_REF] Ursell | Irregular frequencies and the motion of floating bodies[END_REF] found that the use of a set of singularities is necessary for the removal of more than the first irregular frequency in two-dimension. Moreover, the optimal choice of the numbers and locations of the singularities is challenging.

Another method called the modified integral-equation method was introduced by Lee and Sclavounos [START_REF] Lee | Removing the irregular frequencies from integral equations in wave-body interactions[END_REF] by extending the work of Burtin and Miller for acoustic wave scattering problem [11]. They proposed a modified integral equation involving the linear superposition of the classical Green function and its normal derivative with respect to the field point. The normal derivative should be multiplied by a purely imaginary constant to remove all irregular frequencies. It is shown that an optimal value of this imaginary constant exists. It can be determined by minimizing the condition number of the modified integral equation at the first irregular frequency of the Green function. A sphere and a truncated vertical cylinder in three-dimension were tested to illustrate the effectiveness of this approach.

In 1994, Zhu [START_REF] Zhu | Irregular frequency removal from the boundary integral equation for the wave-body problem[END_REF] investigated the modified Green function method and the extended boundary condition method. They compared them to the modified integral equation method. They used the WAMIT code [5]. They found that for a circular cylinder and a rectangular barge, the extended boundary method is the most effective.

Method I: The modified integral-equation method

Burtin and Miller [11] suggested the use of a linear combination of Dirichlet and Neumann boundary conditions which leads to an interior eigenproblem accepting trivial solutions at all frequency:

ϕ + α ∂ϕ ∂n = 0 (2.35)
where ϕ is the velocity potential on the interior boundary and α is a complex constant.

There are no irregular frequencies when α has a nonzero imaginary part according to Lee and Sclavounos [START_REF] Lee | Removing the irregular frequencies from integral equations in wave-body interactions[END_REF] (2.38) There are two disadvantages for this method. The first one is the extra complex constant α. The accuracy of the method depends on the selection of the constant. If its value is too large, numerical errors are introduced because of the additional integral equation associated with the normal derivative. On the other hand, if its value is too small, the condition number of the linear system will still be large and the irregular frequencies cannot be removed completely. Moreover, Zhu and Lee [START_REF] Zhu | Removing the irregular frequencies in wave-body interactions[END_REF] found out that a larger number of panels are required to achieve the same accuracy as with the original integral equation. The second disadvantage is the calculation of the extra normal derivative which increases the computational cost.

Method II: Extended boundary condition method

Paulling and Wood [START_REF] Paulling | Talk in the State-of-the-Art Report on Seakeeping of the 16th American Towing Tank Conference[END_REF] suggested that a 'lid' in the body at the level of the undisturbed free-surface on which a rigid wall boundary condition would be applied would remove the eigenfunctions for the interior problem [START_REF] Paulling | Talk in the State-of-the-Art Report on Seakeeping of the 16th American Towing Tank Conference[END_REF]. In 1975, Ohmatsu [START_REF] Ohmatsu | On the irregular frequencies in the theory of oscillating bodies in a free surface[END_REF] proved this method theoretically and numerically.

In this method, the integral equation becomes:

2πϕ(M ) + S B ϕ(P ) ∂G(M, P ) ∂n P dS(P ) + S F i ϕ(P ) ∂G(M, P ) ∂n p dS(P ) = S B ∂ϕ(P ) ∂n P G(M, P )dS(P ) M ∈ S B (2.39a) -4πϕ(M ) + S B ϕ(P ) ∂G(M, P ) ∂n P dS(P ) + S F i ϕ(P ) ∂G(M, P ) ∂n P dS(P ) = S B ∂ϕ(P ) ∂n P G(M, P )dS(P ) M ∈ S F i (2.39b)
This method requires additional panels on the free-surface inside the body. Therefore the computational cost increases although the additional panels can be much smaller than the panels on the body surface. Note that the required number of panels on the free-surface inside the body increases with the order of the irregular frequencies.

The Green function at the free-surface is oscillating and divergent with increasing frequency. Thus, it can introduce large errors for high frequencies. Therefore a variant of this method was proposed by Malenica and Chen [START_REF] Melenica | On the irregular frequencies appearing in wave diffraction-radiation solutions[END_REF]. It consists in moving down the 'lid' to avoid the oscillation and divergence. The only difference is that the -4π must be changed into -2π in (2.39b). A drawback of this approach is that it can not remove all the irregular frequencies.

Method III: Modified Green function method

In this method, sources are added at points M i inside the body such as ϕ (M i ) = 0 on the interior surface S F i .

∆ϕ (M ) = 0 M ∈ (τ i ) ϕ z (M ) -k 0 ϕ (M ) = 0 M ∈ S F i ϕ (M ) = - - → X (1) t • n M ∈ S B ϕ (M i ) = 0 M i ∈ S F i (2.40)
It leads to the system below which is a overdetermined system.

2πϕ(M ) + S B ϕ(P ) ∂G(M, P ) ∂n P dS(P ) = S B ∂ϕ(P ) ∂n P G(P, M )dS(P ) M ∈ S B (2.41a) S B ϕ(P ) ∂G( s, P ) ∂n P dS(P ) = S B ∂ϕ(P ) ∂n P G( s, P )dS(P ) s ∈ S F i (2.41b)
Challenges of this approach are that s should not be on the nodal points of the eigenmodes, that it involves an overdetermined system which is less efficient than solving a square matrix, and that choosing the sources point and the number of the sources point is difficult for different mesh and different geometries.

An extension is to combine (2.41a) and (2.41b).

2πϕ(M ) + S B ϕ(P ) ∂ Ĝ(M, P ) ∂n P dS(P ) = S B ∂ϕ(P ) ∂n P Ĝ(M, P )dS(P ) M ∈ S B
(2.42) Following this approach, Ogilvie and Shin [START_REF] Ogilvie | Integral-equation solutions for time-dependent free-surface problems[END_REF] added a source at the origin of the interior free-surface. They defined the modified Green function Ĝ in twodimensions as shown in the following: For symmetric problems

Ĝ(M, P ) = G(M, P ) + αG(M, s)G( s, P ) (2.43)
For asymmetric problems

Ĝ(M, P ) = G(M, P )+α 1 G(M, s)G( s, P )+α 2 G(M, s)( ∂G( s, P ) ∂s x )+α 3 G(M, s)( ∂G( s, P ) ∂s y ) (2.44) A more general modified Green function is: Ĝ(M, P ) = G(M, P ) + N k=1 α k G(M, s k )G( s k , P ) (2.45)
Zhu and Lee [START_REF] Zhu | Removing the irregular frequencies in wave-body interactions[END_REF] tested the modified Green function method for bodies with two planes of symmetry with s at the origin. Their numerical results indicate that the irregular frequencies can be removed but the value α depends on the mode of motion. Also the range of α is very restricted if an iterative solver is used due to the poor conditioning of the linear system. For arbitrary bodies and arbitrary motions, the optimum location of the point s and the number of the points is far from being obvious.

The numerical results of removal of the irregular frequencies will be discussed in Chapter 6.

Chapter 3

Comparison of existing methods for the calculation of the infinite water depth free-surface Green function in the frequency domain for the wave-structure interaction problem

Résumé

Dans ce chapitre, les expressions mathématiques et les méthodes numériques pour la fonction de Green à surface libre en profondeur inifnie et dans le domaine fréquentiel sont examinées. Douze expressions différentes sont passées en revue et analysées. Toutes ces expressions sont des solutions mathématiques exactes pour la propagation d'ondes à partir d'une source pulsante située dans le domaine des fluides. Cependant, leur évaluation numérique est difficile. Des méthodes numériques dédiées ont été développées. Elles comprennent des développements en série, des polyn ômes, des interpolations dans une table, des développements multipolaires, des approximations avec des fonctions élémentaires, etc. Dans ce travail, quatre méthodes ont été implémentées: la méthode de Newman [START_REF] Newman | Algorithms for the free-surface Green function[END_REF], la méthode de Delhommeau [27], la méthode de Telste-Noblesse [START_REF] Telste | Numerical evaluation of the Green function of water-wave radiation and diffraction[END_REF] et la méthode de Wu et al. [START_REF] Wu | A global approximation to the Green function for diffraction radiation of water waves[END_REF]. On trouve que le temps de calcul moyen de la méthode de Newman est de 5.745 × 10 -7 . Il en vaut 5. 782 ×10 -8 pour la méthode de Delhommeau. Pour les méthodes de Telste-Noblesse et de Wu et al., elles sont respectivement 4, 642×10 -8 et 1.721×10 -7 . La précision est respectivement de 6D (6 décimales), 5D et 3D pour la méthode de Newman, la méthode de Telste-Noblesse et la méthode de Wu et al. Pour la méthode de Delhommeau, elle est de 3D sauf lorsque les coordonnées verticales sont proches de 0. La précision de la méthode de Delhommeau peut être augmentée de manière significative en affinant la discrétisation des variables d'espace pour les fonctions tabulées et en utilisant des méthodes d'interpolation plus élevées, au prix d'un temps de calcul accru.

In this chapter, the mathematical expressions and numerical methods for the freesurface Green function with zero speed in deep water and in the frequency domain are studied. Twelve different expressions are reviewed. Several existing algorithms are implemented and compared according to the CPU time and accuracies. This work is published in Applied Ocean Research and attached below. There are three errors needed to be corrected.

• R in Eq. ( 44) should be replaced by R 1 .

• z M in Eq. ( 15) should be replaced by z M + z P . The same error was made in Wehausen's book [START_REF] Wehausen | Surface waves[END_REF] and the correct expression can be found in Eq. (5.8b) in Noblesse's paper [START_REF] Noblesse | The Green function in the theory of radiation and diffraction of regular water waves by a body[END_REF].

• Equations ( 30) and ( 31) need to be corrected as:

I(X, Y ) = 1 X ∞ n=0 (-1) n X -2n n! 1 2 • 3 2 • 5 2 ... (2n -1) 2 I 2n (-Y ) (3.1) ∂I(X, Y ) ∂X = - 1 X 2 ∞ n=0 (-1) n X -2n n! 1 2 • 3 2 • 5 2 ... (2n -1) 2 I 2n (-Y ) - 2 X ∞ n=0 (-1) n n X -2n n!X 1 2 • 3 2 • 5 2 ... (2n -1) 2 I 2n (-Y ) (3.2)
• The CPU time of Wu's method is not correct because of wrong detective way.

The Fig. 7(d) and Table 5 should be replaced by the figure and table below. The corrected CPU time is similar as other's result [START_REF] Wu | Wave component in the green function for diffraction radiation of regular water waves[END_REF]. Although Wu's method is not the fastest one, it will benefit its advantage in parallel computation.

Introduction

The diffraction and radiation of water waves by floating bodies with zero mean forward speed and in deep water are phenomena of uttermost importance in ocean engineering. They must be taken into account in the design of marine structures such as oil and gas platforms or marine renewable energy converters. In the industry, boundary element methods (BEM)-based codes are usually used to calculate these effects and determine the pressure fields and wave-induced forces acting on the structures. BEM codes rely on the linear free surface potential flow theory which itself relies on the free-surface Green function and its derivatives.

Two distinct numerical problems must be overcome in the implementation of BEM codes [1]. Firstly, the discretization of the body surface by a large number of panels leads to the construction of a dense linear system of equations which must be solved by suitable algorithms. Secondly, the evaluation of the Green function and its derivatives is challenging because of their singular behaviour at the origin. Typically, the numerical complexity of BEM codes is proportional to O(N 2 ) or O (N 3 ) with N the number of unknowns. The Green function is computed O(N 2 ) times to set up the linear system. The linear system is solved either by an iterative method with O(N 2 ) complexity or by Gauss elimination with O(N 3 ) complexity. The complexity can be reduced down to O(Nlog N) using acceleration algorithms such as the precorrected Fast Fourier Transform(pFFT) [5]. Although the computational time for the evaluation of the Green function and its derivatives is not fully representative of the computational time needed by a BEM code to evaluate the diffraction and radication by a floating body, it is still considered as one of the challenges for efficient three-dimensional computation of hydrodynamic coefficients and forces.

In the frequency domain, this

= < = = = G x x y y z z z G gG z G O R R ( ) ( ) ( ) in 0 / 0 on 0 (1/ ) as M P M P M P z 2 2 (1)
where G ∞ is the Green function in deep water, δ is the Dirac function, R is the distance between the source point M(x M , y M , z M ) and the field point P(x P , y P , z P ). ω is the wave frequency and g is the gravity. In Eq.

(1), the first equation is the Laplace equation; the second equation derives from the linearized free-surface conditions and the last equation is the radiation condition at infinity. This boundary value problem was studied extensively during the 1940s and 1950s and exact mathematical expressions for its solution were obtained and those expressions were reviewed by Wehausen and Laitone [6]. Later, Ursell [7] developed an expression involving a series of spherical harmonics for a heaving semi-sphere. A modified Green function from Haskind's representation was then proposed by Kim [8]. It was implicitly provided by Havelock in [9]. This modified expression was also re-derived by Hearn [10]. Finally, an integral representation in terms of an exponential integral was introduced independently by Guevel [11], Martin [12] and Noblesse [13].

The numerical evaluation of these expressions (and the evaluation of their derivatives) is challenging because of the complexity of the involved mathematical expressions, the mathematical singularity in M = P and the associated computational time [1] [3]. This is an issue for BEM codes because the Green function (and its derivatives) must be evaluated many times to solve the wave-structure interaction problem for each considered frequency. Thus, with the rapid development of computing resources, the focus shifted in the 1980s from the derivation of mathematical expressions for the Green function towards the algorithms and numerical methods for its efficient numerical computation. Noblesse [13] [14] proposed two complementary near-field and farfield single integral representations in terms of the exponential integral. An asymptotic expansion and convergent ascending series are used to calculate the Green function for large and small distances between the source and field point, respectively. Two complementary Taylor series expansions are also provided when the non-dimensional spatial coordinates approach zero. In 1986, Telste and Noblesse [3] published a numerical code for the evaluation of the Green function and its derivatives based on Noblesse's previous study [14]. The computational domain was divided into five sub-domains.

In 1984, Newman [15] used the Romberg quadrature to evaluate the finite integral of the free-surface Green function with double precision accuracy. In the same year, he proposed a new series expansion [16]. It is computationally efficient at small and moderate radial distances between the source and field points. A year later, he [1] proposed an algorithm for the evaluation of Green function and its derivatives in infinite and constant finite depth. The computational domain is divided into several domains for which polynomial approximations are provided including the series expansion from his previous paper [16]. This algorithm is implemented as a standard subroutine named as "FINGREEN" in the boundary element method (BEM)-based code WAMIT [1,17].

In 1991, a new algorithm based on a polynomial series approximation was proposed by Chen et al. [18][19][20]. Double Chebyshev polynomials approximations with special functions are used to evaluate the Green function for infinite and finite water depth. Chen's algorithm is implemented in the BEM code HYDROSTAR. A similar algorithm with double Chebyshev polynomials is also proposed by Wang in [21].

To reduce the computational time, Delhommeau developed a technique involving four tabulated functions and Lagrange interpolations. This technique is implemented in the BEM code AQUADYN [2] and in the open source code NEMOH [22].

Other expressions and numerical methods are available. In a number of research works, the Green function is decomposed into three parts: a free-space singularity, a non-oscillatory local flow and waves. Ba et al. [23] and Ponizy et al. [24] provided approximation methods for the non-oscillatory local flow term based on the use of a coordinates-transformation and a function-transformation. Proper coordinates and function transformations reduce the problem of approximating singular functions for unbounded domains into that of approximating smoothly varying functions over finite domains. Linear table interpolation can thus be used. Peter and Meylan proposed an eigenfunction expansion representation [25]. In 2011, a semi-analytical method was developed with the Haskind-Havelock kernel calculated by a singularity subtractive technique [26]. The same year, the multipole expansion was extended for the infinite water depth free surface Green function [27].

Recently, Wu et al. [4] proposed a global approximation of the Green function and its derivatives based on Noblesse's paper [14,3]. A simple approximation involving elementary functions is given for the local flow component. It does not require dividing the computational domain into multiple sub-domains.

A completely different approach was proposed by Clément [28]. It is based on a second order differential equations for the Green function. However, a challenge with this method is its initialization.

It should be noted that the numerical errors associated to BEMbased codes stem from several sources as discussed in [4]. They include the approximation of the body surface by a large number of panels; the approximation of the variations (piecewise constant, linear, quadratic, or higher-order) of the densities of the singularities (source, dipole) distributions over the surface panels; the numerical integration of the Green function and its gradient over the panels; and the numerical evaluation of the Green function and its gradient. Thus, the effect of numerical errors in the numerical evaluation of the Green function and its derivatives for practical computations is difficult to estimate. Liang et al. [29] investigated the accuracy of linear and second order wave loads for a hemisphere and a freely Floating Production Storage and Offloading (FPSO) unit. They used the Wu et al.'s method. They showed that using this 3D accuracy method for the Green function and its gradient doesn't make much difference for the hydrodynamic coefficients and forces when compared to results obtained using highly accurate methods. Howevever, the conclusion may not be same for other practical cases.

In this paper, various expressions for the free-surface Green function are reviewed and presented. Several numerical methods have been proposed among which four methods were selected and compared with respect to accuracy and computational time. The two first ones are the Newman's and Telste-Noblesse's methods. They were selected because they are widely known in the industry and academia. The Newman's method is implemented in the industry standard BEM code WAMIT. The Delhommeau's method was selected because it is implemented in the open source code NEMOH. The fourth method is that of Wu et al.. It was selected because it provides a global approximation which is easy to implement and doesn't need subdivisions. Thus it may be more suitable for parallel computation.

Mathematical expressions for the free-surface Green function in deep water

In this study, the time factor of the complex potential is taken as e -iωt . The coordinates and variables are depicted in Fig. 1. The mean free surface level is located at the plane z = 0. The vertical axis z points upwards. The source point P(x P , y P , z P ) and the field point M(x M , y M , z M ) are both lying on or under the free surface (z P ≤ 0, z M ≤ 0). The image source point P′(x P , y P , -z P ) is the mirror of the source point P with respect to the mean free surface. The horizontal distance between the source point P and the field point M is denoted by r. The vertical distance between the image source point P′ and the field point M is -Z. The distance between the source and field points is denoted by R and the distance between the image source point and the field point is

= + R r Z 1 2
2 . The angle θ is defined by cosθ =-Z/R 1 . The relations between the coordinates are given by:

= + = + = + + r x x y y Z z z R x x y y z z ( ) ( ) ( ) ( ) (
)

M P M P M P M P M P M P 2 2 2 2 2
(2)

Expressions of the first type

We define the expressions of the Green function of the first type as the expressions that can be written as:

= + + G r Z R R g r Z i k e J k r 4 ( , , ) 1 1 ( , , ) 2 ( )
k Z 1 0 0 0 0 (3)
g -(r, Z, ω) is called the free-surface term whereas the first two terms are the Rankine source (1/R) and the image source (1/R 1 ) contributions, respectively. One can note that these two source distributions cancel each other when the field point is located on the mean free surface z = 0. k 0 is the wave number defined by k 0 = ω 2 /g and J 0 (k 0 r) is the zeroth order Bessel function of the first kind.

Expression 1

The classical expression for the free-surface term reads [28]:

= g r Z k k k e J ( , , ) 2PV (kr)dk 0 0 kZ 0 (4)
where PV∫ represents the Cauchy principal integral. For ω = 0 and ω → ∞, one can show that:

= + = g r Z g r Z ( , , 0) ( , , ) 0 R 2 1
(5)

Expression 2

The expression used in the BEM codes NEMOH and AQUADYN [2], [22] is:

= g r Z k d Re J ( , , ) 2 [ ( ) 1 ] 0 2 2 (6)
where ζ = k 0 (Z + ir cosθ) and

J(ζ) = e ζ [E 1 (ζ) + iπ]. E 1 (ζ)
is the complex exponential function. In this expression, the free-surface term only involves a finite integral of the complex exponential integral whereas it is an infinite integral in the classical expression (4). Therefore, a direct numerical integration may be used to evaluate it. However, it should be noted that the integration kernel tends to infinity with a logarithmic behavior when both r and Z tend to 0. Note that in the BEM codes NEMOH and AQUADYN, the evaluation of the Bessel function in the imaginary part of the Green function is replaced by the evaluation of the following expression:

= G r Z k e d Im Re [ 4 ( , , )] 2 0 2 2 (7)

Second type

We define the expressions of the Green function of the second type as:

= + + G r Z R g r Z i k e J k r 4 ( , , ) 1 ( , , ) 2 ( ) k Z 0 0 0 0 0 ( 8
)
By comparing to the expressions of the first type, one can see that the contribution of the image source is included in the free-surface term g 0 (r, Z, ω).

Expression 3

The free-surface term g 0 is given in Wehausen and Laitone's book [6], equation (13. 17″). It reads:

= + g r Z k k k k e J ( , , ) PV (kr)dk 0 0 0 0 kZ 0 (9)
It can be shown that it is equivalent to expression 1 by using the Lipschitz's integral:

= + = R r Z e J 1/
1/ (kr)dk . Notably, expression 3 is used in Newman's paper [1].

Expression 4

Another expression of g 0 can be found in [9], [25], [30],:

= + + g r Z k k kk K k k k e Y k r ( , , ) 2 (( )cos(kZ) 2 sin(kZ)) (kr) dk 2 ( ) k Z 0 0 2 0 2 0 0 2 0 2 0 0 0 0 ( 10
)
where K 0 (kr) is the modified Bessel function of the second kind and Y 0 (k 0 r) is the Bessel function of the second kind. This expression is known to be difficult to evaluate numerically due to the singular behavior of K 0 (kr). It has a logarithmic singular behavior close to k = 0 and a slowly decaying behavior of the integrand when r is small.

Third type

We define the expressions of the Green function of the third type as:

= + + + + G r Z R R g r Z i k e J k r 4 ( , , ) 1 1 ( , , ) 2 ( ) k Z 1 0 0 0 0 (11)
For the expressions of this type, it can be noted that the sum of the contributions of the Rankine and image source term satisfies the homogeneous Neumann boundary condition on the mean free surface

+ = n R n R ( / (1/ ) / (1/ ) 0) 1 .

Expression 5

Applying Lipschitz's integral into Eq. ( 4) or ( 9), the free-surface term can be written by:

= + g r Z k k k e J ( , , ) 2 PV 1 (kr)dk 0 0 0 kZ 0 (12)
Note that it is equivalent to equation (13.15) in Wehausen and Laitone's book [6]. When r = 0, the free-surface term given can be written:

= + g Z k e k Z E (0, , ) 2 ( )
k Z i 0 0 0 ( 13
)
where E i (x) is the exponential integral. When r = 0 and Z → 0, g + (r, Z, ω) has a logarithmic behavior.

Expression 6

In [9], [30], the following expression was given by Havelock:

= + + g r Z k k k K k k k e Y k r ( , , ) 4 ( cos(kZ) sin(kZ)) (kr) dk 2 ( ) k Z 0 0 0 0 2 0 2 0 0 0 0 (14)
This expression corresponds to equation 3 in Peter and Meylan's paper [25] and to equation (13. 17 ) in Wehausen's book [6]. One may note the similarity with equation (10).

Expression 7

Expression 7 was proposed by Haskind [31]. It corresponds to the equation (13. 17′) in Wehausen's Book [6]:

= + g r Z k e e R k e Y k r ( , , ) 2 dz 2 ( ) k z z k z k Z 0 1 0 0 0 M M 0 0 0 (15)

Expression 8

A modified Haskind Green function was studied by several researchers [8][9][10]32]. It is called the Haskind-Havelock representation in D'elía et al.'s paper [26]. It reads:

= + + + + g r Z k e k r Y k r e t r H ( , , ) ( ) ( ) 2 dt k Z Z k t 0 0 0 0 0 0 2 2 0 0 (16)
where H 0 (k 0 r) is the Struve function as defined in [33]. This expression is equivalent to Eq. (3b) in Newman's paper [1] and equation ( 12) in Liapis's paper [34] where it is called as a Havelock's finite integral.

It is well suited for evaluating the Green function when Z is small. When the source and field points are both located on the free-surface (z = 0), it can be simplified to:

= + + + g r k k r Y k r i k J k r H ( , 0, ) [ ( ) ( )] 2 ( ) 0 0 0 0 0 0 0 0 (17)

Expression 9

The Bessel function of the second kind in equation ( 16) is singular when kr = 0. Kim [8] introduced a regularized Bessel function N 0 to remove this singularity. This regularized Bessel function reads:

= N k r Y k r k r ( ) ( ) 1 2 ln( ) 0 0 0 0 0 (18)
Substituting (18) into equation ( 16) leads to:

= + + + + + g r Z k e k r N k r k R Z e t r H ( , , ) ( ) ( ) 2 ln[ ( )] 2 1 dt k Z Z k t 0 0 0 0 0 0 1 0 2 2 0 0 (19)

Expression 10

The near-field expression was introduced by Noblesse [14], equation (5.11). It reads:

= + + g r Z k e k r k e d H Re E ( , , ) 2 ( ) 4 (
)

k Z 0 0 0 0 0 /2 1 0 (20)
When the source point and the field point have the same horizontal coordinates (i.e. r = 0), the near-field expression simplifies to:

= = + + g r Z k e k Z i Re E ( 0, , ) 2 [ ( 0)] k Z 0 1 0 0 (21)
The near-field expression includes special functions and a finite integral which can be evaluated by direct numerical integration techniques.

Expression 11

The far-field expression is given by equation (5.21) in [14]. It reads:

= + + g r Z k e Y k r k e d Im E ( , , ) 2 ( ) 4 ( )sec k Z 0 0 0 0 0 /2 1 0 (22)
Similarly to the near-field expression given in Eq. ( 21), the far-field expression includes special functions and a finite integral which can be evaluated by direct numerical integration techniques.

Fourth type 2.4.1. Expression 12

An eigenfunction representation of the Green function was proposed by Peter and Meylan [25]:

= + + × + + { }{ } G r Z e H k r k k k K kz k k kz k k kz 4 ( , , ) ik 2 ( ) 1 dk (kr) cos( ) sin(kz ) cos( ) sin( ) k Z M M P P 0 0 (1) 0 2 0 2 2 0 2 0 0 0 0 (23)
H 0 (1) is the Hankel function of the first kind of first order. There is no singularity in the integrand. Numerical quadrature is used to evaluate the integral in their work.

Summary of the analytical expressions

Expressions 1 to 12 of free-surface Green function are summarized in Table 1. In practice, it is convenient to use non-dimensionalized coordinates rather than dimensional ones since it reduces the integral to a function of only two variables. Therefore, let us define the non-dimensional variables X = k 0 r, Y = k 0 Z, Y M = k 0 z M and Y P = k 0 z P . Expressions 1 to 12 can be rewritten as function of these non-dimensional variables given in Table 2.

Numerical methods for the evaluation of the Green function and its derivatives

In this section, the numerical methods for the evaluation of the Green function and its derivatives developed by Newman [1], Delhommeau [2], Telste and Noblesse [3] and Wu et al. [4] are reviewed.

Newman's method

Newman's method proposed in 1985 is based on the non-dimensional versions of ( 9) and (16). For the efficient evaluation of the integrals in these expressions, Newman divided the quadrant (X, Y) into six sub-domains. Corresponding sub-domains are depicted in Fig. 2. The derivatives of the Green function can be calculated by using the chain rule:

= + + + = + + + + + ( ) ( ) k k G X Y ik e J X k k G X Y ik e J X 4 ( ( , ) 2 ( )) 4 ( ( , ) 2 ( )) G r r R R X Y G Z Z R R Y Y 1 1 0 0 0 0 1 1 0 0 0 0 1 1 (24)
where G + (X, Y) is the non-dimensionalized free-surface term. According to Eq. ( 16), ∂G + (X, Y)/∂Y =-G + (X, Y) -2/R 1 . Therefore, only algorithms for G + (X, Y) and ∂G + (X, Y)/∂X are required to evaluate the Green function and its derivatives. For the sake of convenience, the finite integral in equation ( 16) is written as I(X, Y)

= + = + + I X Y e e t r e X t ( , ) dt dt k Z Z k t Y t Y 0 2 2 0 2 2 0 0 (25) 3.1.1. Sub-domain 1 Sub-domain 1 is defined by X > 8 and -Y > 20.
There, the non-dimensional version of equation ( 16) is used. The integral is approximated by Legendre polynomials:

= + I X Y n P d ( , ) ! ( ) n n n 0 4 ( 1)
(26)

+ + = + I X Y X n d d P n P X ( , ) ! XY ( ) ( 1) ( ) n n n n 0 4 3 ( 27
)
where P n is the Legendre polynomial and (d, θ) denote the polar coordinates such that X = d cosθ and -Y = d sinθ and d = k 0 R 1 , α = sinθ. According to [1], the truncation order of 4 leads to a 6 decimals accuracy in this sub-domain.

Sub-domain 2

For -X/Y < 0.5, Eq. ( 9) is used. The free-surface term G + (X, Y) can be rewritten as:

= + = = G X Y X n m Y e Y E ( , ) 2 ( /4) ( !) ( 1)! ( ) ( ) n n m n m Y i 0 2 2 1 2 (28) = + = = X G X Y n X X n m Y e Y E ( , ) 2 2 ( /4) ( !) ( 1)! ( ) ( ) n n m n m Y i 0 2 2 1 2 (29)
The calculation of E i (-Y) depends on the definition in [33]. The series are truncated at n = 9 to achieve a 6 decimals accuracy.

Sub-domain 3

For X > 3.7, -X/Y > 4, the integral is approximated by:

= = { } I X Y I Y ( , ) ( 1)
• • . .. ( ) X n n X n n n 1 0 ! 1 2 3 2 5 2 ( 1) 2 2 n 2 (30)
= = = { } { } I Y n I Y ( 1) • • . .. ( ) ( 1) • • . .. ( ) I X Y X X n n X n n n X n n X n X n n ( , ) 1 0
Y n Y t Y n 0 2 0 2 (32)
The integral I 2n is solution of the recursion relation:

= + I Y n Y n n I ( ) 2 ( ) 2 (2 1
)

n n n n 2 2 2 1 2 2 (33)
The series in equations ( 30) and ( 31) are truncated at n = 3 to achieve a 6 decimals accuracy.

Sub-domain 4

For 0 < X < 3.7, 0 < -Y < 2, G + (X, Y) is approximated using a series expansion [16]:

= + + + + + + + + = = G X Y e J X Y X Y X Y X X X Y X Y C X Y H ( , ) 2 [ ( )log{ / (1 / ) } ( ) ( )( ) ( ) ( ) ] Y X m n m n 0 2 2 1/2 2 0 2 0 2 2 1/2 2 2 1/2 0 1 mn 2 (34)
= + + + + + + + + + + + + + + + = = = = ( ) ( ) G X Y e J X Y X Y X J X Y X H X X Y H X X Y H X X Y X X Y C X Y X Y mC X Y ( , ) 2 [ ( )log{ / (1 / ) } ( ) 1 1
( ) ( )( ) ( )( ) ( )( ) ( ) ( )
( ) 2 ( ) ] X Y Y X Y X Y X Y X Y X X X m n m n m n m n 1 2 2 1/2 0 1/2 1 1/2 2 1 2 1 0 2 2 1/2 2 1 2 2 1/2 2 0 2 2 1/2 2 2 1/2 0 1 mn 2 2 2 1/2 0 1 mn 2 1 2 2 2 2 3 2 2 2 (35)
= + G g r Z 4 ( , , ) R R 1 1 1 + 2iπk 0 e k0 ZJ 0 (k 0 r) = g r Z e J ( , , ) 2PV (kr)dk k k k 0 0 kZ 0 Classical [28] = g r Z d Re J ( , , ) [ ( ) ] k 2 0 2 2 1 Guével [2], [22] -Delhommeau = + G g r Z 4 ( , , ) R 1 0
+ 2iπk 0 e k0 ZJ 0 (k 0 r)

= + g r Z e J ( , , ) PV (kr)dk k k k k 0 0 0 0 kZ 0 Wehausen [1], [6] -Laitone/1 = + + g r Z k k kk ( , , ) {( )cos(kZ) 2 sin(kZ)} dk K k k 0 2 0 2 0 2 0 0 (kr) 2 0 2
Havelock [30], [9],

-2πk 0 e k0 ZY 0 (k 0 r)

[25] = + + + G g r Z 4 ( , , ) R R 1 1 1 + 2iπk 0 e k0 ZJ 0 (k 0 r) = + g r Z k e J ( , , ) 2 PV (kr)dk k k 0 0 1 0 kZ 0 Wehausen [6]
-Laitone/2

= + + g r Z k k ( , , ) { cos(kZ) sin(kZ)} dk k K k k 4 0 0 0 0 (kr) 2 0 2
Havelock [30], [9],

-2πk 0 e k0 ZY 0 (k 0 r) [6], [25]

= + g r Z k e k e Y k r ( , , ) 2 dz 2 ( ) k z M z M e k z R k Z 0 0 0 1 0 0 0 0 Haskind [31], [6] = + + + + g r Z k e k r Y k r H ( , , ) ( ) ( ) dt k Z Z e k t t r 0 0 0 0 0 0 2 0 0 2 2 Modified [9], [8], [10],
Haskind [26], [32]

= + + + g r Z k e k r N k r k R Z H ( , , ) ( ) ( ) ln[ ( )] k Z 0 0 0 0 0 0 2 0 1 Kim [8] + + dt Z e k t t r 2 0 0 1 2 2 = + + g r Z k e k r e d H Re E ( , , ) 2 ( ) ( )
k Z k 0 0 0 0 4 0 0 /2 1 Near-field [14] = + + g r Z k e Y k r e d Im E ( , , ) 2 ( ) ( )sec k
Z k 0 0 0 0 4 0 0 /2 1 Far-field [14] -4πG ∞ = g ∑ = + + { } g r Z e H k r kz ( , , ) ( ) cos( ) sin(kz
) k Z M k k M ik0 2 0 0 (1) 0 1 2 0 0 Peter-Meylan [25] × + + kz kz K (cos( ) sin( )) (kr)dk k k k P k k P 2 2 0 2 0 0
where the coefficient C mn is defined by:

= + + = + + + C n n C C [( 1)( 1) !] ( ) n n n m n 0 1 mn 2 1 1, 2 (36)
The truncation of the double summation series is made both for m and n. According to Newman's paper [1], the 6 decimals accuracy is achieved with a total of 33 terms. In the present work, it was found that 43 terms are needed to achieve the 6 decimals accuracy for the derivatives of the Green function.

Sub-domain 5

In sub-domain 5, the finite integral in equation ( 16) is approximated by:

= + + I X Y X Y e X Y X Y X Y ( , ) 1 ( ) ( , ) Y 2 2 2 2 3/2 (37)
= + + + + + X Y ( , ) I X Y X X X Y e X X Y Y X Y X Y X ( , ) ( ) 3XY ( ) ( ) ( , ) Y 2 2 3/2 2 2 2 5/2 2 2 3/2 (38)
where X Y ( , ) is the slowly-varying residual functions in the given region.

X Y ( , ) and X Y X ( , )/ are approximated using Chebyshev polynomials [35]. To calculate the coefficients of the Chebyshev polynomials, values of I are obtained by direct numerical integration. The Double Chebyshev polynomial expansions are truncated to the desired accuracy by neglecting coefficients smaller than 10 -6 .

Sub-domain 6

The approximation for the remaining sub-domain is not defined in [1]. After numerical tests, the algorithm for sub-domain 3 with n = 7 is used.

Delhommeau's method

In 1989, Delhommeau introduced an algorithm based on interpolations of the Green function from a pre-calculated table [2]. The

Table 2

Expressions of the Green function as function of non dimensional coordinates.

Green function

Real part of the free-surface term Name Reference

= + G k G X Y 4 ( , ) R R 1 1 1 0 + 2iπk 0 e k0 ZJ 0 (k 0 r) = G X Y e J ( , ) 2PV (KX)dK K K 0 1 kY 0 Classical [28] = + + + + G X Y e Y i d Re E ( , ) [ { ( iX cos ) } ] Y Y 2 2 2 iX cos 1 1 iX cos Guével [2], [22] -Delhommeau = + G k G X Y 4 ( , ) R 1 0 0 + 2iπk 0 e k0 ZJ 0 (k 0 r) = + G X Y e J ( , ) PV (KX)dK K K 0 0 1 1 KY 0 Wehausen [1], [6] -Laitone/1 = + + G X Y K K ( , )
{( 1)cos(KY) 2 sin(KY)} dK

K K 0 2 0 2 0 (KX) 2 1
Havelock [30], [9],

-2πe Y Y 0 (X) [25] = + + + G k G X Y 4 ( , ) R R 1 1 1 0 + 2iπk 0 e k0 ZJ 0 (k 0 r) = + G X Y e J ( , ) 2PV (KX)dK K 0 1 1 KY 0 Wehausen [6] -Laitone/2 = + + G X Y K e Y X ( , ) (cos(KY) sin(KY)) dK 2 ( ) K K Y 4 0 0 (KX) 2 1 0
Havelock [30], [9], [6], [25] =

+ + G X Y Y e e Y X ( , , ) 2 dT 2 ( ) M Y M Y M e T X Y Y 2 2 0
Haskind [31], [6]

= + + + + G X Y e X Y X H ( , ) ( ) ( ) dT Y Y e T T X 0 0 2 0 2 2 Modified [9], [8], [10],
Haskind [26], [32]

= + + + + G X Y e X N X X Y Y H ( , ) [ ( ) ( ) ln[( )] Y 0 0 2 2 2 Kim [8] + + dT] Y e T T X 2 0 1 2 2 = + + + + G X Y e X e Y d H Re E ( , ) 2 ( )
( iX cos ) Y Y 0 4 0 /2 iX cos 1
Near-field [14]

= + + + + G X Y e Y X e Y d Im E ( , ) 2 ( )
( iX cos )sec Y Y 0 4 0 /2 iX cos 1 Far-field [14] -4πG ∞ = k 0 G ∑ = + + G X Y Y e H X Y Y ( , , ) ( ) (cos(
) sin( ))

M P i Y M K M 2 0 (1) 1 2 0 1 Peter-Meylan [25] × + + Y Y K X (cos( ) sin( )) ( )dK K K P K P 2 2 1 1 0
Fig. 2. The six sub-domains in Newman's method.

Delhommeau's method relies on equation (6). The derivatives are given by:

= + + + G r r R R k i d ik i e d Re J Re 1 1 2 cos [ ( ) 1 ] 2 cos 1 0 2 /2 /2 0 2 /2 /2 (39)
= + + + + G Z Z R R k d k R ik e d Re J Re 1 1 2 ( ) 2 1 0 2 /2 /2 0 1 0 2 /2 /2
(40)

Elementary functions

In Delhommeau's method, the expressions of the Green function and its derivatives in equations can be expressed as the function of these four elementary integrals:

= = = = = = D X Y i d d D X Y i e d e d Z X Y d Z X Y e d Re J Im J Re Im Re J Re ( , ) ( cos )[ ( ) ] cos [ ( ) ] ( , ) [ ( cos ) ] [ cos ] ( , ) [ ( ) ] ( , ) [ ] 1 /2 /2 1 /2 /2 1 2 /2 /2 /2 /2 1 /2 /2 2 /2 /2 (41)
The algorithm for the calculation of the complex exponential integral function can be found in appendix 3 in [36].

Tabulation method for near and moderate field

In Delhommeau's method, the four elementary integrals in equation ( 41) are interpolated from tabulated data at selected interpolation nodes (X i , Y j ). Lagrange polynomials of the fourth order are used for the interpolation. Direct numerical integration is used for evaluation of the elementary integrals at interpolations nodes. The integration interval [-π/2, π/2] is divided into 5001 points. For the numerical integration, the Simpson's rule is used. The by-default tabulation domain in AQUADYN 2.1 and NEMOH is X i ∈ [0, 100] and Y j ∈ [-16, -1.58 × 10 -6 ]. Since the integrals D 1 et Z 1 go to infinity with a logarithmic behavior when X = 0 and Y tends to 0 -, the interpolation nodes should be more dense close to zero. Thus, the interpolation nodes (X i , Y j ) are defined by:

= = = + = X i i i i 0 for 1
This discretization was found to be sufficient for simple bodies [2]. For complex geometries, refinement of the grid of interpolation nodes may be required.

Asymptotic formulas for the far field

For X > 100 or -Y > 16, the four elementary integrals are approximated by: (44)

+ + { } { } ( ) ( ) ( ) ( ) ( ) ( ) D e X X D e X X Z e X Z e X cos sin sin cos sin cos X k R Y X X Y X X Y k R Y X k R Y X

Telste-Noblesse's method

In the algorithm of Telste-Noblesse, three representations are used: The modified Haskind representation (equation ( 16)), the near-field representation (equation ( 21)) and the far-field representation (equation (22)). The derivatives of Green function are expressed as the following:

= + = + + + + + + ( ) ( ) { } k ik e J X k G X Y ik e J X 2 ( ) 2 ( , ) ( ) G r Z r r R R G X Y X Y G r Z Z Z R R k R Y ( , , ) 1 1 0 2 ( , ) 0 2 1 ( , , ) 1 1 0 2 1 0 2 0 1 1 0 1 (45)
Depending on the representation, the derivatives with respect to X of real part of the free-surface term can be written:

= + = = + + + + + + X Y X e X e X e e d Y X e e d H H Im E Re E [ ( ) ( )] 2 dt 2 [ ( )] [ ( ) ]cos 2 ( ) [ ( ) ]sec G X Y X Y Y t Y t X G X Y X Y G X Y X Y ( , ) 2 1 1 0 ( ) 1 ( , ) 2 1 4 0 /2 1 1 ( , ) 1
4 0 /2 1 1 2 2 2 3 (46)
In the algorithm of Telste-Noblesse, the computational domain is also split into sub-domains. The sub-domains are shown in Fig. 3.

Sub-domain 1

For large to moderate values of d and X ≠ 0, an asymptotic expansion is applied for the far-field representation. For 0 ≥ Y ≥ Y t (Y t being a user-defined variable set to -14.5 in the Telste-Noblesse's algorithm), the free-surface term is expanded as follows: The derivatives with respect to X are given by:

+ = + G X Y X e p d Y ( , ) 2 ( ) 2 (
+ + = + G X Y X X e X q d Y ( , ) 2 ( ) 2 (
)

Y n n n 1 0 3 ( 49
)
where the polynomials q n are defined by:

= = = + + q q q n q n n q n ( ) 1 ( ) 3 ( ) (2 1) ( ) ( 1)( 1) ( ) for 2 n n n 0 1 1 2 (50)
When Y t ≥ Y, the wave term and its derivatives are approximated as follows:

+ = + G X Y p d Y ( , ) 2
( ) for

n n n 0 1 (51) + = + G X Y X X q d Y ( , ) 2 ( )
for

n n n 0 3 (52)
The truncation order of the series is n = 10 outside the ellipse X 2 / 12 2 + Y 2 /15 2 = 1 and n = 5 outside the ellipse X 2 /16 2 + Y 2 /20 2 = 1. The absolute errors are reported to be smaller than 10 -6 in [3].

Sub-domain 2

The sub-domain 2 corresponds to 0 ≤ X ≤ 1.2 -0.15Y, -2 ≤ Y ≤ 0. A convergent ascending series is used. The ascending series is based on the near-field representation. The free-surface term is expressed as:

= + + + + G X Y d Y F Y J X S X Y e T X Y ( , ) 2[{ ln( ) ( )} ( ) 2 ( , )] 2 ( , ) Y 0 0 0 0 (53)
where F 0 (Y), S 0 (X, Y), T 0 (X, Y) are defined as follows:

= = = = = = = { } ( ) ( ) F Y S X Y T X Y F Y ( ) ln 2 ( , ) ( 1)
( , ) ( )

n Y n n n k n k n k k n k X n n n X n 0 1 ( ) • ! 0 1 0 ( 2 ) ! ( ) ! 2 0 1 4 n k 2 2 2 (54)
In these last equations, γ = 0.5772. . . is the Euler's constant. n′ is the largest integer which does not exceed (n -1)/2 and σ = X/(d -Y). F n (Y) is defined as follows:

= + = = F Y n n m k Y k ( ) 1 ! ! 1 ! n k m n k 0 1 2 (55)
A corresponding ascending series for the derivatives with respect to X is given by:

= + + X d Y F Y J X S X Y e T X Y 2 2{ ln( ) ( )} ( ) 2 ( , ) 2 ( , ) G X Y X J X d d Y Y ( , ) ( ) ( ) 0 1 1 1 0 (56)
with

= + = = = = ( ) ( ) S X Y T X Y nF Y ( , ) ( 1)
( , ) ( )

X d n k n k n k n n k Y d k n k X n X n n X n 1 2 0 2 / ( 2 ) / ! ( ) ! 2 1 1 2 1 4 1 2 (57)
5 terms are required for the series T 0 and T 1 and 10 terms are required for the series S 0 and S 1 to achieve 10 -6 accuracy. For X < 0.55 -2Y, less terms are needed.

Sub-domain 3

A Taylor series expansion of the near-field integral representation is used in the neighborhood of the vertical axis X = 0:

= + + + = G X Y J X e Y i r Y X Y Re E ( , ) 2 ( ){ [ ( 0)]} 2 ( ) 4 Y n n n 0 1 1 2 2 (58)
r n (Y) is given by:

= = r Y n n k Y ( ) [1/( !) ] ( 2 1
) !( )

n k n k 2 0 2 1 (59)
It is also solution of the recursion relation:

= + r Y n n Y n Y n r Y n ( ) (2 2) ! (2 1 )( !) ( / ) ( ), 2 n n 2 2 2 1 (60)
Similarly, a Taylor series expansion can be obtained for the derivatives:

= + + = G X Y X J X e Y i X Y n r Y X Y Re E ( , ) 2 ( ){ [ ( 0)]} • ( ) 4 Y n n n 1 1 2 1 2 2 1 (61)
Rational approximations are used for the calculation of

+ e Y i Re E { [ ( 0)]} Y 1
as in [3]. To achieve the accuracy of 10 -6 , five, ten and fifteen terms are respectively retained in the series depending on whether X ≤ min(-0.3Y, 1.8), X ≤ min(-0.5Y, 5.4) or X ≤ min (-0.7Y, 9.1) with 2 ≤ -Y ≤ 15.

Sub-domain 4

In the neighborhood of the horizontal plane Y = 0, a Taylor series expansion of the Haskind integral representation is used:

= + + = G X Y X Y X e s X Y X H ( , ) { ( ) ( )} 2 ( ) Y n n n 0 0 1 (62)
The polynomials s n are solution of the following recursion relation:

= = = + + + s X s X s X n s X Xs X C n n ( ) 1
( ) ( ), for 1
( ) { ( ) }/(2 1), for 1 n X n n n n n 1 2 2 2 1 2 1 2 (63)
with

= C n n ( 1) {1•3•5 (2 1)}/{2•4•6 (2 )} n n ( 64
)
For the derivative, the Taylor series expansion is:

= + = { } G X Y X X Y X e X t X Y X H ( , ) 2 ( ) ( ) 2 ( )
Y n n n 1 1 1 (65)
with

= = = + + + t X t X Xt X n n t X Xt X n C n ( ) 1 ( ) ( )/2 , for 1 ( ) ( )/( 2 1)
, for 1

n n n n n 1 2 2 1 2 1 2 (66)
In the region 1.2 ≤ X ≤ 12, the truncation order is 10 if Y ≥ max (-0.4X + 0.15, -1.5) or 20 if Y ≥ -0.4X.

Sub-domain 5

The sub-domain 5 corresponds to the region where 1 < d < 10 along with 0.4 < -Y/X < 10/7. In this sub-domain, the modified Haskind representation is used. The free-surface term and its derivatives are given by:

= + + + G X Y X Y X e F X I X e H ( , ) { ( ) ( )} 2[ ( ) ( )] Y Y 0 0 0 0 cX (67)
= + + + + { } G X Y X X Y X e F X I X e X H ( , ) 2 ( ) ( ) 2[ ( ) ( )] / Y Y 1 1 1 1 cX ( 68
)
where c is a special constant and F 0 (X), F 1 (X), I 0 (X, Y) and I 1 (X, Y) are defined by the integrals:

= + F X e d ( ) (1 ) c X c 0 0 ( ) 2 1/2 (69) = + F X e d ( ) (1 )
c X c 1 0 ( ) 2 3/2
(70)

= + I X Y e d ( , ) (1 ) c Y X X c 0 / ( ) 2 1/2 (71) = + I X Y e d ( , ) (1 ) c
Y X X c 1 / ( ) 2 3/2 (72)
These integrals involve the integrands + (1 ) 2 1/2 and + (1 ) 2 3/2 . Those terms can be approximated by polynomials. Then the integrals can be integrated analytically which leads to series expansion forms of the solution. The special constant c is set to 0.4, 0.65, 0.85 and 1.1 for the regions 0.4 < -Y/X ≤ 0.65, 0.65 < -Y/X ≤ 0.85, 0.85 < -Y/ X ≤ 1.1 and 1.1 < -Y/X ≤ 10/7, respectively.

Wu et al.'s method

A simple approximation of the Green function was proposed by Wu et al. [4] which valid in the whole domain. It derives from the Noblesse's near-field representation [13]. The free-surface term is split into a local flow term L(X, Y) and an out-going free-surface term W(X, Y):

= + G X Y L X Y W X Y ( , ) { ( , ) ( , )} 0 (73) with = = L X Y e d W X Y X e Re E H ( , ) [ ( ) ] ( , ) 2 ( )
k R Y 1 4 0 /2 1 0 0 1 (74)
where ζ has the same definition as in section 2.1.2. The derivatives of Green function are:

= + = + ( ) ( ) k k W ( ) ( ) G r r R L X W X G Z Z R L Y 1 0 1 0 (75)
with

= + = = + X Y L X Y X Y X e X Y L X Y H ( , ) ( , ) ( , ) 2 [2/
( )] ( , ) * ( , ) L Y Y k R k R W X Y L X X k R 1 1 0 3 1 3 0 1 0 3 1 3 (76)
= + + + L k R L k R L 1 2 1 2 (1 )
a P R 0 1 0 3 1 3 3 (78)
with

= + + = + + L e k R k R Y L A B D (log 2 ) (1 ) (1 )
P Y k R Y R C 2 0 2 1 2 0 2 1 2 1 6 (1 ) 0 1 (79)
where

α =-Z/R 1 , β = r/R 1 , ρ = k 0 R 1 /(1 + k 0 R 1 )
The other term L * a is:

= + + L L k R L L * 2 * 1 4 * 2 (1 ) * a P Q R 0 3 1 3 3 (81)
with

= + = + = + + + L e k R X L e L A B C * 2 2 * (1 ) (1 )
* * (1 ) * (1 ) (1 2 ) * P X k R Y Y Q k R k R k R R 0 1 1 0 1 0 1 0 1 0 3 1 3 (82)
The polynomials A * , B * and C * are defined by: (83)

= + + + + = + + + + + = + + A B C * 2.

Summary of numerical schemes

The mathematical formulas for the calculation of the free-surface term of the Green function and its derivative for the four methods are summarized in Tables 3 and4 .

Comparison of the four methods

We implemented the algorithms of Newman, Delhommeau, Telste and Noblesse and Wu et al. For the Telste-Noblesse's method, the source code is given in their paper [3]. However, it should be noted that we changed the precision from single float to double float. For the Delhommeau's method, the code follows the NEMOH code (which is open source) with a few minor modifications in the tabulated file. For the Newman's and Wu et al.'s methods, we developed the source code based on their papers. It should be noted that the computational time may vary depending on the coding technique.

All the calculations were performed on a PC with intel(R) Core(TM) i7-6700 CPU @ 3.40 GHZ. The codes were written in FORTRAN90 and compiled with intel fortran 14.0.1. The evaluation of the Struve function is based on the algorithm of Newman [38]. For the Bessel functions, the functions available in the FORTRAN compiler were used.

The accuracy and computational time for the four methods were compared. To facilitate the comparisons, only results for G -and its derivative ∂G -/∂X are shown. For the methods that provide G + and G 0 , the image source contribution were added to obtain G -by using the relations

= + + G X Y G X Y X Y ( , ) ( , ) 1/ 0 2 2 , = G X Y ( , ) + + + G X Y X Y ( , ) 2/ 2 2
. For the comparisons, the non-dimensionalized variables (X, Y) are uniformly discretized over the domain X ∈ [0.001, 22], Y ∈ [0.001, 15]. This domain is expected to cover most of the practical cases. The case of the source and field points located at the same position is not considered since it corresponds to the singularities. Analytical expressions are available for Y = 0 (equation ( 16))and for X = 0 (equation ( 13)). Reference values for the real part of the non-dimensionalized free-surface term of the Green function (G -) and its derivative (G -X = ∂G -/∂X) are shown in Fig. 4. They were obtained by direct numerical integration using an adaptive Gaussian 3-point quadrature. The absolute error for the reference values is less than 10 -8 .

The absolute error for the Green function between the values computed from the direct integration and four methods are plotted in Fig. 5. As expected, a 6 decimals accuracy is obtained with the Newman's method over the whole computational domain, see Fig. 5(a). The error is greater in the sub-domain 1, which corresponds to the far-field, than in the other sub-domains. The absolute errors computed with the algorithm of Delhommeau are shown in Fig. 5(b). The absolute error is greater than the other methods especially in the near field. The error for the Telste-Noblesse's method is shown in Fig. 5(c). The error is the greatest in the middle of the domain. Nevertheless, the method achieves a 5 decimals accuracy. The error for the algorithm of Wu et al. is shown in Fig. 5(d). The accuracy is 3 decimals.

The absolute error for ∂G -/∂X for the four methods are shown in Fig. 6. The results are similar to that for G -: the accuracy is 6D, 5D, 3D

Table 3

The expressions of the free-surface term of the Green function and their derivatives which are useful for numerical calculation (1) Expression number Expressions:

G + (X, Y) and + X Y ( , ) G X
Region of fast convergence Reference 8.1

+ + = + G e X Y X n P d H [ ( ) ( )] 2 ! ( ) Y n n n 0 0 0 4 ( 1)
Large X and Large Y

[1] + e X Y X H [ ( ) ( )] G X Y 2 1 1 + + + = + { } n P X 2 ( 1) ( ) n n d n d Pn n 0 4 ! 3 XY ( ) 3.1 = + = = { } G e Y E 2 ( ) n X n n m n m Y m Y i 0 ( 2 / 4) ( !) 2 1 2 ( 1) ! ( )
Small X [1], [14],

= + = = { } e Y E 2 ( ) G X n n X X n n m n m Y m Y i 0 2 ( 2 / 4) ( !) 2 1 2 ( 1) ! ( ) [37] 8.2 = + + = { } G e X Y X I Y H [ ( ) ( )] ( 1) • • . .. ( ) Y X n n X n n n n 0 0 2 0 2 ! 1 2 3 2 5 2 ( 1) 2 2 Large X/Y [1] = + + = { } e X Y X H [ ( ) ( )] ( 1) • • . .. G X Y X n n X n n n 2 1 1 2 2 0 2 ! 1 2 3 2 5 2 ( 1) 2 × + = { } I Y n I Y ( ) ( 1) • • . .. ( ) n X n n X n n X n n 2 4 0 2 ! 1 2 3 2 5 2 ( 1) 2 2 8.3 = + + + + G e J X Y X Y X Y X 2 [ ( )log{ / (1 / ) } ( ) Y 0 2 2 1/2 2 0
Uniformly convergent [1], [16]

+ + + + = = X X Y X Y C X Y H ( )( ) ( ) ( ) ] X m n m n 2 0 2 2 1/2 2 2 1/2 0 1 mn 2 = + + + e J X Y X Y X 2 [ ( )log{ / (1 / ) } G X Y 1 2 2 1/2 + + + + J X ( ) 1 1 Y X Y X Y X Y X Y X 0 2 2 1/2 1 2 2 3 2 2 1/2 + + Y X X X Y X X Y H H ( ) ( )( ) ( )( ) X X 2 1 2 1 2 0 2 2 1/2 2 1 2 2 1/2 + + + + = = X X Y X X Y C X Y H ( )( ) ( ) ( ) m n m n 2 0 2 2 1/2 2 2 1/2 0 1 mn 2 + + = = X Y mC X Y ( ) 2 ( ) ] m n m n 2 2 1/2 0 1 mn 2 1 8.4 = + + + + G e X Y X X Y H [ ( ) ( )] 2 ( , ) Y X Y e Y X Y X Y 0 0 1 2 2 ( 2 2 ) 3/2 Moderate X, Y [1], [15], [35] = + + + e X Y X H [ ( ) ( )] 2 G X Y X X Y e Y X 2 1 1 ( 2 2 ) 3/2 2 + + + X Y ( , ) X Y Y X Y X Y X 3XY ( 2 2 ) 5/2 ( 2 2 ) 3/2 ( , )

Table 4

The expressions of the free-surface term of the Green function and their derivatives which are useful for numerical calculation (2)

Expression number Expressions:

G + (X, Y) and + X Y ( , ) G X
Region of fast convergence Reference 11.1

+ = + G Y X e 2 ( ) 2
Y n Pn d n 0 0 ( ) 1
Large X and Large Y

[3], [14]

+ + = + Y X e X 2 ( ) 2 G X Y n Qn d n 1 0 ( ) 3
10.1

= + + + + G d Y F Y J X S X Y e T X Y 2[{ ln( ) ( )} ( ) 2 ( , )] 2 ( , ) Y 0 0 0 0 Small X and Small Y [3] = + + X d Y F Y J X 2 2{ ln( ) ( )} ( ) G X J X d d Y 0 ( ) ( ) 0 1 -2S 1 (X, Y)e Y -2T 1 (X, Y) 10.2 = + + + = G J X e Y i r Y Re E 2 ( ){ [ ( 0)]} 2 ( ) Y n n X Y n 0 1 1 2 4 2 Small X [3], [14] = + + = J X e Y i n r Y Re E 2 ( ){ [ ( 0)]} • ( ) G X Y X Y n n X Y n 1 1 2 1 2 4 2 1 8.5 = + + = ( ) G X Y X e s X H { ( ) ( )} 2 ( ) Y n n Y X n 0 0 1 Small Y [3], [14] = + = ( ) { } X Y X e t X H ( ) ( ) ( ) G X Y X n n Y X n 2 1 1 2 1 8.6 = + + + G X Y X Y X e F X I X e H ( , ) { ( ) ( )} 2[ ( ) ( )] Y Y 0 0 0 0 cX Moderate X, Y [3] = + + + + { } G X Y X Y X e F X I X e X H ( , ) ( ) ( ) 2[ ( ) ( )] / X Y Y 2 1 1 1 1 cX
for the algorithms of Newman, Telste-Noblesse's and Wu et al., respectively. For Delhommeau's method, the accuracy is 3D almost everywhere in the considered domain except when (X, Y) is close to (0, 0). However, the value of the Green function is large in that situation which makes the relative error small, in the order of 1 to 2%.

Note that for some cases of (X, Y), the accuracy of the derivative can be less than the function as can be seen by close inspection of Fig. 5 and Fig. 6, depending on the method.

The computational time for calculating the Green functions and its derivative over the set of discretized (X, Y) values for the different Fig. 4. The free-surface term of the Green function G -and its derivative G -X using the direct integral method. algorithms is shown in Table 5. As one can expect, the computational time is the greatest for the methods that are the most accurate (direct numerical integration and Newman's method). The direct numerical integration with tolerance of 10 -8 needs 20 times more computational time than the Newman's method. The Delhommeau and Telste-Noblesse's methods are an order of magnitude faster than the Newman's method. However, the accuracy is respectively about three orders of magnitude and slightly less than the Newman's method. With respect to computational time, the Wu et al.'s method is the fastest one.

Overall, it appears that the best method is the Newman's method and Telste-Noblesse's method with respect to accuracy. With respect to computational time, the Wu et al.'s method is best.

As in Newman's and Telste-Noblesse's method, the evaluations algorithms are different in each sub-domains. The computational time for each (X, Y) computational points as function of the method is shown in Fig. 7. It can be seen that the computational time depends not only on the method but also on the sub-domains. The local computational time for the Newman's method is shown in Fig. 7(a). The average computational time for each sub-domain is shown in Table 6. In the sub-domain 2, the computational time increases with increasing Y. The computational time for the Delhommeau's method is shown in Fig. 7(b). Since it is based on the fourth order Lagrange polynomial interpolation, the computational time does not vary much over the interpolation domain. For large X and Y, the asymptotic formulation is used to calculate the Green function which explains why the computational time for Y ≥ 16 is much shorter than for the rest of the domain.

The compuational time for the Telste-Noblesse's method in shown in Fig. 7(c). The computational time varies between O(10 -7 ) and O (10 -8 ). The average computational time for each sub-domain is shown in Table 7. For any given X and Y, the computational time appears to be less than that of the Newman's method.

The computational time for the Wu et al.'s method is shown in Fig. 7(d). It is shown that the computational time suddenly increases for X ≥ 3. It is because despite the fact that the flow term is approximated with a simple polynomials with the elementary functions, the Wu et al.'s method still requires a subroutine to calculate the Struve function (H 0 (k 0 r)). This subroutine is based on the algorithm by Newman [39], which uses different approaches depending on whether X is greater or

Table 5

Computational time for calculating the Green function and its derivatives a function of the different algorithms. smaller than 3. Nevertheless, a key feature of the Wu et al.'s method is that it is not based on subdivisions of the flow domain. Thus, it may be more suited for parallel computations compared to the other methods.

Conclusions

In this paper, we reviewed the available mathematical expressions for the deep water free-surface Green function in frequency domain and four different numerical methods which were developed for its numerical evaluation. The four methods are the Newman's method, the Delhommeau's method, the Telste-Noblesse's method and the Wu et al.'s method. The computational time and accuracy with each method are compared. The Newman's method is the most accurate providing a 6 decimal accuracy. However, it is also the slowest. The Telste and

Table 6

Averaged computation time of Newman's algorithm. Noblesse's method is an order of magnitude faster than the Newman's method, but it is slightly less accurate. The method by Delhommeau lie in-between the Newman and Teslte-Noblesse's method with respect to computational time but it is less accurate than the Telste-Noblesse's method. This makes this latter method preferrable to the Delhommeau's method. The Wu et al's method is also less accurate than the Newman and Teslte-Noblesse's method, but it is also the fastest method and the simplest with respect to numerical implementation (as they do not use different expressions depending on the position in the computational domain). Finally, a limitation of this work is that only the case of deep water was considered despite the case of finite depth is also very important for practical applications. This is left for future work.

Chapter 4 Ordinary differential equations for the infinite water depth free-surface Green functions and its derivatives Résumé

Dans le chapitre précédent, les fonctions de Green à surface libre dans le domaine fréquentiel sont toutes exprimées sous forme intégrale. Dans ce chapitre, nous montrons que la fonction de Green et ses dérivées sont solutions d'équations différentielles par rapport au temps ou à la fréquence. Par ailleurs, de nouvelles équations différentielles ordinaires sont établies pour la fonction de Green par rapport à la variable d'espace, à la fois dans le domaine fréquentiel et dans le domaine temporel.

In the previous chapter, the free-surface Green functions in the frequency domain are all expressed in integral forms. In this chapter, we derive several ordinary differential equations (ODEs) for the free-surface Green function and its derivatives in the frequency domain. This part of work are extension of Clement's work [24], [23]. The ODEs with respect to time and frequency are derived in a new way. The ODEs with respect to spatial variables in the time domain and frequency domain are firstly proposed in this thesis.

Time Domain Green function

Boundary value problem in time domain

The time-domain Green function of the free-surface hydrodynamic problem is defined as the fundamental solution of the following boundary value problem. The variables are nondimensional in terms of reference length [L] and reference time [L/g].

• Continuity or Laplace equation

∆ M F ∞ (r, Z; t) = δ(|P M |)δ(t -t 0 ) ∀z M ≤ 0; ∀t ≥ 0 (4.1)
• Free surface conditions

∂ 2 F ∞ (r, Z; t) ∂t 2 + ∂F ∞ (r, Z; t) ∂z M = 0 z M = 0; ∀t ≥ 0 (4.2)
• Condition at infinity

∇ M F ∞ (r, Z; t) → 0 z M → -∞; ∀t ≥ 0 (4.3) ∇ M F ∞ (r, Z; t) → 0 r → ∞; ∀t ≥ 0 (4.4)
t 0 is assumed to be 0. The solution to this initial-boundary value problem has been derived by Haskind [38] and Brard [10] for deep water. It was generalized by Finkelstein [29] for the finite water depth cases. According to these references, a mathematical expression of time-domain Green function in deep water reads,

F ∞ (r, Z; t) = - 1 4π [δ(t)F 0 (r, Z) + H(t)F (r, Z; t)] (4.5)
with

F 0 (r, Z) = 1 R - 1 R 1 (4.6)
and

F (r, Z; t) = 2 ∞ 0 J 0 (Kr)e KZ √ Ksin( √ Kt)dK (4.7)
H(t) is the Heaviside unit step function, and δ(t) the Dirac distribution and J 0 is the Bessel function of the first kind of order 0. F 0 (r, Z) is often referred to the instantaneous or impulsive part of the Green function while F (r, Z; t) is called the memory part.

The Green function in time domain as function of the natural variables can be obtained from the Green function in time domain with initial variables by a change of variables. It reduces the number of effective variables from three to two. Jami [START_REF] Jami | Etude théorique et numérique de phénomènes transitoires en hydrodynamique navale[END_REF] showed that the memory part of the Green function can be expressed as a function of two variables (µ, τ ).

F (r, Z; t) = 2R -3/2 1 F (µ, τ ) (4.8)
with

F (µ, τ ) = ∞ 0 J 0 (λ 1 -µ 2 )e -λµ √ λ sin( √ λτ )dλ (4.9)
by using the change of variables (r, Z, t) → (R 1 , µ, τ ).

         r Z t ⇒            R 1 = (r 2 + Z 2 ) 1/2 µ = - Z (r 2 + Z 2 ) 1/2 τ = t (r 2 + Z 2 ) 1/4 (4.10)          R 1 µ τ ⇒          r = R 1 1 -µ 2 Z = -µR 1 t = τ R 1 (4.11)

Time-domain differential equations for the Green function

Using Eq. (4.7), one can show that the time-domain Green function is one solution of three differential equations, ∀Z ≤ 0 and ∀t > 0.

∂ 2 F ∂Z 2 + 1 r ∂F ∂r + ∂ 2 F ∂r 2 = 0 (4.12a) ∂ 2 F ∂t 2 + ∂F ∂Z = 0 (4.12b) ∂ 4 F ∂t 4 - ∂ 2 F ∂Z 2 = 0 (4.12c)
Eq. (4.12a) expresses the Laplace equation (4.1). Eq. (4.12b) is similar to (4.2) except that it is valid in the whole fluid domain, not only on the free-surface (z = 0).

Using the change of variables defined by Eqs. (4.10) and (4.11), we define the Jacobian Matrix (see [24])

∂r ∂R 1 ∂r ∂µ ∂r ∂τ ∂Z ∂R 1 ∂Z ∂µ ∂Z ∂τ ∂t ∂R 1 ∂t ∂µ ∂t ∂τ = 1 -µ -R 1 µ 1 -µ 2 0 -µ -R 1 0 τ 2 √ R 1 0 R 1 = -R 1 R 1 1 -µ 2 (4.13)
It vanishes when and only when the source point and the field point are located in the same place. It is a trival case which is a priori excluded from the domain of the present study.

The change of variables gives, at the first level of derivative

∂ ∂r = 1 -µ 2 ∂ ∂R 1 - µ 1 -µ 2 R 1 ∂ ∂µ - τ 1 -µ 2 2R 1 ∂ ∂τ ∂ ∂Z = -µ ∂ ∂R 1 - 1 -µ 2 R 1 ∂ ∂µ + 1 2 µτ R 1 ∂ ∂τ ∂ ∂t = 1 √ R 1 ∂ ∂τ (4.14)
and for the second derivatives (diagonal terms):

∂ 2 ∂r 2 = µ 2 R 1 ∂ ∂R 1 + µ(2 -3µ 2 ) R 2 1 ∂ ∂µ + τ (3 -5µ 2 ) 4R 2 1 ∂ ∂τ + (1 -µ 2 )[ ∂ 2 ∂R 2 1 - 2µ R 1 ∂ 2 ∂µ∂R 1 - τ R 1 ∂ 2 ∂τ ∂R 1 ] + (1 -µ 2 ) R 2 1 [µ 2 ∂ 2 ∂µ 2 + µτ ∂ 2 ∂τ ∂µ + τ 2 4 ∂ 2 ∂τ 2 ] ∂ 2 ∂Z 2 = (1 -µ 2 ) R 1 ∂ ∂R 1 -3µ (1 -µ 2 ) R 2 1 ∂ ∂µ + τ 4 (5µ 2 -2) R 2 1 ∂ ∂τ + µ 2 ∂ 2 ∂R 2 1 + 2µ (1 -µ 2 ) R 1 ∂ 2 ∂µ∂R 1 - µ 2 τ R 1 ∂ 2 ∂τ ∂R 1 + [ 1 -µ 2 R 1 ] 2 ∂ 2 ∂µ 2 - µτ R 1 (1 -µ 2 ) R 1 ∂ 2 ∂τ ∂µ + µ 2 τ 2 4R 2 1 ∂ 2 ∂τ 2 ∂ 2 ∂t 2 = 1 R 1 ∂ 2 ∂τ 2 (4.15)
With the change of variables, Eq. (4.12) can be transformed into three similar partial equations. The first step of the derivation is to perform the change of variables in Eq. (4.10), then using Eqs. (4.14), (4.15) together with a change of the dependent variable F → F through Eq. (4.8). Switching to the indicative notation for the derivatives, we get:

τ 2
Fττ + 5τ Fτ + 4(1 -µ 2 ) Fµµ -8µ Fµ + 3 F = 0 (4.16a)

2 Fττ + µτ Fτ -2(1 -µ 2 ) Fµ + 3µ F = 0 (4.16b) 2 Fττττ - µ 2 τ 2 2 Fττ + τ (1 - 11 2 µ 2 ) Fτ + 2µτ (1 -µ 2 ) Fτµ -2(1 -µ 2 ) 2 Fµµ + 12µ(1 -µ 2 ) Fµ + (3 - 21 2 µ 2 ) F = 0 (4.16c)

Ordinary differential equations in time domain

A fourth order ODE for F with respect to the time variable τ

From the set of partial differential equations 4.16 one can easily recover the fourth order ordinary differential equation fot the time-domain Green function uncovered by Clément in 1998 [24].

Eq. (4.16b) can be rewritten as follows:

F 2k (µ, 0) = 0, F 2k+1 (µ, 0) = (-1) k (k + 1)!P k+1 (µ); k = 0, 1, ... (4.24)
where P k+1 denotes the Legendre polynomial of order k + 1.

A fourth order ODE for F with respect to the time variable t

Returning to the initial variables (r, Z; t), one can show from (4.23) that the memory part of the Green function F is solution of the ordinary differential equation:

(r 2 + Z 2 )F tttt -ZtF ttt + ( 1 4 t 2 -4Z)F tt + 7 4 tF t + 9 4 F = 0 (4.25)
where the subscript t means the time derivatives. The corresponding initial conditions are,

F (r, Z; 0) = 0 F t (r, Z; 0) = -2 Z (r 2 + Z 2 ) 3/2 F tt (r, Z; 0) = 0 F ttt (r, Z; 0) = 2 r 2 -2Z 2 (r 2 + Z 2 ) 5/2 (4.26)
The ordinary differential equations for F and F and corresponding initial conditions were first obtained by [24] by introducing a general fourth order ODE.

A new fourth order ODE for F with respect to the space variable µ

Re-starting from Eq. (4.16), we will now eliminate the time derivatives to keep only the spatial derivatives with respect to µ. Multiplying Eq. (4.16a)by µ and (4.16b) by -5 and summing leads to:

(µτ 2 -10) Fττ = -4µ(1 -µ 2 ) Fµµ -2(5 -9µ 2 ) Fµ + 12µ F (4.27)
Multiplying now Eq. (4.16a) by 2 and (4.16b) by -τ 2 and summing gives

-τ (µτ 2 -10) Fτ = -8(1 -µ 2 ) Fµµ + 2(µ 2 τ 2 + 8µ -τ 2 ) Fµ -3(2 -µτ 2 ) F (4.28)
By differentiating these two equations with respect to µ, one can obtain expressions Fµτ , Fµµτ and Fµττ as differential forms of successive derivatives of F with respect to the single variable µ: F , Fµ , Fµµ , Fµµµ , Fµµµµ .

Then, let us substitute these expressions into the time variable ODE Eq. (4.23) and express the fourth derivative using Eq. (4.19), we obtain:

c 4 Fµµµµ + c 3 Fµµµ + c 2 Fµµ + c 1 Fµ + c 0 F = 0 (4.29)
with

c 0 = -12µ(µ 2 + 1)τ 6 + 24(14µ 2 + 9)τ 4 -2880µτ 2 -2400 c 1 =µ(µ 2 + 1)τ 8 -2(39µ 4 + 8µ 2 + 9)τ 6 + 12µ(167µ 2 -23)τ 4 -160(87µ 2 -35)τ 2 + 2400µ c 2 =µ 2 (µ 2 -1)τ 8 -8µ(µ 2 -1)(7µ 2 + 4)τ 6 + 4(297µ 4 -266µ 2 -63)τ 4 - 320µ(15µ 2 -23)τ 2 -400(45µ 2 -13) c 3 = -8µ 2 (µ 2 -1) 2 τ 6 + 64µ(µ 2 -1)(µ 2 -3)τ 4 + 160(µ 2 -1)(9µ 2 + 7)τ 2 -12800µ(µ 2 -1) c 4 = -16(µτ 2 -10) 2 (µ 2 -1) 2 (4.30)
Using the original integral formula (4.9), one can show that the initial conditions are:

F (0, τ ) = + ∞ 0 λ 1 2 J 0 (λ) sin( √ λτ )dλ Fµ (0, τ ) = - ∞ 0 λ 3 2 J 0 (λ) sin( √ λτ )dλ Fµµ (0, τ ) = + ∞ 0 λ 5 2 J 0 (λ) sin( √ λτ )dλ + ∞ 0 λ 3 2 J 1 (λ) sin( √ λτ )dλ Fµµµ (0, τ ) = - ∞ 0 λ 7 2 J 0 (λ) sin( √ λτ )dλ - ∞ 0 λ 5 2 J 1 (λ) sin( √ λτ )dλ (4.31)
Following [24], [34] (P.609), let us define the auxiliary function K p,q (R, τ ):

K p,q (R, τ ) = ∞ 0 λ -(1/2)+p J q (λR) sin( √ λτ )dλ (4.32)
Thus:

F (0, τ ) = +K 1,0 (1, τ ) Fµ (0, τ ) = -K 2,0 (1, τ ) Fµµ (0, τ ) = +K 3,0 (1, τ ) + K 2,1 (1, τ ) Fµµµ (0, τ ) = -K 4,0 (1, τ ) -3K 3,1 (1, τ ) (4.33)
The RHSs of Eqs. (4.33) can be computed by recurrence from the results in [34] (22.20, P.609).

K 0,0 (R, τ ) = π 2 √ 2 τ R J (1/4) ( τ 2 8R )J (-1/4) ( τ 2 8R ) (4.34)
The first recurrence relation is obtained by differentiating Eq. (4.32) twice with respect to τ , leading to

K p+1,q (R, τ ) = - ∂ 2 ∂τ 2 K p,q (R, τ ) (4.35)
The second one derives from the well-known recurrence relations for the derivatives of the Bessel functions:

K p+1,q+1 (R, τ ) = q R K p,q (R, τ ) - ∂ ∂R K p,q (R, τ ) (4.36)
Then, starting from Eq. (4.34), and using successively Eqs. (4.35) and (4.36) with at least R = 1, one can show that the initial conditions for the ODE in µ of the Green function F are:

F (0, τ ) = π 16 √ 2 τ 3 J (1/4) ( τ 2 8 )J (-1/4) ( τ 2 8 ) + J (3/4) ( τ 2 8 )J (-3/4) ( τ 2 8 ) Fµ (0, τ ) = - π 64 √ 2 τ           (τ 4 -8)J (1/4) ( τ 2 8 )J (-1/4) ( τ 2 8 ) -6τ 2 J (-1/4) ( τ 2 8 )J (-3/4) ( τ 2 8 ) +6τ 2 J (1/4) ( τ 2 8 )J (3/4) ( τ 2 8 ) +τ 4 J (-3/4) ( τ 2 8 )J (3/4) ( τ 2 8 )           Fµµ (0, τ ) = π 256 √ 2 τ 3           (τ 4 -60)J (1/4) ( τ 2 8 )J (-1/4) ( τ 2 8 ) -16τ 2 J (-1/4) ( τ 2 8 )J (-3/4) ( τ 2 8 ) +16τ 2 J (1/4) ( τ 2 8 )J (3/4) ( τ 2 8 ) +(τ 4 -12)J (3/4) ( τ 2 8 )J (-3/4) ( τ 2 8 )           Fµµµ (0, τ ) = - π 1024 √ 2 τ           (τ 8 -220τ 4 -160)J (1/4) ( τ 2 8 )J (-1/4) ( τ 2 8
)

-(30τ 6 -40τ 2 )J (-1/4) ( τ 2 8 )J (-3/4) ( τ 2 8 )
+(30τ 6 -40τ 2 )J (1/4) ( τ 2 8 )J (3/4) ( τ 2 8 )
+(τ 8 -116τ 4 )J (-3/4) (

τ 2 8 )J (3/4) ( τ 2 8 )           (4.37)
This ODE are original work and may be useful for the calculation of influence coefficient in BEM without the convolution operator.

Validation of the ODEs with respect to µ of F

In order to validate this ODE with respect to µ, we compare the results of evaluating the function F (µ, τ ) by integrating the new ODE in µ 4.29 with the ODE for F with respect to τ which has been validated by Clément [24], [22].

Eq. (4.23) and Eq. (4.29) were integrated using the RK4 method with initial conditions given by Eq. (4.24) and Eq. (4.37), respectively. Fig. 4-1 shows the results for F obtained by integrating the ODE with respect to τ with ∆τ = 0.01; while the Fig. 4-2 shows the results for the ODE with respect to µ with ∆µ = 0.01, 0.001, 0.0001.

The order of the ODE in Eq. (4.29) decreases when µτ 2 -10 = 0 or µ 2 -1 = 0. The first case is shown by the black line in Fig. 4-2 (the other case is µ = 1. A

An alternative derivation of ordinary differential equations in the time domain

In 1998, fourth order differential equations were proposed by Clément [24] for the free-surface Green function and its derivatives in the time domain. They were derived using another method. A general fourth order ordinary differential equation was derived for a class of functions including the time-domain Green function of linearized free-surface hydrodynamics and all its spatial derivatives. Bingham compared the accuracy and efficiency of several method for computing the timedomain Green function including the ODE method based on Clément's work. [9] The ODE for F is given above (Eq. 4.25). We recall in the following the ODEs for the vertical and horizontal derivatives:

ODE for F r (r 2 + Z 2 )F rtttt -ZtF rttt + ( 1 4 t 2 -6Z)F rtt + 11 4 tF rt + 21 4 F r = 0 (4.38)
The initial conditions are

F r (r, Z; 0) = 0 F rt (r, Z; 0) = 6 rZ (r 2 + Z 2 ) 5/2 F rtt (r, Z; 0) = 0 F rttt (r, Z; 0) = 6 r(-r 2 + 4Z 2 ) (r 2 + Z 2 ) 7/2 (4.39) ODE for F Z (r 2 + Z 2 )F Ztttt -ZtF Zttt + ( 1 4 t 2 -6Z)F Ztt + 11 4 tF Zt + 25 4 F Z = 0 (4.40)
The initial conditions are

F Z (r, Z; 0) = 0 F Zt (r, Z; 0) = rZ (-2r 2 + 4Z 2 ) 5/2 F Ztt (r, Z; 0) = 0 F Zttt (r, Z; 0) = -6 Z(3r 2 -2Z 2 ) (r 2 + Z 2 ) 7/2 (4.41)

Frequency-Domain Green function

Boundary value problem in the frequency domain

The frequency-domain Green function of the free-surface hydrodynamic problem is the solution of the following boundary-value problem.

∆ M G ∞ (r, Z; ω) = δ(|P M |) ∀z M ≤ 0 (4.42) -k 0 G ∞ (r, Z; ω) + ∂G ∞ (r, Z; ω) ∂z M = 0 z M = 0; ∀t ≥ 0 (4.43)
Conditions at infinity

∇ M G ∞ (r, Z; ω) → 0 z M → -∞; ∀t ≥ 0 (4.44) lim r→∞ √ r ∂G ∞ (r, Z; ω) ∂r -ik 0 G ∞ (r, Z; ω) = 0 (4.45)

Partial differential equations in the frequency domain

Using Eq.4.16a, Eq.4.16b and Eq. 4.23, let us now derive the PDEs in the frequency domain. Multiplying both sides of Eq.4.16a by H(τ ):

τ 2 H Fττ + 5τ H Fτ + 4(1 -µ 2 )H Fµµ -8µH Fµ + 3H F = 0 (4.46)
Let us introduce the auxiliary function S(µ, τ ) = H(τ ) F (µ, τ ). Taking into account the differential relation between the Heaviside and Dirac functions, i.e δ(τ ) = dH(τ )/dτ , one can show:

S = H F Sτ = δ F + H Fτ Sττ = δ τ F + 2δ Fτ + H Fττ Sτττ = δ τ τ F + 3δ τ Fτ + 3δ Fττ + H Fτττ S τ τ τ τ = δ τ τ τ F + 4δ τ τ Fτ + 6δ τ Fττ + 4δ Fτττ + H Fττττ (4.47)
Using Eq. (4.47) in Eq. ( 4.46), one can show:

τ 2 Sττ + 5τ Sτ + 4(1 -µ 2 ) Sµµ -8µ Sµ + 3 S = τ 2 (δ (1) F + 2δ F (1) ) + 5τ δ F (4.48)
Let us take its Fourier transform seen in Appendix a. Note the Fourier transform of S(τ ) is the frequency dependent part G0 ( ) with R 1 G0 ( ) = G( ). For the left hand side:

L = F(τ 2 Sττ + 5τ Sτ + 4(1 -µ 2 ) Sµµ -8µ Sµ + 3 S) (4.49) F(τ 2 Sττ ) = 2 G + 4 G + 2 G F(5τ Sτ ) = -5( G + G) F(3 S) = 3 G F(4(1 -µ 2 ) Sµµ -8µ Sµ ) = 4(1 -µ 2 ) Gµµ -8µ Gµ (4.50) L = 2 G -G + 4(1 -µ 2 ) Gµµ -8µ Gµ (4.51)
For the right hand side:

R = F(τ 2 (δ (1) F + 2δ F (1) ) + 5τ δ F ) = 0 (4.52)
Finally, it remains:

2 G -G + 4(1 -µ 2 ) Gµµ -8µ Gµ = 0 (4.53)
For Eq. (4.16b), we use the same derivation. Multiplying the both sides by H(τ ):

2H Fττ + µτ H Fτ -2(1 -µ 2 )H Fµ + 3µH F = 0 (4.54)
Using Eq. (4.47), one can show:

2 Sττ + µτ Sτ -2(1 -µ 2 ) Sµ + 3µ S = 2(δ (1) F + 2δ F (1) ) + µτ δ F (4.55)
For the left hand side:

L = F(2 Sττ + µτ Sτ -2(1 -µ 2 ) Sµ + 3µ S) (4.56)
F(2 Sττ ) = 2(i ) 2 G = -2 2 G F(µτ Sτ ) = -µ( Gω + G) F(3µ S) = 3µ G F(-2(1 -µ 2 ) Sµ ) = -2(1 -µ 2 ) Gµ (4.57) L = -µ G -(2 2 -2µ) G -2(1 -µ 2 ) Gµ (4.58)
For the right hand side:

R = F(2(δ (1) F + 2δ F (1) + µτ δ F ) (4.59)
F(2δ (1) F ) = 2(-F (1) (0) + i F (0)) = 2(-µ + 0) = -2µ F(4δ F (1) ) = 4 F (1) (0) = 4µ F(µτ δ F ) = 0 (4.60) R = 4µ -2µ = 2µ (4.61)
Finally:

-µ G -(2 2 -2µ) G -2(1 -µ 2 ) Gµ = 2µ (4.62)
which shows that the derivatives of frequency-domain Green function can be obtained directly without other ODEs.

The time domain homogeneous ODE with natural variables (4.23) obviously still holds when both sides are multiplied by H(τ ), leading to:

H Fττττ + µτ H Fτττ + ( τ 2 4 + 4µ)H Fττ + 7 4 τ H Fτ + 9 4 H F = 0 (4.63)
Using Eq. (4.47), one can show:

Sττττ + µτ Sτττ + ( τ 2 4 + 4µ) Sττ + 7 4 τ Sτ + 9 4 S = + δ(4 Fτττ + µτ 3 Fττ + ( 1 4 τ 2 + 4µ)2 Fτ + 7 4 τ F ) + δ τ (6 Fττ + µτ 3 Fτ + ( 1 4 τ 2 + 4µ) F ) + δ τ τ (4 F (1) + µτ F ) + δ τ τ τ F (4.64)
For the left hand side:

L = F[ Sττττ + µτ Sτττ + ( 1 4 τ 2 + 4µ) Sττ + 7 4 τ Sτ + 9 4 S] (4.65)

Ordinary differential equations in frequency domain

A second order for G with respect to the frequency variable

The second order ODE for G with respect to the frequency variables has been shown in (4.70).

A second order for G with respect to the frequency variable ω

Returning back to the original variables (r, Z, ω), this ODE reads:

ω 2 4 G ωω -ω(ω 2 Z + 3 4 )G ω + (ω 4 (r 2 + Z 2 ) + ω 2 Z + 1)G = 2(1 + Zω 2 ) √ r 2 + Z 2 (4.71)
The initial conditions are:

G(r, Z; 0) = 2 R 1 G ω (r, Z; 0) = 0
(µ 3 - 1 2 ) G + ( 4 -µ 2 + 1) G -(1 -µ 2 ) Gµµ + 2µ Gµ = 1 -2 µ (4.73)
From Eq. (4.62) and Eq. ( 4.73), one can eliminate G . Thus, we obtain a second order ODE for G with respect to µ.

µ(1 -µ 2 ) Gµµ + (2µ 2 -2µ 3 2 -1 -µ 2 ) Gµ + ( 4 µ -2 µ 2 -2 ) G = -µ 2 2 (4.74)
To establish the initial values of G(µ = 0), Gµ (µ = 0), Gµµ (µ = 0), we start from the modified Haskind formulation of the Green function [START_REF] Kim | On the harmonic oscillations of a rigid body on a free surface[END_REF]:

G ∞ (r, Z; ω) = 1 R + 1 R 1 -πk 0 e k 0 Z [H 0 (k 0 r) + Y 0 (k 0 r) + 2 π 0 Z e -k 0 t √ t 2 + r 2 dt] + 2iπk 0 e k 0 Z J 0 (k 0 r) (4.75)
Follow the definition above,

G ∞ (r, Z; w) = 1 R - 1 R 1 + G(r, Z; iw) (4.76) G(r, Z; ω) = 2 R 1 -πk 0 e k 0 Z [H 0 (k 0 r) + Y 0 (k 0 r) + 2 π 0 Z e -k 0 t √ t 2 + r 2 dt] + 2iπk 0 e k 0 Z J 0 (k 0 r) (4.77)
Through the change of variables (r, Z, w) → (R 1 , µ, ) and omitting the imaginary part, one can show:

G(µ, ) = 2 R 1 -π 2 R 1 e -µ 2 [H 0 ( 2 1 -µ) + Y 0 ( 2 1 -µ 2 )+ 2 π 0 -µ e -2 x x 2 + 1 -µ 2 dx] (4.78) Using R 1 G0 ( ) = G( ), we have G = R 1 2 G. Thus: G(µ, ) =(2 -π 2 e -µ 2 [H 0 ( 1 -µ) + Y 0 ( 2 1 -µ 2 ) + 2 π 0 -µ e -x x 2 + 1 -µ 2 dx])/2 (4.79)
From this last equation:

G(µ = 0) = (2 -π 2 H 0 ( 2 ) -π 2 Y 0 ( 2 ))/2 (4.80)
After differentiation of Eq. (4.78):

Gµ (µ = 0) = -2 (2 -π 2 H 0 ( 2 ) -π 2 Y 0 ( 2 ))/2 (4.81)
After a second differentiation:

Gµµ (µ = 0) = -4 (π 2 H 0 ( 2 ) + π 2 Y 0 ( 2 ) + πH 1 ( 2 ) + πY 1 ( 2 ) -4)/2 (4.82)
Note that Eq. (4.80) agrees with Eq 5.9 in [START_REF] Noblesse | The Green function in the theory of radiation and diffraction of regular water waves by a body[END_REF], and that Eq. (4.81) agrees with Eqn. [START_REF] Taylor | Wave diffraction theory-some developments in linear and nonlinear theory[END_REF] with µ = 0.

Validation of the ODE of G with respect to µ

The ODE for G with respect to has been validated by Clément [23], [START_REF] Xie | Use of clement's odes for the speedup of computation of the green function and its derivatives for floating or submerged bodies in deep water[END_REF], [START_REF] Shen | Ordinary differential equation algorithms for a frequency-domain water wave green's function[END_REF]. Here, we validate the ODE with respect to µ (Eq. (4.74)).

The ODE is integrated using the RK4 method with the initial conditions given by (4.80), (4.81), (4.82). The results for G are shown in Fig. 4-4(a) and compared to the results with the direct numerical method. One can see that the ODE of G with respect to µ is well validated 4-4(b).

An alternative derivation of ordinary differential equations in the frequency domain

Eq. (4.71) for G(r, Z; ω) can also be obtained from (4.25) by using the Fourier Transform [23].

ODE for G r

Following [23], the ODE for the spatial derivative G r can be obtained as follows.

We start from the time domain domain homogeneous ODE for G r , we multiply by the Heaviside step function H(t) and we use Eq. (4.47). It leads to:

(r 2 + Z 2 )S (4) -ZtS (3) + ( 1 4 t 2 -6Z)S (2) + 11 4 tS (1) + 21 4 S = + δ(4(r 2 + Z 2 )F (3) r -3ZtF (2) r + 2( 1 4 t 2 -6Z)F (1) r + 11 4 tF r ) + δ (1) (6(r 2 + Z 2 )F (2) r -3ZtF r (1) + ( 1 4 t 2 -6Z)F r ) + δ (2) (4(r 2 + Z 2 )F (1) r -ZtF r ) + δ (3) ((r 2 + Z 2 )F r ) (4.83)
Let us take the Fourier transform of the left hand side of Eq. (4.83):

L = F((r 2 + Z 2 )S (4) -ZtS (3) + ( 1 4 t 2 -6Z)S (2) + 11 4 tS (1) + 21 4 S) (4.84)
By using the rules of Fourier transform for derivatives

F f (n) (t) = (iω) n f (iω)
and for the product by a polynomial F{t m f (t)} = i m d m dω m f (iω) .

L 1 = (r 2 + Z 2 )(iω) 4 Ŝ = ω 4 (r 2 + Z 2 ) Ŝ (4.85) L 2 = -Z(i d dω ((iω) 3 Ŝ)) = -Z(ω 3 S ω + 3ω 2 Ŝ) (4.86) L 3 = 1 4 i 2 d 2 dω 2 ((iω) 2 Ŝ) -6Z(iω) 2 Ŝ =
R 2 = -(6(r 2 + Z 2 )F (2) r -3ZtF r (1) (0) + ( 1 4 t 2 -6Z)F r ) (1) (0) + iω(6(r 2 + Z 2 )F (2) r -3ZtF r (1) + ( 1 4 t 2 -6Z)F r )(0) = -(6(r 2 + Z 2 )F (3) r -9ZF (1) r )(0) (4.94) R 3 = (4(r 2 + Z 2 )F (1)
r -ZtF r ) (2) (0) -2iω((4(r 2 + Z 2 )F (1) r -ZtF r )) (1) (0) -ω 2 (4(r 2 + Z 2 )F (1) r -ZtF r )(0) = 4(r 2 + Z 2 )F (3) r -ω 2 (4(r 2 + Z 2 )F (1) r ) -2ZF (1) r (4.95)

R 4 = (r 2 + Z 2 )(-F (3) r (0) + 3iωF (2) r (0) + 3ω 2 F (1) r (0) -iω 3 F r (0)) = (r 2 + Z 2 )F (3) r (0) + 3(r 2 + Z 2 )ω 2 F (1) r (0) (4.96)
The right hand side is then:

R = (r 2 + Z 2 )F (3) r (0) -ω 2 (r 2 + Z 2 )F (1) r (0) -5ZF (1) r (0) = -6r(1 + Zω 2 ) (r 2 + Z 2 ) 3/2 (4.97)
And the ODE for G r is:

1 4 ω 2 G rωω -(Zω 3 + 7 4 ω)G rω + (ω 4 (r 2 + Z 2 ) + 3ω 2 Z + 3)G r = -6r(1 + Zω 2 ) (r 2 + Z 2 ) 3/2 (4.98)
with initial conditions:

G(r, Z; 0) = -2r R 3 1 G rω (r, Z; 0) = 0 (4.99) ODE for G Z
The difference between G r and G Z for the left hand side is the coefficient of the last term, L 5 = 25 4 Ŝ. Therefore:

L = 1 4 ω 2 Ŝωω -(Zω 3 + 7 4 ω) Ŝω + (ω 4 (r 2 + Z 2 ) + 3ω 2 Z + 4) Ŝ (4.100)
The difference for the right hand term are the initial conditions :

R = (r 2 + Z 2 )F (3) r (0) -ω 2 (r 2 + Z 2 )F (1) r (0) -5ZF (1) r (0) = -4Z(2 + ω 2 Z) + 2r 2 ω 2 (r 2 + Z 2 ) 3/2 (4.101)
Thus, the ODE for G Z is:

1 4 ω 2 G Z ωω -(Zω 3 + 7 4 ω)G Z ω + (ω 4 (r 2 + Z 2 ) + 3ω 2 Z + 4)G Z = -4Z(2 + ω 2 Z) + 2r 2 ω 2 (r 2 + Z 2 ) 3/2

Résumé

Les équations différentiels ordinaires dans le domaine temporel dérivés dans le dernier chapitre ont été considérés comme un moyen efficace d'évaluer la fonction de Green et son gradient [22] et peuvent être utilisés pour éviter les intégrales de convolution de l'équation intégrale dans le domaine temporel [19], [20], [21]. Nous pensons que les équations différentielles ordinaires dans le domaine fréquentiel peuvent être utiles pour évaluer la fonction de Green et avoir d'autres applications. Cependant, il est difficile de résoudre ces équations différentielles dans le domaine de fréquence depuis l'origine en raison de la singularité. Dans ce chapitre, une nouvelle équation différentielle ordinaire exempte de singularité est dérivée et est utilisée pour évaluer la fonction de Green dans le domaine fréquentiel et son gradient. Une méthode d'expansion efficace pour obtenir la fonction de Green est également proposée pour les petites fréquences. Les comparaisons avec des évaluations de la fonction de Green utilisant la méthode classique d'intégration directe sont fournies. Elles montrent que l'équation différentielle modifée peut fournir des estimations précises de la fonction de Green.

A new ordinary differential equation for the evaluation of the frequency-domain Green function

The ODEs in time domain derived in the last chapter have been considered as an efficient way to evaluate the Green function and its gradient [22] and can be used to avoid the convolution integrals of time domain boundary integral equation [19], [20], [21]. We believe that ODEs in the frequency domain found in Chapter 4 and [23] may be also useful to evaluate the Green function and have other applications. However, it is challenging to solve the ODEs of Green function and its gradient in the frequency domain at the origin due to the singularity. In this chapter, modified ordinary differential equations free of singularity is derived and is used to evaluated the frequency-domain Green function and its gradient. The details can be seen in attached article which is accepted by Applied Ocean Research. The computational time for calculating the Green function is not discussed in attached article. Here we add this information. The computational time for calculating the Green function in the attached article for G(µ, ) and G(r, Z; ω) are 4.52 × 10 -8 , 4.39 × 10 -8 . We can see that ODE-based method is faster than other existing methods shown in Chapter 3 but this advantage is not so obvious.

Introduction

Boundary element method (BEM)-based codes are widely used in the industry and in academia to investigate wave structure interactions effects on marine structures. BEM codes rely on the frequency domain linear free-surface potential flow theory which involves the free-surface Green function and its gradient. An important part of the numerical burden in BEM codes relates to the numerous numerical evaluations of the free-surface Green function. Indeed, they must be evaluated for a large set of geometrical configurations and for a wide range of frequencies to assess the structure response across all possible incident wave conditions.

Various analytical formulations and algorithms for the calculation of the free-surface Green function and its derivatives have been proposed. Pioneering work was performed by Noblesse [1][2][3] and Newman [4][5][6] using sub-domains methods with series expansions and polynomial approximations. Double Chebyshev polynomials approximations methods with special functions were proposed by Chen [7,8] and Wang [9]. Other efficient methods include tabulated functions with Lagrange interpolations [10], approximation methods with coordinates-transformation [11,12], eigen-function expansion [13], semianalytical method based on a singularity subtraction technique [14], multipole expansion method [15]. Recently, Wu et al. [16,17] proposed a global approximation based on Noblesse's method [2,3]. It uses a simple approximation involving elementary functions for the local flow component. It does not required dividing the computational domain into multiple sub-domains. All those methods were recently reviewed in [18].

The analytical formulations typically include a source term, an image term, and a singular integral. The algorithms implemented in most common BEM softwares (WAMIT [19], HydroStar [20] or NEMOH [21] for example) are based on integral formulations.

A different approach was introduced by Clément [22] who showed that the frequency domain Green function is the solution of a second order ordinary differential equation (ODE) of the frequency variable. Thus, the Green function can be evaluated by integrating numerically this ODE, provided that initial conditions are available. Since then, similar ODEs for the gradient of the frequency domain Green function have been established (see [23] and [24]).

The remaining challenge for the evaluation of the frequency domain Green function and its gradient using ODEs is the evaluation of the initial conditions. Indeed, despite the values of the Green function and its first derivatives are known for ω = 0, the method cannot be started from that point because the ODEs degenerate for ω = 0. Thus, in [24], the ODEs are initialized at the non-zero frequency ω = 1. A powerseries expansion method is used to solve the ODEs for ω < 1 while a trigonometrically fitted block Numerov-type method is used for ω ≥ 1. Note that it is not explained in [24] how the initial conditions are obtained for ω = 1. In [23], initial conditions for arbitrary values of ω were obtained using methods relying on the numerical approximations of the singular integral.

In this paper, we present an alternative method which enables using ω = 0 for the initial conditions. The ODEs for the Green function and its derivatives are recalled in Section 2. The singular behaviour for ω = 0 is highlighted. In Section 3, the singularity is extracted and a modified ODE is introduced. Numerical results and comparisons with direct numerical integration are shown in Section 4. Finally, a conclusion is given in Section 5.

ODEs for the Green function and its derivatives

In this study, only the infinite water depth free-surface Green function is considered. The coordinates and variables are depicted in Fig. 1. The mean free surface level is located at the plane z = 0. The vertical axis z points upwards. The gravity constant is denoted g. The source point P(x P , y P , z P ) and the field point M(x M , y M , z M ) are both lying on or under the free surface (z P ≤ 0, z M ≤ 0). The image source point P′(x P , y P , -z P ) is the mirror of the source point P with respect to the mean free surface. The horizontal distance between the source point P and the field point M is denoted by r. The vertical distance between the image source point P′ and the field point M is -Z. The distance between the source and field points is denoted by R and the distance between the image source point and the field point is

= + R r Z 1 2
2 . The angle θ is defined by cos θ = -Z/R 1 and sin θ = r/R 1 . The relations between the coordinates are given by:

= + = + = + + r x x y y Z z z R x x y y z z ( ) ( ) ( ) ( ) ( ) M P M P M P M P M P M P 2 2 2 2 2
(1)

ODEs with nondimensional original variables

The wave frequency is denoted ω and the time dependent factor of the complex potential is e -iωt . The frequency domain Green function G ∞ can be written:

= + G r Z R R G r Z 4 ( , , ) 1 1 ( , , ) 1 (2)
with The ODEs for the nondimensional Green function and its derivatives can be written [23]:

= + G r Z k k k e J kr dk i k e J k r ( , , ) 2PV ( ) 2
= R R L ¯/ , = R R L ¯/ 1 1 , = r r L ¯/ , = Z Z L ¯/ , = = = k k L L g ¯/ 0 0 2 2 , = k ¯kL, = G ¯GL
+ + + + + = + + G Z G r Z Z G Z r Z 4 ¯¯( ¯¯3 4 ) ¯( ¯(¯¯) ¯¯1) ¯2(1 ¯¯) (¯¯) 2 ¯¯2 ¯4 2 2 2 2 2 2 1/2 (6)
+ + + + + = + + G Z G r Z Z G r Z r Z 4 ¯¯( ¯¯7 4 ) ¯( ¯(¯¯) 3 ¯¯3) 6¯(1 ¯¯) (¯¯) r r r 2 ¯¯¯2 ¯¯4 2 2 2 2 2 2 3/2 (7)
+ + + + + = + + + G Z G r Z Z G Z Z r r Z 4 ¯¯( ¯¯7 4 ) ¯( ¯(¯¯) 3 ¯¯4) 4 ¯(2 ¯¯) 2¯( ¯¯) Z Z Z 2 ¯¯¯2 ¯¯4 2 2 2 2 2 2 2 2 3/2 (8)
For = ¯0, the values of the Green function and its first gradient can be obtained from equation (3). They are given by:

= = G R G ¯2 ¯, ¯0 1 ¯(9a) = = G r R G ¯2¯, ¯0 r r ¯13 ¯¯(9b) = = G Z R G ¯2 ¯, ¯0 Z Z ¯13 ¯¯(9c)
It can be observed that for = ¯0, equations ( 6),( 7) and ( 8) degenerate into zero order equations. Thus, they cannot be used to obtain the second derivatives of the Green function for that frequency. This singularity prevents using = ¯0 as the initial condition for the ODEs. The ODEs apply both to the real and imaginary part of the Green function. However, the imaginary part can be expressed analytically [18], [6], [2]. Thus, in the following, we shall focus only on the evaluation of the real part of the Green function. Moreover, the . ¯on top of the variables will be omitted as hereafter we only discuss the nondimensional Green function.

The singularity of the Green function for ω = 0

The difficulty to use ω = 0 as the initial condition for the ODEs arises not only from the fact that they degenerate for that frequency but also because of the singular behaviour of the Green function for ω = 0.

Let us recall the modified Haskind representation of the Green function [25], [26]:

= + + + G r Z R e r Y r e t r H ( , , ) 2 [ ( ) ( ) 2 dt] Z Z t 1 2 0 2 0 2 0 2 2 2 2 (10)
where H 0 (•) denotes the Struve function as defined by Abramowitz and Stegun [27]. By using a Taylor expansion, the Green function can be rewritten:

Thus, equations ( 23) and ( 24) provide a way to evaluate the horizontal derivative of the Green function from the Green function itself without having to solve an additional ODE.

A new ODE free of the singularity at the origin

Derivation

The ODEs ((6, 7), (8,14), (15,16)) can be written in a general form as follows:

+ + + + + = + x y c x c xy c x c x c y c c x ( ) ( ) 2
where c 21 , c 01 , c 40 , c 20 , c 00 , c 0 , c 2 are parameters independent on the variable x, x ≥ 0. Note that x is the frequency variable, which is either ω or ϖ depending on the considered ODE. y is the non-dimensional Green function or one of its derivatives as function of the non-dimensional or natural variables. The parameters are given for each equation in A.

The initial conditions are:

= y c c (0) / 0 00 (26a) = y (0) 0 (26b)
According to (11), the asymptotic form for y when x is close to 0 is:

= + + + y x c c a x x b x o x ( ) / log ( ) 0 00 1 2 1 2 2 (27)
with

= = + + a b r Z R G 4, 2( log( /(2( )))), for 1 1 2 1 (28a) = = + a b r r R Z R G 0, 4( 1 /(2 ( ))), for r 1 1 1 1 (28b) = = + a b µ G 2, log 2 log(1 ), for 1 1 (28c) = = + a b µ µ R µ G 0, 1 /( ( 1
)), for ˜r

The asymptotic form for the first derivative of y is:

= + + + y x a x x b x o x ( ) (2 log 1) 2 ( ) 1 1 (29)
Let us define the new function z(x) as the difference of the original function y minus the two leading terms in (27)

: = z x y x c c a x x ( ) ( ) / log 0 00 1 2 (30)
Thus, y(x) and its derivative are given by:

= + + y x z c c a x x ( ) / log 0 00 1 2 (31a) = + + y x z a x x a x ( ) 2 log 1 1 (31b)
By introducing (31) into (25), one can show that z(x) is the solution of a the following ODE:

+ + = x z c x xz d x z x f x ( ) ( ) ( ) 2 2 (32)
with

= + c x c x c ( ) 21 2 01 (33a) = + + d x c x c x c ( ) 40 4 20 2 00
(33b) For the case = y G ˜, the coefficients involved in the last bracket of f (x) are c 01 = -3, a 1 = -2, c 2 = -4μ, c 20 = -4μ, c 0 = 4 and c 00 = 4 according to (14). Thus, 2 + 2c 01 + c 00 = 0 and -(3 + c 01 ) a 1 + c 2 -c 20 c 0 /c 00 = 0. It can be shown that it is also the case for y = G, y = G r and G ˜r using ( 6), ( 7) and (15). f(x) can then be rewritten:

= + + f x c a x x c c a x c a c c c x ( ) log [ (2
= + + f x c a x x c c a x c a c c c x ( ) log [ (2 ) log /
For x → 0, c(x) → c 01 and d(x) → c 00 according to equation (33) while f(x) → 0 according to equation (34).

Moreover, it can be noted that for x → 0:

= + z x b x o x ( ) ( ) 1 2 2 (35a) = + z x b x o x ( ) 2 ( ) 1 (35b)
according to the asymptotic expression (30) of y(x). Thus, z(x) → 0 and z′(x) → 0 for x → 0.

In comparison with the second-order ODE (25) for y(x), the second order ODE (32) for z(x) is of the same form on the left hand side. However, it has a non-homogeneous term on the right hand side whose leading order is x 4 log x for x → 0. Therefore, in contrast to y″(x), z″(x) tends to a finite value for x → 0 (that is 0 according to equation ( 34)). The ODE for z equation ( 32) is free of the singularity for x → 0.

Practical representation of the new ODE for its numerical integration

The direct application of a numerical integration scheme (e.g. Runge-kutta (RK4)) for the ODE (32) remains difficult in practice because of the coefficient x 2 associated with z″. Hereafter, we present a method to deal with this issue.

Let us assume that x ≠ 0. The ODE (32) can be re-written as:

+ + = z c x z x d x z x f x ( ) / ( ) / ( ) 2 (36)
Let us define the new variables = u v z z x ( , ) ( , / ). Their derivatives are: = u v z u x v x ( , ) ( , / / ). By introducing the new variables in the second-order ODE (36), one can obtain the system of differential equations:

= x x Y F Y ( , / ) (37)
with Y and F defined by

= = { } u v x x f x c x d x x Y F Y Y and ( , / ) ( ) 0
( ) ( ) 1 1 / (38)
The initial values for x → 0 are

= = = = + { } { } { } b c c b Y Y F 0 0 ˆ| 2 1 ( 2 ) 1 x x Y 0 0 0 1 0 00 01 1 (39)
in which we used (35) to obtain and (33) for c(0) = c 01 , d(0) = c 00 and f (0) = 0 to obtain F 0 defined by (38).

Both the initial values of and F 0 being finite, usual numerical integration schemes can be used to integrate (37) in order to obtain values of the Green function from the initial conditions at x = ω = ϖ = 0. The coefficients in equations (38) and (39) and relation between y(x) and z(x) for G, G r , G ˜, G ˜r are given in Appendix A.

Results

In this section, the results of the evaluation of the Green function and its horizontal derivative by using the ODE are compared to the results obtained using the direct integration method as described in [23].

The Green function as function of the original and natural variables are evaluated by solving the system (37) by using the Runge-Kutta 4 method with a constant step size Δϖ = 0.001.

The initial conditions are for ϖ = 0 from Eq. ( 39). The results are Fig. 2. Results for the evaluation of G ˜with the ODE and error with the direct integration method. The ODE is solved using the RK4 method with Δϖ = 0.001 ϖ ∈ (0, 8], and μ ∈ [0, 1).

shown in Fig. 2 for ϖ ∈ (0, 8], μ ∈ [0, 1). The error is defined as the difference between the values of the Green function obtained by integrating the ODE and the direct integration method with 8 decimal accuracy [18]. One can see that a 6 decimals (6D) accuracy is obtained. Thus, the method using the ODE is as accurate as the most accurate of the other numerical methods for evaluating the Green function [18], that are based on the approximations of integral formulations. The horizontal derivative of the Green function is obtained from the Green function itself using (23). The values and errors are shown in Fig. 3 for ϖ ∈ (0, 8], and μ ∈ [0, 1). A satisfactory absolute accuracy of 4 decimals is obtained.

Results for the Green function as function of the original variables and its horizontal derivative are shown in Fig. 4 and5 . A 6D accuracy is obtained.

Therefore, the ODE-based method can predict the Green function as accurately (6D accuracy) as existing alternative methods. Nevertheless, Fig. 3. Results for the evaluation of G ˜r and error with the direct integration method. The ODE is solved using the RK4 method with Δϖ = 0.001 ϖ ∈ (0, 8], and μ ∈ [0, 1). Fig. 4. Results for the evaluation of G and error with the direct integration method. The ODE is solved using the RK4 method with Δω = 0.001. Fig. 5. Results for the evaluation of G r and error with the direct integration method. The ODE is solved using the RK4 method with Δω = 0.001.

A n = [2n(2n -1) + 2nc 01 + c 00 ] a n + [2(n -1)c 21 + c 20 ] a n-1 + c 40 a n-2 = 0 B n = [2n(2n -1) + 2nc 01 + c 00 ] b n + [2(n -1)c 21 + c 20 ] b n-1 + c 40 b n-2 + (4n -1 + c 01 )a n + c 21 a n-1 = 0 (5.5)
for n ≥ 3.

These recurrence relations can be used to evaluate z, z for small values of x. According to (31), the original function y and its derivative y can be then obtained.

It should be noted that, according to the parameters c 00 and c 01 in equations A.1, A.5, the coefficient of a n and b n will not be zeros for n ≥ 3 for G and G.

Results of evaluation of the Green function for small frequencies using the expansion method

The values of the Green function and its first derivative can be calculated using the expansion method for < 1 or ω < 1. Figures 5-1 and 5-2 shows the absolute errors when using this approach. 11 terms were used for G while 14 terms where used G, which are enough to achieve a 6D accuracy.

Results of boundary element method with ODEbased Green function evaluation

In this thesis, a new boundary element method code has been developed. The programming language is FORTRAN. The geometry is represented by flat quadrilaterals with normal vectors oriented towards the fluid and the mixed of sources and normal dipoles distributions are adopted with constant density on each panel.

The singular Rankine source and image source components 1/R+1/R 1 and its gradient are analytically integrated over a flat quadrangle [START_REF] Hess | Calculation of non-lifting potential flow about arbitrary three-dimensional bodies[END_REF], [START_REF] Newman | Distributions of sources and normal dipoles over a quadrilateral panel[END_REF]. The free-surface term of the Green function and its gradient in Chapter 5 are numerically integrated via Gaussian quadrature rule with 4 Gaussian points. The complex linear matrix is solved by LAPACK routine zgesv.

The extended boundary condition method with a 'lid' on the free-surface is used because it can remove all the irregular frequencies without extra selection of constant, surface or points.

For the calculation of the Green function, the new code implements the direct integral method and the ODE-based method.

In this section, we present validation tests for this new code for a hemisphere, a boxbarge and the KRISO Container ship (KCS).

Note that the reference length it set to L = 1m, the density of water is defined as ρ = 1025kg/m 3 , the gravity is g = 9.81m/s 2 .

Hemisphere

The analytical solution of wave radiation by a floating hemisphere was given by Hulme [START_REF] Hulme | The wave forces acting on a floating hemisphere undergoing forced periodic oscillations[END_REF]. The hemisphere is discretized with 500, 2000, 4500 panels shown in Figs. 6-1(a), 6-1(b), 6-1(c). To remove the irregular frequencies, the interior freesurface was discretized as shown in Fig 6-1(d). The additional number of panels is 225. The added mass coefficient A 11 and A 33 for surge and heave are nondimensionalized by ρV and the corresponding wave-damping coefficients B 11 and B 33 are nondimensionalized by ρωV where ρ is the water density, V the corresponding numerical volume and ω the wave frequency. The center of gravity is defined as (0, 0, 0). The radius of gyration is K xx = 16, K yy = 16, K zz = 8.

The results for the hemisphere are also compared to the analytical solution given by Hulme [START_REF] Hulme | The wave forces acting on a floating hemisphere undergoing forced periodic oscillations[END_REF] and these of Hydrostar, see Figure 6-3. An excellent agree-ment can be observed. The irregular frequency removal is then validated. Figure 6-2 shows the hydrodynamic coefficients with and without irregular frequencies removal. It can be seen in Figure 6-2(a) that for the considered hemisphere, the irregular frequencies are located around 1.4 and 1.8. With the extended boundary condition method, the irregular frequencies are well removed.

Boxbarge

The second test case is a box barge (boxbarge) . In contrast to the hemisphere which has a smooth surface, the boxbarge has sharp corners.

The main dimension and hydrostatic characteristics of the boxbarge are shown in Table . 6.1. The boxbarge is discretized into 720, 1940, 3710 panels and the interior free-surface is discretized with 200 panels as shown Fig. 6-5. The linear roll viscous damping follows that in [START_REF]HydroStar for experts user manual[END_REF].

The added mass coefficient of surge, sway and heave are nondimensionalized by ρV . For for roll, pitch and yaw, they are nondimensionalized by ρV B, ρV L, ρV H, respectively. The comparison of added mass and damping coefficients with or without irregular frequencies are shown in Fig. 6-6. The irregular frequencies are well removed. The added mass and damping coefficients of our inhouse code are also compared with these of Hydrostar in Fig. 67. A good agreement can be found. It can be observed that the results with 1940 panels are converged. The added mass and damping coefficients calculated with the direct integral method and the ODE-based method for the evaluation of the Green function are compared. For ω ∈ [0.1, 1.8] with ∆ω = 0.0017, the maximum absolute difference of nondimensional added mass coefficients and damping coefficients between Direct Integral method and ODE-based method are listed in Table. 6.2 and Table. 6.3, respectively. It is found that the maximum difference is 8.39e -08. This small difference further validates the usability of ODE-based method.

RAOs with ODE-based Green function method associated to boxbarge at an incident angle β = 180 • are presented in Fig. 678. The average time for the calculation of the hydrodynamic coefficients for one frequency are shown in Fig. 6789. In comparison to the direct integral method, the ODE-based method appears to be very efficient. Indeed, the CPU time for the calculation of matrix of influence is not important as that of solving the linear matrix of the boundary element method, specially, when the panels are large. It should be noted that the ODE-based method is efficient only for multi-frequencies and for very large number of panels, it will be difficult to store all the values of the Green function and its derivatives. Compared to direct integral method, the ODE-based method is very efficient method. We can also find that the CPU time for the calculation of matrix of influence is not important as that of solving the linear matrix of the boundary element method, specially, when the panel number is large. It should be noted that the ODE-based method is efficient only for multi-frequencies and for very large number of panels, it will be difficult to store all the values of the Green function and its derivatives. Conclusions and perspectives

Conclusions

In this thesis, wave-structure interaction problems were considered in the framework of the linear potential flow theory. The boundary element method with mixed distribution of sources and dipoles, constant panels, free-surface Green function was used. Only the case of infinite water depth with zero speed was studied. The exact mathematical expressions for the Green function were collected and analyzed for better understanding of its numerical behavior. Several existing algorithms for the evaluation of the frequency-domain Green function based on these mathematical expressions were implemented and compared with respect to their accuracies and computational time. A different way to understand the timedomain and frequency-domain Green functions was proposed by introducing a fourth and a second order ordinary differential equations for the time-domain and frequency-domain Green functions, respectively. The solution of ODEs of Green function in frequency domain was shown to be an alternative way to evaluate the Green function. However, it is challenging to solve the second order ODEs of the Green function and its gradient with usual initial conditions at the original point of frequency variable because of the degenerated properties and the singularities. A new way to remove the difficulty by modifying the ODEs associated with new functions free of singularity was introduced in this thesis. The new ODEs were then transformed into their canonic form by using a novel definition of vector functions. This canonic form can be then solved with initial conditions at the original point since all involved terms are finite. An expansion method was also given to find the analytical solution of Green function for small frequencies. The ODE-based Green function was validated by calculating the hydrodynamics coef-ficient (added mass and damping coefficients) and RAOs of a floating hemisphere, a floating boxbarge and a floating KCS.

A series of ODEs of the Green function and its gradient with respect to spatial variables were also introduced in this thesis. They give another perspective for the calculation of surface integrals or volume integrals of the Green function.

Perspectives

The ODEs for the time-domain Green function and frequency-domain Green function can be shown to be correct by comparing Green function with classical formulations of integral form. ODEs can be used for evaluation of Green function. However, the usefulness or beauty should not limited to develop an alternative way to evaluate Green function. The ODEs should give us more or new insights into the wave diffraction and radiations. Lots of work have been done on the way to evaluate time-domain Green function by using ODEs. One work by Clément concerning time-domain Green function and showed the use of ODEs can avoid the convolution integral involved in the time-domain BEM solver at high-frequency [20], [21]. The frequency-domain Green function is just to be used to evaluate Green function in Shen et al [START_REF] Shen | Ordinary differential equation algorithms for a frequency-domain water wave green's function[END_REF] and in this thesis, no work on application of ODE in frequency domain other than Green function evaluation itself. It should be more interesting to find the applications of ODEs in shedding other light on physical or mathematical features of the problem, in addition to just evaluate Green function for which classical formations exist since a long time and used in practice. Unlike ODEs in the frequency domain(F-ODE) which give Green function for all frequencies at a fixed space points, ODEs with spatial variables (S-ODE) give Green function for all space points for fixed frequency. S-ODE is closer to classical formulation of Green function, at least, in its use in BEM. Both S-ODE and F-ODE have a drawback in that we have to use a very small step in space coordinates or frequency in order to achieve a good accuracy. However, the practical problems by BEM concern only largely-spaced discrete space coordinates and frequencies. Thus most results in the solution of S-ODE or F-ODE are useless in frequency-domain problems, unlike time-domain problems, all time-step Green function could be useful.

Recently, a multi-domain method (MDM) to solve ship seakeeping with forward speed [16], [START_REF] Liang | A new multi-domain method based on an analytical control surface for linear and second-order mean drift wave loads on floating bodies[END_REF], [17], [START_REF] Ten | New formulations of the ship-motion green function[END_REF]. They divide the fluid domain into an internal and external domains by an analytical control surface. The Green function method is applied in the external domain limited from the control surface to infinity. The classical formulations of Green function have been used in the evaluation of matrix Title : An efficient method for the calculation of the free-surface Green function using ordinary differential equations Keywords : potential flow, Green function, ordinary differential equation, frequency domain, boundary element method Abstract: The boundary element method (BEM) with constant panels is a common approach for wave-structure interaction problems. It is based on the linear potential-flow theory. It relies on the frequency-domain free-surface Green function, which is the focus of this thesis. First, the mathematical expressions and numerical methods for the frequency-domain free-surface Green function are investigated. Twelve different expressions are reviewed and analyzed. Several existing numerical methods are compared including their computational time and accuracies. Then, a series of ordinary differential equations (ODEs) for the time-domain and frequency-domain free-surface Green functions and their derivatives are derived. These ODEs can be used to better understand the properties of the Green function and can be an alternative way to calculate the Green functions and their derivatives. However, it is challenging to solve the ODEs for the frequency-domain Green function with initial conditions at the origin due to the singularity. This difficulty is removed by modifying the ODEs by using new functions free of singularity. The new ODEs are then transformed in their canonic form by using a novel definition of the vector functions. The canonic form can be solved with the initial conditions at the origin since all involved terms are finite. An expansion method based on series of logarithmic function together with ordinary polynomials which is very efficient for low frequency problems is also developed to obtain analytical solutions.

Finally, the ODE-based method to calculate the Green function is implemented and an efficient BEM solver is obtained. The removal of irregular frequencies is included. The new solver is validated by comparison of hydrodynamic coefficients to analytical solutions for a heaving and surging hemisphere, and to numerical results obtained with a commercial solver for a box barge and the KCS container ship.

2- 1 4 - 3 4 - 4 5 - 1 6 - 1

 143445161 Schematic of immersed part of an arbitrary geometry . . . . . . . . . 2-2 Definition of source point (P), field point (M) and other notations . . 3-1 Wu et al's method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-1 The result F of ODE with respect to τ with ∆τ = 0.01 . . . . . . . . . 4-2 The result F of ODE with respect to µ with ∆µ = 0.01, 0.001, 0.0001 and the black lines are the degenerated order lines. . . . . . . . . . . Same as Fig. 4-2(a) with a modified legend to highlight the limitation between the oscillating and non oscillating zones. . . . . . . . . . The result G of ODE with respect to µ with ∆µ = 0.0001 and the error by using ODE with ∆µ = 0.0001 compared to direct numerical integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Absolute errors for G and Gr evaluated using the expansion method with 11 terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-2 Absolute errors for G and G r evaluated using the expansion method with 14 terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The visualization of meshes of a hemisphere discretized with 500, 2000, 4500 panels and a interior free-surface discretized with 225 panels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6- 4

 4 The RAOs of hemisphere discretized by 2000 panels by using inhouse code and Hydrostar . . . . . . . . . . . . . . . . . . . . . . . . . 6-5 The visualization of meshes of a Boxbarge discretized with 720, 1940, 3710 panels and a interior free-surface discretized with 200 panels . . 6-6 The nondimensional added mass and damping coefficients of Boxbarge with (dashed lines with circle markers) and without (solid lines) irregular frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-7 The nondimensional added mass and damping coefficients of Boxbarge of inhouse code (solid lines), Hydrostar (green red and blue circle markers) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-8 The RAOs of Boxbarge discretized by 1940 panels by using inhouse code and Hydrostar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-9 The average CPU time for boxbarge . . . . . . . . . . . . . . . . . . . 6-10 The visualization of meshes of a KCS discretized with 400, 1940, 3840 panels and a interior free-surface discretized with 650 panels . . . . . 6-11 The nondimensional added mass and damping coefficients of KCS with (dashed lines with circle markers) and without (solid lines) irregular frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-12 The nondimensional added mass and damping coefficients of KCS of inhouse code (solid lines), Hydrostar (green red and blue circle markers) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-13 The RAOs of KCS discretized by 1280 by using inhouse code and Hydrostar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-14 The average CPU time for KCS . . . . . . . . . . . . . . . . . . . . . .

  a surface at the infinity. (τ e ) stands for the fluid domain shown in Figure 2-1. It is limited by S F , S B , (C).

Figure 2 - 1 :

 21 Figure 2-1: Schematic of immersed part of an arbitrary geometry

Figure 2 - 2 :

 22 Figure 2-2: Definition of source point (P), field point (M) and other notations

Figure 3 - 1 :

 31 Figure 3-1: Wu et al's method

T

Fig. 1 .

 1 Fig. 1. Source and field point.

Fig. 3 .

 3 Fig. 3. Sub-domains used in the Telste-Noblesse's method.

  . developed a simple approximation for the local flow term and its derivative: L ≈ L a , L L * * a . The approximated local flow term L a is:

Fig. 5 .

 5 Fig. 5. The errors of the free-surface term of the Green function G -of Newman's, Delhommeau's, Telste and Noblesse's, Wu et al.'s methods compared to direct integral method.

Fig. 6 .

 6 Fig. 6. The errors of the derivative of free-surface term of the Green function G -X of Newman's, Delhommeau's, Telste and Noblesse's, Wu et al. methods compared to direct integral method.

Algorithm

  Total computation time (s) (N X = 220, N Y = 150)

Fig. 7 .

 7 Fig. 7. Local computation time to calculate the Green function and its derivatives based on Newman's, Delhommeau's, Telste-Noblesse's, Wu et al's and direct integral method.

  (a) ∆µ = 0.01 (b) ∆µ = 0.001 (c) ∆µ = 0.0001

Figure 4 - 2 :

 42 Figure 4-2: The result F of ODE with respect to µ with ∆µ = 0.01, 0.001, 0.0001 and the black lines are the degenerated order lines.

Figure 4 - 3 :

 43 Figure 4-3: Same as Fig. 4-2(a) with a modified legend to highlight the limitation between the oscillating and non oscillating zones.

( 4 . 72 )

 472 A second order for G with respect to the space variable µ Let use now combine the three PDEs for G Eqs. (4.53), (4.62), (4.70) to obtain a new ODE with respect to the space variable µ. Multiplying Eq. (4.53) by -1 4 and then summing with Eq. (4.70), we obtain:

Figure 4 - 4 :

 44 Figure 4-4: The result G of ODE with respect to µ with ∆µ = 0.0001 and the error by using ODE with ∆µ = 0.0001 compared to direct numerical integration

1 G4. 3 Summary

 13 Zω (r, Z; 0) = 0 (4.103) In this chapter, new derivations of ODEs of time-domain free-surface Green function and its derivatives associated with the time variable and ODEs of frequencydomain free-surface Green function and its derivatives associated with the frequency variable are introduced. A new ODE of time-domain free-surface Green function associated with the spatial variable and a new ODE of frequency-domain free-surface Green function associated with spatial variables are proposed and validated. Further work on the uniqueness and existence of ODEs need to be done. The properties of ODEs of time-domain Green function associated with spatial variables are net clear and more theory on ODEs are needed.

  number and J 0 (•) is the Bessel function of the first kind and zero order. Let us consider the nondimensional Green function = G G L ¯. Let:

  where L denotes a reference length. The non-dimensional Green function has the same form as the Green function with dimensional variables:

Fig. 1 .

 1 Fig. 1. Definition of source point (P), field point (M) and other notations.

  (a) G error (b) Gr error

Figure 5 - 1 :Figure 5 - 2 :

 5152 Figure 5-1: Absolute errors for G and Gr evaluated using the expansion method with 11 terms

  The response amplitude operators (RAOs) for the response of the hemisphere to incident waves are presented in Fig6-4. They are compared to Hydrostar's results. Again, an excellent agreement is obtained.

Figure 6 - 1 :

 61 Figure 6-1: The visualization of meshes of a hemisphere discretized with 500, 2000, 4500 panels and a interior free-surface discretized with 225 panels

Figure 6 - 2 :Figure 6 - 3 :Figure 6 - 4 :

 626364 Figure 6-2: The nondimensional added mass and damping coefficients of a heaving and surging hemisphere with (dashed lines with circle markers) and without (solid lines) irregular frequencies

Figure 6 - 5 : 24 Figure 6 - 6 :

 652466 Figure 6-5: The visualization of meshes of a Boxbarge discretized with 720, 1940, 3710 panels and a interior free-surface discretized with 200 panels

Figure 6 - 7 :

 67 Figure 6-7: The nondimensional added mass and damping coefficients of Boxbarge of inhouse code (solid lines), Hydrostar (green red and blue circle markers)

  RAOs associated to KCS with ODE-based Green function at an incident angle β = 180 • are shown in Fig. 6-13. The average time for KCS of one frequency are shown in Fig. 6-14.

Titre:

  Accélération du calcul des efforts hydrodynamiques par utilisation des propriétés différentielles des fonctions de Green de à surface libre Mots clés : écoulement potentiel, la fonction de Green, équation différentielle ordinaire, domaine fréquentiel, la méthode des singularités Résumé : Le calcul des efforts hydrodynamiques de premier ordre sur un ou plusieurs corps perçant la surface libre est aujourd'hui bien maîtrisé, et plusieurs codes de calcul implémentant la méthode des singularités (dite BEM ou méthode d'élément frontière) ont été développés. Le cadre est la théorie linéarisée des écoulements potentiels à une surface libre. Dans ces codes BEM, les singularités utilisées ont la propriété intrinsèque de satisfaire à la fois l'équation de Laplace dans le domaine fluide ainsi que la condition linéarisée de surface libre. Ces singularités, dites fonctions de Green à surface libre, dans le domaine fréquentiel en profondeur infinie et sans vitesse d'avance constituent le point focal de cette thèse. Tout d'abord, les expressions mathématiques existantes pour la fonction de Green de surface libre sont examinées. Douze expressions différentes sont passées en revue et analysées. Plusieurs méthodes numériques existantes sont comparées par rapport à leur temps de calcul et leur précision. Ensuite, une série d'équations différentielles ordinaires (ODEs) pour les fonctions de Green de surface libre dans le domaine temporel et le domaine fréquentiel et leur gradient est établie. Ces ODEs peuvent être utilisées pour mieux comprendre les propriétés de la fonction de Green et peuvent constituer un moyen alternatif de calculer ces fonctions de Green et leurs dérivées. Cependant, il est difficile de résoudre numériquement ces ODEs à cause de l'existence d'une singularité à l'origine. Cette difficulté est éliminée en modifiant les ODEs par l'utilisation de nouvelles fonctions sans singularité. Les nouvelles ODEs sont ensuite écrites sous forme canonique en utilisant une nouvelle définition de la fonction vectorielle. La forme canonique peut être résolue avec les conditions initiales à l'origine puisque tous les termes impliqués sont finis. Une méthode d'expansion basée sur une série de fonctions logarithmiques et de polynômes ordinaires, très efficace pour les problèmes de basse fréquence, a également été développée pour obtenir des solutions analytiques. Enfin, la méthode basée sur les ODE pour calculer la fonction de Green est implémentée et un nouveau solveur BEM est obtenu. L'élimination des fréquences irrégulières est incluse. Le nouveau solveur est validé par comparaison des coefficients hydrodynamiques à des solutions analytiques pour une hémisphère, ainsi qu'à des résultats numériques obtenus avec un solveur commercial pour un chaland parallèlépipédique et le porte-conteneurs KCS.

  3.1 Computation time for calculating the Green function and its derivatives for different algorithms . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Main dimension and hydrostatic characteristics of the boxbarge . . . 6.2 The maximum absolute difference of added mass coefficient with 8 digital accuracy between Direct Integral method and ODE-based method for Boxbarge case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017 . . . . . 6.3 The maximum absolute difference of damping coefficient with 8

  . They extended this Dirichlet and Neumann boundary condition to the exterior problem.

	in which M ∈ S B .
		A simpler integral equation exists for the sum of the incident and diffraction
	potential ϕ ID .
				2πϕ ID (M ) +	S B	ϕ ID (P )	∂G(M, P ) ∂n P	dS(P ) = 4πϕ 0 (M )	(2.37)
	After having applied the Dirichlet and Neumann boundary condition, one has:
	2πϕ ID (M ) +		S B	ϕ ID (P )	∂G(M, P ) ∂n P	dS(P ) + α	∂ ∂n M	S B	ϕ ID (P )	∂G(M, P ) ∂n M	dS(M )
	= 4π[ϕ 0 (M ) + α	∂ϕ 0 (M ) ∂n M	]
	2πϕ R (M ) +	S B	ϕ R (P )	∂G(M, P ) ∂n P	dS(P ) + α	∂ ∂n M	S B	ϕ R (P )	∂G(M, P ) ∂n P	dS(P )
											(2.36a)
	=	S B	∂ϕ R (P ) ∂n P	G(M, P )dS(P ) -2πα	∂ϕ R (M ) ∂n(M )	+ α	∂ ∂n M	S B	∂ϕ R (P ) ∂n P	G(M, P )dS(P )
											(2.36b)

  Green function is solution of the following boundary value problem:

https://doi.org/10.1016/j.apor.2018.10.007 Received 20 July 2018; Received in revised form 25 September 2018; Accepted 9 October 2018 ⁎ Corresponding author. E-mail address: aurelien.babarit@ec-nantes.fr (A. Babarit). Applied Ocean Research 81 (2018) 150-163 Available online 24 October 2018 0141-1187/ © 2018 Elsevier Ltd. All rights reserved.

Table 1

 1 Summary table of the expressions of the Green function

	Type	Real part of the free-surface term	Name	Reference

Table 7

 7 Averaged computation time of Telste-Noblesse's algorithm.

	# of sub-	Corresponding (X, Y)	Averaged computation
	domain		time (s)
	1	(X 2 /12 2 + Y 2 /15 2 ≥ 1)	3.607 × 10 -8
	2	(0 ≤ X ≤ 1.2 -0.15Y,	4.110 × 10 -8
		-2 ≤ Y ≤ 0)
	3	(X ≤ min(0.7Y, 9.1))	7.454 × 10 -8
	4	(X ≥ -0.4X)	4.094 × 10 -8
	5	(Therestpart)	6.699 × 10 -8

Table 6 .

 6

	1: Main dimension and hydrostatic characteristics of the boxbarge
	Length (m)	100
	Breath (m)	50
	Draft (m)	20
	Reference point (m)	(0, 0, 0)
	Roll radius of gyration R XX (m)	20
	Pitch radius of gyration R Y Y (m)	20
	Yaw radius of gyrations R ZZ (m)	20
	Linear roll viscous damping (%)

Table 6 .

 6 2:The maximum absolute difference of added mass coefficient with 8 digital accuracy between Direct Integral method and ODE-based method for Boxbarge case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017For ω ∈ [0.1, 0.18] with ∆ω = 0.0017, the maximum absolute difference for the nondimensional added mass coefficients and nondimensional damping coefficients between Direct Integral method and ODE-based method are listed in Table.6.5 and Table. 6.6, respectively. The small difference further validates the usability of the ODE-based method.

	A 11	8.39e -08	A 12	5.46e -17	A 13	3.13e -
	A 14	1.79e -17	A 15	4.88e -09	A 16	3.59e -
	A 21	7.74e -17	A 22	1.95e -08	A 23	1.68e -
	A 24	3.90e -09	A 25	6.68e -17	A 26	5.98e -
	A 31	7.46e -16	A 32	7.48e -16	A 33	5.85e -
	A 34	2.44e -16	A 35	1.35e -15	A 36	1.26e -
	A 41	6.26e -17	A 42	3.90e -09	A 43	5.34e -
	A 44	7.80e -10	A 45	1.10e -16	A 46	7.20e -
	A 51	1.95e -09	A 52	7.45e -17	A 53	1.02e -
	A 54	6.55e -17	A 55	1.95e -09	A 56	2.47e -
	A 61	3.67e -16	A 62	5.49e -16	A 63	8.63e -
	A 64	2.28e -16	A 65	1.78e -16	A 66	1.95e -

Table 6 .

 6 5: The maximum absolute difference of added mass coefficient with 8 digital accuracy between Direct Integral method and ODE-based method for KCS case for ω ∈ [0.1, 1.8] with ∆ω = 0.0017

	A 11	1.30e -07	A 12	1.43e -08	A 13	4.38e -07
	A 14	7.43e -09	A 15	2.99e -07	A 16	2.12e -07
	A 21	3.74e -08	A 22	4.96e -06	A 23	2.11e -07
	A 24	7.19e -07	A 25	7.47e -08	A 26	2.31e -05
	A 31	3.72e -07	A 32	7.26e -08	A 33	8.47e -06
	A 34	4.93e -08	A 35	5.52e -07	A 36	8.11e -07
	A 41	3.89e -09	A 42	2.62e -07	A 43	3.02e -08
	A 44	8.11e -08	A 45	9.07e -09	A 46	4.62e -06
	A 51	2.58e -07	A 52	2.18e -08	A 53	5.29e -07
	A 54	1.34e -08	A 55	5.14e -07	A 56	2.99e -07
	A 61	2.30e -07	A 62	2.28e -05	A 63	1.56e -06
	A 64	4.00e -06	A 65	4.42e -07	A 66	2.98e -04

Acknowledgments

Acknowledgements

The first author is supported by the Chinese Scholarship Council (CSC).

Acknowledgements

The first author is supported by Chinese Scholarship Council(CSC).

2(1 -µ 2 ) Fµ = 2 Fττ + µτ Fτ + 3µ F (4.17) It shows that the derivatives of the Green function can in fact be computed directly from the solution to the Green function ODE, removing the need to solve additional ODEs which is similar as in Li et al. [START_REF] Li | A precise computation method of transient free surface green function[END_REF]. We also found the similar result in frequency domain.

Then, by differentiating once with respect to τ : 2(1 -µ 2 ) Fµτ = 2 Fτττ + µτ Fττ + 4µ Fτ (4.18) Differentiating one more time leads to:

2(1 -µ 2 ) Fµττ = 2 Fττττ + µτ Fτττ + 5µ Fττ (4. 19) Now Eq. (4.18) can be substituted into Eq. (4.16c) to eliminate the Fµτ terms. We get:

The Fµ term in the RHS is now eliminated by using Eq. (4.17) above, leading to

Now eliminating Fµ in Eqs. (4.16a) and (4.16b), we can express Fµµ in terms of F and its successive derivatives with respect to τ :

Finally, using Eq. (4.22) to express the RHS of Eq. (4.21) leads, after simplifications, to sought ODE for F :

The initial conditions are derived from Eq.(4.9)and equation 6.624.6 pp.734 in [34]: limited case which is excluded here). When µ is smaller than 10 τ 2 , one can see that there is an excellent agreement between the results with ODE with respect to τ and the ODE with respect to µ. However, when µ becomes greater than 10 τ 2 , the results for the ODE with respect to µ diverges as instabilities appear in the integration of the ODE.

From Fig. 4-2, one may wonder whether the degenerated order line defines the limit of the oscillating zone in (µ, τ ) plane, as the knowledge of such a limit would be beneficial for the numerical integration of the function. Fig. 4-3, in which the legend has been modified to highlight this limitation, shows the degenerated order line dose not actually separate the oscillating and non-oscillating zones.

Thus:

For the right hand side, using the initial conditions for F of Eq. (4.24):

Thus:

Rearranging the terms, we finally get:

For developing the Fourier transform of the right hand side of Eq. (4.83), let us recall first the following relations derived from the fundamental property of the Dirac delta function and through integration by parts:

We now take the Fourier transform of the right hand side:

where γ is the Euler's constant.

Let us now differentiate twice this last equation to obtain the second derivative of the Green function as function of ω. One can show that for ω → 0, ∂ 2 G/∂ω 2 → -8 ln ω. Thus, the second derivative of the Green function tends to infinity for ω = 0.

ODEs with nondimensional natural variables

In [23], ODEs as function of the natural variables (μ, ϖ) were introduced. The natural variables are defined from the original variables through the change of variables (r, Z, ω) ↔ (R 1 , μ, ϖ):

Note that as Z has values in [-∞ , 0] and r has values in [0, ∞], the range of possible values for μ is limited to [0, 1].

The set of ODEs as function of the natural variables are [23]:

(2 ) (1 )

The initial conditions for ϖ = 0 as function of the natural variables can be obtained according to (9) are:

As for the case of the original variables, it can be observed that the ODEs degenerate for ϖ = 0. Thus, the ODEs as function of the natural variables cannot be integrated from ϖ = 0 either.

Relations between the Green function and its spatial derivatives

The vertical derivative of the Green function can be obtained directly from the Green function itself [18], [6], [2]. In the following, we show that the same applies to the horizontal derivative.

Let us recall the nondimensional time-domain free-surface Green function, F ∞ (r, Z, t) [27], [28]:

with

where t is the time variable, δ(•) is the Dirac delta function and H(•) is the Heaviside unit step function. The terms associated to the Dirac delta function are often referred to the impulsive part of the Green function while the term with the Heaviside unit step function is called the memory part. F is solution of the following partial differential equation:

Using the change of variables ( 13) with = t R / 1 , one can show:

. Applying the Fourier transform as in [22], one can show that the Green function as function of the natural variables is a solution of the partial differential equation:

Thus, the derivative of the Green function G ˜µ can be obtained from the knowledge of G ˜and G ˜. From G ˜µ and G ˜, the horizontal derivative of G ˜can be obtained using:

, the horizontal derivative of the Green function as function of the original variables can be obtained by using:

it can be noted that 4D accuracy is expected to be sufficient for practical applications according to [16,17,29]

Conclusion

In this paper, we present a new ordinary differential equation for the evaluation of the frequency domain Green function. In contrast to original Clément's ODE, this modified ODE allows for the zero frequency to be used for the initial conditions. Comparisons with the direct integration method show that 6D accurate estimates of the Green function can be obtained with this ODE. This method may be used for calculating the Green function and its gradient. It will be investigated in future work whether it is more efficient than conventional methods for the calculation of hydrodynamic coefficients in BEM codes. G According to equation ( 6),

G r According to equation ( 7),

G ˜According to equation ( 14),

An expansion method of Green function in small frequencies

An efficient expansion method to obtain the Green function in small frequencies is also proposed for small frequencies. The computational time for calculating the Green function by using expansion method for G(µ, ) and G(r, Z; ω) are 3.36 × 10 -8 , 2.71 × 10 -8 , respectively.

An expansion method for the initial conditions for the ODEs of G and G for small x

An expansion method for the initial conditions for the ODEs of G and G is introduced. The gradient of Green function then can be derived by using ( 23) and ( 24). The ODE (32) involves only terms of the form x 2n with n ≥ 1 and x 2n log x with n ≥ 2. Thus, let us look for a solution of the ODE of the form:

The derivatives read:

Introducing (5.1), (5.2) into (32) and after identifications and simplifications, one gets:

Note that a 1 and b 1 are defined according to equation (28).

KCS

The KRISO Container Ship (KCS) is further calculated to validate the solver for a complex geometry.

The main dimension and hydrostatic characteristics of the KCS are shown in Table . 6.4. The KCS is discretized into 400, 1940, 3840 panels and the interior freesurface is discretized with 650 panels as shown in Fig. 6-10. The added mass and the damping coefficients are nondimensionalized as explained in the last subsection. Figure 6-11 shows the removal of irregular frequencies works well. Then, the added mass and damping coefficients are compared to . A good agreement can be found.

Summary

In this section, the extended boundary condition method and ODE-based Green function are used in the boundary element method code. The resulting inhouse code is validated by comparing to the analytical solution and numerical results of Hydrostar for hemisphere, boxbarge and KCS. The ODE-based method to calculate the Green function is validated to be an accurate and efficient method suitable for multi-frequencies calculation.

coefficients. Green function is not explicitly evaluated but its spatial integrations. S-ODE may be interesting in MDM.

The extensions of ODEs of Green function in the finite depth of water and Green function in forward speed are expected.

Appendix Appendix a

The Green function of frequency domain G(r, Z; ω) can be derived as the Fourier transform of the time-domain Green function with the following convention and notations: G(r, Z; ω) = Ŝ(r, Z; ω) = F(S(r, Z; ω)) = where S(r, Z; ω) = H(t)F (r, Z; t).

The basic rules of Fourier transform for a derivative and for the product by a polynomial are shown as below.

Now we will discuss the functional relationships of Fourier transform that we will use with the natural variables. We assumed = ω R 1 and using the natural variables (µ, τ ). with S(r, Z; ω) = H(t) F (r, Z; t).

F( S(µ, τ )) = F( S(µ, 1

It should be noted G(µ, ) = R 1 G0 (µ, ) and G0 (µ, ) is the corresponding value of Fourier transform of S(r, Z; ω). τ m S(µ, τ ) = R

Sn (µ, τ ) = R

From equation 8 and 10, we have

119