Numerical coupling strategy using
HOS-OpenFOAM-MoorDyn for OC3 Hywind SPAR
type platform

Sithik Aliyar, Guillaume Ducrozet, Benjamin Bouscasse, Félicien Bonnefoy, V.

Sriram, Pierre Ferrant

► To cite this version:

Sithik Aliyar, Guillaume Ducrozet, Benjamin Bouscasse, Félicien Bonnefoy, V. Sriram, et al.. Numerical coupling strategy using HOS-OpenFOAM-MoorDyn for OC3 Hywind SPAR type platform. Ocean Engineering, 2022, 263, pp.112206. 10.1016/j.oceaneng.2022.112206. hal-04491026

HAL Id: hal-04491026

https://hal.science/hal-04491026v1

Submitted on 27 May 2025

**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
# Numerical coupling strategy using HOS-OpenFOAM-MoorDyn for OC3 Hywind SPAR type platform

Sithik Aliyar<sup>a,b</sup>, Guillaume Ducrozet<sup>b</sup>, Benjamin Bouscasse<sup>b</sup>, Félicien Bonnefoy<sup>b</sup>,
Sriram Venkatachalam<sup>a</sup> and Pierre Ferrant<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Numerical modelling

CFD

OpenFOAM

Floating bodies

SPAR

Moorings

FOWT

MoorDyn

nonlinear wave structure interactions

## ABSTRACT

The coupled (Potential theory and Navier Stokes) solver is applied to simulate the interaction of sea waves with substructure of floating offshore wind turbines (FOWT) platform, notably similar to the OC3 Hywind SPAR structure. The intention is to develop a numerical tool that allows the study of the survivability of floating structures in extreme sea states. In this study, the moorings are modelled in two ways. One is by considering the mooring lines as a linear spring with defined spring stiffness, and another is by coupling the solver(foamStar) with a lumped-mass mooring dynamics model (MoorDyn). MoorDyn represents mooring line behaviour subject to axial elasticity, hydrodynamic forces, and vertical contact forces with the seabed. The coupled model has been validated against the experiments carried out as part of the SOFTWIND project. The numerical model results of free surface elevation, floating body motions and mooring tensions are compared with the experiments. For wave cases with mild and moderate amplitudes, mooring in the form of a stiffness matrix is sufficient. However, dynamic mooring simulation (MoorDyn) is required for the extreme sea state conditions.

## 1. Introduction

Floating Offshore Wind Turbines (FOWT) technology is a new trend that allows wind turbines to be installed in larger, deeper offshore locations for increased wind potential. But the accurate prediction of floater behavior in such a dynamic environment is still a key challenge (Pinguet, 2021; Huang et al., 2021; Zhang and Kim, 2018). Wave-structure interaction problems can be assessed by empirical formulations, model tests in wave tanks, and numerical methods. The most accurate method is model testing, which provides a direct comparison to physical reality within the limitations of experimental error and wave tank fidelity. However, it is rarely possible to employ a real-scale test model and accurately calculate critical loads for the complex multi-physics problem. There is a restriction in choosing the scaling factors based on the experimental facilities available in model testing (Aliyar et al., 2021). These scaling effects cannot be ignored in predicting the elemental forces acting over the wind turbine models (Make, 2014). Also, moored FOWT physical models add multiple uncertainties like the shifts in the equilibrium position, mooring dynamics, motion effects, etc., which are difficult to control during the experiment. Numerical techniques without scaling limitations, on the other hand, may improve model reliability by lowering uncertainties. Based on whether the fluid viscosity is considered or not, the numerical wave tank (NWT) simulations can be grouped into potential theory-based NWT (PNWT) and Navier Stokes (NS) equation based CFD-NWT (CNWT). Historically, wave-structure interaction was modelled using linear potential flow theory with the boundary element method (BEM), when a smooth geometry and negligible viscous effects prevail (Black et al., 1971; Sun, 2007). However, when flow separation occurs around the body and nonlinear behaviour of a floating structure is the point of interest, assumptions of irrotational and inviscid flow in PNWT will not represent a realistic scenario for reproducing model testing. Multiphase CNWT, together with specific wave generation and absorption techniques (Jacobsen et al., 2012; Higuera et al., 2013), can reproduce open sea and experimental conditions (Elhanafi et al., 2017; Gu et al., 2018; Devolder et al., 2018; Guilcher et al., 2007; Miquel et al., 2018). A large computational domain is required to develop nonlinear wave evolution to avoid boundary effects

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

<footnote><sup>a</sup>Dept of Ocean engineering, IIT Madras, Chennai, 600036, Tamil Nadu, India</footnote>

<footnote><sup>b</sup>Nantes Université, École Centrale de Nantes, CNRS, LHEEA, UMR 6598, Nantes, F-44000, France</footnote>

<footnote>* This document is the results of the research project funded by the Joint doctoral fellowship by IITM-ECN</footnote>

<footnote>Sithik.Aliyar@ec-nantes.fr (S. Aliyar); guillaume.ducrozet@ec-nantes.fr (G. Ducrozet); benjamin.bouscasse@ec-nantes.fr (B. Bouscasse); felicien.bonnefoy@ec-nantes.fr (F. Bonnefoy); vsriram@iitm.ac.in (S. Venkatachalam); pierre.ferrant@ec-nantes.fr (P. Ferrant) ORCID(s): 0000-0003-0188-4613 (S. Aliyar)</footnote>
resulting from waves interacting with any fixed or floating structures in the domain. Also, high temporal and spatial discretization is mandatory in CFD solvers to reproduce the exact events and avoid excessive numerical dissipation during wave propagation. These factors will significantly increase the computational costs, resulting in the conventional CFD approaches being impractical.

In recent years, many hybrid coupling methods have appeared and grown in popularity to overcome the weakness of both PNWT and CNWT in the name of domain decomposition (DD). Generally, the complex interaction (viscous effects and violent free-surface deformation) only occurs near the structure, and in the far field, these effects can often be neglected. For this reason, it is viable to divide the computational domain into a viscous inner sub-domain and an inviscid outer sub-domain. Potential flow models solve the flows in the outer subdomain, and CFD models solve the inner subdomain. Information between the models is shared in a common region, either in one-way or two-way interaction. Many papers have published the latest advancements in domain decomposition, devoted to one way (Hildebrandt et al., 2013; Lachaume et al., 2003; Biausser et al., 2004; Paulsen et al., 2014; Choi et al., 2020; Agarwal et al., 2021b)and two-way coupling the (Sriram et al., 2014; Lu et al., 2017; Kemper et al., 2019). For more details about the different hydrodynamic coupling, (Di Paolo et al., 2020) detailed review can be referred to. In addition, some articles portray the robustness of DD for coastal engineering applications (Lachaume et al., 2003; Biausser et al., 2004; Sriram et al., 2014) and bottom fixed support structures (Hildebrandt et al., 2013; Paulsen et al., 2014; Li et al., 2021; Sriram et al., 2021; Agarwal et al., 2021a), but its application to the floating structures is the recent research interest among the researchers (Xu et al., 2021; Li and Bachynski, 2021).

Several studies have used CFD to model the floating structures in single or multiple degrees of freedom (DoF) (Elhanafi et al., 2017; Simonsen et al., 2013; Dunbar et al., 2015). The assessment of moored floating structures needs interdisciplinary expertise in hydrodynamics, floating body dynamics, mooring, and their interaction. Three types of mooring analysis have been commonly adopted, static, quasi-static, and dynamic analysis. Static analysis is often used in preliminary design to compute the static response of moored vessels (Karnoski and Palo, 1988). The offshore sector has a tradition of designing the mooring system through quasi-static analysis, justified by the low responsiveness of the large masses and corresponding low velocities. In operational conditions, this assumption might be sufficient (Thomsen et al., 2017). As the quasi-static model bypasses hydrodynamic drag and inertial forces on mooring lines, it underpredicts line tensions under extreme conditions. Also, in deep waters, a mooring line's response may be delayed, causing it to violate the catenary profile assumptions of a quasi-static model (Vickers, 2012). Thus, employing a dynamic mooring model to reliably predict loads acting on mooring lines and the floating structure may be advantageous. MDD (Dewey, 1999), SEAWAY (Journee and Adegeest, 2003), OrcaFlex (Randolph and Quiggin, 2009), Moody (Ferri and Palm, 2015), and MoorDyn (Hall, 2015) are some of the mooring dynamics models developed in recent years. But there are only a few papers in which CNWT is used in conjunction with the above library to solve mooring line tensions (Nematbakhsh et al., 2015; Benitz et al., 2015; Palm et al., 2016; Jiang et al., 2021; Lee et al., 2021). In addition, most of the research so far has used a floating box, buoy, or cylinder as a floating structure for validation, and there is a lack of validation for real-time floating systems like SPAR, TLP, or Semi-submersible for industrial applicability.

In this study, we have used *foamStar* (Choi et al., 2020), an in-house code co-developed by Bureau Veritas and Ecole Centrale de Nantes in the OpenFOAM framework for solving wave-structure interaction problems. It is based on the standard multi-phase solver *interDymFoam* in the OpenFOAM package with two special modules. One is wave generation, similar to waves2Foam (Jacobsen et al., 2012) but it is further developed with a faster implementation and to allow the coupling with the nonlinear potential open-source solvers HOS-Ocean and HOS-NWT. Another is the mechanical solver designed to solve the 6 DoF dynamics of rigid and flexible floating/fixed bodies in waves. In-build modules can support the simulation for the decay test, regular and irregular waves interaction with the structure, following one way DD approach. The solver has been developed and validated for regular and irregular wave propagation (Choi et al., 2020; Li et al., 2021) and other applications which can be referred at (Kim, 2021; Choi, 2019; Li, 2018).

This paper aims to implement a static and dynamic mooring model (MoorDyn (Hall, 2015)) into *foamStar* to conduct a fully coupled analysis of moored floating structures replicating the model test in NWT. MoorDyn is a lumped-mass mooring model that has been associated with a potential flow solver and successfully validated against measurements in a wave tank (Hall and Goupee, 2015). The equation governing dynamics of mooring lines is coupled to the six DoF nonlinear rigid body motions equations in *foamStar*. MoorDyn calculates the tensions of the mooring lines, which are added as an external force to compute the resulting response and motion of the moored floating structures. Based on our knowledge, this is the first example of coupling between a DD-based solver (here *foamStar*), floating
body dynamics, and a dynamic mooring model with a systematic validation of the structural motions and mooring tensions with the experiment.

Validation of the two-way coupling between the *foamStar* and *MoorDyn* is carried out with the experimental results on the floating SPAR with a delta mooring system. The moored decay tests are examined to ensure the computational model performance in comparison with the experimental model in terms of its natural period, mooring tension, and damping properties. The SPAR is then subjected to regular and irregular waves in the following study, and simulations are compared to experimental data, verifying the developed coupling model's ability to effectively forecast both body motions and mooring line tensions. The paper has been organized in the following manner. Section 2 discusses the governing equations underpinning the *foamStar* and *MoorDyn* models, as well as the framework of coupling between the solvers. Section 3 shows the general setup and parametric studies used in this study. Section 4 shows moored floating structures validation with both regular and irregular waves, and the paper concludes with a summary of the hybrid model's capability and the performance.

## 2. Numerical Methods

In the present work, the one-way domain decomposition approach (Di Paolo et al., 2020) using *foamStar* (Choi et al., 2020) as the CFD solver and CN-Stream, HOS-NWT as the potential solvers has been adopted as shown in the Figure 1. In simple terms, the potential solver's velocity, pressure, and surface elevation will be transferred into the viscous solver through the inlet relaxation zone, and all parameters from the viscous flow domain will be blended to zero in the outlet relaxation zone. This section provides the governing equations behind both solvers, followed by the adopted coupling approach. Next, the mooring models adopted in the solver are presented, and its coupling methodology has been briefly discussed.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** General description of the Numerical Wave Tank (NWT) under one way Domain decomposition strategy between potential and viscous model</figcaption></figure>

### 2.1. *foamStar*

The open-source solver interDymFoam (of OpenFOAM) is used to solve the multiphase problem by coupling the Navier-Stokes equations with a Volume of Fluid (VoF) method (Hirt and Nichols, 1981). The computational domain is subdivided into a finite number of control volumes(CV), which can be unstructured. The integrals over each CV are numerically approximated using the midpoint rule. The pressure Poisson equation is solved based on a PIMPLE algorithm, which combines the Semi-Implicit Pressure Linked Equations (SIMPLE) algorithm with the Pressure Implicit with Splitting of Operators (PISO) algorithm. The unknown variables at the center of a cell face are determined by combining a central differencing scheme (CDS) with an upwind differencing scheme (UDS). The *k-omega* SST turbulence model was used for the turbulence requirement, however, it was discovered after the case study that it had the least impact on the cases picked, hence it is not included in this study.

