# Second Order Difference- and Sum-Frequency Wave Loads in the Open-Source Potential Flow Solver NEMOH

Ruddy Kurnia, Guillaume Ducrozet, Jean-Christophe Gilloteaux

► To cite this version:

Ruddy Kurnia, Guillaume Ducrozet, Jean-Christophe Gilloteaux. Second Order Difference- and Sum-Frequency Wave Loads in the Open-Source Potential Flow Solver NEMOH. 41st International Conference on Ocean, Offshore and Arctic Engineering, ASME, Jun 2022, Hamburg, Germany. 10.1115/OMAE2022-79163. [hal-03846967](https://hal.science/hal-03846967)

HAL Id: hal-03846967

https://hal.science/hal-03846967

Submitted on 20 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SECOND ORDER DIFFERENCE- AND SUM-FREQUENCY WAVE LOADS
IN THE OPEN-SOURCE POTENTIAL FLOW SOLVER NEMOH

R. Kurnia

Ecole Centrale Nantes
LHEEA Lab. (ECN and CNRS)
Nantes 44000, FRANCE

Email: ruddy.kurnia@ec-nantes.fr

G. Ducrozet

Ecole Centrale Nantes
LHEEA Lab. (ECN and CNRS)
Nantes 44000, FRANCE

J-C. Gilloteaux

Ecole Centrale Nantes
LHEEA Lab. (ECN and CNRS)
Nantes 44000, FRANCE

ABSTRACT

Theoretical and numerical aspects of the open-source po-
tential flow boundary element solver, NEMOH, for the first or-
der hydrodynamic coefficients computations in the frequency do-
main are described in [Babarit, A. and Delhommeau, G., 2015].
[Philippe, M. et al., 2015] described the implementation and
verifications of the second order difference-frequency quadratic
transfer functions (QTFs) in the NEMOH code. In the latter pa-
per, the QTFs are verified for standard cases, a bottom-mounted
cylinder and a hemisphere. The present study reports the imple-
mentation and verification of the complete QTFs, for difference-
and sum-frequency loads. The QTFs are composed of quadratic
and potential parts. The quadratic part depending on the first
order hydrodynamic quantities, is implemented using the near-
field approach. The potential part, which depends on the second
order potential, is solved using the indirect method. Verifica-
tion is achieved by comparing the NEMOH result with a com-
mercial software HYDROSTAR for a hemisphere and the OC5-
DeepCwind semisubmersible.

INTRODUCTION

The open source potential flow boundary element solver, NEMOH<fnref n="1" />, is developed at Ecole centrale de Nantes [1, 2]. The code has been extensively validated and used in several ocean engineering applications, especially for the design and optimization of wave energy converters [3, 4].

Complementary, offshore wind energy is nowadays seen as
an alternative green and renewable energy with a huge potential
that will be a key element in the energetic transition to come.

Potential flow codes are now often used for studying offshore
floating wind turbine platforms. In [5], a commercial potential
flow code WAMIT [6] was used to study the second order differ-
ence wave frequency loads on the individual components of the
OC-5 DeepCwind semisubmersible under bichromatic waves.

The difference low-frequency wave, as a product of the pair
frequencies of nonlinear waves, excites large wave loads on a
structure if the difference wave frequency is close to the natural
frequency of the structure. On the other hand, the sum-frequency
wave can also excite large heave force and pitch moment on cer-
tain types of structures, such as the Tension Leg Platform (TLP)
as in [7]. Therefore, investigating second order wave excitation
forces at the difference and sum frequencies of paramount im-
portance.

Methods for solving second order diffraction wave problem
and to calculate second order wave loads have been developed
in past years by many researchers [8, 9, 10]. The second order
wave loads are composed of the quadratic part and the poten-
tial part. Several formulations are available for calculating the
quadratic part: the near-field formulation [9], the far-field for-
mulation [11], the middle field formulation [12] and the lagally
formulation [13]. The potential part can be solved directly by
solving the second order diffraction problem as in [10] or indi-
rectly using the green formulation as in [8, 9].

In the following sections, we describe the second order
wave loads formulations used and implemented in NEMOH. The
quadratic part is based on the near-field formulation and the po-
tential part is based on the indirect method. Then, NEMOH
full QTF implementation is validated through the comparison to
HYDROSTAR commercial software [14]. Four geometries are
tested: a simple hemisphere and the three individual components
of OC5-semisubmersible.

<footnote n="1">¹https://lheea.ec-nantes.fr/valorisation/
logiciels-et-brevets/nemoh-presentation</footnote>
## NOTATIONS

We start by defining notations used throughout this paper. The fixed earth Cartesian coordinate system is denoted by **x** = (x̄, z) with x̄ = (x, y) the horizontal coordinates perpendicular to the z axis in the opposite direction of gravity, g. The free surface elevation at time, t, is expressed as η(x̄, t) = 0 with respect to the mean water level at z = 0 and the fluid potential is denoted as Φ(**x**, t). The interior fluid domain is denoted by V<sub>Ω</sub>, the points located at the free surface fluid domain is denoted by S<sub>F</sub> and the ones located at the body boundary is S<sub>B</sub>. Water mass density is denoted by ρ.

Floating body has 6 degrees of freedom (DOF) determined at center of gravity (COG), ξ = (X, θ) where the positions, X = (X, Y, Z) and the orientations, θ = (θ₁, θ₂, θ₃). Displacement of points at the hull are specified by a body of vector r with respect to the COG as X = X + R(r). R is a rotation operator where R(r) ≈ θ × r. The velocity of the points at the hull is expressed as X.

We express the wetted part on the body hull S<sub>B</sub> as a function z = ζ(x,y,t). We use the normalized normal vector directed toward the fluid domain, **n** = −N/|N| with N = (−∇₂ζ, 1) where ∇₂ is the two dimensional gradient in x̄. Then the six-dimensional generalized normal vector is defined as **v** = (n, r × n)<sup>T</sup>, with ( )<sup>T</sup> is the matrix transpose operator. The normal vectors on the body at the mean position, S<sub>B₀</sub> are expressed as **v**₀ and **n**₀.

## HYDRODYNAMIC PROBLEM

In this section we describe the potential flow problem.

Fluid is assumed to be inviscid and the flow is irrotational (∇ × **U** = 0) allowing to write the fluid velocity **U** as **U** = ∇Φ. Fluid is also considered incompressible so that the fluid potential in the interior domain is described by Laplace equation.

The fluid potential has to satisfy a set of boundary conditions. First, the constant pressure condition on the free surface. Second, the diffraction and radiation conditions on the body hull. Third, the impermeable condition at bottom and the last, radiation wave condition at the far field.

Taylor series development is applied to approximate the fluid potential on the free surface by the potential at the mean water level. The hydrodynamic problem can then be reformulated at the first order, n = 1, and at the second order nonlinear, n = 2, by applying the perturbation series, as

$$ egin{cases} 
\nabla^2 \Phi^{(n)}(\mathbf{x}, t) = 0 & \mathbf{x} \in V_{\Omega} \\
\partial_t^2 \Phi^{(n)} + g \partial_z \Phi^{(n)} = Q_F^{(n)} & \mathbf{x} \in S_F \\
\partial_n \Phi^{(n)} = \mathbf{\dot{X}}_B^{(n)} \cdot \mathbf{n}_0 + Q_B^{(n)} & \mathbf{x} \in S_B \\
\partial_z \Phi^{(n)}(\mathbf{x}) = 0 & \text{at } z = -D \\
\lim(\Phi^{(n)} - \Phi_I^{(n)}) = 0 & \text{for } x^2 + y^2 \to \infty,
\end{cases} 
\quad (1) $$

where the incoming wave potential, Φ<sub>I</sub><sup>(n)</sup>, will be discussed in the following sections. The free surface forcing term, Q<sub>F</sub><sup>(n)</sup>, and the body forcing term, Q<sub>B</sub><sup>(n)</sup>, are zero for the first order terms and for the second order terms are given by

$$ Q_F^{(2)} = -2\nabla\Phi^{(1)} \cdot \partial_t\nabla\Phi^{(1)} + \frac{1}{g}\partial_t\Phi^{(1)}\partial_z\left[\partial_{tt}\Phi^{(1)} + g\partial_z\Phi^{(1)}\right] \Big|_{z=0} \quad (2) $$

$$ Q_B^{(2)} = [\mathbf{\dot{X}}^{(1)} - \nabla\Phi^{(1)}] \cdot R^{(1)}(\mathbf{n}_0) - \left[ (\mathbf{\dot{X}}^{(1)} \cdot \nabla) \nabla\Phi^{(1)} \right] \cdot \mathbf{n}_0 \Big|_{S_{B_0}} \quad (3) $$

In NEMOH, the first order hydrodynamic problem is formulated in the source distribution boundary integral equations and solved using a panel method as described in [1].

The second order hydrodynamic problem is not solved directly in NEMOH but the second order force for sum- and difference-frequency can be calculated using the indirect method as in [8, 9] and described in the following sections.

## SECOND ORDER HYDRODYNAMIC FORCES

The hydrodynamic force is obtained from integration of the hydrodynamic pressure over the body hull. Following [9], the integration over the body hull, S<sub>B</sub>, is composed of the integration over the mean wetted surface position, S<sub>B₀</sub>, and of the integration over the perturbed body surface, εS. Taylor expansion is applied to the potential over the body hull to be expressed as the potential at the mean body hull position. The hydrodynamic pressure, P<sub>H</sub>, is obtained by subtracting the hydrostatic part in the Bernoulli pressure. Then, by applying a perturbation series up to and including second order term, ε², the truncated hydrodynamic pressure is defined as

$$ 
\begin{aligned}
P_H^{(1)} &= -\rho \partial_t \Phi^{(1)} \\
P_H^{(2)} &= -\rho \left[ \partial_t \Phi^{(2)} + \frac{1}{2} |\nabla \Phi^{(1)}|^2 + \mathbf{\dot{X}}^{(1)} \cdot \partial_t \nabla \Phi^{(1)} \right].
\end{aligned} 
$$

The hydrodynamic force is then given as

$$ 
\begin{aligned}
\mathbf{F}_H(t) &= -\iint_{S_{B_0} \cup \varepsilon S} P_H \mathbf{v} \, dS = \varepsilon \mathbf{F}_H^{(1)} + \varepsilon^2 \mathbf{F}_H^{(2)}, \\
\mathbf{F}_H^{(1)} &= -\iint_{S_{B_0}} P_H^{(1)} \mathbf{v}_0 dS \\
\mathbf{F}_H^{(2)} &= R^{(1)} \mathbf{F}_H^{(1)} - \iint_{S_{B_0}} P_H^{(2)} \mathbf{v}_0 dS - \iint_{\varepsilon S} P_H^{(1)} \mathbf{v}_0 dS.
\end{aligned} 
\quad (4) $$
The second order hydrodynamic force in Eq. (4) can be composed as a sum of the quadratic terms of the first order quantities, $\mathbf{F}_{H_1}^{(2)}$, and a second order potential term, $\mathbf{F}_{H_2}^{(2)}$.

Similarly, the second order excitation force, $\mathbf{F}_{exc}^{(2)}$, is also separated in the quadratic term and the potential term. In the quadratic term of the excitation force, there is also a contribution from the second order hydrostatic force as

$$ \mathbf{F}_{exc1}^{(2)} = R^{(1)}(\mathbf{F}_{S}^{(1)}) + \rho g \iint_{\mathcal{E}S} z^{(1)} \mathbf{v}_{0}dS + \mathbf{F}_{H_1}^{(2)}. $$

Note that the sum of the first order hydrostatic force and the first order hydrodynamic force is the inertia force, $\mathbf{F}_I^{(1)}$. The perturbed surface, $\mathcal{E}S$, is located between an instantaneous waterline, $\zeta_{wl}$ and the free surface along the waterline. The integral then can be performed horizontally along the waterline, $\Gamma$ and vertically between $\zeta_{wl}$ and $\eta$. Let $\Theta$ be the angle of hull with respect to vertical axis along the waterline, so $dS = d\Gamma dz/\cos\Theta$. With the first order free surface is defined as $\eta^{(1)} = -\partial_t \Phi^{(1)}/g$, then the integral over the perturbed surface in the excitation force can be calculated. The excitation force is then expressed, with $\Gamma_0$ is the waterline of the body at the mean position, as

$$ \mathbf{F}_{exc1}^{(2)} = R^{(1)}(\mathbf{F}_{I}^{(1)}) - \frac{\rho g}{2} \int_{\Gamma_0} [\eta^{(1)} - \zeta_{wl}^{(1)}]^2 \mathbf{v}_0 d\Gamma \quad (5) \\ + \iint_{S_{B_0}} \rho \left[ \frac{1}{2} |\nabla \Phi^{(1)}|^2 + \mathcal{X}^{(1)} \cdot \partial_t \nabla \Phi^{(1)} \right] \mathbf{v}_0 dS $$

This quadratic term of the second order excitation force, Eq. (5), has been verified and the results shown in [2].

The potential term of the second hydrodynamic force is

$$ \mathbf{F}_{H_2}^{(2)} = \iint_{S_{B_0}} \rho \partial_t \Phi^{(2)} \mathbf{v}_0 dS. $$

As in the first order potential, the second order total potential, $\Phi^{(2)}$, can be composed as a sum of the second order incoming potential, $\Phi_I^{(2)}$, and the perturbed potential, $\Phi_P^{(2)}$. The perturbed potential is $\Phi_P^{(2)} = \Phi_D^{(2)} + \Phi_R^{(2)}$, where $\Phi_D^{(2)}$ is the diffraction potential and $\Phi_R^{(2)}$ is the radiation potential. The second order radiation potential is due to the second order motion, $\Phi_R^{(2)} = \dot{\xi}^{(2)} \cdot \mathbf{\psi}$, where $\mathbf{\psi}$ is the normalized radiation potential. Then the force is composed as a sum of the potential part of the second order excitation force and the second order radiation force. The potential part of the second order excitation force, $\mathbf{F}_{exc2}^{(2)}$, is a sum of the second order Froude-Krylov force, $\mathbf{F}_{H_I}^{(2)}$, and the second order diffraction force, $\mathbf{F}_{H_D}^{(2)}$, $\mathbf{F}_{exc2}^{(2)} = \mathbf{F}_{H_I}^{(2)} + \mathbf{F}_{H_D}^{(2)}$.

The second order potential is assumed to be an harmonic function, which frequency is the sum- or difference-frequency modes, with $\omega^+ = \omega_1 + \omega_2$ and $\omega^- = \omega_1 - \omega_2$. The potential can then be expressed as $\Phi^{(2)\pm}(\mathbf{x}, t) = \text{Re}\{\Phi_1^{(2)}(\mathbf{x})e^{-i\omega_1^\pm t}\}$. Correspondingly, the excitation force is expressed as

$$ \begin{aligned} \mathbf{F}_{exc2}^{(2)} &= \text{Re} \left\{ \left[ \mathbf{F}_{H_I}^{(2)} + \mathbf{F}_{H_D}^{(2)} \right] e^{-i\omega_1^\pm t} \right\}, \quad \text{where} \tag{6} \\ \mathbf{F}_{H_I}^{(2)} &= -i\omega_1^\pm \rho \iint_{S_{B_0}} \Phi_{I_1}^{(2)} \mathbf{v}_0 dS \\ \mathbf{F}_{H_D}^{(2)} &= -i\omega_1^\pm \rho \iint_{S_{B_0}} \Phi_{D_1}^{(2)} \mathbf{v}_0 dS. \end{aligned} $$

The Froude-Krylov force, $\mathbf{F}_{H_I}^{(2)}$, in Eq. (6), can be calculated by given the second order incoming potential as discussed in the next section.

Instead of solving the second order diffraction problem, the diffraction force is obtained with the indirect method using Green formulation and an assisting function that is a normalized radiation potential [8]. As also described in [9], the normalized radiation potential, $\mathbf{\psi}$, is the first order solution of the radiation potential that satisfies the boundary value problem (BVP) similarly as in the first order problem of Eq. (1) but restricted to the radiation problem. The BVP, Eq. (1), has also to be satisfied by the second order diffraction potential. Then the BVP is reformulated in the boundary integral equations using the Green-Gauss theorem. The integrals over the bottom surface, $S_D$, and over the far field boundary vanish, the boundary integral of the second order diffraction potential can then be expressed as

$$ \iint_{S_{B_0}} \Phi_{D1}^{(2)} \mathbf{v}_0 dS = - \iint_{S_{B_0}} (\partial_{\mathbf{n}} \Phi_I^{(2)} - Q_B^{(2)\pm}) \mathbf{\psi}^{\pm} dS \\ - \frac{1}{g} \iint_{S_F} Q_{FD}^{(2)\pm} \mathbf{\psi}^{\pm} dS \quad (7) $$

where the body forcing term $Q_B^{(2)\pm}$ as in Eq. (3) but now for sum and difference frequencies. The diffraction free surface forcing term is obtained as, $Q_{FD}^{(2)\pm} = Q_F^{(2)\pm} - Q_{FI}^{(2)\pm}$ where $Q_F^{(2)\pm}$ and $Q_{FI}^{(2)\pm}$ are obtained as in Eq. (2) but with the sum- and difference-frequency total potential and the incoming potential, respectively. The potential part of the second order diffraction force, $\mathbf{F}_{H_D}^{(2)}$, Eq. (6), can now be obtained by the boundary integral in Eq. (7).

## DIFFERENCE AND SUM FREQUENCY FORCES

In this section, we give an explicit expression for the second order excitation force on a structure excited by uni-directional bichromatic incoming wave.
Following identity of product of two bi-harmonic functions will be used. Suppose a harmonic function with coefficient, $a, \Psi_{lj} = a_{lj}e^{-i\omega_j t}$ and a biharmonic function $\varphi_l = Re\{\Psi_{l1} + \Psi_{l2}\} = Re\{\Psi_l\} = \frac{\Psi_{l1} + \Psi_l^*}{2}$, where * is a complex conjugate, then

$$ \varphi_1 \cdot \varphi_2 = \frac{1}{2} \begin{bmatrix} \mathrm{Re}\{\Psi_{11}\Psi_{21}\} + \mathrm{Re}\{\Psi_{12}\Psi_{22}\} \\ + \mathrm{Re}\{\Psi_{11}\Psi_{21}^*\} + \mathrm{Re}\{\Psi_{12}\Psi_{22}^*\} \\ + \mathrm{Re}\{\Psi_{11}\Psi_{22}^* + \Psi_{12}^*\Psi_{21}\} \\ + \mathrm{Re}\{\Psi_{11}\Psi_{22} + \Psi_{12}\Psi_{21}\} \end{bmatrix} . \quad (8) $$

The product of the two biharmonic functions, Eq. (8), is composed of several terms; the double frequency terms $e^{-i2\omega_1 t}$, $e^{-i2\omega_2 t}$, the constant terms, the difference-frequency, $e^{-i(\omega_1 - \omega_2)t}$, and the sum-frequency, $e^{-i(\omega_1 + \omega_2)t}$. Here we focus only on the difference- and sum-frequency products in Eq. (8) and rewrite the products, with a complex conjugate operator applied to a complex variable, i.e. $\gamma$ as $\gamma^c = \gamma^*$ and $\gamma^{c+} = \gamma$, as

$$ [\varphi_1 \cdot \varphi_2]^\pm = \frac{1}{2} \mathrm{Re} \left\{ \Psi_{11} \Psi_{22}^{c\pm} + \Psi_{12}^{c\pm} \Psi_{21} \right\}. \quad (9) $$

This product is used to evaluate the second order incident potential and the second order excitation force in the next subsections.

## Second order Froude-Krylov force

Incident Airy potential for bichromatic wave propagation above constant bottom with two radial frequencies $\omega_j$, with $j = 1$ and 2, wave number vectors $\vec{k}_j = (k_j \cos \beta, k_j \sin \beta)$, $\beta$ is angle from the positive x-axis and the wave number $k$ is related with $\omega$ in the dispersion relation $\omega = \Omega(k, D) = \sqrt{gk \tanh kD}$, is defined as

$$ \begin{aligned} \Phi_I(\mathbf{x}, t) &= \mathrm{Re} \left\{ \Phi_{I_1}(\mathbf{x}) e^{-i\omega_1 t} + \Phi_{I_2}(\mathbf{x}) e^{-i\omega_2 t} \right\} \\ \Phi_{I_j}(\mathbf{x}) &= -i \frac{a_j g}{\omega} \mathcal{Z}(k_j, D, z) e^{i\vec{k}_j \cdot \vec{x}} \end{aligned} \quad (10) $$

where the Airy profile $\mathcal{Z}(k, D, z) = \frac{\cosh k(D+z)}{\cosh kD}$.

Using the product identity in Eq. (9) to evaluate the free surface forcing term in Eq. (1) with the biharmonic incident potential, Eq. (10), the difference- and sum-frequency incident potential is expressed, with $\omega^\pm = \omega_1 \pm \omega_2$, $k^\pm = k_1 \pm k_2$, as

$$ \begin{aligned} \Phi_I^{(2)\pm}(\mathbf{x}, t) &= \mathrm{Re} \left\{ \Phi_{I_1}^{(2)\pm}(\mathbf{x}) e^{-i(\omega_1 \pm \omega_2)t} \right\} \\ \Phi_{I_1}^{(2)\pm}(\mathbf{x}) &= \frac{ia_1 a_2 g^2 e^{i\vec{k}_1 \cdot \vec{x}} \mathcal{Z}(k^\pm, D, z)}{-\omega^{(\pm)} + \Omega^2(k^\pm, D)} \\ & \quad \begin{bmatrix} \frac{\omega^\pm}{\omega_1} k_1 k_2 (1 \mp \tanh k_1 D \tanh k_2 D) \\ + \frac{1}{2} \left[ \frac{k_1^2}{\omega_1^2} \sech^2 k_1 D \pm \frac{k_2^2}{\omega_2^2} \sech^2 k_2 D \right] \end{bmatrix}. \end{aligned} \quad (11) $$

The Froude-Krylov force, $F_{Hl}^\pm$, Eq. (6), with the incoming potential is given in Eq. (11), can then be calculated directly. This force contributes in the potential part of the second order excitation force.

## Quadratic part

The quadratic part of the difference- and sum-frequency excitation forces, with the 6-DOF quadratic transfer function (QTF), $T_{FQ}$, is given as follow

$$ \begin{aligned} \mathbf{F}_{exc1}^{(2)} &= \mathrm{Re} \left\{ \mathbf{T}_{FQ}^{-} a_1 a_2^* e^{-i(\omega_1 - \omega_2)t} \right\} + \mathrm{Re} \left\{ \mathbf{T}_{FQ}^{+} a_1 a_2^* e^{-i(\omega_1 + \omega_2)t} \right\}, \\ \mathbf{T}_{FQ}^{\pm} &= [\mathbf{F}_{11}^{\pm} + \mathbf{F}_{12}^{\pm} + \mathbf{F}_{13}^{\pm} + \mathbf{F}_{14}^{\pm}] / a_1 a_2^*. \end{aligned} \quad (12) $$

By applying the product identity in Eq. (9) to evaluate all the products in Eq. (5), the $F_{11}^\pm, F_{12}^\pm, F_{13}^\pm$ and $F_{14}^\pm$ are derived as follows

$$ F_{11}^{\pm} = -\frac{\rho g}{2} \int_{\Gamma_0} \left[ \eta_1^{(1)} - \zeta_{wl_1}^{(1)} \right] \left[ \eta_2^{(1)} - \zeta_{wl_2}^{(1)} \right]^{c\pm} \mathbf{v}_0 d\Gamma \quad (13) $$

$$ F_{12}^{\pm} = \frac{\rho}{2} \iint_{S_{B_0}} \left[ \nabla \Phi_1^{(1)} \cdot \nabla \Phi_2^{(1)C_{\pm}} \right] \mathbf{v}_0 dS $$

$$ F_{13}^{\pm} = \frac{\rho}{2} \iint_{S_{B_0}} \left[ \mathcal{X}_1^{(1)} \cdot \left( -i\omega_2 \nabla \Phi_2^{(1)} \right)^{c\pm} - \mathcal{X}_2^{(1)C_{\pm}} \cdot i\omega_1 \nabla \Phi_1^{(1)} \right] \mathbf{v}_0 dS $$

$$ F_{14}^{\pm} = \frac{1}{2} \left[ R_1^{(1)} (\mathbf{F}_{I_2}^{(1)C_{\pm}}) + R_2^{(1)C_{\pm}} (\mathbf{F}_{I_1}^{(1)}) \right]. $$

## Potential part

The potential part of the difference- and sum-frequency excitation forces, with the 6-DOF quadratic transfer function (QTF), $T_{FP}$, is given as follows

$$ \begin{aligned} \mathbf{F}_{exc2}^{(2)} &= \mathrm{Re} \left\{ \mathbf{T}_{FP}^{-} a_1 a_2^* e^{-i(\omega_1 - \omega_2)t} \right\} + \mathrm{Re} \left\{ \mathbf{T}_{FP}^{+} a_1 a_2^* e^{-i(\omega_1 + \omega_2)t} \right\}, \\ \mathbf{T}_{FP}^{\pm} &= \left[ \mathbf{F}_{Hl}^{\pm} + \mathbf{F}_{HDB1}^{\pm} + \mathbf{F}_{HDF1}^{\pm} + \mathbf{F}_{HDF2}^{\pm} \right] / a_1 a_2^*, \end{aligned} \quad (14) $$

where $F_{Hl}^\pm$ is the Froude-Krylov force, the diffraction force is composed of several terms: the body forcing term $F_{HDB1}^\pm$, the free surface forcing term in the finite domain, $F_{HDF1}^\pm$, and the asymptotic free surface forcing term in the infinite domain, $F_{HDF2}^\pm$.

**Body forcing term** Using the product identity in Eq. (9) to evaluate the body forcing term in Eq. (3), the diffraction force is calculated and composed of two terms, the term with only first spatial derivatives in $F_{HDB1}^\pm$ and only the second spatial derivative
$F_{HDB_{12}}^{\pm}$. The first force term is

$$F_{HDB_1}^{(2)\pm} = i\omega^{\pm}\rho \iint_{S_{B_0}} (\partial_n \Phi_I^{(2)\pm} - Q_{B_1}^{(2)\pm}) \psi^{\pm} dS \quad (15)$$

$$Q_{B_1}^{(2)\pm} = \frac{1}{2} \left[ \begin{aligned} & (\dot{\mathbf{X}}_1^{(1)} - \nabla \Phi_1^{(1)}) \cdot R_2^{(1)\mathcal{C}_{\pm}}(\mathbf{n}_0) \\ & + (\dot{\mathbf{X}}_2^{(1)\mathcal{C}_{\pm}} - \nabla \Phi_2^{(1)\mathcal{C}_{\pm}}) \cdot R_1^{(1)}(\mathbf{n}_0) \end{aligned} \right]_{S_{B_0}}.$$

The normal derivative of the incident potential can be calculated for given the incident potential in Eq. (11).

The second force term in $F_{HDB_{12}}^{(2)\pm}$ is calculated by applying the Stokes theorem such that the force is only expressed in first spatial derivative as

$$\mathbf{F}_{HDB_2}^{(2)\pm} = \frac{i\omega^{\pm}\rho}{2} \iint_{S_{B_0}} \left[ \begin{aligned} & \left( \nabla\Phi_2^{(1)\mathcal{C}_{\pm}} \cdot \nabla \right) \left( \psi^{\pm}\mathcal{X}_1^{(1)} \right) \\ & + \left( \nabla\Phi_1^{(1)} \cdot \nabla \right) \left( \psi^{\pm}\mathcal{X}_2^{(1)\mathcal{C}_{\pm}} \right) \\ & - \left( \nabla \cdot \left( \psi^{\pm}\mathcal{X}_1^{(1)} \right) \right) \nabla\Phi_2^{(1)\mathcal{C}_{\pm}} \\ & - \left( \nabla \cdot \left( \psi^{\pm}\mathcal{X}_2^{(1)\mathcal{C}_{\pm}} \right) \right) \nabla\Phi_1^{(1)} \end{aligned} \right] \cdot \mathbf{n}_0 \, dS \\ - \frac{i\omega^{\pm}\rho}{2} \int_{\Gamma} \left[ \begin{aligned} & \left( \psi^{\pm}\mathcal{X}_1^{(1)} \right) \times \nabla\Phi_2^{(1)\mathcal{C}_{\pm}} \\ & + \left( \psi^{\pm}\mathcal{X}_2^{(1)\mathcal{C}_{\pm}} \right) \times \nabla\Phi_1^{(1)} \end{aligned} \right] \cdot d\Gamma. \quad (16)$$

**Free surface forcing term** The diffraction force due to the free-surface forcing in the finite domain $S_{F_1}$, $F_{HDF_1}^{(2)\pm}$, is calculated at a radius $R_{F_1}$ between the body waterline $R_{wl}$ and a finite distance $R_e$ from the origin, $R_{wl} < R_{F_1} \le R_e$. The calculation is grouped for terms with first order spatial derivative, $F_{HDF_{11}}^{(2)}$, and for terms with second order spatial derivative, $F_{HDF_{12}}^{(2)}$. Using the product identity in Eq. (9) to evaluate the first terms diffraction free surface forcing body forcing $Q_{FD_1}^{(2)}$, the force is given as

$$F_{HDF_{11}}^{(2)\pm} = \frac{i\omega^{\pm}\rho}{g} \iint_{S_{F_1}} Q_{FD_1}^{(2)\pm} \psi^{\pm} dS \quad (17)$$

$$Q_{FD_1}^{(2)\pm} = \left[ \begin{array}{l} i(\omega_1 \pm \omega_2) \left[ \nabla \Phi_1^{(1)} \cdot \nabla \Phi_{P_2}^{(1)\mathcal{C}_{\pm}} + \nabla \Phi_{P_1}^{(1)} \cdot \nabla \Phi_{I_2}^{(1)\mathcal{C}_{\pm}} \right] \\ - \frac{i\omega_1}{2g} \left[ \Phi_1^{(1)}(-\omega_2^2 \partial_z \Phi_{P_2}^{(1)\mathcal{C}_{\pm}}) + \Phi_{P_1}^{(1)}(-\omega_2^2 \partial_z \Phi_{I_2}^{(1)\mathcal{C}_{\pm}}) \right] \\ \mp \frac{i\omega_2}{2g} \left[ \Phi_2^{(1)\mathcal{C}_{\pm}}(-\omega_1^2 \partial_z \Phi_{P_1}^{(1)}) + \Phi_{P_2}^{(1)\mathcal{C}_{\pm}}(-\omega_1^2 \partial_z \Phi_{I_1}^{(1)}) \right] \end{array} \right]_{z=0}.$$

The second diffraction force terms due to the free-surface

forcing in the finite domain is given as

$$F_{HDF_{12}}^{(2)\pm} = \frac{i\omega^{\pm}\rho}{g} \iint_{S_{F_1}} Q_{FD_2}^{(2)\pm} \psi^{\pm} dS \quad (18)$$

$$Q_{FD_2}^{(2)\pm} = \frac{1}{2} \begin{bmatrix} - \left( i\omega_1 \Phi_1^{(1)} \right) k_2^2 \Phi_{P_2}^{(1)\mathcal{C}_{\pm}} \\ - \left( i\omega_1 \Phi_{P_1}^{(1)} \right) k_2^2 \Phi_{I_2}^{(1)\mathcal{C}_{\pm}} \\ \mp \left( i\omega_2 \Phi_2^{(1)\mathcal{C}_{\pm}} \right) k_1^2 \Phi_{P_1}^{(1)} \\ \mp \left( i\omega_2 \Phi_{P_2}^{(1)\mathcal{C}_{\pm}} \right) k_1^2 \Phi_{I_1}^{(1)} \end{bmatrix}.$$

Note that here the second derivative in $z$ is evaluated directly as $\partial_z \Phi_I = k^2 \Phi_I$ and similarly for the perturbed potential, that satisfy dispersion relation in the free surface, $\partial_{zz} \Phi_P = k^2 \Phi_P$. This approach is different from original work of [9] that used the Green formulation to express the second $z$ derivative such that it becomes only the first derivative. The Green formulation leads to an integral with a large error, particularly in the high frequency case, due to the first derivatives of the functions, i.e. the perturbed potential and the radiation potential, are not continuous on the body boundary.

The asymptotic diffraction force due to the free-surface forcing in the infinite domain $S_{F_2}$, $F_{HDF_2}^{(2)\pm}$, is calculated at a radius $R_{F_2}$ between the truncated radius $R_e$ to the infinity, $R_e \le R_{F_2} \le \infty$. Following [9], the perturbed potential on the free surface is expressed in an asymptotic form a cylindrical coordinate system $(r, \vartheta, z)$, with the Kochin function, $\mathcal{H}_P(\vartheta)$, defined as

$$\mathcal{H}_P(\vartheta) = \sum_{l=0}^{\infty} C_{IP} \cos l\vartheta + S_{IP} \sin l\vartheta \quad (19)$$

$$C_{IP} = \frac{-1}{4\pi} \iint_{S_{B_0}} \sigma_P(\mathbf{x}') \mathcal{Z}(k, D, z') \varepsilon_l(-i)^l J_l(kr') \cos l\alpha' dS$$

$$S_{IP} = \frac{-1}{4\pi} \iint_{S_{B_0}} \sigma_P(\mathbf{x}') \mathcal{Z}(k, D, z') \varepsilon_l(-i)^l J_l(kr') \sin l\alpha' dS$$

where $\varepsilon_0 = 1$, $\varepsilon_l = 2$ for $l \ge 1$, $r' = \sqrt{x'^2 + y'^2}$, $\alpha' = \tan^{-1}(y/x')$, $\vartheta = \tan^{-1}(y/x)$ and $J_l$ is the first kind Bessel function order $l$. $\sigma_P(\mathbf{x}')$ is the perturbed source distribution, $\mathbf{x}'$ is location of points on the body hull. For the radiation potential, $C_{IR}^{\pm}$ and $S_{IR}^{\pm}$ are defined similar as in Eq. (19) but now with the radiation source distribution function. The source distributions are obtained from the first order NEMOH hydrodynamic solution. Then, the asymptotic radiation, perturbed and incoming potentials are defined re-
spectively as

$$
\begin{align}
\Psi_{R1}^{\pm} &= \frac{\sqrt{8\pi k^{\pm}}}{\sqrt{r}} \mathcal{F}(k^{\pm}, D, z) e^{i(k^{\pm}r + \pi/4)} \sum_{l=0}^{\infty} \mathbf{C}_{lR}^{\pm} \cos l\vartheta + \mathbf{S}_{lR}^{\pm} \sin l\vartheta \nonumber \\
\Phi_{P_{1,2}} &= \frac{\sqrt{8\pi k_{1,2}}}{\sqrt{r}} \frac{\mathcal{X}(k_{1,2}, D, z) e^{i(k_{1,2}r + \pi/4)} \mathcal{H}_P(\vartheta)}{k_{1,2}D(1 - \tanh^2 k_{1,2}D) + \tanh k_{1,2}D} \nonumber \\
\Phi_{I_{1,2}} &= \frac{-ia_{1,2}g}{\omega_{1,2}} \sum_{l=0}^{\infty} \varepsilon_l i^l J_l(k_{1,2}r) \cos l(\vartheta - \beta). \tag{20}
\end{align}
$$

The force, as expressed in Eq. (17) and Eq. (18) but now for the domain $S_{F_2}$ with the potentials expressed as in Eq. 20, is given as

$$
\mathbf{F}_{H_{DF2}}^{(2)\pm} = \frac{i\omega^{\pm}\rho}{g} \left[ \mathcal{H}_{1}^{\pm}(k_1, k_2, \omega_1, \omega_2) \left[ \mathbf{I}_{DF11}^{\pm} + \mathbf{I}_{DF12}^{\pm} \right] + \mathcal{H}_{2}^{\pm}(k_1, k_2, \omega_1, \omega_2) \left[ \mathbf{I}_{DF21}^{\pm} + \mathbf{I}_{DF22}^{\pm} \right] \right] \quad (21)
$$

where

$$
\begin{align*}
\mathcal{H}_1^{\pm}(k_1, k_2, \omega_1, \omega_2) &= \mp i\omega^{\pm}k_1 k_2 \\
\mathcal{H}_2^{\pm}(k_1, k_2, \omega_1, \omega_2) &= i\omega^{\pm} \frac{\Omega^2(k_1, D) \Omega^2(k_2, D)}{g} \\
&\quad \mp \frac{i\omega_1 \omega_2}{2} \left( \frac{k_1^2}{\omega_1 \cosh^2(k_1 D)} \pm \frac{k_2^2}{\omega_2 \cosh^2(k_2 D)} \right) \\
\mathbf{I}_{DF11}^{\pm} &= \int_0^{2\pi} \int_{R_e}^{\infty} \cos(\vartheta - \beta) \Phi_{I_1}^{(1)} \Phi_{I_2}^{(1)\mathcal{C}_{\pm}} \Psi^{\pm} r dr d\vartheta \\
\mathbf{I}_{DF12}^{\pm} &= \int_0^{2\pi} \int_{R_e}^{\infty} \cos(\vartheta - \beta) \Phi_{P_1}^{(1)} \Phi_{I_2}^{(1)\mathcal{C}_{\pm}} \Psi^{\pm} r dr d\vartheta \\
\mathbf{I}_{DF21}^{\pm} &= \int_0^{2\pi} \int_{R_e}^{\infty} \Phi_{I_1}^{(1)} \Phi_{P_2}^{(1)\mathcal{C}_{\pm}} \Psi^{\pm} r dr d\vartheta \\
\mathbf{I}_{DF22}^{\pm} &= \int_0^{2\pi} \int_{R_e}^{\infty} \Phi_{P_1}^{(1)} \Phi_{I_2}^{(1)\mathcal{C}_{\pm}} \Psi^{\pm} r dr d\vartheta.
\end{align*}
$$

The integral over $r$ is conducted semi analytically.

RESULTS

In this section, we perform an extensive validation of the NEMOH QTF through the computation of different test-cases. First, a floating hemisphere is studied as a standard and simple configuration. Then, the OC5-DeepCwind semi-submersible platform is modelled. The NEMOH results are compared with the commercial hydrodynamic software, HYDROSTAR [14]. For the comparison, the same meshes are used both in NEMOH and HYDROSTAR. Note that in both NEMOH and HYDROSTAR softwares, the near-field formulation is used for the QTF computation. In addition, both codes have a specific

<figure><img src="image_2.png" /><figcaption>FIGURE 1. HEMISPHERE MESH WITH 514 PANELS.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. COMPARISON OF ADDED MASS (TOP), DAMPING COEFFICIENT (MIDDLE), EXCITATION FORCE (BOTTOM) OF SURGE (LEFT) AND HEAVE (RIGHT) BETWEEN HYDROSTAR AND NEMOH WITH DIFFERENT NUMBER OF PANELS.</figcaption></figure>

treatment of irregular frequencies. In NEMOH the irregular fre-
quencies are removed by extended the boundary integral equa-
tions with additional meshes on the lid of the structures [1]. In
HYDROSTAR, similar approach is implemented, except that the
lid meshes are automatically generated. All the meshes used here
are generated using an open source mesh generator software,
GMSH [15]. The convergence study of the first order quantities
of the hydrodynamic coefficients are conducted and compared
for different number of meshes. The first order results are basis
for calculating the QTFs. All the QTFs results are normalized by
ρg, where ρ = 1025 kg/m³ and g = 9.81 m/s².

Hemisphere

This test-case considers a hemisphere with a radius 10 m in
a fluid domain with 30 m water depth. Fig. 1 shows the half
<figure><img src="image_4.png" /><figcaption>**FIGURE 3.** CONTOUR PLOTS OF THE NORMALIZED SURGE DIFFERENCE-FREQUENCY TOTAL QTF MAGNITUDE FOR THE FLOATING HEMISPHERE. THE LEFT COLUMN: HYDROSTAR, THE MIDDLE COLUMN: NEMOH, THE RIGHT COLUMN: THE DIFFERENCE.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 4.** CONTOUR PLOTS OF THE NORMALIZED SURGE SUM-FREQUENCY QUADRATIC QTF MAGNITUDE FOR THE FLOATING HEMISPHERE. THE LEFT COLUMN: THE DATA IN KIM AND YUE [10], THE MIDDLE COLUMN: NEMOH, THE RIGHT COLUMN: THE DIFFERENCE.</figcaption></figure>

space of the hemisphere that is discretized with 514 quadrilateral panels. The free surface is discretized by a circular mesh in the radius $R_{F1}$, $10 \le R_{F1} \le 600$ m, with $N_r = 600$ points in the radial and $N_\theta = 50$ points in the angular. Those meshes have been chosen after a convergence study on the results, as shown in Fig. 2 for the first order results. The QTFs results are shown with the wave frequencies in the interval $f = (0.14, 0.22)$ Hz and $\Delta f = 0.008$ Hz.

In a previous study [2], the NEMOH results for the quadratic terms of the difference-frequency QTFs for a fixed and floating hemisphere have been shown and verified against the numerical data proposed in [16]. Fig. 3 shows the contour plots of the surge difference-frequency total QTF, $|T_F^-| = |T_{FQ}^- + T_{FP}^-|$ Eqs. (12) to (14), that are obtained with HYDROSTAR and NEMOH. In addition, the difference between the two solvers results is shown in Fig. 3 to evaluate their similarities. It is observed that NEMOH and HYDROSTAR are in good agreement, with larger (but still limited) differences observed for the higher frequencies $f > 0.2$ Hz.

The surge sum-frequency quadratic QTF, $|T_{FQ}^+|$ Eq. (12), result obtained by NEMOH is compared with the corresponding data in [10]. The result show that an excellent agreement is achieved between the NEMOH result and the data in [10], see in Fig. 4. In NEMOH the irregular frequency removal method

is not optimal yet, as can be seen in Fig. 2, that leads to a large error in the potential part of the sum-frequency QTF, Eq. (14). We compared also the NEMOH result against HYDROSTAR, the agreement is achieved only in the frequency interval less than 0.15 Hz where the sum-frequency is not close to an irregular frequency at 0.3 Hz, but the result is not shown here.

This validate the full QTF calculations within NEMOH in a simple configuration and that the next section is devoted to a more complex geometry.

## OC5-semisubmersible

Similar to the study [5], in this section we examine the second order wave loads on the OC5-DeepCWind semisubmersible platform at first on the individual elements and then on the interconnected complete structure. The individual components are i) an isolated cylinder and ii) an isolated cylinder with a heave plate, that is labeled as ‘Assembly’. The interconnected components is a combination of three ‘Assemblies’, which is labeled as ‘Configuration’, resembling the full OC5-semisubmersible platform but without the smaller central columns and pontoon/braces. The dimensions of the components are: i) cylinder with a 6 m radius and 14 m draft and ii) same cylinder dimensions with a heave plate of 12 m radius and 6 m draft. For the ‘configuration’, the center of the three as-
<figure><img src="image_6.png" /><figcaption>**FIGURE 5.** BODY BOUNDARY MESHES FOR THE CYLINDER, AT THE LEFT TOP PLOT, THE 'ASSEMBLY', AT THE RIGHT TOP, AND THE 'CONFIGURATION', AT THE BOTTOM.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 6.** COMPARISON OF ADDED MASS (TOP), DAMPING COEFFICIENT (MIDDLE), EXCITATION FORCE (BOTTOM) OF SURGE FOR THE CYLINDER (LEFT) AND THE 'CONFIGURATION' BETWEEN HYDROSTAR AND NEMOH.</figcaption></figure>

