Green function in infinite depth

The boundary integral equations relate the potential \(\Phi\) to the normal velocity \(u \cdot n\) via the Green function \(G\). Let us now discuss the evaluation of this function for an infinite water depth. See also [X18] for a review of several expression and evaluation methods for this expression.

Green function’s Hankel form

The infinite depth Green function takes the following form

(61)\[-4\pi G(x, \xi) = \frac{1}{\|x - \xi\|} + k \mathcal{G}\left(k \sqrt{(x_1 - \xi_1)^2 + (x_2 - \xi_2)^2}, k (x_3 + \xi_3) \right)\]

The first term of \(G\) is the usual Green function for the 3D Laplace equation without our specific boundary conditions. The wave term \(\mathcal{G}\) is complex-valued and it is introduced to satisfy the boundary conditions (16).

Property 1 (Symmetries)

The infinite depth Green function \(G\) is scaling-invariant in the sense of

(62)\[\forall x, \xi, \quad G(x, \xi, k) = G(k \xi, k x, 1).\]

(hence the dependency to \(k\) will never be explicitly written in the following formulas)

It is also symmetric in the sense of

(63)\[\forall x, \xi, \quad G(x, \xi) = G(\xi, x).\]

The first term of (61) is invariant under all rotations and translations, whereas the wave term is invariant under isometric transformations that don’t change the vertical coordinate (reflection across a vertical plane, rotation around a vertical axis, translation following an horizontal vector).

Introducing the dimensionless variables \(r = k \sqrt{(\xi_1 - x_1)^2 + (\xi_2 - x_2)^2}\) and \(z = k (x_3 + \xi_3)\), the wave term in its usual Hankel transform form reads:

(64)\[\mathcal{G}(r, z) = \frac{1}{k} \int_0^\infty \frac{\kappa + k}{\kappa - k} \exp \left(\kappa z \right) J_0 \left(\kappa r \right) d \kappa\]

where \(J_0\) is the Bessel function of the first kind. The integrand above is singular for \(\kappa = k\). Avoiding the singularity with a proper integral path avoiding the singularity in the complex plane leads to

(65)\[\mathcal{G}(r, z) = \frac{1}{k} \mathop{PV} \int_0^\infty \frac{\kappa + k}{\kappa - k} \exp \left(\kappa z \right) J_0 \left(\kappa r \right) d \kappa + 2 \pi i e^{z} J_0(r)\]

where \(\mathop{PV}\) stands for the principal value of the integral.

Singularities extraction

Lemma 1 (Lipschitz integral)

Also called Lipschitz-Hankel integral

(66)\[\forall z < 0,\, \forall r \in \mathbb{R}, \quad \frac{1}{\sqrt{r^2 + z^2}} = \int_0^\infty \exp \left(\kappa z \right) J_0 \left(\kappa r \right) d \kappa\]