The governing equations are based on mass and momentum conservation equations written as follows,

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$
$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla \cdot (\rho\mathbf{u}\mathbf{u}^T) - \nabla \cdot (\mu(\nabla\mathbf{u} + \nabla\mathbf{u}^T)) = -\nabla p_d - (\mathbf{g} \cdot \mathbf{x})\nabla\rho \quad (2) $$

where $\mathbf{u}, \mathbf{x}$ and $\mathbf{g} = [0, 0, -g]^T$ are the fluid velocity, position vector and gravitational acceleration vector, respectively. The dynamic pressure $p_d = p - \rho\mathbf{g} \cdot \mathbf{x}$ is as introduced in Rusche (2002). The free surface is identified using the volume of fluid (VoF) method, where the two-phase problem is treated as a single fluid with a volume-fraction parameter $\alpha \in [0, 1]$ subjected to transport equation,

$$ \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{u}) + c_{\alpha} \nabla \cdot (\mathbf{u}_n \alpha (1 - \alpha)) = 0 \quad (3) $$

where the third term on the left-hand-side is the artificial compression to keep the interface sharp. $\mathbf{u}_n$ is the fluid velocity normal to the interface and $c_\alpha$ is the compression coefficient. In each cell, then fluid properties (density and viscosity) are computed as a mixture between air ($\alpha = 0$) and water ($\alpha = 1$):

$$ \rho = \alpha \rho_w + (1 - \alpha) \rho_a \qquad \mu = \alpha \mu_w + (1 - \alpha) \mu_a \quad (4) $$

where index $w$ and $a$ indicate water and air respectively.

## 2.2. Wave models and coupling

The incident wave generation is based on the potential theory approach. With the assumption of an ideal fluid and irrotational flow, it is possible to define a velocity potential ($\phi$) as,

$$ \mathbf{u}_I = \nabla \phi \quad (5) $$

Then the incident wave potential satisfies Laplace's equation in the fluid domain ($\Omega$)

$$ \Delta\phi = 0 \quad x \in \Omega \quad (6) $$

The non linear kinematic and dynamic free surface boundary conditions for incident waves are given as,

$$ \frac{\partial \eta_I}{\partial t} + \frac{\partial \phi}{\partial x} \frac{\partial \eta_I}{\partial x} + \frac{\partial \phi}{\partial y} \frac{\partial \eta_I}{\partial y} + \frac{\partial \phi}{\partial z} = 0 \quad \text{at } z = \eta_I(x, t) \quad (7) $$

$$ \frac{\partial \phi}{\partial t} + \frac{1}{2}(\nabla \phi^2) + gz = 0 \quad \text{at } z = \eta_I(x, t) \quad (8) $$

where $x$ and $y$ are horizontal coordinate, $t$ is time. The flat bottom condition is given as,

$$ \frac{\partial \phi}{\partial z} = 0 \quad \text{at } z = -h \quad (9) $$

where $h$ is water depth. In the present work, spectral methods are used to solve the potential and surface elevation. The dynamics of non-breaking gravity waves are accurately modeled using the potential flow theory. The numerical methods at use in the present study are based on the stream function theory for regular waves and the Higher Order Spectral (HOS) method for irregular waves. Both methods are briefly presented in this section followed by coupling with foamStar.
### 2.2.1. Regular Waves: Stream function wave theory