assemblies are placed at $(x, y, z) = (14.43, 25, 0), (-28.87, 0, 0)$ and $(14.43, -25, 0)$ m.

Fig. 5 shows meshes of the considered structures that are used in the computation. The cylinder is discretized with 634 quadrilateral panels, the 'assembly' with 920 panels and the 'configuration' with 1970 panels. The panels are finer than the ones used in [5] and ensures the convergence of the 1st-order results with NEMOH, as shown in Fig. 6.

<figure><img src="image_8.png" /><figcaption>**FIGURE 7.** COMPARISON OF THE SURGE DIAGONAL QTF (THE REAL AND IMAGINARY PARTS) AT $\Delta f = 0.01$ Hz BETWEEN HYDROSTAR AND NEMOH. THE TOP PLOT IS FOR THE CYLINDER, THE MIDDLE PLOT IS FOR THE 'ASSEMBLY' AND THE BOTTOM PLOT IS FOR THE 'CONFIGURATION'.</figcaption></figure>

For the full QTF computation, including the free surface terms, on the cylinder structure, the free surface is discretized by a circular mesh in the radius $R_{F1}$, $6 \le R_{F1} \le 50$ m, with $N_r = 120$ points in the radial and $N_{\vartheta} = 200$ points in the angular.

The domain is assumed of infinite depth and waves with frequencies $f = (0.01, 0.21)$ Hz and $\Delta f = 0.01$ Hz are considered.