(see e.g. https://mathworld.wolfram.com/LipschitzsIntegral.html)

Using the Lipschitz integral, the wave part of the infinite depth Green function can be rewritten as:

(67)\[\begin{split}\mathcal{G}(r, z) & = \frac{1}{\sqrt{\tilde{r}^2 + \tilde{z}^2}} + \frac{2}{k} \int_0^\infty \frac{k}{\kappa - k} \exp \left(\kappa z \right) J_0 \left(\kappa r \right) d \kappa \\ & = \frac{1}{\sqrt{\tilde{r}^2 + \tilde{z}^2}} + \mathcal{G}^+(r, z)\end{split}\]

that is

(68)\[-4 \pi G(x, \xi) = \frac{1}{|x - \xi|} + \frac{1}{|x - S_0(\xi)|} + k \, \mathcal{G}^+(kr, kz)\]

where \(S_0(\xi)\) is the reflection of \(\xi\) across the free surface \(S_0(\xi) = (\xi_1, \xi_2, -\xi_3)\).

The newly introduced term is a similar to the Rankine term \(\frac{1}{|x - \xi|}\) except one point has been mirrored across the free surface, and is thus referred to as the reflected Rankine term. It is singular when \(x = S_0(\xi)\), which only appends when computing the interaction of a panel horizontal on the free surface with itself. Note that for horizontal panels on the free surface, there is another logarithmic singularity in \(\mathcal{G}^+(kr, kz)\) which needs to be handled.

Note

For convenience, the reflected Rankine term \(\frac{1}{|x - S_0(\xi)|}\) can also be written as \(\frac{1}{|S_0(x) - \xi|}\). This is especially useful when computing the integral of this term \(\int_{\Gamma} \frac{1}{|x - S_0(\xi)|} d\xi\).

Alternatively, the Lipschitz integral can also be used to rewrite (64) as

(69)\[\begin{split}\mathcal{G}(r, z) & = - \frac{1}{\sqrt{\tilde{r}^2 + \tilde{z}^2}} + \frac{2}{k} \int_0^\infty \frac{\kappa}{\kappa - k} \exp \left(\kappa z \right) J_0 \left(\kappa r \right) d \kappa \\ & = - \frac{1}{\sqrt{\tilde{r}^2 + \tilde{z}^2}} + \mathcal{G}^-(r, z)\end{split}\]

that is

(70)\[-4 \pi G(x, \xi) = \frac{1}{|x - \xi|} - \frac{1}{|x - S_0(\xi)|} + k \, \mathcal{G}^-(kr, kz)\]

The notation \(\mathcal{G}^-\) and \(\mathcal{G}^+\) are based on the one from [X18]. In Capytaine, the variants above are referred to as the low-frequency variant for (67) and the high-frequency variant for (69), due to the asymptotic behavior of the Green function:

(71)\[-4 \pi G(x, \xi) \rightarrow_{k \rightarrow 0} \frac{1}{|x - \xi|} + \frac{1}{|x - s(\xi)|}\]

that is

(72)\[k \mathcal{G}^+(r, z) \rightarrow_{k \rightarrow 0} 0\]

and

(73)\[-4 \pi G(x, \xi) \rightarrow_{k \rightarrow \infty} \frac{1}{|x - \xi|} - \frac{1}{|x - s(\xi)|}\]

that is

(74)\[k \mathcal{G}^-(r, z) \rightarrow_{k \rightarrow \infty} 0\]

As discussed in [A24], working with the low-frequency variant is usually more accurate and is thus the default in recent version of Capytaine.

Property 2

From the definitions of \(\mathcal G^+\) and \(\mathcal G^-\), we have

(75)\[\frac{\partial \mathcal G^+}{ \partial \tilde z} = \mathcal G ^-(r, z) = \mathcal{G}^+(r, z) + \frac{2}{\sqrt{ r^2 + z^2}}\]

Guével-Delhommeau formulation

Lemma 2 (Integral form of the Bessel function)

(76)\[J_0(x) = \frac{2}{\pi} \Re \int_0^{\pi/2} \exp\left(i x \cos \theta \right) d\theta\]

(see e.g. https://en.wikipedia.org/wiki/Bessel_function#Bessel’s_integrals):

Inputing the integral form of the Bessel function into (67) gives:

(77)\[\begin{split}\mathcal{G}^+(r, z) & = \frac{4}{\pi} \mathop{PV} \int_0^\infty \Re \int_0^{\pi/2} \frac{\exp \left(\kappa \left(z + i r \cos\theta \right) \right)}{\kappa - k} d\theta d\kappa \\ & \qquad \qquad \qquad \qquad + 4 i \Re \int_0^{\pi/2} \exp( z + i r \cos \theta) d\theta\end{split}\]

which can be rewritten as

(78)\[\begin{split}\mathcal{G}^+(r, z) & = \frac{4}{\pi} \Re \left( \int_0^{\pi/2} e^{\zeta(r, z, \theta)} \left[ E_1(\zeta(r, z, \theta)) + i \pi \right] ) \, \mathrm{d} \theta \right) \\ & \qquad \qquad \qquad \qquad + 4 i \Re \left( \int_{0}^{-\pi/2} e^{\zeta (r, z, \theta)} \, \mathrm{d} \theta \right)\end{split}\]

where

(79)\[\zeta (r, z, \theta) = z + i r \cos \theta.\]

and \(E_1\) is the first exponential integral, defined as

(80)\[E_1(\zeta) = \int_\zeta^\infty \frac{e^{-t}}{t} \mathrm{d} t.\]

where we used the following property

(81)\[\begin{split}\int_0^\infty \frac{e^{-k z}}{k - a} dk = \begin{cases} e^{-a z} (E_1(az) + i \pi) & \Im(z) \ge 0 \\ e^{-a z} (E_1(az) - i \pi) & \Im(z) < 0 \end{cases}\end{split}\]

Note

In [Del87] integrals of the form \(\int_{-\pi/2}^{\pi/2} \ldots d \theta\) are used. Given the parity of the integrand with respect to \(\theta\), we prefer to simplify them as \(2 \int_{0}^{\pi/2} \ldots d \theta\).

Lemma 3

Using the integral form of the Bessel function Lemma 2 into Lemma 1 gives

(82)\[\frac{2}{\pi} \Re \int_{0}^{\pi/2} \frac{1}{\zeta(\theta)} \, \mathrm{d} \theta = - \frac{1}{\sqrt{r^2 + z^2}}.\]

which can also be found in [Del89]

The above lemma allows to retrieve the expression of the Green function found e.g. in [BD15]:

(83)\[\begin{split}\mathcal{G}^-(r, z) & = \frac{4}{\pi} \Re \left( \int_{0}^{\pi/2} \left( e^{\zeta(r, z, \theta)} \left[ E_1(\zeta(r, z, \theta)) + i \pi \right] - \frac{1}{\zeta(r, z, \theta)} \right) \, \mathrm{d} \theta \right) \\ & \qquad \qquad \qquad \qquad + 4 i \Re \left( \int^{\pi/2}_{0} e^{\zeta (r, z, \theta)} \, \mathrm{d} \theta \right)\end{split}\]

Alternative formulations

Lemma 4

The zeroth order Bessel function of the first kind \(J_0\) and the Struve function \(H_0\) are such that

(84)\[\begin{split}J_0(r) & = \frac{2}{\pi} \int_{0}^{\pi/2} \cos(r\cos(\theta)) \, \mathrm{d} \theta \\ H_0(r) & = \frac{2}{\pi} \int_{0}^{\pi/2} \sin(r\cos(\theta)) \, \mathrm{d} \theta \\\end{split}\]

hence

(85)\[2 \int_{0}^{\pi/2} i e^{\zeta} \, \mathrm{d} \theta = \pi e^z \left(- H_0(r) + i J_0(r) \right)\]

The function \(\mathcal{G}\) can also be rewritten as

(86)\[\begin{split}\mathcal{G}(r, z) & = \frac{1}{\sqrt{r^2 + z^2}} + \frac{2}{\pi} \int^{\pi/2}_{-\pi/2} \Re \left( e^\zeta E_1(\zeta) \right) \, \mathrm{d} \theta + 2 \int^{\pi/2}_{-\pi/2} i e^{\zeta (r, z, \theta)} \, \mathrm{d} \theta \\ & = \frac{1}{\sqrt{r^2 + z^2}} + \frac{2}{\pi} \int^{\pi/2}_{-\pi/2} \Re \left( e^\zeta E_1(\zeta) \right) \, \mathrm{d} \theta + 2 \pi e^z \left( - H_0(r) + i J_0(r) \right)\end{split}\]

Noblesse [N82] splits the function \(\mathcal{G}\) into a near field term \(N\) and a wave field \(W\) such that

(87)\[\begin{split}N(r, z) & = \frac{1}{\sqrt{r^2 + z^2}} + \frac{2}{\pi} \int^{\pi/2}_{-\pi/2} \Re \left( e^\zeta E_1(\zeta) \right) \, \mathrm{d} \theta \\ W(r, z) & = 2 \pi e^z \left( - H_0(r) + i J_0(r) \right)\end{split}\]

Note that \(E_1\), \(J_0\) and \(H_0\) are available for instance in the Scipy library.

The literature also contains alternative formulation with another variant of the \(\zeta\) function as seen in the lemma below.

Lemma 5

For any function \(f\), the following two formulations of the integral are equivalent [Del89]:

(88)\[\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} f \left(\zeta(\theta) \right) \mathrm{d} \theta = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} f \left(\tilde{\zeta}(\theta) \right) \mathrm{d} \theta\]