The open-source CN-Stream model (the library is published with GNU General public license and can be downloaded from the URL https://github.com/LHEEA/CN-Stream) is used to describe nonlinear regular waves propagating over a flat bottom. It is based on the stream function theory (Rienecker and Fenton, 1981) that allows an efficient and accurate solution of regular wave trains up to the breaking wave limit. The model assumes periodic boundary conditions and solves the propagation in a moving frame, making the problem steady. Spatial discretization is achieved thanks to a spectral method that decomposes the unknowns (free surface elevation and velocity potential) on a set of basis functions that satisfy intrinsically Laplace equation (Eqn. (6)) as well as the bottom boundary condition (Eqn. (9)). The modal amplitudes are numerically evaluated to satisfy the free surface boundary conditions. They can be used to reconstruct the wave field (velocity and pressure through Bernoulli's equation (Eqn. (8)) at required locations. The reader can find more information in (Rienecker and Fenton, 1981; Ducrozet et al., 2019).

### 2.2.2. Irregular Waves: Higher Order Spectral method

Stream function methodology is limited to regular waves. In the present work, we are addressing another type of unidirectional wave, irregular sea states (refer Section 3.2.5). With the objective of an accurate and computationally efficient numerical model, the High-Order Spectral (HOS) method (Dommermuth et al., 1988; West et al., 1987) has been adopted to describe the complex sea states. This numerical approach relies on a pseudo-spectral discretization scheme that exhibits a high (exponential) convergence rate and that allows the use of Fast Fourier Transforms, which ensures high efficiency.

The original HOS method (Dommermuth et al., 1988; West et al., 1987) was dedicated to the time evolution of an initially prescribed wave field in an open domain with periodic boundary conditions. The present study intends to replicate experiments carried out in a wave tank where the waves are generated thanks to a wavemaker in a finite-size domain. Then, in this work, a dedicated Numerical Wave Tank (NWT) based on HOS is used (Ducrozet et al., 2012). This open-source solver named HOS-NWT has been developed at LHEEA Lab, (Ecole Centrale Nantes (ECN)) and now constitutes a Digital Twin of the ECN wave tank facility. The geometry is the exact physical one, including perfectly reflective sidewalls. The initial conditions are at rest and the waves are generated thanks to a wavemaker. For direct comparison with experiments, the exact same wavemaker motion is used in the numerical and the physical wave tank. In addition, a numerical absorbing zone mimics the physical beach located close to the wall opposite the wavemaker to prevent wave reflection. This NWT has been applied and validated on different wave configurations of various levels of complexity (Ducrozet et al., 2012, 2016; Michel et al., 2020; Canard et al., 2020).

### 2.2.3. Relaxation Zone - Coupling between solvers

Because of the unsteady nature of the water waves and the moving free surface, CFD solvers must have appropriate temporal and spatial discretization to accurately simulate wave-related problems. Furthermore, practical marine and ocean engineering applications require a sufficient simulation time (to eliminate the transient flow regime at the start of the simulations and achieve a periodic condition). Also, the NS solver is known for its energy loss in wave propagation, which dampens the waves in long duration and its prohibitive computational cost. Hence the objective is to couple the viscous solver with the potential solvers to eliminate its shortcomings, and this coupling happens in the relaxation zone (Figure 1). These relaxation zones are used to prevent waves reflected from outlet boundaries from interfering with wavemaker boundaries, as well as to attenuate the waves reflected internally in the computational domain (Jacobsen et al., 2012).

A simple way of generating a wave based on one-way DD (foamStar) is to impose wave velocity and free surface elevation from potential wave models at the inlet of the CFD domain. In foamStar, the wave generation and absorption method are based on an explicit scheme that relaxes the computed solution towards a given target flow field (Jacobsen et al., 2012; Seng, 2012; Mayer et al., 1998; Engsig-Karup, 2006; Choi et al., 2020), and the scheme requires a weight function that varies between 0 and 1 in the relaxation zone. It is achieved as follows:

$$ \zeta = w\zeta_{\text{target}} + (1 - w)\zeta_{\text{CFD}} \quad (10) $$

where $\zeta$ is either the velocity field(u) or the VOF field ($\alpha$), $w$ is the weight function, and $\zeta_{\text{target}}$ will be the given target flow field usually set to be either the incident wavefield or calm-water (to use relaxation as an absorbing zone). Two common choices of weight function are the polynomial weight function (Mayer et al., 1998; Engsig-Karup, 2006; Choi et al., 2020) and the exponential weight function (Mayer et al., 1998; Engsig-Karup, 2006; Jacobsen et al., 2012;
Seng, 2012; Choi et al., 2020). The present study adopts the explicit relaxation scheme using the polynomial weight function:

$$w(x) = -2x^3 + 3x^2 \qquad (11)$$

Where $x \in [0, 1]$, is the normalised coordinate in the relaxation zone. As mentioned, the above technique can be applied in the wave absorption zone as if target values ($\zeta_{target}$) are set to no waves, then the entire wavefield would be damped, and if the target values are set to incident waves, then only the scattered waves are damped. Choi et al. (2020) can be referred for some examples of 3D wave generation of realistic sea spectrum. To summarize, In the present case, the relaxation zone explicitly blends the wave properties from CN-Stream or HOS-NWT within the foamStar domain. In the case of irregular waves (HOS-NWT) wave properties are blended using the Grid2Grid wrapper program, which uses inverse FFTs and interpolation technique to reconstruct the wave amplitude and the velocity field at the inlet for required time and spatial locations (Choi et al., 2017).

## 2.3. Rigid body motions

There are two reference frames used in the present study, the Inertial ($R_0$) and the body reference frame($R_b$) (refer Figure 5). All the frames are based on Cartesian systems and defined according to the right-hand rule. The Inertial reference frame's origin is positioned at the center of the domain coinciding with the body, where the body-fixed frame origin is set at the body's Center of gravity. A positive surge is defined forward, and a positive heave is defined upwards. The rotation angles need to be defined carefully since the transformation between the global reference frame and the body-fixed frame is not unique. The foamStar library supports the Euler angle and quaternion representations among different methods of transformation.

The present work follows the definition of Euler intrinsic angles with zyx convention. The translation and rotational equations of motion are expressed in the body-fixed frame as follows:

$$[m] \ddot{q} = Q_e + Q_v + Q_c - [K]q - [C]\dot{q} \qquad (12)$$

$$[I] \ddot{\theta} = M_e + M_v + M_c - [K]_{rot} \theta - [C]_{rot} \dot{\theta} \qquad (13)$$

Where $m$, $I$, $K$, $C$ are the generalized mass, Inertia, stiffness, and damping matrix respectively; $q$ and $\theta$ represents the translation and rotational motions; $Q_e$ and $M_e$ are the external forces and moments acting on the body; $Q_v$ and $M_v$ are the quadratic velocity forces (Coriolis/centrifugal forces) and moments; $Q_c$ and $M_c$ represents the constraint forces and moments. Gravity, mooring, propulsion, and fluid are some of the components of external forces and moments. Among them, fluid force is evaluated as the surface integration of pressure and viscous stress over the body surface. Here, $n$ is a unit surface normal vector pointing inside the body and $\tau$ indicates the viscous stress tensor.

$$Q_{e,fluid} = \int_{body} (pn - \tau \cdot n) dS \qquad (14)$$

## 2.4. Mesh morphing

In general, any problem in computational fluid dynamics (CFD) implies the use of a grid or mesh to replace a continuous problem domain with a finite field. This model works flawlessly when the problem does not involve any solid body motion in the computational domain. However, many offshore engineering applications require solid body motion in the computational domain, which has a consequence on the fluid flow. Sliding interface, Overset mesh adaptation, and Mesh morphing are the three different ways dynamic mesh motion operates in OpenFOAM. A sliding interface allows one or more internal subdomains to translate or rotate relative to a static background mesh domain, with the interface between the domains sliding relative to each other. The main disadvantage of sliding interfaces is that they are only appropriate for single DoF motions. A complex motion can be simulated using a mesh overset technique, but it may prove very computationally expensive mainly due to the mesh refinement and adaptation, the remapping, and the re-computation. This is not necessary if the mesh is only deformed and the mesh's topology remains unchanged. Mesh morphing is the most commonly used method in CFD floating body simulations. The method deforms meshes into new
meshes by updating the node positions near the floating body while preserving the shape of the target floating body. Because only the nodal positions are updated, the necessary computation time is shorter than that of the sliding interface and overset mesh methods. However, in the case of large deformations, solutions may diverge due to the skewness of the mesh, and meshes that have been optimized for a generated wave signal may also get impacted. Assuming that the amplitude of the floating body motion is relatively small in the present case, the mesh morphing deformation method is adopted.

Different morphing techniques can be applied within the finite volume framework (Jasak and Tuković, 2010) like Spring analogy, Pseudo-solid, Laplacian smoothing, Algebraic interpolation, solid body motion, etc. The mesh update strategy is based on the Laplacian mesh smoothing technique in the present work and the mesh motion of the computational domain is calculated by solving the cell center Laplace smoothing equation (Jasak and Tukovic, 2006). The Laplacian equation is given as,

$$ \nabla \cdot (\gamma(r)\nabla\mathbf{x}) = 0; \quad \gamma(r) = \frac{1}{r^2} \qquad (15) $$

where $\mathbf{x}$ is the point displacement field and $\gamma$ is the diffusivity coefficient. $r$ is the distance to the moving boundary. Several choices for $\gamma$ are presented in (Löhner and Yang, 1996; Jasak and Tuković, 2010). In the present work, the inverse quadratic distance relation has been chosen as the coefficient will be significant, and the gradient will be small, yielding a small relative displacement for the points close to the moving boundary, preserving the cells' quality near the moving boundaries (Jasak and Tukovic, 2006).

## 2.5. Mooring system

Mooring systems are widely used in offshore engineering for station keeping of floating bodies in various sea states. In the present study, two different mooring systems are adopted. One is based on the static mooring condition, and another approach is based on the dynamic mooring condition. Both are discussed briefly.

### 2.5.1. Static mooring model: Stiffness Matrix

The external mooring used to control the motions in six DoF is simplified as the linear spring model with stiffness matrix $[K]$. The restoring force is given as,

$$ [Q_{e,M}] = [K_M][\delta q] \qquad (16) $$

where $Q_{e,M}$ represents external force and moment from mooring, and $\delta q$ represents the change in motion (in six DoF). This force is added to the right-hand side of the equation of motion (Eqn.(13)) that can be interpreted as attaching a spring to the body's center of gravity. The working procedure of the stiffness matrix method inside the foamStar algorithm is depicted in the general flowchart shown in the Figure 3.

### 2.5.2. Dynamic mooring model: MoorDyn

The dynamic mooring solver MoorDyn (Hall, 2015) is a lumped-mass mooring model that lumps all the effects of mass, external forces, and inertial reactions at a finite number of nodes, N, along the mooring line connected by massless spring as shown in the Figure 2. It follows the right-handed inertial reference frame defined with the z-axis being measured positive up from the water plane. The location of each node i is defined by the position vector $\mathbf{r}_i$ which includes $[X_i, Y_i, Z_i]$ components of the vector. The equation of motion for each node i is written as follows:

$$ (\mathbf{m}_i + \mathbf{a}_{pi} + \mathbf{a}_{qi})\ddot{\mathbf{r}}_i = \mathbf{T}_i + \mathbf{C}_i + \mathbf{W}_i + \mathbf{B}_i + \mathbf{D}_{pi} + \mathbf{D}_{qi} \qquad (17) $$

where $\mathbf{m}_i$ is the mass matrix for node i. The mass of the mooring cable is discretized into point masses at each node by allotting each node half the combined mass of the two neighboring cable segments.

The mooring line is represented as a cylinder of diameter $d$, unstretched length $L$, and density $\rho$. For a mooring line with $N$ segments connecting $N + 1$ nodes, the unstretched length of each line segment is $l_0 = L/N$. The forces acting on the nodes are divided into internal and external forces. The internal forces include net buoyancy ($W_i$), tension ($T_i$) and internal damping ($C_i$) forces. The external forces contain transverse and tangential hydrodynamic forces ($D_i$) based on Morison's equation, added mass force ($a_{pi}$ and $a_{qi}$) and the seabed contact force ($B_i$). The tension, $T_{i+1/2}$ in
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Mooring line discretization in MoorDyn with a representation of internal and external forces in inertial reference frame $R_0$ (Hall and Goupee, 2015)</figcaption></figure>

any line segment is calculated as,

$$ T_{i+1/2} = \begin{cases} E \frac{\pi}{4} d^2 \left( \frac{l}{l_0} - 1 \right) & \text{if } l > l_0 \\ 0 & \text{if } l \le l_0 \end{cases} \qquad (18) $$

Where $E$ is Young's modulus. The calculations assume that line tensions are always positive $l > l_0$ (compression line forces are not allowed). By multiplying a unit vector directed along a line segment, the above equation can be rewritten in vector form as

$$ \mathbf{T}_{i+(1/2)} = \begin{cases} E \frac{\pi}{4} d^2 \left( \frac{\|\mathbf{r}_{i+1}-\mathbf{r}_i\|}{l_0} - 1 \right) \left( \frac{\mathbf{r}_{i+1}-\mathbf{r}_i}{\|\mathbf{r}_{i+1}-\mathbf{r}_i\|} \right) & \text{if } \|\mathbf{r}_{i+1}-\mathbf{r}_i\| > l_0 \\ 0 & \text{if } \|\mathbf{r}_{i+1}-\mathbf{r}_i\| \le l_0 \end{cases} \qquad (19) $$

The net buoyancy ($W_i$) at node i assuming the force is evenly divided among two connecting nodes is described in vector form as follows:

$$ \mathbf{W}_i = \hat{\mathbf{e}}_z \frac{W_{i-\frac{1}{2}} + W_{i+\frac{1}{2}}}{2} \qquad (20) $$

Where

$$ W_i = \frac{\pi}{4} d^2 l (\rho_{lw} - \rho) g \qquad (21) $$

Where $\rho_{lw}$ is the density of water and $\rho$ is the density of the mooring line. In order to improve numerically stability, a numerical internal damping force is included:

$$ \mathbf{C}_{i+(1/2)} = C_{\text{inter}} \frac{\pi}{4} d^2 \dot{\mathbf{e}}_{i+(1/2)} \left( \frac{\mathbf{r}_{i+1} - \mathbf{r}_i}{\|\mathbf{r}_{i+1} - \mathbf{r}_i\|} \right) \qquad (22) $$

Where $C_{\text{inter}}$ is the numerical internal damping coefficient and $\dot{\mathbf{e}}_{i+1/2}$ is the strain rate of segment. The linear spring damper model is activated only when the node interacts with the sea bottom ($z_{\text{bot}}$) and it is given as:

$$ \mathbf{B}_i = dl[(z_{\text{bot}} - z_i)k_b - \dot{z}_i c_b] \hat{\mathbf{e}}_z \qquad (23) $$

where $k_b$ is the seabed stiffness coefficient, and $c_b$ is the seabed damping coefficient.
The drag and added mass are calculated at each node i based on Morison's equation. The direction of a line passing between two adjacent nodes is approximated by the tangential direction $\hat{q}_i$ (Figure 2):

$$ \hat{q}_i = \frac{\mathbf{r}_{i+1} - \mathbf{r}_i}{||\mathbf{r}_{i+1} - \mathbf{r}_i||} \qquad (24) $$

Transverse and tangential drag forces are calculated using the transverse drag coefficient $C_{dn}$, and the tangential drag coefficient, $C_{dt}$ respectively. Similarly, transverse and tangential added mass forces are calculated using the transverse added mass coefficient $C_{an}$, and the tangential added mass coefficient, $C_{at}$ respectively. Both the forces are computed as follows:

$$ \mathbf{D}_{pi} = \frac{1}{2} \rho_w C_{dn} dl \left\| (\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i - \ddot{\mathbf{r}} \left\| [(\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i - \ddot{\mathbf{r}}] \right\| \right\| \qquad (25) $$

$$ \mathbf{D}_{qi} = \frac{1}{2} \rho_w C_{dt} dl \left\| (-\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i \right\| \left\| [(-\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i] \right\| \qquad (26) $$

$$ \mathbf{a}_{pi} \ddot{\mathbf{r}}_i = \rho_w C_{an} \frac{\pi}{4} d^2 l [(\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i - \dot{\mathbf{r}}_i] \qquad (27) $$

$$ \mathbf{a}_{qi} \ddot{\mathbf{r}}_i = \rho_w C_{at} \frac{\pi}{4} d^2 l (-\ddot{\mathbf{r}} \cdot \hat{\mathbf{q}}_i) \hat{\mathbf{q}}_i \qquad (28) $$

The second-order system of ordinary differential equations presented in Eqn. (17) can be reduced to a system of first-order differential equations, and then solved using a second-order Runge-Kutta integration scheme with a constant time step

## 2.6. foamStar-MoorDyn coupling

MoorDyn is coupled to the foamStar solver in two way coupling approach where the information between the solvers is passed at each time step, and each solver will influence the solutions of the other. The core functions in MoorDyn are compiled into a shared library to use the dynamic loading of the functions inside the foamStar. Generally, there are two ways of coupling between the solvers that can be adopted, strong and weak coupling. Each solver is solved separately in the weak coupling, and some variables are exchanged between the solvers. Also, it can adopt different time step sizes for each solver. However, in the strong coupling, both the solvers save all the variables in the common solver, and each solver only composes its governing equation. The present work has adopted the weakly coupling between the foamStar and MoorDyn that it implies that mooring lines will influence the platform's motion but not the flow field. The coupling procedure follows three phases: initialization, information sharing, and closing.

**Initialization phase:** Before the time loop begins, the solvers are being set up, and the mooring system is created. An initialization function is then called calculating its initial conditions based on body position specified by the equation of the motion.

**Information sharing phase:** As time progress, foamStar transfers the six degrees of freedom displacement (x), velocities (xv), and current timestep size into the MoorDyn. MoorDyn calculates the fairleads' position, velocities, and accelerations based on the input. Along with that, the mooring system tensile load and restoring forces are calculated and shared with the foamStar as one of the external forces [$Q_{e,M}$]. Then the foamStar calculates the platform motions, accelerations etc., and the loop continues till the last time step.

**Closing phase:** The last stage releases the variables and memories used to calculate the tensile load of the mooring system after the end of calculations.

It should be noted that foamStar uses a body reference system ($R_b$) for its governing equation, and MoorDyn uses an inertial frame of reference ($R_0$) for its governing equation. Hence, body displacements and velocities are transferred to the inertial reference frame at each timestep and mooring tension to the body reference frame. The detailed framework of the coupling procedure is shown in the Figure 3.
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Computational algorithm of coupling between the potential theory wave models, *foamStar*, and *MoorDyn* solver</figcaption></figure>

# 3. General Setup and parametric study

## 3.1. Physical model

Arnal (2020), experimentally tested a 1:40 scale 10 MW SPAR floating wind turbine (FOWT) in the ECN wave basin. The basin (50 m x 30 m x 5 m) is fitted with a 48-hinged flap-wavemaker on its shorter side for wave generation. A beach at the opposite end of the basin absorbs the incoming wave energy and reduces the reflections. The basin can generate regular waves up to 1 meter height and irregular waves up to 0.8 meter height. The experiment's goal is to see how the wind turbine controllers affect the global motions of the FOWT. The experiment's general layout is depicted in the Figure 4. Eight wave probes were used to measure wave heights. Two wave probes, WG1 and WG2, are set 15 and 17 meters from the paddle, respectively, near the sidewall. Wave probes WG3 to WG8 are spaced 15cm apart at a distance of 7.5m from the SPAR. The motions were measured using a MOCAP system based on a Qualisys environment that included multiple cameras, reflecting markers, and QTM software for real-time processing. The processing software takes the 3D positions of various markers as input, smooths them out, and then computes the body's 6 DoF positions. Wave elevation data were collected at a sample rate of 100 HZ, and 6-DoF motions are sampled at the rate of 200 HZ.

The experiment model's components include the floater, mooring system, and the tower. SPAR was chosen as a floater with proportions based on the existing OC3 Hywind SPAR platform (Jonkman, 2010) and *Hywind Scotland* as a reference system. The SPAR is made up of two parts: an aluminum shell that gives the buoy its outward shape, and a ballast, battery, controller, power units, suspended wires, and other components that make up the buoy's core. A vertically truncated cylinder with the taper and transition piece to the tower forms the SPAR floater. The primary dimensions and mass attributes have been summarized in Table 1. The three mooring cables were catenary chains positioned symmetrically 120° apart for the mooring mechanism. The main mooring line was fastened to the bottom of the wave tank and connected to two "bridles," which are part of the delta connection on each cable. The attachment points for moorings were 0.335m below the waterline. Figure 4 depicts the various pieces that make up the mooring model. The underwater unidirectional load cells (LS1, LS2, and LS3) connect the bottom chains to the bridles and are used to measure mooring tension. Studless stainless steel chains were chosen in order to preserve the chain's mass properties and geometry. Model scale mooring properties are shown in Table 1. Because the goal is to understand the system's hydrodynamics, the details of the upper part of the model, which includes aerodynamics components (actuator assembly, nacelle etc.), are not considered in this study. More details related to experiment can be found in Arnal (2020).
<table><caption>Table 1</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>SPAR Mass (including floater, tower and RNA)</td><td>kg</td><td>329.4</td></tr><tr><td>Centre of Gravity with respect to Inertial Frame</td><td>m</td><td>-0.004, 0, -1.535</td></tr><tr><td>Mass Moment of Inertia with respect of CoG FOWT (I<sub>xx</sub>, I<sub>yy</sub>, I<sub>zz</sub>)</td><td>kg/m<sup>2</sup></td><td>490, 490, 10</td></tr><tr><td>Diameter of floater (SPAR)</td><td>m</td><td>0.45</td></tr><tr><td>Diameter at MSL (after tapering)</td><td>m</td><td>0.28</td></tr><tr><td>Taper height</td><td>m</td><td>0.2</td></tr><tr><td>Draft of the SPAR</td><td>m</td><td>2.285</td></tr><tr><td colspan="3">Mooring (Delta or Crow foot mooring)</td></tr><tr><td>Anchors depth</td><td>m</td><td>5</td></tr><tr><td>chain diameter</td><td>mm</td><td>3.7</td></tr><tr><td>Mass of chain in air</td><td>kg/m</td><td>0.275</td></tr><tr><td>Fairlead depth</td><td>m</td><td>-0.335</td></tr><tr><td>Bridle chain length</td><td>m</td><td>1.2</td></tr><tr><td>Bottom chain length (Line 1, Line 2, Line 3)</td><td>m</td><td>15.69, 15.74, 15.72</td></tr><tr><td>Mass of load cell in water (Clumped mass)</td><td>kg</td><td>0.43</td></tr></tbody></table>

Mass and dimensional properties of the SPAR with its mooring arrangement parameters

<figure><img src="image_6.png" /><figcaption>Figure 4: The general layout of the experiment in side and top view with the mooring cable configuration, (Arnal, 2020). Figures are not drawn to scale</figcaption></figure>

## 3.2. Numerical setup

Regular and irregular wave simulations are used to validate the coupling in the domain decomposition. The decay test is also used to verify the natural period of SPAR and to estimate the damping coefficient along with the mooring tension. In the NWT, simulations were run with and without the floating structure. The mesh convergence for wave propagation was investigated using simulations without the structure. The following section describes the computational setup, mesh convergence analysis and finally, the validation of the mooring study is reported to demonstrate the model's capability.

### 3.2.1. Computational domain

In the present coupling, two frames of reference are used, as indicated in the Figure 5, the Inertial frame($R_0$) and the body reference frame ($R_b$). In $R_0$, HOS-NWT, has its origins in the NWT left corner at the still water level (SWL). In $R_b$ foamStar (CFD), which has its origin at the body's center of gravity. Finally, the MoorDyn (Mooring
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Schematics of the computational domain within the HOS domain with Relaxation zones (RZ) are represented in terms of wavelength ($\lambda$). $R_0$ represents the inertial reference frame, and $R_b$ represents the body reference frame. HOS is based on $R_0$, but its origin is different</figcaption></figure>

<table><caption>**Table 2**
foamStar Boundary conditions</caption><thead><tr><td>Boundary</td><td>Velocity</td><td>Pressure</td><td>Alpha</td></tr></thead><tbody><tr><td>Inlet</td><td>waveVelocity</td><td>fixedFluxPressure</td><td>waveAlpha</td></tr><tr><td>Outlet</td><td>waveVelocity</td><td>fixedFluxPressure</td><td>waveAlpha</td></tr><tr><td>Side</td><td>waveVelocity</td><td>fixedFluxPressure</td><td>waveAlpha</td></tr><tr><td>Atmosphere</td><td>pressureInletOutletVelocity</td><td>totalPressure</td><td>inletOutlet</td></tr><tr><td>SPAR-Buoy</td><td>movingWallAndFrameVelocity</td><td>fixedFluxPressure</td><td>zeroGradient</td></tr><tr><td>Bottom</td><td>noslip</td><td>fixedFluxPressure</td><td>zeroGradient</td></tr></tbody></table>

library) has its origins in an inertial frame of reference ($R_0$) at the SWL in the center of the *foamStar* domain. The generation of the regular incident waves is achieved by using the stream function wave theory using an analytical approach, with the far-field boundary condition. The computational domain for the regular wave study is similar to the *foamStar* domain in the Figure 5. The HOS computational domain reproduces the experimental set-up but without the structure. Only through the inlet relaxation zone, the wave generated has been transferred into the *foamStar* domain and in general, all relaxation zones can suppress the waves generated from inside the viscous domain. The *foamStar* domain is solely confined near the SPAR region, which is 6 times the diameter of the SPAR (excluding RZ) and was determined based on the expected motions of the SPAR in 6 DoF from the experiment. The length of the relaxation zones is one wavelength for inlet and outlet and half of the wavelength for side zones decided based on the numerical tests. The numerical domain water depth of 5m has been maintained in HOS and *foamStar* as similar to that of the experiment. This research aims to evaluate the moored floating SPAR's hydrodynamics. Therefore, RNA assembly in the upper tower is not considered in modelling or analysis. However, while considering the total mass characteristics of FOWT, the mass, inertia, and COG of the ignored parts are also taken into account. The tower's freeboard has been kept at 1.2m to appropriately capture the wave run-up during the wave interaction.

### 3.2.2. Boundary conditions

Boundary conditions must be set at the boundaries of the CFD domain in order to guarantee the uniqueness of the governing equations. Table 2 lists the boundary conditions employed in the present investigation and most of them are standard OpenFOAM boundary conditions. The boundary conditions, WaveVelocity and waveAlpha (free
<table><caption>**Table 3**

Numerical schemes</caption><thead><tr><td>Term</td><td>Scheme</td><td>Order</td></tr></thead><tbody><tr><td>Gradient</td><td>Gauss linear</td><td>Second</td></tr><tr><td>Divergence</td><td>Gauss vanLeer</td><td>Second</td></tr><tr><td>Laplacian</td><td>Gauss linear corrected</td><td>Second</td></tr><tr><td>Time</td><td>CrankNicolson</td><td>Second</td></tr></tbody></table>

surface elevation in the form of VOF parameter $\alpha$) sets the values from respective wave models at the inlet of the CFD domain. At other RZ, the target velocity or VOF parameter is blended to zero with the weight function by these boundary conditions. Thus, representing the simplest one-way domain decomposition approach for wave generation. The pressureInletOutletVelocity boundary condition is a blend of pressureInletVelocity and inletOutlet boundary conditions. In the pressureInletVelocity, velocity is computed from the difference between total and static pressure where the direction is normal to the patch faces. The inletOutlet boundary condition is normally the same as zeroGradient, but it switches to fixedValue if the velocity vector next to the boundary aims inside the domain (reverse flow). The movingWallAndFrameVelocity boundary condition was created in-house and is utilized for translation and rotational DoF motions of the body. This boundary condition corrects the flux due to mesh motion in such a way that the total flux through the moving wall is always zero. The fixedFluxPressure boundary condition sets the pressure to ensure that the flux on the boundary is as defined in the velocity boundary condition and the totalPressure boundary condition sets the pressure to zero at the atmospheric boundary. More detailed descriptions of the boundary conditions can be referred in (OpenFOAM).

### 3.2.3. Solver settings

The Courant-Friedrich-Levy number (Co) is often used in CFD solvers for temporal discretization and it indicates that the how much information (U) travels across a computational grid cell during the unit timestep. If the Courant number is more than one, the information propagates through more than one grid cell at each time step, causing the solution to be inaccurate and potentially resulting in nonphysical results. For adequate resolution and capture of the transient flow field in ocean wave modeling, relatively small Co numbers are necessary. The key advantage of the interFoam solver's PIMPLE algorithm is the possibility to use large Co numbers. When smaller Co numbers are necessary to accurately capture minor variations in the fluid, then the benefit becomes less meaningful. In the present study, the maximum Co number of 0.5 was adopted along with 10 outer (SIMPLE) pressure-momentum corrections, 3 inner (PISO) pressure correction loops have been adopted. The adopted residuals limit in the present study was $10^{-8}$ for all simulations. Table 3 summarises the methods of discretization schemes in space and time used in the present study.

### 3.2.4. SPAR body mesh validation: Hydrostatics

The mesh was created with OpenFoam's SnappyHexMesh tool. Near the SPAR, where boundary layer separation occurs, cells were refined into much smaller ones, and refinement boxes were formed at varying distances from the SPAR to minimize a rapid shift from very fine grids to coarser grids. This is critical as the inadequate meshing will lead to loose body volume, causing the change in the draft and subsequently affecting the hydrodynamic characteristics. Thus, the accuracy of the meshed SPAR volume ($V_{OF}$) is estimated and then compared to the physical volume of the SPAR ($V_p$) for the error ($\frac{V_p-V_{OF}}{V_p}$). The SPAR can be refined in several ways to enhance the accuracy, including refining the background mesh (BlockMesh), surface refinement, and tapering edge refinement as shown in the Figure 6. Figure depicts a typical five-step refinement for various background mesh cell configurations (expressed in terms of SPAR diameter). It can be seen that as expected finer refining reduces the error to zero. To further comprehend the refinement aspect, the new parameter PPSD (Points Per SPAR Diameter) is introduced. The fine mesh configuration has been chosen (it has accuracy of more than 99% for even coarser background cells). Based on these results, the number of points per SPAR diameter is increased, and the result is illustrated in the Figure 7b. These PPSD results are also compared in terms of the number of cells that are created, in addition to the error comparison. As demonstrated in Figure 7b, the PPSD of 35 has been found to be cost-effective (in terms of the number of cells) with a volume loss of less than 0.1% and hence it is adopted in the present study.
<figure><img src="image_8.png" /><figcaption>Figure 6: Mesh improvement using different types of refinement over the SPAR body</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Validation of SPAR underwater volume for different mesh refinements</figcaption></figure>

Hydrostatic characteristics are used to validate the new mass (which must be revised to account for the little volume lost while maintaining the draft measured experimentally) and other mass attributes. The SPAR model has been placed in the NWT, allowing it to settle and an equilibrium position has been established in all six degrees of freedom. A slight mass adjustment has been made to bring all motions (especially heave) within $10^{-6}m$, and the investigation was done for 50s, proving the optimal configuration for further investigations.

### 3.2.5. Mesh validation: Wave propagation

The fundamental prerequisite for wave-structure interaction simulations is that the incident waves have to be of accurate quality. This section validates the incident waves for various spatial and temporal discretization combinations. In total, four-wave test cases are investigated for SPAR interaction, as indicated in Table 4. Regular waves are represented by LC 1, while irregular waves are represented by LC 2 (generated based on the Bretschneider spectrum). The chosen cases in Table 4 are based on the number of cases attempted in the experiment, which is likewise investigated with and without wind factors. In the present cases, LC 1.1 represents a severe situation in regular waves, while LC 1.2 represents a normal sea state. For irregular sea states, LC 2.1 and LC 2.2 can be portrayed as normal and severe sea states, respectively.

The mesh and temporal convergence studies for all four test cases for wave propagation in a 2D NWT are discussed in this subsection. This 2D study's computational domain is depicted in Figure 8. The domain is a two-dimensional plane that is an exact reproduction of the three-dimensional domain illustrated in the Figure 5, but without the structure.
<table><thead><tr><th rowspan="2">Parameters</th><th colspan="4">Value</th></tr><tr><th>LC 1.1</th><th>LC 1.2</th><th>LC 2.1</th><th>LC 2.2</th></tr></thead><tbody><tr><td>Wave period (T or T<sub>p</sub>) [s]</td><td>1.9</td><td>1.58</td><td>1.72</td><td>2.37</td></tr><tr><td>Wave height (H or H<sub>s</sub>) [m]</td><td>0.187</td><td>0.132</td><td>0.145</td><td>0.273</td></tr></tbody></table>

**Table 4**

Progressive wave test case: wave parameters with LC 1 represents regular waves and LC2 represents irregular waves

<figure><img src="image_10.png" /><figcaption>**Figure 8:** Sketch of the HOS and CFD domains for wave propagation validation in 2D NWT and the position $x_{HOS} = 17m$ is identical to $x_0 = 0$</figcaption></figure>

<table><thead><tr><th>Mesh level</th><th>Mesh size (Δx, Δz)</th></tr></thead><tbody><tr><td>Very Coarse</td><td>λ/32, H/2</td></tr><tr><td>Coarse</td><td>λ/64, H/5</td></tr><tr><td>Medium</td><td>λ/100, H/8</td></tr><tr><td>Fine</td><td>λ/150, H/10</td></tr><tr><td>Very Fine</td><td>λ/200, H/15</td></tr></tbody></table>

**Table 5**

Progressive waves test case: Test matrix for convergence study

The regular waves (LC 1.1 and LC 1.2) are generated utilizing the *foamStar*'s built-in Stream function wave theory and periodic boundary conditions. HOS-NWT is used for irregular waves (LC 2.1 and LC 2.2). The HOS-NWT exactly reproduces the experimental conditions paddle motion. As previously mentioned, in the $R_0$ reference system, the origin is different for the HOS domain and CFD domain. The HOS domain ($L_{HOS}$) has a full-scale length of 50 meters and a water depth of 5 meters. The full-scale length of the CFD domain is $1\lambda+6D_{SPAR}+1\lambda$, and its origin is positioned from $x_{HOS} = 17m$ at the inertial frame origin. The relaxation zones length is set to $1\lambda$ for both inlet and outlet, respectively. To be noted, for irregular wave generation cases, the HOS domain must entirely incorporate the CFD computational domain into it.

Upper zone, free surface zone, and lower zone are the three zones that make up the computational domain in the vertical direction (Figure 8). Between $z = -H_s$ and $z = H_s$ is the free-surface refinement region. According to (Kim, 2021), the cell aspect ratio ($\Delta z/\Delta x$) in the free-surface refinement region is maintained as 0.5, and the outer zones are kept at ratio four. To understand the solver's convergence, several discretizations are used to investigate the computing grid and time steps. Table 5 shows the cell length $\Delta x$ and height $\Delta z$ that were examined. Based on the wave period ($T$ or $T_p$) considered, different temporal resolutions for wave propagation are performed for each cell
arrangement (Table 5), such as T/100, T/200, T/400, and T/800. The grid and time steps are configured to modify the cell size while keeping the time step constant, and the grid is fixed while the time step is only varied in the second set of investigations. From Mesh32TH2T100 through Mesh200H15T400, eighteen combinations of mesh and temporal resolutions are employed. Mesh32TH2T100 indicates that there are 32 cells per wavelength in the wave propagation direction, 2 cells per wave height in the vertical direction, and 100 time steps per wave period. With all of the above factors mentioned, the outcomes or validation of the regular wave is described initially followed by the irregular wave.

**Regular wave validation:** For the regular wave validation, *foamStar* computational domain was considered with periodic boundary conditions. In this study, the duration of the simulation is kept at 35 wave periods for investigation. Figures 9a and 9b, presents the results of the simulation in the case of LC1.1. Similar observations are noticed for LC 1.2 and hence results are not reproduced here. Figure 9a displays the peak wave amplitude ($\eta_{max}$) for each combination of mesh and time step measured at nearly $20T_p$ and Figure 9b represents the first harmonic amplitude variation for the whole duration of the wave and its damping behavior. In both the figures, the parameters clearly demonstrates convergence when modifying the mesh and time step. Thus, a cell size of 150 cells per wavelength, 10 cells per wave height, and 400 time steps per wave period seem to be computationally efficient and considered for further analysis.

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Progressive waves test case: Harmonic amplitudes of the free surface elevation at the WG5 for different mesh and time configurations</figcaption></figure>

**Irregular wave validation:**

<figure><img src="image_12.png" /><figcaption>**Figure 10:** Cross-Correlation Coefficient for LC 2.1 and LC 2.2 for different mesh and time combination</figcaption></figure>

The irregular waves for various combinations of spatial and temporal discretization (as in regular waves) measured at WG5 were used in the following analysis. The results from the *foamStar* are cross-correlated with the HOS
simulations to investigate the spatial and temporal convergence (Table 5). The suggested technique has the advantage of being able to assess the time shift between signals as well as the similarity between them (Fernández et al., 2014). In the present case, the cross-correlated time series is normalized with auto-correlation coefficient at zero lag such as,

$$ \overline{\eta_{xy}(t)} = \frac{\eta_{xy}(t)}{\sqrt{\eta_{xx}(0)\eta_{xx}(0)}} \qquad (29) $$

where $\eta_{xy}(t)$ represents two real continuous surface elevation of HOS and any of the foamStar wave case. The maximum cross-correlation coefficient of +1 indicates a perfect match between the foamStar and the HOS wave, but a difference in the coefficient other than +1 reveals the quantification of inaccuracy between the signals. Typical cross-correlation results for LC 2.1 and LC 2.2 is shown in Figure 10. Although the experiment and HOS were conducted for a total of 600 sec, to avoid the large computational cost, the comparisons were carried out for 100s (i.e. around $50T_p$). Eighteen different configurations were carried out for both LC 2.1 and LC 2.2, as indicated by the circles and squares in the Figure 10. The lines connect meshes that are of the same type but with varying temporal discretization. The findings demonstrate that even the coarser combinations have an accuracy of approximately 0.97 and that refining the discretizations will bring the accuracy closer to 1. For the wave structure interaction study, the mesh combination of Mesh150H12T400 which is nearly on and above 0.99 in both load cases is proven to be an affordable combination and hence this combination will be used in the present study.

### 3.2.6. Irregular wave Validation: HOS wave Vs Experiment

<figure><img src="image_13.png" /><figcaption>Figure 11: Comparison between experimental and HOS-NWT numerical time series of wave elevation for the two Irregular wave test cases - LC2.1 and LC 2.2</figcaption></figure>

It is important to remember that the HOS-NWT formulation was effectively validated with studies on various cases during the last decade (Ducrozet et al., 2012, 2006; Canard et al., 2020; Hasan et al., 2019). LC 2.1 and LC 2.2 imply a wide variety of wavemaker motions with large nonlinearities in order to validate our current concept and implementation. The experimental paddle motion has been used in the HOS-NWT to simulate the wave characteristics and compared with the experimental measurements. For this simulation, 325 modes in the flow direction ($N_x$) and 33 modes in the vertical direction ($N_z$) have been used. Throughout the simulation, the HOS order is kept at 3. Figures 11a and Figure 11b shows the comparison of the HOS-NWT wave generation and the experiment at the WG2 probe location, respectively for both the LC. In both figures, we can see that the agreement between the experiment and numerical models is identical. After 50 seconds, however, slight deviations in the peak magnitude of the crest or trough were noted. This error is clearly connected with wave amplitude, indicating high nonlinearity but the cross-correlation coefficient for these simulations is around 0.994, which makes the discrepancy quite small, hence disregarded. It should also be noted that these deviations in the incident wave will have a minor impact on the motion responses in wave-floating body simulations.

### 3.2.7. MoorDyn Validation

MoorDyn is a dynamic mooring solver based on lumped mass approach. The version of MoorDyn used in this paper was v1.01.03C. The solution for the mooring solver is based on a three-step process. First, it generates a Quasi
Static model based on the input conditions. Second, to allow the system to achieve the initial equilibrium, dynamic
relaxation will be carried out, which includes dynamic drag and acceleration forces. Finally, tension at each fairlead as
well as the motion (position, velocity, and acceleration) of each node at each time step will be evaluated.

In the present study, the standalone version of the MoorDyn model was first validated for a static problem, and
compared with other quasi-static models such as MoorPy (Hall et al., 2021), MAP (Masciola, 2018), and AQWA. The
single delta line is the subject of the first validation, as it is one part of the mooring configuration used in the SPAR
interaction study. Masciola et al. (2013) can be used to look up the nomenclature of the single delta line considered. This
line's properties and dimensions are taken from the MAP examples section. This delta line introduces the connecting
node between the lines, the position of which will be iterated by the solvers, and the accuracy of the pretension is
dependent on the final node's position. The static pretension comparison between the solvers for the case is as shown
in Figure 12a. Since Lines 2 and 3 are symmetric, the figure shows a single bar representing both lines. For Line 2/3, the
difference between the models is quite small, but for Line 1, the difference is nearly 2% of the pretension between the
models. It is realized that the original iteration of the connection node position displaces into a new position, resulting
in the Line tension difference. The dynamics assumed by MoorDyn could also be a concern. However, because the
error is so minimal, the MoorDyn model is agreed upon for its robustness.

<figure><img src="image_14.png" /><figcaption>**Figure 12:** Comparison of line tensions along single delta line (left) and complete delta type mooring(right) for different mooring models</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13:** Typical delta mooring model adopted in the present study with its top view(left) and isometric view (right)</figcaption></figure>

For the final validation, the mooring setup used in the present study, the crowfoot model having 3 main lines
connecting 6 bridle lines with the load sensor (LS1 to LS3) as shown in the Figure 4 are tested. This configuration is
validated in a standalone MoorDyn model for a total of 9 mooring lines and 9 nodal locations, and the typical setup is
as illustrated in Figure 13. It is important to validate input characteristics such as connecting nodes, unstretched length,
mass, and stiffness properties from the experiment. The primary mooring properties from the experiment are listed in
the Table 1. In order to obtain proper tension, in the numerical model symmetry adjustment of the anchor and fairlead
positions (120° to each main line) are considered, because of the fact that achieving symmetry while doing physical
tests was difficult. The axial stiffness (AE) and equivalent diameter are obtained from ORCINA recommendations
<table><caption>**Table 6**

Parameters of the mooring lines in prototype scale</caption><thead><tr><th>Parameters<br>(Scale 40)</th><th>Units</th><th>Prototype Values</th></tr></thead><tbody><tr><td>Diameter of chain</td><td>[m]</td><td>0.266</td></tr><tr><td>Mass density</td><td>[kg/m]</td><td>440</td></tr><tr><td>EA</td><td>[N]</td><td>6E9</td></tr><tr><td>Added mass coefficient</td><td></td><td>1.0,0.0</td></tr><tr><td>Drag coefficient</td><td></td><td>1.8, 0.25</td></tr><tr><td>Clumped mass (sensor)</td><td>[kg]</td><td>2.75e4</td></tr><tr><td>UnstrLen (Major and Bridle)</td><td>[m]</td><td>638.2, 48</td></tr></tbody></table>

(ORCINA). In the transverse direction (Eqn. (25)), quadratic damping coefficients are set to 1.8, and in the axial direction (Eqn. (26)), they are set to 0.25 based on the numerical analysis. The added mass coefficient in the transverse direction are set to the value of 1.0. The convergence study for a number of segments in the mooring line was studied, and it was determined that from the 50 segments in the main line and 4 segments in the bridle were the most efficient and converged solution begins, anyhow with a factor of safety of 1.5, 75-6 combination was chosen.

Although MoorDyn is simple in concept, its practical implementation can be difficult, as it is sensitive to the input parameters chosen, causing numerical instability (Paduano et al., 2020). Similar instability is observed in the present study when working with the mooring properties on experimental model values. However, using the prototype values shows that numerical instability can be treated. Hence the mooring properties used in the experiment are modified to prototype values (1:40 scaled model of FOWT) using Froude and Cauchy scaling, and the updated mooring properties are shown in Table 6. Similar setups were also built in various mooring models, and pretension at the static conditions for major lines Line 1 to 3 (Figure 4) are compared as shown in Figure 12b. The pretension aligns well with MAP, and only minor variations were noticed compared to MoorPy and AQWA proving the efficiency of the MoorDyn and input conditions considered. The final system of coupling considered is as shown in Figure 14, wherein the mooring analysis will be carried out in prototype conditions for motions and velocities and the computed forces are scaled-down and coupled with the CFD.

<figure><img src="image_16.png" /><figcaption>**Figure 14:** Flow chart of the short algorithm between the *foamStar* (in model scale) and *MoorDyn* library (in prototype scale).</figcaption></figure>
Finally, the MoorDyn results are validated from the experimental pull-out test for the estimation of stiffness. To characterize the horizontal mooring stiffness for the platform surge DoF, various pull-out tests were done in the experiment. One such test comparison with experiments is shown in Figure 15, wherein the mooring stiffness is nearly $75Nm^{-1}$ based on the slope of the force displacement relationship. The mooring restoring force can be quite similar to a linear force, which is interesting for a catenary system since the total restoring force in DeepCWind campaigns and the related OC4 numerical model (Robertson et al., 2014) tends to be nonlinear at larger offset values. It is important to note that the restoring force stated here refers to the overall restoring of all 9 lines, not just the catenary line.

<figure><img src="image_17.png" /><figcaption>Figure 15: Comparison of the results of the Pullout test in the MoorDyn solver with the findings of the experiment</figcaption></figure>

### 3.2.8. Stiffness Matrix Coefficients

The stiffness coefficient as determined by the pull-out test for surge ($K_{11}$) mooring lines is $75Nm^{-1}$. Similarly, the sway ($K_{22}$) and yaw coefficients ($K_{66}$) were calculated using similar conditions and also modeling the mooring lines in ORCAFLEX (ORCINA) for validation. The values are $80Nm^{-1}$ and $127Nm rad^{-1}$, respectively for surge and yaw. Non-diagonal stiffness in terms of restoring stiffness in heave, roll, or pitch from the mooring appeared to be extremely minimal in this investigation. In the case of solving wave structure interaction with stiffness matrix, the mooring line mass of 11.4 kg ($\pm$1 kg) is added to the SPAR's total mass, which is then fine-tuned to maintain the same draft.

# 4. Results and Discussion

This section discusses the validation of the developed numerical approach with experimental measurements. The first moored decay test is discussed, followed by the interactions of regular and irregular waves with SPAR for mild and severe sea state conditions.

## 4.1. Moored Decay

In the coupled *foamStar*-MoorDyn solver, two independent DoF decay tests are carried out with delta mooring lines. The SPAR is perturbed along the DoF of interest, then released to move freely from that point. Surge and heave are the decay tests performed. The mooring tensions and their accompanying modeling results are discussed.

### 4.1.1. Surge decay

The initial displacement for the Surge in *foamStar* is accomplished using the `moveDynamicMesh` utility, which only conducts the dynamic mesh motion without solving the flow field while monitoring the mesh quality aspects. The SPAR floater with delta mooring lines is moved 0.11m and then released to move freely from its initial position. The COG is set to a displaced position, and the mooring lines are given initial displacement and zero initial velocity and are free to move with the body motions. The surge had a strong influence on the pitch and a weak influence on the heave, making it difficult to estimate the proper surge natural period. Hence only for this testing, only Surge DoF taken into account. The findings of the motion shown in comparison to the experiment (Figure 16a and Figure 16b), as well as the mooring tension in Line 1(ref Figure 4 for Line nomenclature), appear to be quite good. All of the line tensions
<figure><img src="image_18.png" /><figcaption>**Figure 16:** Comparison between experimental and numerical time series of surge decay test with its corresponding DoF and mooring tension</figcaption></figure>

were in good agreement, but *Line1* was the most critical (as it is in every test scenario), therefore it is only shown. When the SPAR is shifted towards *Line1* for initial displacement, there is a loss of up to 8N in pretension at first, and the variation approaches ±10N. Surge's natural period is found to be 18.82s, while the experiment's natural period is 18.6 – 18.7s. The overall damping ratio is calculated using the logarithmic decrement method, and it is found that the percentage changes between 2.2% to 2.3% for the experiment, but from the numerical simulation it is between 2.35% to 2.45%. In the instance of the surge decay test, the overall difference is found to be less than 1%, indicating that the coupling model performs well for both motion and tension estimation in the surge decay problems.

### 4.1.2. Heave decay

<figure><img src="image_19.png" /><figcaption>**Figure 17:** Comparison between experimental and numerical time series of heave decay test with its corresponding DoF and mooring tension</figcaption></figure>

The heave decay follows a similar procedure as that of the surge simulation, with the SPAR being moved in the heave direction for an initial displacement (in this case 0.045m) and then released in calm water. The decay simulation in heave is carried out with all 6 DoF motions active since the heave motion has only slight influences from all other motions. Also note that the mooring in the form of stiffness matrix was examined as well (although it did not have any stiffness in the heave direction) for the present heave decay. The motion recorded in both mooring scenarios is shown in the Figure 17a. The figure has a few noteworthy observations as follows. The mooring in the form of a stiffness matrix component can be thought of as having no resistance in the heave direction ($K_{33}$). As a result, the motion characteristics only depict the body motions for a given mass and displacement. In this situation, the natural period appears to be 4.91s, and the damping ratio appears to be 1.5%. Incorporating the mooring into the problem changes the damping of the problem by adding 0.5 – 0.6% to total damping, resulting in a natural period of 4.8s. The natural period in the experiments was determined to be 4.8 – 4.9s, and the damping was found to be similar to that of the moored simulation. Figure 17b compares the numerical tension of *Line 1* (at LS1) with the experiment, to illustrate
the effect of mooring tension variation. The vertical axis is kept the same as in the surge free decay tension figure to show the reader that the influence in the heave is too small when compared to the surge but cannot be ignored.

## 4.2. Wave Interaction with SPAR

<figure><img src="image_20.png" /><figcaption>**Figure 18:** Typical representation of computational domain in the wave structure interaction study with mooring</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19:** Typical representation of different time instants of the interaction of regular waves with the moored floating SPAR using *foamStar*-MoorDyn coupling for case LC 1.1</figcaption></figure>

In this section, the performance of the *foamStar* and MoorDyn coupling in reproducing the experiments on wave interactions with SPAR are reported. The wave interaction's computational domain and relaxation zones are the same as those described in Section 3. A common Mesh150H12T400 is chosen in order to make a similar discretization for
<figure><img src="image_22.png" /><figcaption>**Figure 20:** Comparison between experimental and numerical time series of wave elevation at WG4 of the Floating SPAR for the two regular wave test cases - LC1.1 and LC 1.2</figcaption></figure>

regular (Mesh150H10T400) and irregular waves studies (Mesh150H12T400), and it will be represented in the results with the prefix 'fs' denoting *foamStar*. The Figure 18 depicts a typical computational domain with mooring lines. The fine refinement in the free surface and near the SPAR can be seen in this figure, along with the typical contour of the free surface. The SPAR body is fine-tuned based on the previous estimates in order to maintain PPSD 35 and a little mass modification to ensure the same drafts. Based on the wave being considered, the wavelength will vary for each load case, influencing the number of cells. In the present study of LC 1.1 to LC 2.2, it varies from 1 to 2.5 million cells, and the code was executed concurrently using MPI on the LIGER supercomputer at Ecole Centrale Nantes. The duration of the simulation as a whole varies from case to case. Only 25T is kept as the simulation duration for regular wave cases (LC 1) because the steady state is reached within 4 to 5 peaks of the simulation, but 50Tₚ has been kept as the simulation duration for irregular wave cases (LC 2). One CPU per 30000 cells is chosen for parallel processing and computational time varies from 10 hours in LC 1 to 36 hours in LC 2. Each LC involves two methods of coupling to handle the mooring line, one in the form of a stiffness matrix and the other with dynamic mooring (MoorDyn). As previously indicated, two regular wave test cases (LC 1.1 and LC 1.2) and two irregular wave test cases (LC 2.1 and LC 2.2) are being carried out and discussed in the following sections.

### 4.2.1. Regular wave interaction

The hydrodynamic responses of the SPAR under regular waves are investigated in this section. All the waves propagate in a positive surge direction. The waves with moderate steepness ($H/\lambda = 0.03$) are simulated and the wave characteristics areas listed in Table 4. At various time instants, a typical interaction between the wave, structure, and mooring is demonstrated in the Figure 19. The wave LC 1.1 ($T = 1.9s$) has an initial ramp of 10s, hence the steady-state time instant of 17.5s to 19s was chosen for the presenting outcome. The magnitude of velocity in the x-direction as the wave progresses is shown in red and blue on the graph. For clarification, the free surface contour is also displayed ($\alpha = 0.5$). The mooring is represented by segment tension, which varies from 18N to 45N and is colored blue to red. The crest of the wave approaches before SPAR (already has its own motion during transient and few wave periods before) at 9.25T. Line 1 has the highest mooring tension and Lines 2 and 3 have slightly lower tension. As the wave progresses, the crest reaches the SPAR (9.5T) with an impact, (9.75T) then the SPAR moves along with the wave orbital, reducing the tension in Line 1 to 30-35N and causing significant tension in Line 2 and Line 3. Line 3 is on the opposite side of the figure, making it difficult to visualize in color. Because Line 3 is symmetric, it can be said that the behavior will be similar to Line 2. When the trough reaches the SPAR at the final instant (10T), it gently returns to its previous motion, which is similar to 9.25T, lessening the tension in Lines 2 and 3 and achieving equilibrium before the next cycle begins.

The comparison between experiments and numerical simulation for wave probe WG4 is shown in Figure 20. Out of all the probes, the probe WG4 which is at the same location along the SPAR is shown for the comparison. The initial ramping in the experiment is different from the numerical ramping leading to a small difference up to 10s. After achieving a steady-state, the wave appears to achieve suitable amplitudes and is in phase with the experiment. For both the LC, the total duration of simulation recorded in the experiment is 30s. Both the LC 1.1 and 1.2 were initially evaluated for 30s, however, to quantify the surge effects (Surge natural period roughly 19s) the LC 1.1 simulations were extended for additional 15 seconds. The numerical and experimental time series for surge, heave, and pitch angles (LC
Numerical coupling strategy using HOS-OpenFOAM-MoorDyn

<figure><img src="image_23.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(g)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(h)</figcaption></figure>

**Figure 21:** Comparison between experimental and numerical time series of the surge, heave, and pitch of the floating SPAR and wave excitation forces for the two regular wave test cases - LC1.1 and LC 1.2
<figure><img src="image_31.png" /><figcaption>**Figure 22:** Comparison between experimental and numerical time series of tensions of the front, line 1(LS1), and rear mooring lines, line 2 and line 3 (LS2 and LS3), for the two regular wave test cases - LC1.1 and LC 1.2</figcaption></figure>

1.1 and LC 1.2) are shown in the Figure 21. To demonstrate the influence of LC 1.1 over LC 1.2, the vertical axis limit is kept constant, indicating that LC 1.1 has larger amplitudes than LC 1.2. Sway, roll, and yaw motions were not shown since they were negligible, both in the numerical simulations as well as in the experiments. The surge motions shown in Figures 21a and 21b, appears to oscillate with increased amplitudes at the ramping time itself before the steady-state due to the difference in the ramping. The rate of change of oscillation amplitude compares well to the surge experimental period. In the case of LC 1.2, the surge motion follows a similar trend as LC 1.1, albeit with a much smaller amplitude. The behavior of the stiffness matrix component and MoorDyn is virtually identical, although MoorDyn has a very slightly higher amplitude. This is due to the fact that MoorDyn includes coupling of all DoF, whereas, in the stiffness matrix approach, we have considered only three DoFs. The heave and pitch motions of the LCs are depicted in the Figures 21c and 21d, and Figures 21e and 21f. The motions were well simulated by both the stiffness matrix approach
<figure><img src="image_32.png" /><figcaption>**Figure 23:** Comparison between experimental and numerical time series of Wave elevation WG4 of the Floating SPAR for the two Irregular wave test cases - LC2.1 and LC 2.2</figcaption></figure>

and MoorDyn and both the motion amplitudes are quite small, with heave amplitudes of nearly 2cm and 1cm and pitch amplitudes of 1.25º and 0.8º. From the results, it proves that the accurate stiffness matrix components are sufficient to handle the problem for this small amplitude of motions (despite LC 1.1 having higher amplitude regular waves). Thus, dynamic mooring analysis may not be required in such cases. However, MoorDyn has an advantage if the requirement is for the estimation of restoring force or the mooring tension assessment. Figure 21g and Figure 21h represent the wave excitation forces over the floating SPAR for the load cases LC 1.1 and LC 1.2 respectively. As stated in section (2.3), the total force acting on the SPAR is influenced by a number of factors like mass, buoyancy, gravity, mooring, fluids (including pressure and viscosity), and so on. Using fluid forces, we can derive solely the wave excitation forces operating over the SPAR by taking into account the normal pressure acting over the SPAR and its tangential viscous component.

Following (Huseby and Grue, 2000), the wave excitation forces (removing the hydrostatic forces), are normalized by $ρgS^{\frac{H}{2}}$, where S represents the cross-sectional area of the SPAR. Because the motions closely match the experiment, the forces calculated by the solver should be fairly accurate, and also the numerical cases (stiffness matrix form and MoorDyn form) have similar motions, only results from one form (MoorDyn) is presented. The obtained force trend is proportional to each case's wave elevation, and the x component of the force ($F_x$) amplitude appears to dominate the vertical force ($F_z$) by closely 1.5 times. Otherwise, the normalized amplitude has only a 5% - 7% difference between LC 1.1 and LC 1.2.

Figure 22 shows the line tensions from both LC. It should be noted that the experiment uses three load sensors at the LS position, therefore only Lines 1, 2, and 3 are considered, and since the bridle lines had no tension measurements, they were ignored. In the MoorDyn, the top most segment tensions in the major lines were used for comparison, which was near to LS in the experiment. The experimental data (sampled at 200 Hz) was subjected to an average filter, whereas numerical results are saved at a sampling rate (as determined by the CFD solver) with no filter. The figure shows good comparisons and it's worth noting that Line 1 is at the receiving end of the wave, at the back of the SPAR. As a result, the tension oscillations are up to four times larger than the oscillations of Lines 2 and 3 at the SPAR's front junction. Another point to note is that LS2 and LS3 are symmetric, therefore the tension and variation are identical in both the lines. Despite the input conditions, LC 1.2 tension amplitudes are small when compared to LC 1.1, it demonstrates that *foamStar*-MoorDyn coupling could solve for both the smaller and higher amplitude regular wave interactions with the floating body. Similar observations for motions and tensions were also seen with the quantitative analysis as presented in Section 4.2.3. Hence concluding that the regular wave generated inside the *foamStar* by the CN-Stream function wave model has been evaluated for floating SPAR interaction with mooring models, and it appears that the coupling is effective in accurately predicting the motions and forces. Despite the accuracy, validation of coupling with HOS-NWT is the next requirement as this way multiple types of waves can be generated, making the model more robust.

### 4.2.2. Irregular Wave structure interaction

Now the SPAR response was investigated for the uni-directional irregular sea state generated based on the Bretschneider spectrum (Table 4) in HOS-NWT. Two different irregular sea state was analyzed, Case LC 2.1 is a
<figure><img src="image_33.png" /><figcaption>**Figure 24:** Comparison between experimental and numerical time series of the surge, heave and pitch of the floating SPAR and wave excitation forces for the two irregular wave test cases - LC2.1 and LC 2.2</figcaption></figure>
relatively mild sea state, and Case LC 2.2 corresponds to severe sea state conditions. HOS-NWT is used to simulate the incidence wave field in the whole wave tank based on the motion of the wavemaker from the experiment. The simulation has been carried out for 500-600s in the experiment and HOS, however in the CFD, it is arranged to keep 50 $T_p$ to restrict the computational cost. The computational domain is comparable to that of the domain used in regular wave interaction study, with the exception that $\lambda_p$ instead of $\lambda$ has been adopted. Based on the peak wavelength and significant wave height ($H_s$), the domain zones were updated to maintain 150 cells per wavelength and 12 cells per wave height, 400 timesteps per peak wave period (fsM150H12T400). The time history comparison between the *foamStar* and experimental measurements for WG4 is shown in Figure 23. Unlike regular waves, here the experimental wave paddle motion has been used, so there is no ramping issue in both the test cases. However, it should be noted that the HOS has few discrepancies in comparison with the experimentally measured wave for some amplitudes and that these discrepancies are also seen in the *foamStar* produced waves. Also, as noted in the parametric study, discrepancies created owing to mesh and time discretization in *foamStar* could also be added to the total discrepancy. However, the difference is very small in the case of a small amplitude case (LC 2.1) and a little higher for a higher amplitude case (LC 2.2). It should also be emphasized that the waves recorded are with the SPAR in place and hence the influence of scattered waves from the SPAR will also contribute to the wave elevation recorded.

The SPAR response in the surge, heave, and pitch for the cases of LC 2.1 and LC 2.2 are shown in the in Figure 24 and mooring tensions from the corresponding cases are shown in Figure 25. The influence of the two mooring approaches is discussed similar to regular wave interaction with SPAR. In the case of LC 2.1, there were very few to no discrepancies between the motions predicted in surge, heave, and pitch with the experiment. Figure 24a and Figure 24b represent surge motion recorded for LC 2.1 and LC 2.2. For LC 2.1, fsM150H12[MoorDyn] has predicted comparable motions using the mooring as similar to the stiffness matrix. However, for the case of LC 2.2., the comparisons were good at the beginning, but when the wave amplitudes are increased, the stiffness matrix approach leads to more deviation modifying the SPAR response. However, the fsM150H12[MoorDyn] prediction compared well with the experimental surge motion. For the heave (Figure 24c, Figure 24d) and pitch (Figure 24e, Figure 24f) motions also, mild case (LC 2.1) compared well with the experiment but for the severe case (LC 2.2) stiffness matrix approach has certain discrepancies. The differences are also identical to the wave amplitude variations reported. The MoorDyn overpredicts the pitch by a few decimals in the amplitudes and this little over prediction may be seen throughout the LC, although it is of modest size. The cause for this might be that the mooring configuration in the experiment is not perfectly symmetric, resulting in small unsettled moments in the pitch direction. Also, with all of the nacelle RNA supports, it is possible that the $I_{yy}$ provided by the experiment is not quite accurate. Overall, the findings showed that the coupled model's motion prediction is precise for both the sea state conditions. The stiffness matrix approach prediction is very good for mild sea states, but for severe sea states, a dynamic mooring model is required for accurate motion prediction similar to that of the experiment. A similar conclusion is quantitatively obtained in the Section 4.2.3. The Figure 24g and Figure 24h represents the normalized hydrodynamic force components obtained over the SPAR for the cases of LC 2.1 and LC 2.2, respectively. As noted in the regular wave cases, the force trend corresponds very well to the wave elevation, and the amplitude ranges are confined between $\pm0.5$ for both the load cases. Overall, the range of forces with respect to wave height demonstrates that regular waves develop nearly 1.5 times higher force range than the irregular waves, whereas irregular waves display peak amplitudes at certain locations exceeding 0.9-1, indicating the strength of its impact.

The abrupt change of tension (Figure 25d to Figure 25f) happening in very short interval of time for the LC 2.2 required robust mooring model to capture tensions accurately. In the case of LC 2.1, the tension forces match quite well in all three lines. In large surge drift offset circumstances, the most substantial tension is evident in the seaward cable. The force amplitudes of the leeward and seaward mooring lines are identical in amplitude when there is little or no drift. For LC 2.2, the predicted mooring tension appears to be in phase with the experiment, however, the amplitudes are slightly higher than expected, the reason has been addressed earlier. The amplitude variation of tension between 15 to $25N$ in a short time frame of less than a second (LC 2.2) is 3 to 4 times higher than any variations reported earlier (particularly LC 2.1). The coupled solver captures all of these tension variations with little to no inconsistencies, showing the efficiency of the coupling. In comparison to all the cases reported, the LS2 and LS3 lines experience more impact tension variations only in LC 2.2 which is also accurately captured by the coupled solver.

### 4.2.3. Quantitative analysis

The qualitative comparison presented in the previous sections is quantitatively compared in this section using a cross-correlation analysis. The coefficient is computed with Eqn. 29, but the motions and the mooring tensions of the
<figure><img src="image_34.png" /><figcaption>**Figure 25:** Comparison between experimental and numerical time series of tensions of the front, line 1 (LS1), and rear mooring lines, line 2 and line 3 (LS2 and LS3), for the two irregular wave test cases - LC2.1 and LC 2.2</figcaption></figure>

numerical model are compared with the experiment instead of the surface elevation. Recalling from Section 3.2.5, the maximum cross-correlation coefficient of +1 indicates a perfect match between the *foamStar* and the experimental results, and the distance to 1 quantifies the difference between the signals. The surge motions in the regular wave were sensitive to the ramp in the wave generation, and the Figures 26a and 26b illustrate this discrepancy between the signals and the experiment for both the stiffness matrix and dynamic mooring cases. In the stiffness matrix cases, the coefficient are 0.89 and 0.93 for LC 1.1 and 1.2, respectively. However, a minor improvement is seen with *MoorDyn* (coefficient of 0.94 in both cases), even though the ramp effect is still present. The stiffness matrix approach achieved a good comparison of 0.98 for the mild sea condition LC 2.1 (Figure 26c) but lost accuracy to 0.86 for LC 2.2 (Figure 26d).
<figure><img src="image_35.png" /><figcaption>(a) Case LC 1.1</figcaption></figure>

Relative Comparison - Case LC 2.1

<figure><img src="image_36.png" /><figcaption>(c) Case LC 2.1</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(b) Case LC 1.2</figcaption></figure>

Relative Comparison - Case LC 2.2

<figure><img src="image_38.png" /><figcaption>(d) Case LC 2.2</figcaption></figure>

**Figure 26:** Relative analysis using Cross Correlation between foamStar cases and Experiment where [K] represents stiffness matrix approach, MoorDyn represents dynamic mooring coupled approach

With values higher than 0.96 in both LC 2.1 and LC 2.2, MoorDyn appeared to perform well, supporting the validity of the discussions from earlier sections. In LC 1.1, LC 1.2, and LC 2.1, heave and pitch motions are well reproduced by the numerical model, with coefficient ranges of 0.99 and higher by both the approaches. While MoorDyn appears to perform well in LC 2.2, reaching a coefficient up to 0.96 and above, the stiffness matrix cases are less accurate, and the coefficient stays at 0.9 to 0.94 (Figure 26d). The mooring tensions appear to match well in all the cases with the coefficients in the range of 0.99 to 0.999 (Figure 26), proving the efficiency of the mooring model and its coupling.

**5. Conclusion**

The paper presents the development of a coupled numerical model that can simulate the station keeping of real-time floating offshore structures in both normal and severe sea state conditions. The efficiency and safety of the floating offshore structures have direct implications from the structure's elasticity, structural dynamics, hydrodynamics and aerodynamics interaction. Hydrodynamics is among the most critical fields since it affects all other disciplines directly or indirectly. In the present study, the hydrodynamic interaction of fully nonlinear regular and irregular waves over a floating moored SPAR was simulated using an OpenFOAM-based *foamStar* solver with the domain decomposition technique. The CN-Stream wave model was used for regular wave generation, while for irregular wave generation, the HOS-NWT was used. The floating structure's station keeping was characterized by a static (as a stiffness matrix) and dynamic mooring model (two-way coupling with MoorDyn). The numerical model results were compared with the experiment carried by (Arnal, 2020) in the ECN wave basin. The physical model consists of the SPAR structure moored to the tank bottom using delta type mooring with nine mooring lines.

Meshing a floating SPAR was a critical task since poor meshing would modify the volume of the body, altering the draft and pushing the motions to be irrelevant. Hence, a new parameter, PPSD (points per SPAR diameter) is introduced, and a value of 35 is recommended for an accurate meshing of the floating structure to avoid volume loss. Different types of mooring setups were investigated to see the efficiency of the dynamic mooring model MoorDyn, and the results were compared with other mooring models (MAP, MoorPy, AQWA) along with experimental data. Since working in the model scale in the MoorDyn caused numerical instability in the delta type mooring case, the MoorDyn
solver was made to operate in the prototype scale. Thus, forces from the mooring line is down-scaled and implemented in the *foamStar* solver. Surge and heave decay tests were numerically simulated, and it's interesting to find that in the heave free decay, mooring driven damping accounts for nearly 25% of overall damping. In the case of surge decay, the natural period, damping, and mooring tension from the coupled model provided acceptable results in the experiments. Further, good agreement between the numerical and experimental motions (surge, pitch, and heave) was observed for the regular and irregular wave interaction over the SPAR. The tensions in the mooring lines also compared well with the measurements. Mooring in the form of a stiffness matrix provides an accurate comparison for wave cases with mild and moderate amplitudes. However, to address mild, moderate, and even severe sea state conditions with both regular and irregular waves, dynamic mooring simulation (MoorDyn) is proven to be vital.

Hence, the developed coupled modeling concept can be used to simulate the more complex wave-floating structure interactions. These capabilities will aid in the design and analysis of the floating structures and the mooring to efficiently survive in extreme wave conditions. Apart from hydrodynamics with rigid body motions, hydroelasticity of the structure, aerodynamics for wind will be incorporated in the developed solver to extend its capability for practical aero-hydro-elastic applications

# 6. Acknowledgment

This work has been performed in the framework of the Chaire Hydrodynamique et Structure Marines CENTRALE NANTES - BUREAU VERITAS. The first author acknowledges IITM and LHEEA, ECN for the financial support for this Ph.D. study.

# References

Agarwal, S., Saincher, S., Sriram, V., Yan, S., Xie, Z., Schlurmann, T., Ma, Q., Yang, X., Wan, D., Gong, J., et al., 2021a. A comparative study on the nonlinear interaction between a focusing wave and cylinder using state-of-the-art solvers: part b. International Journal of Offshore and Polar Engineering 31, 11–18.

Agarwal, S., Sriram, V., Yan, S., Murali, K., 2021b. Improvements in mlpg formulation for 3d wave interaction with fixed structures. Computers & Fluids 218, 104826.

Aliyar, S., Meyer, J., Sriram, V., Hildebrandt, A., 2021. Experimental investigation of offshore crane load during installation of a wind turbine jacket substructure in regular waves. Ocean Engineering 241, 109979.

Arnal, V., 2020. Experimental modelling of a floating wind turbine using a «software-in-the-loop» approach. Ph.D. thesis. École Centrale de Nantes.

Benitz, M.A., Schmidt, D.P., Lackner, M.A., Stewart, G.M., Jonkman, J., Robertson, A., 2015. Validation of hydrodynamic load models using cfd for the oc4-deepwind semisubmersible, in: International Conference on Offshore Mechanics and Arctic Engineering, American Society of Mechanical Engineers. p. V009T09A037.

Biausser, B., Fraunié, P., Grilli, S.T., Mercer, R., 2004. Numerical analysis of the internal kinematics and dynamics of 3-d breaking waves on slopes. International Journal of Offshore and Polar Engineering 14.

Black, J.L., Mei, C.C., Bray, M.C.G., 1971. Radiation and scattering of water waves by rigid bodies. Journal of Fluid Mechanics 46, 151–164.

Canard, M., Ducrozet, G., Bouscasse, B., 2020. Generation of 3-hr long-crested waves of extreme sea states with hos-nwt solver, in: International Conference on Offshore Mechanics and Arctic Engineering, American Society of Mechanical Engineers. p. V06BT06A064.

Choi, Y., 2019. Two-way coupling between potential and viscous flows for a marine application. Ph.D. thesis. École centrale de Nantes.

Choi, Y., Gouin, M., Ducrozet, G., Bouscasse, B., Ferrant, P., 2017. Grid2grid: Hos wrapper program for cfd solvers. arXiv preprint arXiv:1801.00026.

Choi, Y.M., Kim, Y.J., Bouscasse, B., Seng, S., Gentaz, L., Ferrant, P., 2020. Performance of different techniques of generation and absorption of free-surface waves in computational fluid dynamics. Ocean Engineering 214, 107575.

Devolder, B., Stratigaki, V., Troch, P., Rauwoens, P., 2018. Cfd simulations of floating point absorber wave energy converter arrays subjected to regular waves. Energies 11, 641.

Dewey, R.K., 1999. Mooring design & dynamics, a matlab® package for designing and analyzing oceanographic moorings. Marine Models 1, 103–157.

Di Paolo, B., Lara, J.L., Barajas, G., Losada, I.J., 2020. Wave and structure interaction using multi-domain couplings for Navier-Stokes solvers in OpenFOAM®. Part I: Implementation and validation. Coastal Engineering, 103799doi:10.1016/j.coastaleng.2020.103799.

Dommermuth, D.G., Yue, D.K., Lin, W., Rapp, R., Chan, E., Melville, W., 1988. Deep-water plunging breakers: a comparison between potential theory and experiments. Journal of Fluid Mechanics 189, 423–442.

Ducrozet, G., Bonnefoy, F., Le Touzé, D., Ferrant, P., 2006. Implementation and validation of nonlinear wavemaker models in a hos numerical wave tank. International Journal of Offshore and Polar Engineering 16.

Ducrozet, G., Bonnefoy, F., Le Touzé, D., Ferrant, P., 2012. A modified high-order spectral method for wavemaker modeling in a numerical wave tank. European Journal of Mechanics-B/Fluids 34, 19–34.

Ducrozet, G., Bouscasse, B., Gouin, M., Ferrant, P., Bonnefoy, F., 2019. Cn-stream: Open-source library for nonlinear regular waves using stream function theory. arXiv preprint arXiv:1901.10577.

Ducrozet, G., Fink, M., Chabchoub, A., 2016. Time-reversal of nonlinear waves: Applicability and limitations. Physical Review Fluids 1, 054302.
Dunbar, A.J., Craven, B.A., Paterson, E.G., 2015. Development and validation of a tightly coupled cfd/6-dof solver for simulating floating offshore wind turbine platforms. Ocean Engineering 110, 98–105.

Elhanafi, A., Macfarlane, G., Fleming, A., Leong, Z., 2017. Experimental and numerical investigations on the hydrodynamic performance of a floating-moored oscillating water column wave energy converter. Applied energy 205, 369–390.

Engsig-Karup, A., 2006. Unstructured Nodal DG-FEM solution of high-order Boussinesq-type equations. Ph.D. thesis.

Fernández, H., Sriram, V., Schimmels, S., Oumeraci, H., 2014. Extreme wave generation using self correcting method—revisited. Coastal Engineering 93, 15–31.

Ferri, F., Palm, J., 2015. Implementation of a dynamic mooring solver (moody) into a wave to wire model of a simple wec. Department of Civil Engineering, Aalborg University.

Gu, H., Stansby, P., Stallard, T., Moreno, E.C., 2018. Drag, added mass and radiation damping of oscillating vertical cylindrical bodies in heave and surge in still water. Journal of Fluids and Structures 82, 343–356.

Guilcher, P., Ducorzet, G., Alessandrini, B., Ferrant, P., 2007. Water wave propagation using sph models, in: Proceedings 2nd International Spheric Workshop, pp. 119–122.

Hall, M., 2015. Moordyn user's guide. Orono, ME: Department of Mechanical Engineering, University of Maine.

Hall, M., Goupee, A., 2015. Validation of a lumped-mass mooring line model with deepwind semisubmersible model test data. Ocean Engineering 104, 590–603.

Hall, M., Housner, S., Sirnivas, S., Wilson, S., et al., 2021. MoorPy (Quasi-Static Mooring Analysis in Python). Technical Report. National Renewable Energy Lab.(NREL), Golden, CO (United States).

Hasan, S., Sriram, V., Selvam, R.P., 2019. Evaluation of an eddy viscosity type wave breaking model for intermediate water depths. European Journal of Mechanics-B/Fluids 78, 115–138.

Higuera, P., Lara, J.L., Losada, I.J., 2013. Realistic wave generation and active wave absorption for navier-stokes models: Application to openfoam®. Coastal Engineering 71, 102–118.

Hildebrandt, A., Sriram, V., Schlurmann, T., 2013. Simulation of focusing waves and local line forces due to wave impacts on a tripod structure, in: The Twenty-third International Offshore and Polar Engineering Conference, OnePetro.

Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (vof) method for the dynamics of free boundaries. Journal of computational physics 39, 201–225.

Huang, Y., Zhuang, Y., Wan, D., 2021. Hydrodynamic study and performance analysis of the oc4-deepwind platform by cfd method. International Journal of Computational Methods 18, 2050020.

Huseby, M., Grue, J., 2000. An experimental investigation of higher-harmonic wave forces on a vertical cylinder. Journal of fluid Mechanics 414, 75–103.

Jacobsen, N.G., Fuhrman, D.R., Fredsøe, J., 2012. A wave generation toolbox for the open-source CFD library: OpenFoam®. International Journal for Numerical Methods in Fluids 70, 1073–1088. doi:10.1002/fld.2726.

Jasak, H., Tukovic, Z., 2006. Automatic mesh motion for the unstructured finite volume method. Transactions of FAMENA 30, 1–20.

Jasak, H., Tuković, Ž., 2010. Dynamic mesh handling in openfoam applied to fluid-structure interaction simulations, in: Proceedings of the V European Conference on Computational Fluid Dynamics ECCOMAS CFD 2010.

Jiang, C., el Moctar, O., Schellin, T.E., Paredes, G.M., 2021. Comparative study of mathematical models for mooring systems coupled with cfd. Ships and Offshore Structures 16, 942–954.

Jonkman, J., 2010. Definition of the Floating System for Phase IV of OC3. Technical Report. National Renewable Energy Lab.(NREL), Golden, CO (United States).

Journee, J.M., Adegeest, L.J., 2003. Theoretical manual of strip theory program" seaway for windows". TU Delft, Report 1370.

Karnoski, S., Palo, P., 1988. Validation of a static mooring analysis model with full-scale data, in: Offshore Technology Conference, OnePetro.

Kemper, J., Windt, C., Graf, K., Ringwood, J., 2019. Development towards a nested hydrodynamic model for the numerical analysis of ocean wave energy systems, in: European Tidal and Wave Energy Conference Proceedings, EWTEC.

Kim, Y.J., 2021. Numerical improvement and validation of a naval hydrodynamics CFD solver in view of performing fast and accurate simulation of complex ship-wave interaction. Ph.D. thesis. École centrale de Nantes.

Lachaume, C., Biausser, B., Fraunié, P., Grilli, S.T., Guignard, S., 2003. Modeling of breaking and post-breakingwaves on slopes by coupling of bem and vof methods, in: The Thirteenth International Offshore and Polar Engineering Conference, OnePetro.

Lee, S.C., Song, S., Park, S., 2021. Platform motions and mooring system coupled solver for a moored floating platform in a wave. Processes 9, 1393.

Li, H., Bachynski, E.E., 2021. Experimental and numerical investigation of nonlinear diffraction wave loads on a semi-submersible wind turbine. Renewable Energy 171, 709–727.

Li, Z., 2018. Two-phase spectral wave explicit Navier-Stokes equations method for wave-structure interactions. Ph.D. thesis. École centrale de Nantes.

Li, Z., Bouscasse, B., Ducorzet, G., Gentaz, L., Le Touzé, D., Ferrant, P., 2021. Spectral wave explicit navier-stokes equations for wave-structure interactions using two-phase computational fluid dynamics solvers. Ocean Engineering 221, 108513.

Löhner, R., Yang, C., 1996. Improved ale mesh velocities for moving bodies. Communications in numerical methods in engineering 12, 599–608.

Lu, X., Chandar, D.D.J., Chen, Y., Lou, J., 2017. An overlapping domain decomposition based near-far field coupling method for wave structure interaction simulations. Coastal Engineering 126, 37–50.

Make, M., 2014. Predicting scale effects on floating offshore wind turbines.

Masciola, M., 2018. Map++ documentation. NREL: Golden, CO, USA.

Masciola, M., Jonkman, J., Robertson, A., 2013. Implementation of a multisegmented, quasi-static cable model, in: The Twenty-third International Offshore and Polar Engineering Conference, OnePetro.

Mayer, S., Garapon, A., Sørensen, L., 1998. A fractional step method for unsteady free-surface flow with applications to non-linear wave dynamics. International Journal for Numerical Methods in Fluids - INT J NUMER METHOD FLUID 28, 293–315. doi:10.1002/(SICI)
1097-0363(19980815)28:23.0.C0;2-1.

Michel, G., Bonnefoy, F., Ducrozet, G., Prabhudesai, G., Cazaubiel, A., Copie, F., Tikan, A., Suret, P., Randoux, S., Falcon, E., 2020. Emergence of peregrine solitons in integrable turbulence of deep water gravity waves. Physical Review Fluids 5, 082801.

Miquel, A.M., Kamath, A., Alagan Chella, M., Archetti, R., Bihs, H., 2018. Analysis of different methods for wave generation and absorption in a cfd-based numerical wave tank. Journal of Marine Science and Engineering 6, 73.

Nematbakhsh, A., Bachynski, E.E., Gao, Z., Moan, T., 2015. Comparison of wave load effects on a tlp wind turbine by using computational fluid dynamics and potential flow theory approaches. Applied Ocean Research 53, 142-154.

OpenFOAM, OpenFOAM v9 User Guide. URL: https://cfd.direct/openfoam/user-guide.

ORCINA, OrcaFlex—Documentation. URL: https://www.orcina.com/releases/orcaflex-101/.

Paduano, B., Giorgi, G., Gomes, R.P., Pasta, E., Henriques, J.C., Gato, L., Mattiazzo, G., 2020. Experimental validation and comparison of numerical models for the mooring system of a floating wave energy converter. Journal of Marine Science and Engineering 8, 565.

Palm, J., Eskilsson, C., Paredes, G.M., Bergdahl, L., 2016. Coupled mooring analysis for floating wave energy converters using cfd: Formulation and validation. International Journal of Marine Energy 16, 83-99.

Paulsen, B.T., Bredmose, H., Bingham, H.B., 2014. An efficient domain decomposition strategy for wave loads on surface piercing circular cylinders. Coastal Engineering 86, 57-76.

Pinguet, R., 2021. Hydrodynamics of semi-submersible floater for offshore wind turbines in highly nonlinear waves using Computational Fluid Dynamics (CFD), and validation of overset meshing technique in a numerical wave tank. Ph.D. thesis. Ecole Centrale Marseille.

Randolph, M., Quiggin, P., 2009. Non-linear hysteretic seabed model for catenary pipeline contact, in: International Conference on Offshore Mechanics and Arctic Engineering, pp. 145-154.

Rienecker, M.M., Fenton, J.D., 1981. A fourier approximation method for steady water waves. Journal of fluid mechanics 104, 119-137.

Robertson, A., Jonkman, J., Masciola, M., Song, H., Goupee, A., Coulling, A., Luan, C., 2014. Definition of the semisubmersible floating system for phase II of OC4. Technical Report. National Renewable Energy Lab.(NREL), Golden, CO (United States).

Rusche, H., 2002. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. Ph.D. thesis.

Seng, S., 2012. Slamming and whipping analysis of ships. DTU Mechanical Engineering.

Simonsen, C.D., Otzen, J.F., Joncquez, S., Stern, F., 2013. Efd and cfd for kcs heaving and pitching in regular head waves. Journal of Marine Science and Technology 18, 435-459.

Sriram, V., Agarwal, S., Yan, S., Xie, Z., Saincher, S., Schlurmann, T., Ma, Q., Stoesser, T., Zhuang, Y., Han, B., et al., 2021. A comparative study on the nonlinear interaction between a focusing wave and cylinder using state-of-the-art solvers: Part a. International Journal of Offshore and Polar Engineering 31, 1-10.

Sriram, V., Ma, Q., Schlurmann, T., 2014. A hybrid method for modelling two dimensional non-breaking and breaking waves. Journal of computational physics 272, 429-454.

Sun, H., 2007. A boundary element method applied to strongly nonlinear wave-body interaction problems.

Thomsen, J.B., Eskilsson, C., Ferri, F., 2017. Assessment of available numerical tools for dynamic mooring analysis: Wp1. 2 & m1.

Vickers, A., 2012. Improve the understanding of uncertainties in numerical analysis of moored floating wave energy converters. Ph.D. thesis. University of Exeter.

West, B.J., Brueckner, K.A., Janda, R.S., Milder, D.M., Milton, R.L., 1987. A new numerical method for surface hydrodynamics. Journal of Geophysical Research: Oceans 92, 11803-11824.

Xu, H., Zhang, Y., Santo, H., Chua, K.H., Law, Y.Z., Chan, E.S., 2021. Coupling of potential flow and cfd model for fluid and structure interactions, in: International Conference on Offshore Mechanics and Arctic Engineering, American Society of Mechanical Engineers. p. V008T08A009.

Zhang, Y., Kim, B., 2018. A fully coupled computational fluid dynamics method for analysis of semi-submersible floating offshore wind turbines under wind-wave excitation conditions based on oc5 data. Applied sciences 8, 2314.

CRediT authorship contribution statement

**Sithik Aliyar**: Conceptualization of this study, Methodology, Software, Writing - Original draft preparation.

**Guillaume Ducrozet**: Data curation, Writing - Original draft preparation. **Benjamin Bouscasse**: Methodology,

Writing - Original draft preparation. **Félicien Bonnefoy**: Data curation, Writing - Original draft preparation. **Sriram**

**Venkatachalam**: Writing - Original draft preparation. **Pierre Ferrant**: Writing - Original draft preparation.

Author biography without author photo.