The comparison of the surge difference-frequency QTF, $|T_F^-| = |T_{FQ}^- + T_{FP}^-|$ but the $T_{FP}^-$ in Eq. (14) without the free surface forces, for the fixed structures between NEMOH and HYDROSTAR is shown in Fig. 7 for the diagonal QTF line plot and Fig. 8 for the contour plots.

Fig. 7 shows the plots of the real part and the imaginary part of the surge diagonal QTF at the surge natural frequency of the OC5 platform, $\Delta f = 0.01$ Hz. For all the structures, the real part of the QTF has a dominant magnitude compared to the imaginary part. In the diagonal QTF plot for the 'configuration'; a local peak is observed. This local peak of the QTF at the difference-frequency equal to the natural frequency of the structure is important for understanding the low frequency response of the structure.

In the contour plots, Fig. 8, the gradient of the QTF contour lines for the cylinder is increasing for the increasing frequencies, $f_1$. For the 'assembly', the additional heave plate structure causes a larger peak QTF magnitude at the largest difference frequency. In the 'configuration', the additional two columns of the 'assemblies' increase the peak QTF magnitude at the highest difference frequency.

Fig. 9 shows the difference frequency total QTF, $|T_F^-| = |T_{FQ}^- + T_{FP}^-|$ Eqs. (12) to (14), for the fixed cylinder. The effect
<figure><img src="image_9.png" /></figure>