where \(\zeta\) is defined in (79) and \(\tilde{\zeta}\) is defined as

(89)\[\tilde{\zeta} (\theta) = k \left( x_3 + \xi_3 + i \left( (x_1 - \xi_1) \cos\theta + (x_2 - \xi_2) \sin\theta \right) \right).\]

Proof. Let us note that:

\begin{align*} (x_1 - \xi_1) \cos(\theta) + (x_2 - \xi_2) \sin(\theta) & = \Re \left( \left( x_1 - \xi_1 + i (x_2 - \xi_2) \right) e^{-i \theta} \right) \\ & = \Re \left( r e^{i (\alpha - \theta)} \right) \\ & = r \cos \left( \alpha - \theta \right) \\ \end{align*}

where \(r\) and \(\alpha\) are defined by

\[ r e^{i \alpha} = (x_1 - \xi_1) + i (x_2 - \xi_2). \]

Finally note that:

\[ \int_{-\frac{\pi}{2}-\alpha}^{\frac{\pi}{2}-\alpha} f \left(\zeta(\theta) \right) \mathrm{d} \theta = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} f \left(\zeta(\theta) \right) \mathrm{d} \theta \]

Gradient of the Green function

The gradient of the Green function with respect to its first variable (that is \(x\)) can be written as

(93)\[\begin{split}-4 \pi \nabla_1 G(x, \xi) = - \frac{x - \xi}{\|x - \xi\|^3} + k \begin{pmatrix} \frac{\partial r}{\partial x_1} \frac{\partial \mathcal{G}}{\partial r} \\ \frac{\partial r}{\partial x_2} \frac{\partial \mathcal{G}}{\partial r} \\ \frac{\partial z}{\partial x_3} \frac{\partial \mathcal{G}}{\partial z} \end{pmatrix}\end{split}\]

