# Progress in Coupling Potential Wave Models and Two-Phase Solvers With the SWENSE Methodology

Zhaobin Li, Benjamin Bouscasse, L. Gentaz, Guillaume Ducrozet, Pierre Ferrant

► To cite this version:

Zhaobin Li, Benjamin Bouscasse, L. Gentaz, Guillaume Ducrozet, Pierre Ferrant. Progress in Coupling Potential Wave Models and Two-Phase Solvers With the SWENSE Methodology. ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, Jun 2018, Madrid, Spain. 10.1115/OMAE2018-77466. [hal-02885122](https://hal.science/hal-02885122)

HAL Id: hal-02885122

https://hal.science/hal-02885122

Submitted on 21 May 2021

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
PROGRESS IN COUPLING POTENTIAL WAVE MODELS AND TWO-PHASE
SOLVERS WITH THE SWENSE METHODOLOGY

Zhaobin Li

LHEEA Lab

École Centrale de Nantes, France
CNRS UMR 6598

Email: zhaobin.li@ec-nantes.fr

Benjamin Bouscasse

LHEEA Lab

École Centrale de Nantes, France
CNRS UMR 6598

Email: benjamin.bouscasse@ec-nantes.fr

Lionel Gentaz

LHEEA Lab

École Centrale de Nantes, France
CNRS UMR 6598

Email: lionel.gentaz@ec-nantes.fr

Guillaume Ducrozet

LHEEA Lab

École Centrale de Nantes, France
CNRS UMR 6598

Email: guillaume.ducrozet@ec-nantes.fr

Pierre Ferrant

LHEEA Lab

École Centrale de Nantes, France
CNRS UMR 6598

Email: pierre.ferrant@ec-nantes.fr

ABSTRACT

This paper presents the recent developments of the Spectral Wave Explicit Navier-Stokes Equations (SWENSE) method to extend its range of application to two-phase VOF solvers. The SWENSE method solves the wave-structure interaction problem by coupling potential theory and the Navier-Stokes (NS) equations. It evaluates the incident wave solution by wave models based on potential theory in the entire computational domain, leaving only the perturbation caused by the structure and the influence of the viscosity to be solved with CFD. The method was proven in previous studies to be accurate and efficient for wave-structure interaction problems, but it was derived for single-phase NS solvers only. The present study extends the SWENSE method by proposing a novel formulation which is convenient to implement in two-phase NS solvers. A customized SWENSE solver is developed with the open source CFD package OpenFOAM. An improvement in accuracy and stability is observed in wave simulations compared with conventional two-phase VOF solvers. The horizontal force on a vertical cylinder in regular waves is also calculated. First results show a good agreement with the experiment on the first harmonic component.

NOMENCLATURE

u Velocity Field

p Pressure Field

α Volume of Fluid (VOF) Field

χ Variables without subscript represent the total solution

χ₁ Subscript I represents the incident wave solution

χc Subscript C represents the complementary solution

INTRODUCTION

Using Computational Fluid Dynamics (CFD) softwares to simulate wave-structure interaction is of great interest, since CFD codes are able to consider more complex phenomena occurring in the real world, i.e. the viscous effects, the violent free surface deformation, etc., compared with potential theory based methods.

However, simulating wave-structure interaction with a conventional CFD tool is computational resource demanding. It often requires a computational zone large enough both on the upstream of the structure to prevent the disturbed waves traveling back to the inlet boundary and a large damping zone on the
downstream to avoid the wave reflections at the outlet boundary. Although it is possible to reduce the damping zone's size by giving the wave solutions at boundaries with domain decomposition methods [1, 2], fine meshes are always needed to prevent excessive numerical attenuation of incoming waves [3]. All of these requirements make the computational cost of wave-structure interaction problems with CFD tools expensive for practical applications.

To overcome this problem, the Spectral Wave Explicit Navier-Stokes Equations (SWENSE) method was proposed to reduce the computational cost of CFD by coupling the solution from potential theory and the Navier-Stokes (NS) equations [4]. In this formalism, the incident wave solution is evaluated by a dedicated model in the entire computational domain and is subtracted from the NS equations, leaving only the perturbation part caused by the viscosity and the presence of structures to be solved with CFD. Since the incident wave solution is explicitly given, it is not influenced by the CFD's resolution. Coarse meshes can be used in the farfield region. A good mesh quality is only necessary close to the structure, since in most cases a well resolved disturbance field is only required near the structure. The total number of cells to reach a given accuracy is reduced [5].

The original SWENSE method was proved in previous studies to be accurate and efficient in its original single-phase form [6, 7]. Recently, a two-phase SWENSE method was proposed [8]. However, the latter behaves differently from the original SWENSE method: the incident waves are influenced by the mesh quality. This is because it only decomposes the velocity field, instead of subtracting the entire incident solutions from the NS equations.

In the present study, the authors propose a novel two-phase SWENSE method that complies with the original single-phase SWENSE method. Mathematics proves that by using the novel formulations, the incident waves are maintained regardless of CFD's resolution. The method is implemented in a two-phase solver based on the Volume of Fluid (VOF) method using the open-source CFD package OpenFOAM [9]. The method is tested with two cases. The first concerns the propagation of regular waves in a 2D wave tank. The second test deals with the calculation of the wave forces on a vertical cylinder.

The paper is organized as follows. In the first part, we introduce briefly the SWENSE method in its original single-phase form. The extension to two-phase solver is presented in detail in the second part and is followed by its implementation in OpenFOAM. Lastly, the two test cases are shown with discussions.

## SINGLE-PHASE SWENSE METHOD

The SWENSE method decomposes the wave-structure interaction problem into an incident part and a complementary part:

1. the incident part concerns the propagation of the incident waves in the computational domain without structures. This solution is given directly by a dedicated wave model based on potential theory;

2. the complementary part serves as a correction to the incident part due to the disturbance caused by the viscosity and the presence of the structures.

This decomposition is shown by Eqn. 1 where a primitive field of the flow $\chi$ (pressure, velocity, and free surface elevation) is divided into the incident part $\chi_I$ and the complementary part $\chi_C$. $\chi_I$ is explicitly evaluated by the wave model; $\chi_C$ is to be calculated by the CFD solver. The decomposition is also illustrated by Fig. 1.

$$ \chi = \chi_I + \chi_C \quad (1) $$

<figure><img src="image_2.png" /><figcaption>FIGURE 1: SWENSE DECOMPOSITION</figcaption></figure>

The advantages of the SWENSE method are:

1. The accuracy of the incident waves: it is explicitly obtained by a dedicated wave model so it is not influenced by the CFD resolution;

2. The efficiency: for a given accuracy, the SWENSE method can use coarser meshes in the farfield compared with the conventional CFD methods;

3. The farfield boundary conditions are simple even in complex irregular sea states: the complementary fields are forced to vanish at the farfield boundary.

The governing equations of the single-phase SWENSE method are derived via solution decomposition. The Euler equations behind the wave models are subtracted from the Navier-Stokes equations, as follows [4, 6].

## Navier-Stokes Equations

The Navier-Stokes equations for incompressible single-phase flow are:
$$ \nabla \cdot \mathbf{u} = 0 \qquad (2) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \frac{\nabla p}{\rho} + \mathbf{g} + \nu \nabla^2 \mathbf{u} \qquad (3) $$

where $\mathbf{u}$ is the velocity, $p$ is the pressure. $\rho$ stands for the density of water and $\nu$ stands for the kinematic viscosity. $\mathbf{g}$ stands for the gravitational acceleration.

## Euler Equations

The incident wave solution comes from the perfect fluid Euler equations, as follows:

$$ \nabla \cdot \mathbf{u}_I = 0 \qquad (4) $$

$$ \frac{\partial \mathbf{u}_I}{\partial t} + \mathbf{u}_I \cdot \nabla \mathbf{u}_I = - \frac{\nabla p_I}{\rho} + \mathbf{g} \qquad (5) $$

## SWENSE Equations

Subtracting the Euler equations from the Navier-Stokes equations, we obtain:

$$ \nabla \cdot (\mathbf{u} - \mathbf{u}_I) = 0 \qquad (6) $$

$$ \begin{aligned} & \frac{\partial (\mathbf{u} - \mathbf{u}_I)}{\partial t} + (\mathbf{u} - \mathbf{u}_I) \cdot \nabla (\mathbf{u} - \mathbf{u}_I) + (\mathbf{u} - \mathbf{u}_I) \cdot \nabla \mathbf{u}_I + \mathbf{u}_I \cdot \nabla (\mathbf{u} - \mathbf{u}_I) \\ & = - \frac{\nabla (p - p_I)}{\rho} + \nu \nabla^2 (\mathbf{u} - \mathbf{u}_I) + \nu \nabla^2 \mathbf{u}_I \end{aligned} \qquad (7) $$

The last term vanishes as $\nabla^2 \mathbf{u}_I$ equals to zero when $\mathbf{u}_I$ is derived from potential wave theory.

Noting $\mathbf{u}_C = \mathbf{u} - \mathbf{u}_I$ and $p_C = p - p_I$, the governing equations of the single-phase SWENSE method read:

$$ \nabla \cdot \mathbf{u}_C = 0 \qquad (8) $$

$$ \frac{\partial \mathbf{u}_C}{\partial t} + \mathbf{u}_C \cdot \nabla \mathbf{u}_C + \mathbf{u}_C \cdot \nabla \mathbf{u}_I + \mathbf{u}_I \cdot \nabla \mathbf{u}_C = - \frac{\nabla p_C}{\rho} + \nu \nabla^2 \mathbf{u}_C \qquad (9) $$

## TWO-PHASE SWENSE METHOD

This section presents the development of the two-phase SWENSE governing equations. Although the same idea as in the single-phase SWENSE is followed, the extension of the

SWENSE method to two-phase solvers is not that straightforward, because a direct subtraction of the Euler equations from the Navier-Stokes equations generates numerical instabilities. A modification on the Euler equations is necessary and it is the vital key to obtain a numerically stable form of the governing equations so that the method can be used by two-phase VOF based CFD solvers.

## Two-phase Navier-Stokes Equations

The Navier-Stokes equations for two-phase incompressible flow, written with Volume of Fluid (VOF) method are:

$$ \rho = \alpha \rho^w + (1 - \alpha) \rho^a \qquad (10) $$

$$ \mu = \alpha \mu^w + (1 - \alpha) \mu^a \qquad (11) $$

$$ \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0 \qquad (12) $$

$$ \nabla \cdot \mathbf{u} = 0 \qquad (13) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \frac{\nabla p}{\rho} + \mathbf{g} + \nu \nabla^2 \mathbf{u} \qquad (14) $$

where $w$ and $a$ stand for the properties in the water and in the air respectively. The $\alpha$ is the VOF field representing the volume rate of water in computational cells. It equals to 1 for a cell full of water and equals to 0 when a cell is full of air.

## Euler Equations

The SWENSE decomposition requires the incident wave solution above the incident wave free surface, when the actual free surface is above the incident one. Thanks to the spectral methods used by SWENSE, the information above the incident wave free surface can still be calculated and is continuous across the interface. Moreover the results obtained from this treatment verifies the Euler equations mathematically in the whole computational domain. With the assumption that such result can be used as the incident solution of the air phase, this method unifies the governing equations of incident waves in both air and water with the single-phase Euler equations, as follows.

$$ \rho_I = \rho^w \qquad (15) $$

$$ \nabla \cdot \mathbf{u}_I = 0 \qquad (16) $$

$$ \frac{\partial \mathbf{u}_I}{\partial t} + \mathbf{u}_I \cdot \nabla \mathbf{u}_I = - \frac{\nabla p_I}{\rho_I} + \mathbf{g} \qquad (17) $$

However, using directly Eqn.17 to derive the two-phase SWENSE momentum equation causes stability problems. To overcome this difficulty, a modified Euler momentum equation
is proposed by introducing a different incident pressure $p_I^*$:

$$p_I^* = \rho \frac{p_I}{\rho_I} \qquad (18)$$

so

$$\nabla p_I^* = \nabla \left( \frac{\rho}{\rho_I} p_I \right) = \frac{p_I}{\rho_I} \nabla \rho + \frac{\rho}{\rho_I} \nabla p_I \qquad (19)$$

and

$$\frac{\nabla p_I}{\rho_I} = \frac{\nabla p_I^*}{\rho} - \frac{p_I}{\rho_I} \frac{\nabla \rho}{\rho} \qquad (20)$$

Eqn.17 written in its modified version by using $p_I^*$ reads:

$$\frac{\partial \mathbf{u}_I}{\partial t} + \mathbf{u}_I \cdot \nabla \mathbf{u}_I = - \frac{\nabla p_I^*}{\rho} + \frac{p_I}{\rho_I} \frac{\nabla \rho}{\rho} + \mathbf{g} \qquad (21)$$

## SWENSE Equations

The two-phase SWENSE equations are obtained by subtracting the modified Euler equations from the two-phase Navier-Stokes equations. The continuity equation for SWENSE variable is obtained directly by the subtracting Eqn.13 with Eqn.16.

$$\nabla \cdot (\mathbf{u} - \mathbf{u}_I) = 0 \qquad (22)$$

Noting $\mathbf{u}_C = \mathbf{u} - \mathbf{u}_I$, the continuity equation using the complementary variable reads:

$$\nabla \cdot \mathbf{u}_C = 0 \qquad (23)$$

Subtracting Eqn.21 from the two-phase Navier-Stokes equation Eqn.14, and using the notion of $p_C = p - p_I^*$, the two-phase SWENSE momentum equation written with the complementary variables is obtained as:

$$\frac{\partial \mathbf{u}_C}{\partial t} + \mathbf{u}_C \cdot \nabla \mathbf{u}_C + \mathbf{u}_C \cdot \nabla \mathbf{u}_I + \mathbf{u}_I \cdot \nabla \mathbf{u}_C - v \nabla^2 \mathbf{u}_C = - \frac{\nabla p_C}{\rho} - \frac{p_I}{\rho_I} \frac{\nabla \rho}{\rho} \qquad (24)$$

The VOF field is not decomposed as $\alpha = \alpha_I + \alpha_C$, because it is difficult to define the value range of the complementary field $\alpha_C$ so that the boundedness of the total field $\alpha$ is ensured. In the present method, the $\alpha$ field is transported by the total velocity field $\mathbf{u}$ using the advection equation Eqn.12.

A good property of the SWENSE momentum equation is that it keeps the accuracy of the incident waves regardless of the mesh quality. Take a pure incident wave propagation as an example. In CFD solvers, no disturbance of waves is expected if the viscosity effect is negligible, i.e. the result of CFD should equal to the incident wave solution given by potential theory at all time. In the SWENSE formalism, that means the complementary fields remain zero. The following paragraph proves mathematically that the SWENSE momentum equation ensures this characteristics.

Letting the initial value of $\mathbf{u}_C$ and the $p_C$ equal to zero, all the convective and source terms except the last one on the RHS in the Eqn.24 can be canceled out. The last term is a product of $\frac{\nabla \rho}{\rho}$ and $\frac{p_I}{\rho_I}$. The term $\frac{\nabla \rho}{\rho}$ equals to 0 except at the free surface. Coincidentally, the term $\frac{p_I}{\rho_I}$ is 0 at the free surface according to the boundary condition of the potential wave theory: the pressure on the free surface equals to 0. So the last term also equals to 0. It is obvious that the values of $\mathbf{u}_C$ and $p_C$ equal to 0 at next time step. So no disturbance waves will be generated by the CFD solver.

## IMPLEMENTATION IN TWO-PHASE VOF SOLVER

The proposed two-phase SWENSE method is implemented in a two-phase VOF solver using the open-source CFD package OpenFOAM.

OpenFOAM uses finite volume method with unstructured polyhedral meshes. Its native solver for incompressible two-phase flow, interFoam, adopts VOF method and MULES algorithm to keep the boundedness of VOF field and the sharpness of the interface. The PIMPLE algorithm is used to obtain converged results of the velocity-pressure-VOF coupling at each time step.

Based on this solver, various software packages for water wave modeling have been developed by the community of OpenFOAM, including waves2Foam [10], foamStar [11], etc. These packages use interFoam as the core solver and extend it with various wave theories and methods to generate and absorb waves at the inlet and outlet of the computational domain.

The new solver foamStar-SWENSE is developed in the frame of the foamStar package of Bureau Veritas. It replaces the original two-phase Navier-Stokes equations with the two-phase SWENSE equations. It adopts the relaxation zone method in foamStar to attenuate the complementary waves in the farfield.

## TEST CASE ONE: TWO-DIMENSIONAL WAVE TANK

The validity of the proposed method is tested firstly by a regular wave propagation problem in a 2D wave tank. In order to demonstrate the robustness and the accuracy of the method, this test case is intentionally designed to have a large computational domain and last for a sufficiently long time to increase the numerical difficulties.
## Test Case Setup

The test uses steady-state propagating waves with a moderately high wave steepness ($ka = 0.224$). The incident wave information is obtained by the solution method of Rienecker & Fenton [12] based on the stream function wave theory. The wave characteristics and the parameters are listed in Tab.1.

<table><caption>TABLE 1: WAVE CHARACTERISTICS FOR THE CASE OF TWO-DIMENSIONAL WAVE TANK</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Wave length (λ)</td><td>0.8082 m</td></tr><tr><td>Wave period (T)</td><td>0.7017 s</td></tr><tr><td>Wave height (H)</td><td>0.0575 m</td></tr><tr><td>Wave steepness (ka)</td><td>0.2235</td></tr><tr><td>Water depth (h)</td><td>0.6000 m</td></tr><tr><td>Relative water depth (kh)</td><td>4.6645</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>FIGURE 2: COMPUTATIONAL DOMAIN</figcaption></figure>

A two-dimensional rectangular computational domain is used as illustrated by Fig.2. The waves travel from the left to the right. Periodic boundary conditions are applied on the inlet and the outlet boundaries. The origin of the coordinate system locates at the left of the computational domain on the free surface position at rest. The axis x points towards the direction of the wave propagation and the axis y points upward. The length of the computational domain equals to ten times the wave length ($λ$), and the height of the computational domain equals to the water depth (h) plus 2.5 times the wave height (H). The discretization information is summarized in Tab.2.

The initial values of the VOF field, the incident wave velocity, and the incident pressure fields are calculated using the stream function wave theory. The initial conditions for the complementary

velocity field and the complementary pressure field are set to zero.

For comparison, another simulation using a conventional two-phase VOF solver *foamStar* is carried out. It is worth noticing that the discretization used by *foamStar* is twice more refined both in time and in space to stabilize the comparative simulation. Indeed, the computation of *foamStar* always stops due to stability problem if it uses the same discretization as *foamStar-SWENSE*.

<table><caption>TABLE 2: SIMULATION PARAMETERS</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Domain size (L<sub>x</sub> × L<sub>y</sub>)</td><td>10λ × (h + 2.5H)</td></tr><tr><td>Mesh size (Δx, Δy)</td><td>λ/100, H/20</td></tr><tr><td>Simulation time (t)</td><td>40T</td></tr><tr><td>Time step (Δt)</td><td>T/400</td></tr></tbody></table>

## Numerical Results

The free elevation of the two numerical simulations are compared with the reference value from the wave theory. Fig.3 plots the free surface elevation $ζ$ against the distance from the inlet x every ten wave periods. The surface elevation is non-dimensionalized by the incident wave amplitude A and the distance is non-dimensionalized by the wave length $λ$.

It is observed that the free surface's shape is well kept by *foamStar-SWENSE* even after a very long simulation time. The difference between the SWENSE method's results and the reference values is not obvious during the entire 40 wave periods. On the contrary, a discrepancy between the *foamStar*'s result and the reference solution is remarkable after 20 wave periods.

A moving window harmonic analysis is applied to the time history of the free surface elevation at the center of domain. Fig. 4 plots the first and second harmonic amplitudes against the simulation time. The amplitude is non-dimensionalized by the incident wave amplitude and the time is non-dimensionalized by the wave period. For *foamStar-SWENSE*: Fig.4a shows the first harmonic amplitude of the waves which is almost constant for the first 20 wave periods. After 20 periods, the amplitude starts to decrease; at the end of the simulation, the first harmonic amplitude is about 94% of the reference value. Fig.4b shows the second harmonic component amplitude which is also well kept with only 10% maximum difference from the target value. In the contrary, the results of *foamStar* have larger oscillations in both the first and second harmonic amplitudes.
<figure><img src="image_4.png" /><figcaption>FIGURE 3: FREE SURFACE ELEVATION</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 4: HARMONIC AMPLITUDES OF FREE SURFACE ELEVATION AT THE CENTER OF DOMAIN</figcaption></figure>

## Discussion

The above results illustrate the capability of the two-phase SWENSE method for wave propagation problems and its advantages over conventional VOF solvers. Despite the large computational domain and the long simulation time, the two-phase SWENSE method gives quite good results, showing a good accuracy and stability.

However, it is worth noticing that the present study shows a discrepancy between the simulation result and the incident wave solution even though no difference is expected according to SWENSE's theory. This difference is due to the numerical error in evaluating the term $\frac{\partial_t \nabla \rho}{\rho_t \rho}$ in Eqn.24. Although it equals to zero in SWENSE's theory, the numerical evaluation of the term is non-zero on the Finite Volume mesh. The numerical error on that term creates non-zero complementary fields as soon as the simulation starts. Fig.5 shows the complementary velocity $\mathbf{u}_C$ and the complementary pressure field $p_C$ after the first time step

at $t$ equals to $\Delta t$. From the picture it is confirmed that the complementary fields are created only near the free surface with very small values. The maximum $\mathbf{u}_C$ value is about 1% of the maximum incident velocity. This small error is not problematic in most cases. However, it is responsible for the discrepancy observed between the simulation result and the target value after a long simulation.

Regarding *foamStar*, the numerical errors manifest from the 15th wave period approximately. This is because the solver has more difficulties in balancing between the accuracy and stability in such a large-domain and long-time simulation. The dilemma is: for transient problems such as waves, both the temporal and spatial discretization should be at least second-order accurate [13, 14] which is the upper limit for Finite Volume Method codes using unstructured meshes; moreover, the pure second-order schemes create stability problems and are usually blended with first order ones [9, 15]. The latter induces numerical damping and should be avoid if possible. The result is that despite a more refined temporal and spatial discretization, obvious deformations of the wave free surface shape still happens in the result of *foamStar* after 20 wave periods. This is a common difficulty for numerical wave tanks using the second-order finite volume method.

The same difficulty remains for *foamStar-SWENSE*, since the same temporal and spatial discretization schemes. However, its influence is much smaller, as only the complementary part is affected. The incident wave solution is maintained, leading to a much more accurate and stable result. Compared with the conventional two-phase CFD methods, the primary advantage of the SWENSE method is an accurate and stable computation of the incident waves, which is of vital importance in the simulation of wave-structure interaction problems.

## TEST CASE TWO: CYLINDER IN REGULAR WAVES

To validate the SWENSE method for the wave-structure interaction problem, we consider a bottom-mounted surface piercing vertical circular cylinder in regular waves in this section. The setup of the simulations follows the published experimental study [16].

The cylinder has a radius $R$ equals to $0.03m$, being fixed in a water tank with depth $h$ equals to $0.6m$. The incident wave frequency $f$ equals to $1.425Hz$. The wave steepness $ka$ is $0.24$. The parameters of the test case are summarized in Tab.3.

The total simulation time corresponds to 15 wave periods. The inline wave force on the cylinder is recorded and then analyzed with moving window harmonic analysis. The first harmonic amplitude is compared with the reference data. Three different resolutions are used in this study to investigate the convergence property.
<figure><img src="image_6.png" /><figcaption>FIGURE 5: COMPLEMENTARY FIELD AFTER THE FIRST TIME STEP</figcaption></figure>

<table><caption>TABLE 3: WAVE PARAMETERS FOR THE CASE OF CYLINDER IN WAVES</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Cylinder radius (r)</td><td>0.030m</td></tr><tr><td>Water depth (h)</td><td>0.600m</td></tr><tr><td>Wave frequency (f)</td><td>1.425Hz</td></tr><tr><td>Wave number (k)</td><td>8.172</td></tr><tr><td>Wave steepness (ka)</td><td>0.240</td></tr><tr><td>Relative water depth (kh)</td><td>4.903</td></tr><tr><td>Simulation time (t)</td><td>15T</td></tr></tbody></table>

Note: wave length and wave number are calculated with the linear dispersion relation [16]

## Numerical Setup

A 3D cylindrical mesh is used. The origin of the coordinate system locates at the free surface position of the water at rest, axis x points to the incident wave propagation direction and axis z points upward. The domains radius equals to 2 wave lengths ($R = 2\lambda$). The mesh is fine near the cylinder and gradually enlarged along the radius direction. In the farfield, the relaxation zone with a length of $1.5\lambda$ is used to absorb the disturbed waves, leaving a pure CFD zone with one wave length $r \in (-0.5\lambda, 0.5\lambda)$.

In the tangential direction, only half of the domain is modeled by using the symmetric assumption. Along the vertical axis, the computational domain extends from the tank bottom at $z = -0.60m$ until $z = 2.5H$ in the air, with $H = 2a$ being the wave height. In this direction, the meshes are uniformly refined in the zone near the air-water interface. The size of this zone is $z \in (-1.5a, 1.5a)$. Out of this refined zone, the mesh size increases gradually.

The information of the finest discretization is summarized in Tab.4 and the grid is illustrated by Fig.6. The other two discretization are derived from this setup by enlarging the grid size and time step proportionally. The refinement factor is approximately 1.4. The information about them is shown in Tab. 5.

## Results and Discussion

The time history of the horizontal force on the cylinder is firstly recorded and analyzed by the moving window harmonic analysis. The amplitude of the first harmonic component is extracted and compared with the experimental data [16] in Fig.7. The non-dimensional first harmonic amplitude is defined as follows.

$$F'_{1} = \frac{F_{1}}{\rho g R^{3}} \cdot \left(\frac{R}{A}\right) \qquad (25)$$

The first order amplitudes at steady state are extracted. The
<table><caption>TABLE 4: FINEST DISCRETIZATION INFORMATION</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain size (L<sub>r</sub> × L<sub>θ</sub> × L<sub>z</sub>)</td><td>2λ × 180 × (h + 5a)</td></tr><tr><td>Grids in radial direction (n<sub>r</sub>)</td><td>100</td></tr><tr><td>Grids in tangential direction (n<sub>θ</sub>)</td><td>80</td></tr><tr><td>Grids in vertical direction (n<sub>z</sub>)</td><td>136</td></tr><tr><td>Grids in one wave amplitude (a/Δz)</td><td>12</td></tr><tr><td>Time steps in one wave period (T/Δt)</td><td>400</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIGURE 6: GRID FOR CYLINDER IN WAVES</figcaption></figure>

results are 6.09, 6.23, and 6.30 for the discretization number 1, 2 and 3.

In terms of the accuracy, the results show a good agreement with the experimental data. The relative differences compared with the experimental data are -2.71%, -0.48% and 0.64% for the coarsest to the finest resolution.

Regarding the efficiency, the necessary number of cells is relatively small compared with conventional two-phase VOF solvers [11] where the recommended mesh is Cartesian with at least 60 cells per wave length [17]. Fig.8 plots the first harmonic amplitude of the in-line force against the number of cells for the present study and the reference results from [11] using *foam-Star*, a conventional two-phase VOF solver in OpenFOAM. The necessary number of cells to achieve a given accuracy is much smaller for *foamStar-SWENSE* than *foamStar*. The reduction of cell number leads to shorter computational time which is a second key advantage of the SWENSE method.

<table><caption>TABLE 5: DISCRETIZATIONS FOR CONVERGENCE STUDY</caption><thead><tr><th rowspan="2">Index</th><th colspan="2">Representative mesh information</th><th rowspan="2">Time step<br>T/Δt</th></tr><tr><th>Total cells</th><th>Cells per wave amplitude</th></tr></thead><tbody><tr><td>1</td><td>88000</td><td>6</td><td>200</td></tr><tr><td>2</td><td>246512</td><td>9</td><td>285</td></tr><tr><td>3</td><td>668000</td><td>12</td><td>400</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>FIGURE 7: FIRST HARMONIC AMPLITUDES COMPARED WITH THE EXPERIMENTAL DATA</figcaption></figure>

## CONCLUSION

The present paper reported the latest developments of the SWENSE method in coupling wave models and two-phase VOF CFD solvers for wave-structure interaction problems. The two-phase SWENSE equations were mathematically derived. After being implemented in OpenFOAM, the method was tested through two cases. The 2D wave tank case demonstrated the capability of the new method to solve wave propagation problem; an improvement of stability and accuracy was observed. In a first step of validation, the method is tested on calculating wave forces on a vertical cylinder. The necessary mesh number to achieve a given accuracy was proven to be smaller by using the SWENSE method than using conventional two-phase VOF methods. The efficiency of the proposed method to solve wave-structure interaction problems was demonstrated.
<figure><img src="image_9.png" /><figcaption>FIGURE 8: MESH CONVERGENCE STUDY OF THE FIRST ORDER AMPLITUDES OF HORIZONTAL FORCE</figcaption></figure>

## ACKNOWLEDGMENT

This work has been performed in the framework of the Chaire Hydrodynamique et Structure Marines CENTRALE NANTES - BUREAU VERITAS. The authors acknowledge Mr. YoungMyung Choi from Ecole Centrale de Nantes and Dr. Vuko Vukčević from University of Zagreb for their help with the post-processing of the test cases' results. The first author acknowledges the financial support of China Scholarship Council (CSC) for his PhD study.

## REFERENCES

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

[2] Choi, Y., Bouscasse, B., Seng, S., Ducrozet, G., Gentaz, L., and Ferrant, P., 2018. "Generation of regular and irregular waves in navier-stokes CFD solvers by matching with nonlinear potential wave solution at the boundary". In Proceedings of the ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering.

[3] Kim, J., O'Sullivan, J., and Read, A., 2012. "Ringing analysis of a vertical cylinder by Euler overlay method". In Proceedings of the ASME 2012 31st International Conference on Ocean, Offshore and Arctic Engineering, Vol. 10, p. 1.

[4] Ferrant, P., Gentaz, L., Alessandrini, B., and Le Touzé, D., 2003. "A potential/RANSE approach for regular water

wave diffraction about 2-D structures". *Ship Technology Research, 50*(4), pp. 165-171.

[5] Monroy, C., Ducrozet, G., Bonnefoy, F., Babarit, A., Gentaz, L., and Ferrant, P., 2010. "RANS simulations of a calm buoy in regular and irregular seas using the SWENSE method". In *The Twentieth International Offshore and Polar Engineering Conference, International Society of Offshore and Polar Engineers*.

[6] Gentaz, L., Luquet, R., Alessandrini, B., and Ferrant, P., 2004. "Numerical simulation of the 3D viscous flow around a vertical cylinder in non-linear waves using an explicit incident wave model". In *ASME 2004 23rd International Conference on Offshore Mechanics and Arctic Engineering, American Society of Mechanical Engineers*, pp. 157-163.

[7] Li, Z., Gentaz, L., Ducrozet, G., and Ferrant, P., 2017. "Calculation of high-order wave loads on a vertical circular cylinder using the SWENSE method". In *32nd International Workshop on Water Waves and Floating Bodies*.

[8] Vukčević, V., Jasak, H., and Malenica, Š., 2016. "Decomposition model for naval hydrodynamic applications, part I: Computational method". *Ocean Engineering, 121*, pp. 37-46.

[9] OpenCFD Ltd, 2017. *OpenFOAM User Guide*. at https://www.openfoam.com/.

[10] Jacobsen, N. G., Fuhrman, D. R., and Fredsøe, J., 2012. "A Wave Generation Toolbox for the Open-Source CFD Library: OpenFoam®". *Int. J. for Numer. Meth. Fluids, 70*(9), pp. 1073-1088.

[11] Monroy, C., Seng, S., and Malenica, S., 2016. "Développements et validation de l'outil CFD OpenFOAM pour le calcul de tenue à la mer". In *Proceedings of the 15th Journées de l'Hydrodynamique*.

[12] Rienecker, M., and Fenton, J., 1981. "A Fourier approximation method for steady water waves". *Journal of Fluid Mechanics, 104*, pp. 119-137.

[13] Anderson, J. D., and Wendt, J., 1995. *Computational Fluid Dynamics*, Vol. 206. Springer.

[14] Monroy, C., and Seng, S., 2017. "Time-stepping schemes for seakeeping in OpenFOAM". In *Proceedings of the 10th International Workshop on Ship and Marine Hydrodynamics*.

[15] Monroy, C., Seng, S., and Malenica, S., 2017. "On the use of Euler and Crank-Nicolson time-stepping schemes for seakeeping simulations in OpenFOAM". In *Proceeding of VII International Conference on Computational Methods in Marine Engineering*.

[16] Huseby, M., and Grue, J., 2000. "An experimental investigation of higher-harmonic wave forces on a vertical cylinder". *Journal of fluid Mechanics, 414*, pp. 75-103.

[17] Numeca International, 2014. *User Manual FINE/Marine v3*. Numeca International, Brussels, Belgium.