**FIGURE 8.** CONTOUR PLOTS OF THE NORMALIZED SURGE DIFFERENCE FREQUENCY QTF MAGNITUDE (WITHOUT THE FREE SURFACE FORCES) FOR THE FIXED STRUCTURES; ON THE TOP ROW FOR THE CYLINDER, THE MIDDLE ROW FOR THE 'ASSEMBLY' AND ON THE BOTTOM FOR THE 'CONFIGURATION'. HYDROSTAR RESULTS ARE IN THE LEFT COLUMN, NEMOH RESULTS ARE IN THE MIDDLE COLUMN AND THE DIFFERENCE IN THE RIGHT COLUMN.

<figure><img src="image_10.png" /></figure>

**FIGURE 9.** CONTOUR PLOTS OF THE NORMALIZED SURGE DIFFERENCE FREQUENCY TOTAL QTF MAGNITUDE (INCLUDING THE FREE SURFACE FORCES) FOR THE FIXED CYLINDER. HYDROSTAR RESULTS ARE IN THE LEFT COLUMN, NEMOH RESULTS ARE IN THE MIDDLE COLUMN AND THE DIFFERENCE IN THE RIGHT COLUMN.

of the free surface forcing term in Eq. (14) is rather small and concentrated around the largest frequency difference.