with

(94)\[\begin{split}\frac{\partial r}{\partial x_1} & = k^2 \frac{x_1 - \xi_1}{r} \\ \frac{\partial r}{\partial x_2} & = k^2 \frac{x_2 - \xi_2}{r} \\ \frac{\partial z}{\partial x_3} & = k.\end{split}\]

or equivalently using \(\mathcal G^+\) from (78):

(95)\[\begin{split}-4 \pi \nabla_1 G(x, \xi) = - \frac{x - \xi}{\|x - \xi\|^3} - \frac{x - S_0(\xi)}{\|x - S_0(\xi)\|^3} + k \begin{pmatrix} \frac{\partial r}{\partial x_1} \frac{\partial \mathcal{G}^+}{\partial r} \\ \frac{\partial r}{\partial x_2} \frac{\partial \mathcal{G}^+}{\partial r} \\ \frac{\partial z}{\partial x_3} \frac{\partial \mathcal{G}^+}{\partial z} \end{pmatrix}\end{split}\]

The derivative of \(\mathcal G^+\) with respect to \(z\) can be handled with Property 2. Let us now focus on the derivative with respect to \(r\). From the definition of the exponential integral, we have

(96)\[\frac{d}{d\zeta}\left(e^\zeta \left( E_1(\zeta) + i \pi \right)\right) = e^\zeta \left( E_1(\zeta) + i \pi \right) - 1/\zeta\]

hence

