============= 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. .. contents:: Contents Linear boundary value problem ============================= Hypotheses ---------- 1. The fluid is inviscid. 2. The fluid is incompressible: :math:`\nabla \cdot u = 0` with :math:`u` the flow velocity. 3. The flow is irrotational: :math:`\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 :math:`h`. Mathematical problem -------------------- The mass conservation equation can be rewritten as the Laplace equation .. math:: \nabla^2 \phi = 0, :label: laplace where :math:`\phi` is the velocity potential defined as :math:`u = \nabla \phi`. Since the problem is linear, we look for a solution in the frequency domain: .. math:: \phi = \mathrm{Re} \left( \Phi e^{-i \omega t} \right). The partial differential equation :eq:`laplace` is completed with the following boundary conditions: * linearized free surface boundary condition: .. math:: g \frac{\partial \Phi}{\partial z} - \omega^2 \Phi = 0, \qquad \text{on } z = 0. :label: bc_fs * no velocity boundary condition at the (flat) sea bottom: .. math:: \frac{\partial \Phi}{\partial z} = 0, \qquad \text{on } z = -h. :label: bc_bottom * a given velocity :math:`u` on the floating body surface :math:`\Gamma`: .. math:: \nabla \Phi \cdot n = u \cdot n, \qquad \text{on } \Gamma, :label: bc_body where :math:`n` denotes the normal vector at the surface of the floating body. .. * in the far field, .. math:: \sqrt{R} \left( \frac{\partial \Phi}{\partial R} - i k \right) \left( \Phi - Phi_0 \right) \rightarrow 0, \qquad \text{when } R \rightarrow \infty, 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 (:math:`p = i \rho \omega \Phi`) on the floating body (see also Post-processing_). The incoming Airy's wave fields is given by .. math:: \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 :math:`k` is defined by the dispersion relation :math:`\omega^2 = k g \tanh (k h)`, and by .. math:: \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 :math:`k` is defined by :math:`\omega^2 = k g`. In the above equations, :math:`\beta` is the angle of the incoming wave. The angle :math:`\beta = 0` corresponds to waves propagating in the :math:`x` direction from :math:`x=-\infty` to :math:`x=+\infty`. The angle :math:`\beta = \pi/2` corresponds to waves propagating in the :math:`y` direction from :math:`y=-\infty` to :math:`y=+\infty`. Integral problem ---------------- The partial differential equation can be rewritten as a boundary integral problem. Let us introduce the Green function :math:`G(\xi, \cdot)`, which is solution of the partial differential equation: .. math:: \nabla^2 G(x; \xi) = \delta(\xi - x), \qquad \forall x, where the :math:`\nabla` is meant as the derivative with respect to :math:`x`. The above equation is associated with the boundary condition :eq:`bc_fs` and :eq:`bc_bottom`, where :math:`\xi` is a given point in the domain and :math:`\delta` is the Dirac distribution. With the help of this Green function :math:`G`, the potential of the surface of the floating body :math:`\Gamma` can be rewritten as a function of a source distribution :math:`\sigma`: .. math:: \Phi(x) = \iint_\Gamma \sigma(\xi) G(x; \xi) \, \mathrm{dS}(\xi) :label: potential_representation for all point :math:`x` in the fluid or on the hull of the floating body :math:`\Gamma`. .. note:: There is a typo in equation :eq:`potential_representation` in [BD15]_. The integral on the other boundaries of the domain is zero due to the properties of the Green function. The differentiation of :eq:`potential_representation` differs depending whether :math:`x` is in the bulk of the fluid or on the hull. On the hull, one has [Del87]_: .. math:: \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). :label: normal_velocity_on_hull_representation where :math:`x` is a point on :math:`\Gamma` and :math:`n` is the vector normal to :math:`\Gamma` in :math:`x`. For any vector :math:`t` tangential to :math:`\Gamma` at :math:`x`, one has .. math:: \frac{\partial \Phi}{\partial t}(x) = (u \cdot t)(x) = \iint_\Gamma \sigma(\xi) \, (\nabla G(x; \xi) \cdot t) \, \mathrm{dS}(\xi). :label: tangential_velocity_on_hull_representation Finally, for :math:`x` in the bulk of the fluid, one has .. math:: \nabla \Phi(x) = u(x) = \iint_\Gamma \sigma(\xi) \, \nabla G(x; \xi) \, \mathrm{dS}(\xi). :label: velocity_in_bulk_representation .. note:: Dimensional analysis: :math:`\Phi` is in m²·s¯¹. :math:`\sigma` is in m·s¯¹. :math:`G` is in m¯¹. TODO: discuss direct method as in [AN24]_ Discretization ============== The equations :eq:`potential_representation` and :eq:`normal_velocity_on_hull_representation` can be discretized using a collocation method. Considering a mesh of the surface of the floating body :math:`\Gamma = \cup_i \Gamma_i`: .. math:: \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), where for all :math:`i`, :math:`x_i` is the center of the face :math:`\Gamma_i` and :math:`n_i` is its normal vector. Each element of the matrices :math:`S` and :math:`K` can be seen as the interaction between two faces of the mesh. .. note:: :math:`K` should not be confused with the similar matrix :math:`D` defined as: .. math:: 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 :math:`G` is done with respect to a different variable. The matrix :math:`D` is used in the `direct` boundary integral equation, as e.g. in HAMS [Liu19]_. In the mathematical literature, :math:`D` is also referred to as the `double layer operator` and :math:`K` as the `adjoint double layer operator`. The matrices :math:`S` and :math:`K` relates the vectors :math:`\Phi`, :math:`u` and :math:`\sigma` through the following approximations of :eq:`potential_representation` and :eq:`normal_velocity_on_hull_representation`: .. math:: \Phi = S \sigma, \qquad u = K \sigma. :label: discrete_BEM_problem The resolution of the discrete problem with Nemoh consists of two main steps: 1. The evaluation of the coefficients of the complex-valued matrices :math:`S` and :math:`K` 2. The resolution of the complex-valued linear problem :math:`K \sigma = u`. Once :math:`\sigma` has been computed, :math:`\Phi` can be easily deduced. Then other magnitudes such as the Froude-Krylov forces or the added mass can be derived. .. mermaid:: :caption: A simplified flowchart of the internals of Capytaine solver 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 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 :math:`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 :math:`\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* .. math:: \omega_e = \omega - k U \cos (\beta) where :math:`k` is the wavenumber and :math:`\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 :math:`\delta\!r(x, y, z)` reads .. math:: \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 | :math:`\delta \! r` | :math:`\frac{\partial \delta\! r}{\partial x}` | +=======+=====================+================================================+ | Surge | :math:`(1, 0, 0)` | :math:`(0, 0, 0)` | +-------+---------------------+------------------------------------------------+ | Sway | :math:`(0, 1, 0)` | :math:`(0, 0, 0)` | +-------+---------------------+------------------------------------------------+ | Heave | :math:`(0, 0, 1)` | :math:`(0, 0, 0)` | +-------+---------------------+------------------------------------------------+ | Roll | :math:`(0, -z, y)` | :math:`(0, 0, 0)` | +-------+---------------------+------------------------------------------------+ | Pitch | :math:`(z, 0, -x)` | :math:`(0, 0, -1)` | +-------+---------------------+------------------------------------------------+ | Yaw | :math:`(-y, x, 0)` | :math:`(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 .. math:: 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 .. math:: \eta = -\frac{1}{g} \left( -i \omega_e \phi + U \frac{\partial \phi}{\partial x} \right) The computation of :math:`\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. .. mermaid:: :caption: A simplified flowchart of the internals of Capytaine solver **with forward speed**, where red boxes are the supplementary steps introduced by forward speed. 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 Post-processing =============== Forces on body surfaces ----------------------- Forces acting on body surfaces are computed by integration of the pressure field. .. math:: 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 :math:`p = i \omega \rho \Phi` stands for the complex-valued pressure fields in frequency-domain, :math:`n` is the normal vector on the hull :math:`\Gamma` (oriented towards the fluid in Capytaine, see :doc:`../user_manual/conventions`) and :math:`\delta\!r_j` is the local displacement of the hull of the degree of freedom :math:`j`. For a single rigid body, the degrees of freedom reads: +---------------+-----------------------------------------------------------------+ | Dof :math:`j` | Local hull displacement | +===============+=================================================================+ | Surge | :math:`\delta\!r_j(x) = (1, 0, 0)` | +---------------+-----------------------------------------------------------------+ | Sway | :math:`\delta\!r_j(x) = (0, 1, 0)` | +---------------+-----------------------------------------------------------------+ | Heave | :math:`\delta\!r_j(x) = (0, 0, 1)` | +---------------+-----------------------------------------------------------------+ | Roll | :math:`\delta\!r_j(x) = (1, 0, 0) \times (x-x_0, y-y_0, z-z_0)` | +---------------+-----------------------------------------------------------------+ | Pitch | :math:`\delta\!r_j(x) = (0, 1, 0) \times (x-x_0, y-y_0, z-z_0)` | +---------------+-----------------------------------------------------------------+ | Yaw | :math:`\delta\!r_j(x) = (0, 0, 1) \times (x-x_0, y-y_0, z-z_0)` | +---------------+-----------------------------------------------------------------+ where :math:`(x_0, y_0, z_0)` is the rotation center and :math:`\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 :math:`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 :func:`~capytaine.bem.airy_waves.airy_wave_pressure` function. 2. The diffraction forces :math:`F_{D}`, from the integration of the diffracted wave field (all bodies held fixed). 3. The radiation forces :math:`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 :math:`j` of the radiation force :math:`F_{R}` is further rewritten as .. math:: F_{R, j} = \sum_k \left[\omega^2 A_{jk} + i \omega B_{jk}\right] X_k where :math:`A_{jk}` is the added mass matrix, :math:`B_{jk}` is the radiation damping matrix and :math:`X_k` is the amplitude of the motion of the body along the degree of freedom :math:`k`. In other words, one has .. math:: 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 and .. math:: 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 where :math:`\Phi_k` is the potential field computed with the normal velocity on the hull :math:`\frac{\partial \Phi_k}{\partial n} = -i \omega \delta \! r_k \cdot n`. In Capytaine's wording, the degree of freedom :math:`k` defining the normal velocity on the hull is called ``radiating_dof``, while the degree of freedom :math:`j` used in the integration of the force is the ``influenced_dof``. .. note:: From Green second identity .. math:: \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, .. math:: \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 :math:`\frac{\partial \Phi_k}{\partial n} = -i \omega \delta \! r_k \cdot n`, some software such as the version 1 of Capytaine use :math:`\frac{\partial \tilde \Phi_k}{\partial n} = \delta \! r_k \cdot n`, that is :math:`\tilde \Phi_k = \frac{\Phi_k}{-i \omega}`. It leads to the following definition of the added mass and radiation damping .. math:: 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 and .. math:: 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 This form is convenient since the all the :math:`\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 :math:`\tilde \Phi` in version 1 of Capytaine was not consistent with the use of :math:`\Phi` for diffraction problem and it was easy to forget the missing :math:`-i\omega` for some post-processing of :math:`\tilde \Phi` for radiation problems. In version 2.0 of Capytaine, :math:`\Phi` is used everywhere instead of :math:`\tilde \Phi`. Since version 2.1, another method has been implemented to take into account the cancelling of the :math:`\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 .. math:: 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 :math:`x(t) = \mathrm{Re} \left( X e^{-i \omega t} \right)`, .. math:: \left[-\omega^2 M_{jk} - i \omega C_{jk} + K_{jk}\right] X_k = F_j, where :math:`M_{jk}` is the inertia matrix, accounting for the mass distribution, :math:`C_{jk}` is the mechanical damping matrix, :math:`K_{jk}` is the stiffness matrix which comprises mechanical and hydrostatic effects, and :math:`F_j` are generic external forces. .. note:: The hydrostatic contribution to matrix :math:`K_{jk}` accounts for a variation of hydrostatic force in direction :math:`j` due to a unit motion in direction :math:`k`. It is a geometric property of the body. As seen above, forces :math:`F_j` can be decomposed as .. math:: F_j = F_{FK, j} + F_{D, j} + F_{R, j} The full system becomes .. math:: \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 .. math:: H X = F_{ex} where :math:`H` denotes the following transfer function matrix .. math:: H_{jk} = \left[-\omega^2 (M_{jk} + A_{jk}) - i \omega (C_{jk} + B_{jk}) + K_{jk}\right] and :math:`F_{ex}` denotes the excitation force. .. math:: F_{ex, j} = F_{FK, j} + F_{D, j}. The oscillation amplitude is obtained by solving the complex-valued linear system. .. note:: Matrices :math:`A_{jk}` and :math:`B_{jk}` depend on :math:`\omega`, and so does :math:`H_{jk}` and :math:`X_k`. Free surface elevation ---------------------- The potential at the reference surface :math:`z = 0` can be connected to the free surface elevation by the dynamic condition .. math:: \dfrac{\partial \phi}{\partial t} = - g \eta which, in frequency domain, is .. math:: \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 .. math:: \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 .. math:: \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 :math:`\rho(x)` and :math:`\theta(x)` are the cylindrical coordinates of :math:`x` and :math:`H` is the so-called Kochin function which can be computed as .. math:: 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 .. prf:proof:: Let us first note, following [Del87]_ (p.79), that for large :math:`\rho(x)`: .. math:: :nowrap: \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 :doc:`inf_depth_green_function`) reads .. math:: - 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 .. math:: - 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 :eq:`potential_representation`, we have .. math:: \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 .. math:: 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: .. math:: \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 where :math:`\beta` is the wave direction, :math:`H` the Kochin function, :math:`k` the wavenumber, :math:`h` the water depth and :math:`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: .. math:: \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: .. math:: H(\theta) = e^{i\frac{\pi}{2}} ( H_D(\theta) + \sum_{k=1}^6 X_k H_{R_k} (\theta)) where :math:`X_k` is the motion RAO of the body corresponding to the degree of freedom :math:`k`, :math:`H_{R_k}` is the Kochin function associated with the radiated potential of degree of freedom :math:`k` and :math:`H_{D}` is the Kochin function associated with the diffracted potential.