The NEMOH results for the surge difference frequency QTF of the fixed structures are in good agreement with the HYDROSTAR results. The results also agree with the results re-

ported in [5]. This validate the full QTF calculations within NEMOH in the complex geometry.
## CONCLUSION

The full second order wave loads have been implemented in the open source potential flow code NEMOH. The full second order waves loads composed of the quadratic part and the potential part. The near-field method as in [9] is applied for calculating the quadratic part. The indirect-method as in [8,9], that using the Green-Stokes formulation and the assisting normalized radiation potential, is implemented for calculating the potential part.

Verification has been performed against an existing commercial code HYDROSTAR [14] and the good agreement is achieved for the simple geometry, the floating hemisphere, and more complex geometry, the OC5-DeepCwind semisubmersible platform.

In the floating hemisphere case, we verified the difference-frequency full QTF and the quadratic sum-frequency QTF. The quadratic sum-frequency QTF Nemoh results is verified with the data presented in [10]. The full sum-frequency QTF was not shown because the results were influenced by the irregular frequency. An improvement of the irregular frequency removal method in NEMOH is needed.

In the OC5-DeepCwind semisubmersible platform case, we verified the difference-frequency QTF for the three different components: the cylinder, the ‘assembly’ and the ‘configuration’, of the fixed platform. We also verified the difference-frequency full QTF, including the free surface forces, for the cylinder.