(97)\[\begin{split}\frac{\partial \mathcal{G}^+}{\partial r} = & \frac{4}{\pi} \Re \left( \int_{0}^{\pi/2} \frac{\partial \zeta}{\partial r} \left( e^\zeta \left( E_1(\zeta) + i \pi \right) - \frac{1}{\zeta} \right) \, \mathrm{d}\theta \right) \\ & \qquad \qquad \qquad \qquad + 4 i \Re \left( \int^{\pi/2}_{0} \frac{\partial \zeta}{\partial r} e^{\zeta} \, \mathrm{d} \theta \right)\end{split}\]

that is

(98)\[\begin{split}\frac{\partial \mathcal{G}^+}{\partial r} = & \frac{4}{\pi} \Re \left( \int_{0}^{\pi/2} i \cos(\theta) \left( e^\zeta \left( E_1(\zeta) + i \pi \right) - \frac{1}{\zeta} \right) \, \mathrm{d}\theta \right) \\ & \qquad \qquad \qquad \qquad + 4 i \Re \left( \int^{\pi/2}_{0} i \cos(\theta) e^{\zeta} \, \mathrm{d} \theta \right)\end{split}\]

Note

The derivative of \(G\) with respect to \(x_1\) and \(x_2\) are antisymmetric in the sense of

\[ \frac{\partial G}{\partial x_1} (\xi, x) = - \frac{\partial G}{\partial x_1}(x, \xi). \]

Its derivative with respect to \(x_3\) is symmetric in infinite depth.

In finite depth, some terms of the derivative with respect to \(x_3\) are symmetric and some are antisymmetric.

Higher order derivatives

From Property 2, one has

(100)\[\begin{split}\frac{\partial \mathcal{G}}{\partial z} &= \mathcal{G}(r, z) + \left( 1 + \frac{\partial}{\partial z} \right) \frac{1}{\sqrt{r^2 + z^2}} \\ \frac{\partial^2 \mathcal{G}}{\partial z \partial r} &= \frac{\partial \mathcal{G}}{\partial r} + \left( \frac{\partial}{\partial r} + \frac{\partial^2}{\partial z \partial r} \right) \frac{1}{\sqrt{r^2 + z^2}}\end{split}\]

and

(101)\[\begin{split}\frac{\partial^2 \mathcal{G}}{\partial z^2} &= \mathcal{G}(r, z) + \left( 1 + 2 \frac{\partial}{\partial z} + \frac{\partial^2}{\partial z^2} \right)\frac{1}{\sqrt{r^2 + z^2}} \\ &= \mathcal{G}(r, z) + \frac{1}{\sqrt{r^2 + z^2}} - 2 \frac{z}{(r^2 + z^2)^{3/2}} - \frac{r^2 - 2 z^2}{(r^2 + z^2)^{5/2}}\end{split}\]

Since the Green function is solution of the Laplace equation, it follows that

(102)\[\frac{\partial^2 \mathcal{G}}{\partial r^2} + \frac{1}{r} \frac{\partial \mathcal{G}}{\partial r} + \frac{\partial^2 \mathcal{G}}{\partial z^2} = 0\]

then

(103)\[\begin{split}\frac{\partial^2 \mathcal{G}}{\partial r^2} = - \frac{1}{r} \frac{\partial \mathcal{G}}{\partial r} - \mathcal{G} - \left( 1 + 2 \frac{\partial}{\partial z} + \frac{\partial^2}{\partial z^2} \right)\frac{1}{\sqrt{r^2 + z^2}} \\\end{split}\]

All higher order derivative can be expressed with the help of \(\mathcal{G}\) and \(\frac{\partial \mathcal{G}}{\partial r}\).

Note

The same derivation is done in e.g. [N20] using instead the function \(F = \mathcal{G} - \frac{1}{\sqrt{r^2 + z^2}}\) for which the expressions are slightly simpler.

Delhommeau’s method for evaluation

The current version of Capytaine can use either the low-frequency variant (78) or high-frequency variant (83) when evaluating the Green function and its integral over a panel. For this purpose, the following values needs to be computed:

(104)\[\begin{split}I_1(r, z) & = \frac{4}{\pi} \Re \left( \int^{\pi/2}_{0} e^\zeta \left( E_1(\zeta) + i \pi \right) \, \mathrm{d} \theta \right) \\ I_2(r, z) & = \frac{4}{\pi} \Re \left( \int^{\pi/2}_{0} \left( e^\zeta \left( E_1(\zeta) + i \pi \right) - \frac{1}{\zeta} \right) \, \mathrm{d} \theta \right) \\ I_3(r, z) & = 4 \Re \left( \int^{\pi/2}_{0} e^{\zeta} \, \mathrm{d} \theta \right) \\ I_4(r, z) & = \frac{4}{\pi} \Re \left( \int^{\pi/2}_{0} i \cos(\theta) \left( e^\zeta \left( E_1(\zeta) + i \pi \right) - \frac{1}{\zeta} \right) \, \mathrm{d} \theta \right) \\ I_5(r, z) & = 4 \Re \left( \int^{\pi/2}_{0} i \cos(\theta) e^{\zeta} \, \mathrm{d} \theta \right)\end{split}\]

The integral are computed using Simpson’s rule for around 1000 points between \(0\) and \(\pi/2\).

Then (78) and (98) reads

(105)\[\begin{split}\mathcal{G}^+(r, z) & = I_1(r, z) + i I_3(r, z) \\ \frac{\partial \mathcal{G}^+}{\partial r} & = I_4(r, z) + i I_5(r, z).\end{split}\]

and (83) reads (still using (98) for the derivative):

(106)\[\begin{split}\mathcal{G}^-(r, z) & = I_2(r, z) + i I_3(r, z) \\ \frac{\partial \mathcal{G}^+}{\partial r} & = I_4(r, z) + i I_5(r, z).\end{split}\]

To limit the computational cost of the evaluation of these integrals, they are precomputed for selected values of \(r\) and \(z\) and stored in a table. When evaluating the Green function, the values of the integrals are retrieved by interpolating the values in the tables.

For large values of \(r\) and \(z\), these integrals are asymptotically approximated by the following expressions:

(107)\[\begin{split}I_1(r, z) & \simeq - 2 e^z \sqrt{\frac{2\pi}{r}} \sin(r - \pi/4) + \frac{2 z}{(r^2 + z^2)^{3/2}} - \frac{2}{\sqrt{r^2 + z^2}} \\ I_2(r, z) & \simeq - 2 e^z \sqrt{\frac{2\pi}{r}} \sin(r - \pi/4) + \frac{2 z}{(r^2 + z^2)^{3/2}} \\ I_3(r, z) & \simeq 2 e^z \sqrt{\frac{2\pi}{r}} \cos(r - \pi/4) \\ I_4(r, z) & \simeq - 2 e^z \sqrt{\frac{2\pi}{r}} \left( \cos(r - \pi/4) - \frac{1}{2r} \sin(r-\pi/4) \right) + \frac{2 r}{(r^2 + z^2)^{3/2}} \\ I_5(r, z) & \simeq - 2 e^z \sqrt{\frac{2\pi}{r}} \left( \sin(r - \pi/4) + \frac{1}{2r} \cos(r - \pi/4) \right)\end{split}\]

Incorporating these asymptotic approximation in the expression of the Green function, one gets:

(108)\[\begin{split}\mathcal{G}(r, z) \simeq & -\frac{1}{\sqrt{r^2 + z^2}} - 2 e^z \sqrt{\frac{2\pi}{r}} \left(\sin(r - \pi/4) - i\cos(r - \pi/4)\right) \\ & \qquad\qquad\qquad\qquad + 2 \frac{z}{(r^2 + z^2)^{3/2}}\end{split}\]

Integration

TODO

As seen in Property 2, new reflected-Rankine-type terms might appear in the derivative of the Green wave term. By default, they are integrated with the same method used for the same numerical quadrature method as the rest of the wave term. The setting gf_singularities="low_freq_with_rankine_term" is an attempt to integrate them exactly using the same code as the main reflected Rankine term.