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

1. The fluid is inviscid.

2. The fluid is incompressible: \(\nabla \cdot u = 0\) with \(u\) the flow velocity.

3. The flow is irrotational: \(\nabla \times u = 0\).

4. The wave amplitude is small with respect to the wavelength.

5. The amplitude of the body motion is small with respect to its dimension.

6. The sea bottom is flat. The water depth is denoted \(h\).

Mathematical problem

The mass conservation equation can be rewritten as the Laplace equation

(14)\[\nabla^2 \phi = 0,\]

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:

(15)\[\phi = \mathrm{Re} \left( \Phi e^{-i \omega t} \right).\]

The partial differential equation (14) is completed with the following boundary conditions:

• linearized free surface boundary condition:

(16)\[g \frac{\partial \Phi}{\partial z} - \omega^2 \Phi = 0, \qquad \text{on } z = 0.\]

• no velocity boundary condition at the (flat) sea bottom:

(17)\[\frac{\partial \Phi}{\partial z} = 0, \qquad \text{on } z = -h.\]

• a given velocity \(u\) on the floating body surface \(\Gamma\):

(18)\[\nabla \Phi \cdot n = u \cdot n, \qquad \text{on } \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

(19)\[\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

(20)\[\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:

(21)\[\nabla^2 G(x; \xi) = \delta(\xi - x), \qquad \forall x,\]

where the \(\nabla\) is meant as the derivative with respect to \(x\).

The above equation is associated with the boundary condition (16) and (17), 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\):

(22)\[\Phi(x) = \iint_\Gamma \sigma(\xi) G(x; \xi) \, \mathrm{dS}(\xi)\]

for all point \(x\) in the fluid or on the hull of the floating body \(\Gamma\).

Note

There is a typo in equation (22) in [BD15].

The integral on the other boundaries of the domain is zero due to the properties of the Green function.

The differentiation of (22) differs depending whether \(x\) is in the bulk of the fluid or on the hull.

On the hull, one has [Del87]:

(23)\[\frac{\partial \Phi}{\partial n}(x) = (u \cdot n)(x) = \frac{\sigma(x)}{2} + \iint_\Gamma \sigma(\xi) \, (\nabla G(x; \xi) \cdot n) \, \mathrm{dS}(\xi).\]

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

(24)\[\frac{\partial \Phi}{\partial t}(x) = (u \cdot t)(x) = \iint_\Gamma \sigma(\xi) \, (\nabla G(x; \xi) \cdot t) \, \mathrm{dS}(\xi).\]

Finally, for \(x\) in the bulk of the fluid, one has

(25)\[\nabla \Phi(x) = u(x) = \iint_\Gamma \sigma(\xi) \, \nabla G(x; \xi) \, \mathrm{dS}(\xi).\]

Note

Dimensional analysis:

\(\Phi\) is in m²·s¯¹.

\(\sigma\) is in m·s¯¹.

\(G\) is in m¯¹.

TODO: discuss direct method as in [AN24]

Discretization

The equations (22) and (23) can be discretized using a collocation method. Considering a mesh of the surface of the floating body \(\Gamma = \cup_i \Gamma_i\):

(26)\[\begin{split}\Phi_i & = \Phi(x_i), \\ \sigma_i & = \sigma(x_i), \\ u_i & = (u \cdot n)(x_i) \\ S_{ij} & = \iint_{\Gamma_j} G(x_i, \xi) \mathrm{dS}(\xi), \\ K_{ij} & = \frac{\delta_{ij}}{2} + \iint_{\Gamma_j} \nabla_x G(x_i; \xi) \cdot n_i \, \mathrm{dS}(\xi),\end{split}\]

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:

(27)\[D_{ij} = \frac{\delta_{ij}}{2} + \iint_{\Gamma_j} \nabla_\xi G(x_i; \xi) \cdot n_j \, \mathrm{dS}(\xi).\]

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 (22) and (23):

(28)\[\Phi = S \sigma, \qquad u = K \sigma.\]

The resolution of the discrete problem with Nemoh consists of two main steps:

1. The evaluation of the coefficients of the complex-valued matrices \(S\) and \(K\)

2. 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. 2 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:

1. The magnitude \(U\) of the forward speed is small.

2. 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:

1. Doppler shift: The frequency used in the computation is replaced by the encounter frequency

(29)\[\omega_e = \omega - k U \cos (\beta)\]

where \(k\) is the wavenumber and \(\beta\) is the wave direction. For this purpose, the wave_direction parameter can be passed to radiation problem.

2. Normal velocity on hull: The boundary condition on the body radiating with a dof defined by the displacement \(\delta\!r(x, y, z)\) reads

(30)\[\frac{\partial \phi}{\partial n} = - i \omega_e \delta\!r \cdot n - U \frac{\partial \delta\! r}{\partial x} \cdot n\]

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.

3. Gradient of potential in pressure: The equation relating the potential to the pressure is updated as follows

(31)\[p = -\rho \left( -i \omega_e \phi + U \frac{\partial \phi}{\partial x} \right)\]

Similarly the relationship between the potential and the free surface elevation reads

(32)\[\eta = -\frac{1}{g} \left( -i \omega_e \phi + U \frac{\partial \phi}{\partial x} \right)\]

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. 3 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.

(33)\[F_j = - \int_\Gamma p(x) \, n(x) \cdot \delta\!r_j(x) \, dx = - i \omega \rho \int_\Gamma \Phi(x) \, n(x) \cdot \delta\!r_j(x) \, dx\]

where \(p = i \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_j\) is the local displacement of the hull of the degree of freedom \(j\).

For a single rigid body, the degrees of freedom reads:

Dof \(j\)	Local hull displacement
Surge	\(\delta\!r_j(x) = (1, 0, 0)\)
Sway	\(\delta\!r_j(x) = (0, 1, 0)\)
Heave	\(\delta\!r_j(x) = (0, 0, 1)\)
Roll	\(\delta\!r_j(x) = (1, 0, 0) \times (x-x_0, y-y_0, z-z_0)\)
Pitch	\(\delta\!r_j(x) = (0, 1, 0) \times (x-x_0, y-y_0, z-z_0)\)
Yaw	\(\delta\!r_j(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:

1. 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.

2. The diffraction forces \(F_{D}\), from the integration of the diffracted wave field (all bodies held fixed).

3. 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 \(j\) of the radiation force \(F_{R}\) is further rewritten as

(34)\[F_{R, j} = \sum_k \left[\omega^2 A_{jk} + i \omega B_{jk}\right] X_k\]

where \(A_{jk}\) is the added mass matrix, \(B_{jk}\) 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

(35)\[\begin{split}A_{jk} & = \frac{1}{\omega^2} \Re \left[ i \omega \rho \int_\Gamma \Phi_k(x) \, n(x) \cdot \delta \! r_j(x) \, dx \right] \\ & = - \frac{\rho}{\omega} \int_\Gamma \Im [\Phi_k(x)] \, n(x) \cdot \delta \! r_j(x) \, dx\end{split}\]

and

(36)\[\begin{split}B_{jk} & = \frac{1}{\omega} \Im \left[ i \omega \rho \int_\Gamma \Phi_k(x) \, n(x) \cdot \delta \! r_j(x) \, dx \right] \\ & = \rho \int_\Gamma \Re [\Phi_k(x)] \, n(x) \cdot \delta \! r_j(x) \, dx\end{split}\]

where \(\Phi_k\) is the potential field computed with the normal velocity on the hull \(\frac{\partial \Phi_k}{\partial n} = -i \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 \(j\) used in the integration of the force is the influenced_dof.

Note

From Green second identity

(37)\[\int_\Gamma \left[ \Phi_j \frac{\partial \Phi_k}{\partial n} - \Phi_k \frac{\partial \Phi_j}{\partial n}\right] dx = 0\]

one has, when using the definition of the normal velocity of the radiation problem above,

(38)\[\iint_{\Gamma} \Phi_j \; \delta\!r_k \cdot n = \iint_{\Gamma} \Phi_k \; \delta\!r_j \cdot n\]

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} = -i \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}{-i \omega}\).

It leads to the following definition of the added mass and radiation damping

(39)\[\begin{split}A_{jk} & = \frac{1}{\omega^2} \Re \left[ i \omega \rho \int_\Gamma (- i \omega \tilde \Phi_k(x)) \, n(x) \cdot \delta \! r_j(x) \, dx \right] \\ & = \rho \int_\Gamma \Re [\tilde \Phi_k(x)] \, n(x) \cdot \delta \! r_j(x) \, dx\end{split}\]

and

(40)\[\begin{split}B_{jk} & = \frac{1}{\omega} \Im \left[ i \omega \rho \int_\Gamma (- i \omega \tilde \Phi_k(x)) \, n(x) \cdot \delta \! r_j(x) \, dx \right] \\ & = \rho \omega \int_\Gamma \Im [\tilde \Phi_k(x)] \, n(x) \cdot \delta \! r_j(x) \, dx\end{split}\]

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 \(-i\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

(41)\[M_{jk} \ddot{x}_k + C_{jk} \dot{x}_k + K_{jk} x_k = F_j,\]

and in frequency domain, with the assumed time dependence \(x(t) = \mathrm{Re} \left( X e^{-i \omega t} \right)\),

(42)\[\left[-\omega^2 M_{jk} - i \omega C_{jk} + K_{jk}\right] X_k = F_j,\]

where \(M_{jk}\) is the inertia matrix, accounting for the mass distribution, \(C_{jk}\) is the mechanical damping matrix, \(K_{jk}\) is the stiffness matrix which comprises mechanical and hydrostatic effects, and \(F_j\) are generic external forces.

Note

The hydrostatic contribution to matrix \(K_{jk}\) accounts for a variation of hydrostatic force in direction \(j\) due to a unit motion in direction \(k\). It is a geometric property of the body.

As seen above, forces \(F_j\) can be decomposed as

(43)\[F_j = F_{FK, j} + F_{D, j} + F_{R, j}\]

The full system becomes

(44)\[\left[-\omega^2 (M_{jk} + A_{jk}) - i \omega (C_{jk} + B_{jk}) + K_{jk}\right] X_k = F_{FK, j} + F_{D, j}\]

that is

(45)\[H X = F_{ex}\]

where \(H\) denotes the following transfer function matrix

(46)\[H_{jk} = \left[-\omega^2 (M_{jk} + A_{jk}) - i \omega (C_{jk} + B_{jk}) + K_{jk}\right]\]

and \(F_{ex}\) denotes the excitation force.

(47)\[F_{ex, j} = F_{FK, j} + F_{D, j}.\]

The oscillation amplitude is obtained by solving the complex-valued linear system.

Note

Matrices \(A_{jk}\) and \(B_{jk}\) depend on \(\omega\), and so does \(H_{jk}\) and \(X_k\).

Free surface elevation

The potential at the reference surface \(z = 0\) can be connected to the free surface elevation by the dynamic condition

(48)\[\dfrac{\partial \phi}{\partial t} = - g \eta\]

which, in frequency domain, is

(49)\[\eta = \dfrac{j \omega}{g} \Phi\]

For a fully coupled problem (bodies free to oscillate, i.e. diffraction and radiation combined), the free surface elevation can be computed as

(50)\[\eta = \eta_{\text{incident}} + \eta_{\text{diffracted}} + \sum_k \eta_{\text{radiated}, k} X_k.\]

Far-field coefficients

Far from the floating body, the potential of an infinite depth problem can be written as

(51)\[\Phi(x) \sim - \sqrt{\frac{2 \pi k}{\rho(x)}} e^{k x_3} H(\theta(x)) e^{i k \rho(x)} e^{i \pi/4} + O\left(\frac{1}{\rho(x)}\right)\]

where \(\rho(x)\) and \(\theta(x)\) are the cylindrical coordinates of \(x\) and \(H\) is the so-called Kochin function which can be computed as

(52)\[H(\theta) = \frac{1}{4 \pi} \int_\Gamma \sigma(\xi) e^{k \xi_3} e^{- i k (\xi_1 \cos \theta + \xi_2 \sin \theta)} \, d \xi\]

Proof. Let us first note, following [Del87] (p.79), that for large \(\rho(x)\):

\begin{align*} r & = \sqrt{(x_1 - \xi_1)^2 + (x_2 - \xi_2)^2} \\ & = \sqrt{\rho(x)^2 + \rho(\xi)^2 - 2 \rho(x) \rho(\xi) \cos(\theta(x) - \theta(\xi))} \\ & \sim \rho(x) - \rho(\xi) \cos(\theta(x) - \theta(\xi)) \\ & \sim \rho(x) - (\xi_1 \cos\theta(x) + \xi_2 \sin\theta(x)) \end{align*}

Then the asymptotic behavior of the infinite depth Green function (see Green function in infinite depth) reads

(54)\[- 4 \pi G(x, \xi) \sim 2 k \sqrt{\frac{\pi}{k r}} e^{k z + i k r}(1+i) + O\left(\frac{1}{r}\right)\]

which can be rewritten as

(55)\[- 4 \pi G(x, \xi) \sim \sqrt{\frac{2 \pi k}{\rho(x)}} e^{k z + i k (\rho(x) - (\xi_1 \cos\theta(x) + \xi_2 \sin\theta(x))} e^{i\pi/4} + O\left(\frac{1}{\rho(x)}\right)\]

Incorporating into the source representation (22), we have

(56)\[\Phi \sim - e^{k x_3} e^{ik\rho(x)} e^{i \pi/4} \sqrt{\frac{2\pi k}{\rho(x)}} \underbrace{\frac{1}{4\pi} \int_\Gamma \sigma(\xi) e^{k \xi_3} e^{-ik(\xi_1 \cos \theta + \xi_2 \sin \theta))} d\xi}_{H(\theta)}\]

Capytaine computes the Kochin function in finite depth as

(57)\[H(\theta) = \frac{1}{4 \pi} \int_\Gamma \sigma(\xi) \frac{\cosh (k (\xi_3+h))}{\cosh (k h)} e^{- i k (\xi_1 \cos \theta + \xi_2 \sin \theta)} \, d \xi\]

Note that other work might define the Kochin function slightly differently, with a different normalization.

Mean drift force

Far field formulation:

The expression of the mean drift force for the degrees of freedom Surge and Sway is the following:

(58)\[\begin{split}\left\langle{\begin{matrix} F_x \\ F_y \end{matrix}}\right\rangle = -2 \pi \rho \omega\binom{\cos \beta}{\sin \beta} \Im( H(\beta)) -2 \pi \rho \frac{k\left(k_0 h\right)^2}{h\left[\left(k h\right)^2-\left(k_0 h\right)^2+k_0 h\right]} \int_0^{2 \pi}|H(\theta)|^2\binom{\cos \theta}{\sin \theta} d \theta\end{split}\]

where \(\beta\) is the wave direction, \(H\) the Kochin function, \(k\) the wavenumber, \(h\) the water depth and \(k_0\) the deep water wavenumber.

Note

The coefficient in front of the integral above can be very large at high frequency and can make the result very sensitive to small numerical inaccuracies in the Kochin function. Unfortunately, numerical incurracies can be common at high frequency, that is when the wavelength is small with respect to the mesh resolution. In other words, mesh convergence should be checked carefully for the mean drift force at high frequency.

Here is the expression for the Yaw moment:

(59)\[\left\langle{M_z}\right\rangle = 2 \pi \frac{\rho \omega}{k}\Re (\dot H(\beta)) - \frac{2 \pi \rho (k_0h)^2}{h[(kh)^2 - (k_0h)^2 + k_0h]}\Im (\int_0^{2 \pi} H(\theta)^* \dot H(\theta) \mathrm{d} \theta )\]

The Kochin function has to be rebuild from the contributions of all the radiation problems and the diffraction problem:

(60)\[H(\theta) = e^{i\frac{\pi}{2}} ( H_D(\theta) + \sum_{k=1}^6 X_k H_{R_k} (\theta))\]

where \(X_k\) is the motion RAO of the body corresponding to the degree of freedom \(k\), \(H_{R_k}\) is the Kochin function associated with the radiated potential of degree of freedom \(k\) and \(H_{D}\) is the Kochin function associated with the diffracted potential.