In this paper we compared and showed only the QTF in the surge axis. We conducted also comparison of the QTF for heave and pitch axes between NEMOH and HYDROSTAR. The same conclusion as for the surge QTF, the good agreement between the two codes is achieved.

The NEMOH with the QTF module will be soon released publicly, which will be, to our knowledge, the sole and only open-source software that provides the second order module.

## ACKNOWLEDGMENT

This work was done within the framework of the FLOAT-ECH project. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 101007142.

## REFERENCES

[1] Babarit, A., and Delhommeau, G., 2015. "Theoretical and numerical aspects of the open source BEM solver NEMOH". In Proceedings of the 11th European Wave and Tidal Energy Conference.

[2] Philippe, M., Adrien Combourieu, C. P., Robaux, F., Delhommeau, G., and Babarit, A., 2015. "Introducing Second Order Low Frequency Loads in the Open-Source Boundary Element Method Code Nemoh". In Proceedings of the 11th European Wave and Tidal Energy Conference.

[3] Schmitt, P., Windt, C., Nicholson, J., and Elsässer, B., 2016. "Development and validation of a procedure for numerical vibration analysis of an oscillating wave surge converter". *European Journal of Mechanics - B/Fluids, 58*, pp. 9-19.

[4] Penalba, M., Kelly, T., and Ringwood, J. V., 2017. "Using nemoh for modelling wave energy converters: A comparative study with wamit". In The 12th European Wave and Tidal Energy Conference At Cork, Ireland.

