Theory manual¶
This document presents several aspects of the theoretical background for Nemoh and Capytaine. We refer also to [Del87] [Del89] [Del93] [BD15] and [AD18].
Warning
This document is a work in progress. It is incomplete and might contain errors.
Linear boundary value problem¶
Hypotheses¶
The fluid is inviscid.
The fluid is incompressible: \(\nabla \cdot u = 0\) with \(u\) the flow velocity.
The flow is irrotational: \(\nabla \times u = 0\).
The wave amplitude is small with respect to the wavelength.
The amplitude of the body motion is small with respect to its dimension.
The sea bottom is flat. The water depth is denoted \(h\).
Mathematical problem¶
The mass conservation equation can be rewritten as the Laplace equation
where \(\phi\) is the velocity potential defined as \(u = \nabla \phi\).
Since the problem is linear, we look for a solution in the frequency domain:
The partial differential equation (13) is completed with the following boundary conditions:
linearized free surface boundary condition:
no velocity boundary condition at the (flat) sea bottom:
a given velocity \(u\) on the floating body surface \(\Gamma\):
where \(n\) denotes the normal vector at the surface of the floating body.
The normal velocity on the floating body surface is the input of the problem. It depends on the type of problem:
- Radiation problem:
For the radiation problem, the normal velocity on the body surface corresponds to the motion of the body along one of its degrees of freedom. The resolution of the Laplace problem allows to derive the added mass and the radiation damping associated with this degree of freedom (see also Post-processing).
- Diffraction problem:
For the diffraction problem, the velocity on the floating body is given by the velocity of Airy’s wave field. Once the problem has been solved, the linear Froude-Krylov force is computed by the integration of the pressure (\(p = i \rho \omega \Phi\)) on the floating body (see also Post-processing).
The incoming Airy’s wave fields is given by
(18)¶\[\Phi_0 = - i \frac{g}{\omega} \frac{\cosh (k (z+h))}{\cosh (k h)} e^{i k (x \cos \beta + y \sin \beta)}\]in finite depth, where the wave number \(k\) is defined by the dispersion relation \(\omega^2 = k g \tanh (k h)\), and by
(19)¶\[\Phi_0 = - i \frac{g}{\omega} e^{k z} e^{i k (x \cos \beta + y \sin \beta)}\]in infinite depth, where the wave number \(k\) is defined by \(\omega^2 = k g\).
In the above equations, \(\beta\) is the angle of the incoming wave. The angle \(\beta = 0\) corresponds to waves propagating in the \(x\) direction from \(x=-\infty\) to \(x=+\infty\). The angle \(\beta = \pi/2\) corresponds to waves propagating in the \(y\) direction from \(y=-\infty\) to \(y=+\infty\).
Integral problem¶
The partial differential equation can be rewritten as a boundary integral problem. Let us introduce the Green function \(G(\xi, \cdot)\), which is solution of the partial differential equation:
where the \(\nabla\) is meant as the derivative with respect to \(x\).
The above equation is associated with the boundary condition (15) and (16), where \(\xi\) is a given point in the domain and \(\delta\) is the Dirac distribution.
With the help of this Green function \(G\), the potential of the surface of the floating body \(\Gamma\) can be rewritten as a function of a source distribution \(\sigma\):
for all point \(x\) in the fluid or on the hull of the floating body \(\Gamma\).
The integral on the other boundaries of the domain is zero due to the properties of the Green function.
The differentiation of (21) differs depending whether \(x\) is in the bulk of the fluid or on the hull.
On the hull, one has [Del87]:
where \(x\) is a point on \(\Gamma\) and \(n\) is the vector normal to \(\Gamma\) in \(x\). For any vector \(t\) tangential to \(\Gamma\) at \(x\), one has
Finally, for \(x\) in the bulk of the fluid, one has
Note
Dimensional analysis:
\(\Phi\) is in m²·s¯¹.
\(\sigma\) is in m·s¯¹.
\(G\) is in m¯¹.
Expression of the Green function¶
In infinite depth¶
The integral problem above relates the potential \(\Phi\) to the normal velocity \(u \cdot n\) via the Green function \(G\). Let us know discuss the evaluation of this function for an infinite water depth. See also [X18].
Green function¶
The infinite depth Green function takes the following form
The function \(G\) is symmetric in the sense of
The first term of \(G\) is the usual Green function for the 3D Laplace equation without our specific boundary conditions. The \(\mathcal{G}\) term is complex-valued and it is introduced to satisfy the boundary conditions (15).
Introducing the dimensionless variables \(r = k \sqrt{(\xi_1 - x_1)^2 + (\xi_2 - x_2)^2}\) and \(z = k (x_3 + \xi_3)\), this term reads
where
where \(E_1\) is the first exponential integral, defined as
and
The first term of (27) is actually a Rankine-type singularity similar to the first term of (25), except that one of the point has been reflected through the free surface.
Variants of the formulation¶
The above lemma allows to retrieve the expression of the Green function found e.g. in [BD15]:
(Note the minus sign in front of the first term.)
The zeroth order Bessel function of the first kind \(J_0\) and the Struve function \(H_0\) are such that
hence
The function \(\mathcal{G}\) can also be rewritten as
Noblesse [N82] splits the function \(\mathcal{G}\) into a near field term \(N\) and a wave field \(W\) such that
Note that \(E_1\), \(J_0\) and \(H_0\) are available for instance in the Scipy library.
For any function \(f\), the following two formulations of the integral are equivalent [Del89]:
where \(\zeta\) is defined in (30) and \(\tilde{\zeta}\) is defined as
where \(r\) and \(\alpha\) are defined by
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 can be written as
with
and, using the identity \(J'(\zeta) = J(\zeta) - 1/\zeta\),
and
that is, using Lemma {number}
Note
The derivative of \(G\) with respect to \(x_1\) and \(x_2\) are antisymmetric in the sense of
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 derivative¶
From (46), one has
and
Since the Green function is solution of the Laplace equation, it follows that
then
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 computation¶
The current version of Capytaine can use either the low-frequency variant (27) or high-frequency variant (32) when evaluating the Green function and its integral over a panel. For this purpose, the following values needs to be computed:
then (27) reads
and (32) reads
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:
Incorporating these asymptotic approximation in the expression of the Green function, one gets:
Integration¶
TODO
As seen in (46), 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.
In finite depth¶
Green function¶
TODO
Gradient of the Green function¶
TODO
Symmetries¶
The first term of (25) is invariant under all rotations and translations, whereas the other terms are 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).
Discretization¶
The equations (21) and (22) can be discretized using a collocation method. Considering a mesh of the surface of the floating body \(\Gamma = \cup_i \Gamma_i\):
where for all \(i\), \(x_i\) is the center of the face \(\Gamma_i\) and \(n_i\) is its normal vector. Each element of the matrices \(S\) and \(K\) can be seen as the interaction between two faces of the mesh.
Note
\(K\) should not be confused with the similar matrix \(D\) defined as:
Note that the derivation of \(G\) is done with respect to a different variable.
The matrix \(D\) is used in the direct boundary integral equation, as e.g. in HAMS [Liu19]. In the mathematical literature, \(D\) is also referred to as the double layer operator and \(K\) as the adjoint double layer operator.
The matrices \(S\) and \(K\) relates the vectors \(\Phi\), \(u\) and \(\sigma\) through the following approximations of (21) and (22):
The resolution of the discrete problem with Nemoh consists of two main steps:
The evaluation of the coefficients of the complex-valued matrices \(S\) and \(K\)
The resolution of the complex-valued linear problem \(K \sigma = u\).
Once \(\sigma\) has been computed, \(\Phi\) can be easily deduced. Then other magnitudes such as the Froude-Krylov forces or the added mass can be derived.
flowchart TD;
h[Water depth] --> gf(Assembling matrices);
ω[Wave frequency ω] --> gf(Assembling matrices);
m[Mesh] --> gf;
gf -- K matrix --> ls(Linear solver);
un[Normal velocity on hull] --> ls;
gf -- S matrix --> mvp(Matrix vector product);
ls -- sources distribution σ --> mvp;
mvp -- potential distribution Φ --> int("Integrate on mesh");
m --> int;
int --> f["Hydrodynamic forces\n(aka added mass and radiation damping)"]
classDef input fill:#DDDDDD,color:#333333,stroke:#444444
classDef step fill:#88BBBB,color:#003333,stroke:#226666
classDef output fill:#FFE3AA,color:#553900,stroke:#AA8439
class ω,m,un,h input
class gf,ls,mvp,int step
class f output
Fig. 1 A simplified flowchart of the internals of Capytaine solver¶
Problem with forward speed¶
We refer to [D22] for a detailed description of the theory behind the approximate forward speed model used in Capytaine.
It relies on the following hypotheses:
The magnitude \(U\) of the forward speed is small.
The body is thin enough, such that the flow around the body assuming a rigid free surface (also called double-body flow) can be approximated by \(\overrightarrow{u} = (-U, 0, 0)\) in the reference frame of the body.
Then, the following modification are done to the solver to take forward speed into account:
Doppler shift: The frequency used in the computation is replaced by the encounter frequency
where \(k\) is the wavenumber and \(\beta\) is the wave direction.
For this purpose, the wave_direction parameter can be passed to radiation problem.
Normal velocity on hull: The boundary condition on the body radiating with a dof defined by the displacement \(\delta\!r(x, y, z)\) reads
The above relationship has currently only been implemented for the six dofs of single rigid bodies, as follows
Dof |
\(\delta \! r\) |
\(\frac{\partial \delta\! r}{\partial x}\) |
|---|---|---|
Surge |
\((1, 0, 0)\) |
\((0, 0, 0)\) |
Sway |
\((0, 1, 0)\) |
\((0, 0, 0)\) |
Heave |
\((0, 0, 1)\) |
\((0, 0, 0)\) |
Roll |
\((0, -z, y)\) |
\((0, 0, 0)\) |
Pitch |
\((-z, 0, x)\) |
\((0, 0, 1)\) |
Yaw |
\((y, -x, 0)\) |
\((0, -1, 0)\) |
In other words, the supplementary term is zero except for pitch and yaw.
Gradient of potential in pressure: The equation relating the potential to the pressure is updated as follows
Similarly the relationship between the potential and the free surface elevation reads
The computation of \(\frac{\partial \phi}{\partial x}\) makes the problems with forward speed typically 50% slower that problems without.
The overall workflow with forward speed thus looks as follows.
graph TD
h[Water depth] --> gf(Assembling matrices);
omega[Wave frequency ω] --> doppler(Doppler shift);
fs[Forward speed U] --> doppler(Doppler shift);
fs --> un;
doppler -- Encounter frequency --> gf(Assembling matrices)
m[Mesh] --> gf;
gf -- K matrix --> ls(Linear solver);
dof[Degree of freedom] --> un(Normal velocity on hull);
un -- RHS of linear problem --> ls;
gf -- S matrix --> mvp(Matrix vector product);
ls -- sources distribution σ --> mvp;
ls -- sources distribution σ --> grad;
gf -- extended K matrix --> grad;
mvp -- potential distribution Φ --> int("Integrate pressure on mesh");
grad(Matrix vector product) -- gradient of Φ --> int;
m --> int;
fs --> int;
int --> f["Hydrodynamic forces"]
classDef input fill:#DDDDDD,color:#333333,stroke:#444444
classDef step fill:#88BBBB,color:#003333,stroke:#226666
classDef newstep fill:#FFAAAA,color:#550000,stroke:#113939
classDef output fill:#FFE3AA,color:#553900,stroke:#AA8439
class fs,omega,m,h,dof input
class doppler,un,grad newstep
class gf,ls,mvp,int step
class f output
Fig. 2 A simplified flowchart of the internals of Capytaine solver with forward speed, where red boxes are the supplementary steps introduced by forward speed.¶
Post-processing¶
Forces on body surfaces¶
Forces acting on body surfaces are computed by integration of the pressure field.
where \(p = j \omega \rho \Phi\) stands for the complex-valued pressure fields in frequency-domain, \(n\) is the normal vector on the hull \(\Gamma\) (oriented towards the fluid in Capytaine, see Conventions and differences to other codes) and \(\delta\!r_i\) is the local displacement of the hull of the degree of freedom \(i\).
For a single rigid body, the degrees of freedom reads:
Dof |
Local hull displacement |
|---|---|
Surge |
\(\delta\!r(x) = (1, 0, 0)\) |
Sway |
\(\delta\!r(x) = (0, 1, 0)\) |
Heave |
\(\delta\!r(x) = (0, 0, 1)\) |
Roll |
\(\delta\!r(x) = (1, 0, 0) \times (x-x_0, y-y_0, z-z_0)\) |
Pitch |
\(\delta\!r(x) = (0, 1, 0) \times (x-x_0, y-y_0, z-z_0)\) |
Yaw |
\(\delta\!r(x) = (0, 0, 1) \times (x-x_0, y-y_0, z-z_0)\) |
where \((x_0, y_0, z_0)\) is the rotation center and \(\times\) denotes the cross product.
The potential field can be decomposed into three contributions, and so does the resulting force:
The Froude-Krylov forces \(F_{FK}\), from the integration of the incident wave field pressure (incoming plane waves). In Capytaine, the incident wave pressure can be retrieved with the
airy_wave_pressure()function.The diffraction forces \(F_{D}\), from the integration of the diffracted wave field (all bodies held fixed).
The radiation forces \(F_{R}\), which is itself a linear combination of the forces exerted by the fluid on the body in response to a motion of each degree of freedom.
The component \(i\) of the radiation force \(F_{R}\) is further rewritten as
where \(A_{ik}\) is the added mass matrix, \(B_{ik}\) is the radiation damping matrix and \(X_k\) is the amplitude of the motion of the body along the degree of freedom \(k\).
In other words, one has
and
where \(\Phi_k\) is the potential field computed with the normal velocity on the hull \(\frac{\partial \Phi_k}{\partial n} = -j \omega \delta \! r_k \cdot n\).
In Capytaine’s wording, the degree of freedom \(k\) defining the normal velocity on the hull is called radiating_dof, while the degree of freedom \(i\) used in the integration of the force is the influenced_dof.
Note
From Green second identity
one has, when using the definition of the normal velocity of the radiation problem above,
from which we can deduce the symmetry of the added mass matrix and the radiation dampings matrix.
Note
As an alternative to \(\frac{\partial \Phi_k}{\partial n} = -j \omega \delta \! r_k \cdot n\), some software such as the version 1 of Capytaine use \(\frac{\partial \tilde \Phi_k}{\partial n} = \delta \! r_k \cdot n\), that is \(\tilde \Phi_k = \frac{\Phi_k}{-j \omega}\).
It leads to the following definition of the added mass and radiation damping
and
This form is convenient since the all the \(\omega\) in the expression of the added mass disappears, which make it possible to compute the value of the added mass at frequency such as zero or infinity.
However, the implementation of \(\tilde \Phi\) in version 1 of Capytaine was not consistent with the use of \(\Phi\) for diffraction problem and it was easy to forget the missing \(-j\omega\) for some post-processing of \(\tilde \Phi\) for radiation problems.
In version 2.0 of Capytaine, \(\Phi\) is used everywhere instead of \(\tilde \Phi\). Since version 2.1, another method has been implemented to take into account the cancelling of the \(\omega\) in the expression of the added mass allowing to compute the added mass at zero and infinite frequency.
Dynamic coupling and impedance¶
Consider a body or a system of bodies. The general linear equation of motion can be expressed in time domain as
and in frequency domain, with the assumed time dependence \(x(t) = \mathrm{Re} \left( X e^{-j \omega t} \right)\),
where \(M_{ij}\) is the inertia matrix, accounting for the mass distribution, \(C_{ij}\) is the mechanical damping matrix, \(K_{ij}\) is the stiffness matrix which comprises mechanical and hydrostatic effects, and \(F_i\) are generic external forces.
Note
The hydrostatic contribution to matrix \(K_{ij}\) accounts for a variation of hydrostatic force in direction \(i\) due to a unit motion in direction \(j\). It is a geometric property of the body.
As seen above, forces \(F_i\) can be decomposed as
The full system becomes
that is
where \(H\) denotes the following transfer function matrix
and \(F_{ex}\) denotes the excitation force.
The oscillation amplitude is obtained by solving the complex-valued linear system.
Note
Matrices \(A_{ij}\) and \(B_{ij}\) depend on \(\omega\), and so does \(H_{ij}\) and \(X_j\).
Free surface elevation¶
The potential at the reference surface \(z = 0\) can be connected to the free surface elevation by the dynamic condition
which, in frequency domain, is
For a fully coupled problem (bodies free to oscillate, i.e. diffraction and radiation combined), the free surface elevation can be computed as
Far-field coefficients¶
TODO