[5] Tom, N., Robertson, A., Jonkman, J., Wendt, F., and Böhm, M., 2019. "Bichromatic Wave Selection for Validation of the Difference-Frequency Transfer Function for the OC6 Validation Campaign". In International Conference on Offshore Mechanics and Arctic Engineering, Vol. ASME 2019 2nd International Offshore Wind Technical Conference. V001T01A022.

[6] WAMIT, Inc, 2011. WAMIT SOFTWARE v6.1.

[7] Kim, M., 1991. "Second-order sum-frequency wave loads on large-volume structures". *Applied Ocean Research, 13*(6), pp. 287-296.

[8] Molin, B., 1979. "Second-order diffraction loads upon three-dimensional bodies". *Applied Ocean Research, 1*(4), pp. 197-202.

[9] Chen, X.-B., 1988. "Etudes des réponses du second ordre d'une structure soumise à une houle aléatoire". PhD Thesis, Université de Nantes.

[10] Kim, M.-H., and Yue, D. K. P., 1990. "The complete second-order diffraction solution for an axisymmetric body part 2. bichromatic incident waves and body motions". *Journal of Fluid Mechanics, 211*, p. 557-593.

[11] Maruo, H., 1960. "The drift of a body floating on waves". *journal of Ship Research*.

[12] Chen, X.-B., 2007. "Middle-field formulation for the computation of wave-drift loads". *Journal of Engineering Mathematics, 59*, pp. 61-82.

[13] Rouault, M., Molin, B., and Ledoux, A., 2007. "A Lagally formulation of the second order slowly varying drift forces". In 11ème Journées de l'Hydrodynamique, pp. ISSN 1161-1847.

[14] Bureau Veritas, 2020. HYDROSTAR SOFTWARE v8.14.

[15] Geuzaine, C., and Remacle, J.-F., 2021. GMSH SOFTWARE v4.8.4.

[16] Kim, M.-H., and Yue, D. K. P., 1990. "The complete second-order diffraction solution for an axisymmetric body part 2. bichromatic incident waves and body motions". *Journal of Fluid Mechanics, 211*, p. 557-593.