# Numerical Study on the Temporal Discretization Schemes in Two-Phase Wave Simulation

Young Jun Kim, Benjamin Bouscasse, Sopheak Seng, David Le Touzé

► To cite this version:

Young Jun Kim, Benjamin Bouscasse, Sopheak Seng, David Le Touzé. Numerical Study on the Temporal Discretization Schemes in Two-Phase Wave Simulation. International Conference on Offshore Mechanics and Arctic Engineering (OMAE 19), Jun 2019, Glasgow, United Kingdom. 10.1115/OMAE2019-96278. hal-02885112

HAL Id: hal-02885112

https://hal.science/hal-02885112

Submitted on 8 Oct 2020

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 STUDY ON THE TEMPORAL DISCRETIZATION SCHEMES IN
TWO-PHASE WAVE SIMULATION

**Young Jun Kim**

École Centrale de Nantes, France
CNRS UMR 6598
Bureau Veritas, France
Email: young-jun.kim@ec-nantes.fr

**Sopheak Seng**

Bureau Veritas,
Paris La Defense, France
sopheak.seng@bureauveritas.com

**Benjamin Bouscasse**

École Centrale de Nantes, France
CNRS UMR 6598
benjamin.bouscasse@ec-nantes.fr

**David Le Touze**

École Centrale de Nantes, France
CNRS UMR 6598
david.letouze@ec-nantes.fr

## ABSTRACT

The generation and propagation of waves in a viscous flow solver are indispensable part of naval computational fluid dynamic (CFD) applications. This paper presents numerical simulations of two-dimensional wave propagation in the framework of two-phase finite volume method (FVM) with different temporal discretization schemes. Implicit Euler, Crank-Nicolson(CN) and second-order backward temporal discretization schemes are compared by using viscous flow solver based on the open source library OpenFOAM. The combinations of each temporal discretization scheme and explicit limiter are used for the formulation of the Volume Of Fluid (VOF) field convection equation. A new formulation using the second-order backward temporal discretization scheme with explicit limiter are investigated. Two-dimensional periodic domains are considered to compare different time-stepping methods. Also, five different refinement levels of meshes are used to study the convergence properties of each method. The non-linear wave is generated with stream function wave theory using 'foamStar', which is a specialized OpenFOAM library package developed by Bureau Veritas in collaboration with École Centrale de Nantes.

## INTRODUCTION

The number of applications of computational fluid dynamic (CFD) in naval and offshore field has increased significantly for decades. Accurate wave generation and propagation is one of the most important components of naval and offshore CFD applications. However, with a free-surface URANS (Unsteady Reynolds-Averaged Navier-Stokes) solver, it is sometimes difficult to generate, propagate and maintain the desired wave for long simulation runs. Moreover, the order of convergence obtained in space and time in a wave propagation problem often less than the expected theoretical value.

Previous study [1] focused on the order of convergence of the interface treatment and spatial schemes. As a continuous research, the present study aims to analyze the effect of different temporal discretization algorithms. Three type of time-stepping methods are investigated. These are implicit Euler (Euler) scheme, Crank-Nicolson (CN) scheme and second-order backward finite difference (BFD2) scheme. An implementation of these schemes exists already in OpenFOAM. In addition, the study on Euler and CN schemes with MULES (Multidimensional Universal Limiter for Explicit Solver [2]) is already exist [3]. However, there is no formulations for BFD2 in MULES. There-
fore, a new formulation for bounded convection of Volume Of Fluid (VOF) field using BFD2 time-stepping method with explicit limiter are suggested in this paper. All works done in this paper are based on open-source c++ code package OpenFOAM [4] with wave generation library [5].

## Governing equations

The typical flow model for incompressible two-phase Newtonian flow can be derived into the following form of momentum and continuity equations:

$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla\cdot(\rho\mathbf{u}\mathbf{u}) - \nabla\cdot(\mu(\nabla\mathbf{u})) - \nabla\mathbf{u}\cdot\nabla\mu = -\nabla p_{rgh} - (\mathbf{g}\cdot\mathbf{x})\nabla\rho + \sigma_{\kappa}\nabla\alpha \quad (1) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad (2) $$

where $\mathbf{u}$ is the continuous velocity field; $\rho$ is the averaged density field ($\rho_{air} = 1.0 \text{ Kg/m}^2$, $\rho_{water} = 1000 \text{ Kg/m}^2$); $\mu$ is the averaged dynamic viscosity ($\mu_{air} = 1.0 \times 10^{-5} \text{ Ns/m}^2$, $\mu_{air} = 1.0 \times 10^{-3} \text{ Ns/m}^2$); $\mathbf{g}$ is the gravitational field vector ((0, 0, -9.81) m/s²); $\mathbf{x}$ is the position vector from reference to cell center; $p_{rgh}$ is the dynamic pressure. The surface tension is considered negligible ($\sigma_{\kappa} = 0$) in this study. The velocity and pressure coupling is resolved transiently using PISO algorithm [6, 7].

The VOF method with phase indicator function $\alpha$ is used. Thus, $\alpha = 1$ at fully submerged cell and $\alpha = 0$ in the air and the averaged density is computed as $\rho = \alpha\rho_{water} + (1 - \alpha)\rho_{air}$. The VOF convection equation with artificial compression term is:

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

where, $\mathbf{u}_r$ is the relative velocity between air and water at the interface. The third term on the left-hand-side is the compression term which artificially compress the VOF field to minimize the smearing [8, 9].

A careful treatment is required for the discretization of Eq. (3) due to the discontinuous property of the VOF field [1]. Eq. (4) shows the finite volume (FV) representation of Eq. (3).

$$ \frac{d(V_P\alpha)}{dt} + \sum_f [\phi\alpha_f + \phi_r\alpha_f(1-\alpha_f)]_{UB} = 0 \quad (4) \\ F(\phi, \phi_r, \alpha_P, \text{'scheme'}) = [\phi\alpha_f + \phi_r\alpha_f(1-\alpha_f)] $$

<figure><img src="image_2.png" /><figcaption>FIGURE 1. COMPUTATIONAL ALGORITHM FOR VOF CONVECTION</figcaption></figure>

The subscript $P$ is cell-averaged value and subscript $f$ is face-averaged value. The subscripts $BD$ and $UB$ indicate the bounded and unbounded value, respectively. The definition of bounded/unbounded are to be determined from the spatial discretization schemes. The word 'scheme' in Eq.(4) indicates the dependency of the cell-to-face interpolation scheme.

The computational algorithm to solve the Eq.(4) is illustrated in Fig. 1. To solve the VOF field accurately while satisfying the boundedness criteria, two steps are required. These are named *Predictor* step and *Corrector* step. In the *Predictor* step, bounded intermediate VOF field is computed with bounded(upwind) flux, which gives minimum spatial order and accuracy. In the *Corrector* step, the correction flux is obtained by subtracting the bounded flux from the higher-order flux. Then, MULES limit the correction flux and maximize the available flux using boundedness criteria and the information of neighbor cells. Finally, the intermediate VOF and bounded flux is corrected by limited correction flux to increase the order of accuracy. The *Corrector* step has a iteration named *MULES loop* (Fig. 1) which help to find the middle ground between boundedness and order of accuracy. Since the formulations of VOF convection algorithm are differ from the time integration scheme, the *Predictor* step and *Corrector* step for each temporal scheme are given in details.
## MULES with Euler scheme

By applying Implicit Euler temporal discretization, Eq.(4) can be transform into Eq.(5), where h is the current time-step.

$$ \frac{(V_P \alpha_P)^{t+h} - (V_P \alpha_P)^t}{h} + \sum_f [\phi \alpha_f + \phi_r \alpha_f (1 - \alpha_f)]_{UB}^{t+h} = 0 \quad (5) $$

An intermediate VOF-field $\alpha^*$ is introduced and the Eq.(5) is separated into Eq.(6) (Predictor) and Eq.(7) (Corrector).

$$ \frac{(V_P^{t+h} \alpha_P^*) - (V_P \alpha_P)^t}{h} + \sum_f [\phi^{t+h} \alpha_f^*]_{BD} = 0 \quad (6) $$ $$ F_{BD}^{t+h} = [\phi^{t+h} \alpha_f^*]_{BD} $$

Eq.(6) solves intermediate VOF field with bounded (Upwind) discretization scheme. Note that the artificial compression term is not included for the Predictor step.

$$ \frac{(V_P^{t+h} \alpha_P^{**}) - (V_P \alpha_P^*)}{h} + \sum_f \lambda(F_{UB}^{t+h} - F_{BD}^{t+h}) = 0 \quad (7) $$ $$ F_{UB}^{t+h} = [\phi^{t+h} \alpha_f^* + \phi_r^{t+h} \alpha_f^*(1 - \alpha_f^*)]_{UB} $$

Eq.(7) correct the intermediate VOF field with MULES limier $\lambda$. It is solved iteratively and stop at a user-defined number of iterations. The variables $\alpha_P^{**}$, $\alpha_P^*$ and $F_{BD}^{t+h}$ are the iterative variables inside the MULES loop. At the end of iteration, the VOF field update to $\alpha_P^* := \alpha_P^{**}$ and the bounded flux update to Eq.(8).

$$ F_{BD}^{t+h} := F_{BD}^{t+h} + \lambda(F_{UB}^{t+h} - F_{BD}^{t+h}) \quad (8) $$

## MULES with Crank-Nicolson scheme

The formulation of Crank-Nicolson method with MULES limier is discussed in this section. The application of the Crank-Nicolson scheme to Eq.(4) yields Eq.(9).

$$ \frac{(V_P^{t+h} \alpha_P^{t+h}) - (V_P \alpha_P)^t}{h} + \sum_f [C_{CN} F_{UB}^{t+h} + (1 - C_{CN}) F_{UB}^t] = 0 \quad (9) $$ $$ F_{UB}^{t+h} = [\phi \alpha_f + \phi_r \alpha_f (1 - \alpha_f)]_{UB}^{t+h} $$

The coefficient $C_{CN}$ is the Crank-Nicolson coefficient where $C_{CN} = 0.5$ yields the classical Crank-Nicolson scheme and

$C_{CN} = 1.0$ yields the Implicit Euler scheme. Similar to the formulation of Euler time integration scheme, Eq.(9) is separated into Predictor step (Eq.(10)) and Corrector step (Eq.(11)).

$$ \frac{(V_P^{t+h} \alpha_P^*) - (V_P \alpha_P)^t}{h} + \sum_f [\phi^{CN} \alpha_f^*]_{BD} = 0 \quad (10) $$ $$ \phi^{CN} = (C_{CN} \phi^{t+h} + (1 - C_{CN}) \phi^t) $$

The Predictor equation (Eq.(10)) solves the intermediate VOF field $\alpha^*$ with bounded discretization scheme where the flux $\phi^{CN}$ is a time blended flux.

$$ \frac{(V_P^{t+h} \alpha_P^{**}) - (V_P^{t+h} \alpha_P^*)}{h} + \sum_f \lambda(F_{UB}^{CN} - F_{BD}^{CN}) = 0 \quad (11) $$

$$ F_{BD}^{CN} = [\phi^{CN} \alpha_f^*]_{BD} $$ $$ F_{UB}^{CN} = [C_{CN} F_{UB}^{t+h} + (1 - C_{CN}) F_{UB}^t] \quad (12) $$ $$ F_{UB}^{t+h} = [\phi^{t+h} \alpha_f^* + \phi_r^{t+h} \alpha_f^*(1 - \alpha_f^*)] $$ $$ F_{UB}^t = [\phi^t \alpha_f^* + \phi_r^t \alpha_f^*(1 - \alpha_f^*)] $$

The Corrector equations (Eq.(11) and Eq. (12)) correct the VOF field and the VOF flux($F_{BD}^{CN}$) in the same way that the Euler scheme did.

$$ F_{BD}^{CN} := F_{BD}^{CN} + \lambda(F_{UB}^{CN} - F_{BD}^{CN}) \quad (13) $$

Since the corrected VOF flux is time blended value, a time extrapolation from $F_{BD}^{CN}$ to $F^{t+h}$ is required. Eqn.(14) gives the extrapolation of VOF flux where, $C_{oc}$ is off-centering coefficient ($0 \le C_{oc} \le 1$) and $F^t$ is a VOF flux at the previous time step.

$$ F^{t+h} = (1 + C_{oc}) F_{BD}^{CN} - C_{oc} F^t \quad (14) $$ $$ C_{oc} = (1 - C_{CN}) / C_{CN} $$

## MULES with second-order backward scheme

The formulation of second-order backward finite difference (BFD2) time-stepping method with explicit limier is discussed in this section. The application of BFD2 scheme to Eq.(4) yields Eq.(15) where, $h_0$ is the previous time-step.
$$ \frac{1}{h} [c(V_P^{t+h} \alpha_P^{t+h}) - c_0(V_P \alpha_P)^t + c_{00}(V_P \alpha_P)^{t-h_0}] \\ + \sum_f [\phi \alpha_f + \phi_r \alpha_f (1 - \alpha_f)]_{UB}^{t+h} = 0 \quad (15) $$

$$ c = 1 + \frac{h}{h+h_0} \qquad c_{00} = \frac{h^2}{h_0(h+h_0)} \qquad c_0 = c + c_{00} $$

The Eq.(15) divided into Eq.(16) (Predictor) and Eq.(17) (Corrector) using $(1/c)$ scaling.

$$ \frac{(V_P^{t+h} \alpha_P^*) - (V_P \alpha_P)^t}{h} + \sum_f [\phi^{t+h} \alpha_f^*]_{BD} = 0 \quad (16) $$ $$ F_{BD}^{t+h} = [\phi^{t+h} \alpha_f^*]_{BD} $$

The Predictor equation does not solve $\alpha_f^*$ with BFD2 scheme but used Euler like formulation. As a reaction, the remaining higher-order components moves to the Corrector step. Therefore, some source terms $S_n$ are included in Eq.(17).

$$ \frac{(V_P^{t+h} \alpha_P^{**}) - (V_P^{t+h} \alpha_P^*)}{h} + \sum_f \lambda [F_{UB}^n - F_{BD}^{t+h}] = S_n \quad (17) $$ $$ F_{UB}^n = [\phi^{t+h} \alpha_f^* + \phi_r^{t+h} \alpha_f^* (1 - \alpha_f^*)]_{UB} \quad n = 1, 2, 3... $$

$$ S_1 = \frac{1}{c} [(c-1)F_{UB}^* + (c_0-c) \frac{(V_P \alpha_P)^t}{h} - c_{00} \frac{(V_P \alpha_P)^{t-h_0}}{h}] \quad (18) $$ $$ S_n = \frac{1}{c} [(c-1)F_{UB}^n - F_{UB}^{n-1}] \quad n = 2, 3... $$

Similar to different temporal discretization schemes, Eq.(17) is solved multiple times to correct VOF field and VOF flux. The source term $S_n$ must change for every the *MULES loop*. In the Eq.(18), $S_1$ refers to the source of the *MULES loop* and $S_n$ refers to the source terms at $n^{th}$ *MULES* iteration. Every end of *corrector* iteration, the intermediate VOF is updated and Eq.(19) updates the VOF flux.

$$ F_{BD}^{t+h} := F_{BD}^{t+h} + \lambda (F_{UB}^n - F_{BD}^{t+h}) \quad (19) $$

## Simulation set-up

Three different temporal discretization formulation for VOF convection equation are investigated in this section. For convenience, the name of each solver is present in Table 1.

<table><caption>TABLE 1. FOUR SOLVERS WITH DIFFERENT TEMPORAL DISCRETIZATION SCHEME</caption><tr><td>Solver</td><td>Temporal discretization scheme</td></tr><tr><td>Euler</td><td>Implicit Euler</td></tr><tr><td>CN</td><td>Crank-Nicolson</td></tr><tr><td>Backward</td><td>BFD2</td></tr></table>

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

To investigate the ability of each temporal discretization method in naval application, two-dimensional wave propagation test case is considered. The computational domain is exactly one wave length in x-direction. Fig. 2 shows the computational domain, boundary conditions and initial VOF field. A periodic boundary condition is applied to 'Inlet' and 'Outlet' boundaries. For the 'Bottom' patch, slip boundary condition is applied.

To measure the wave in space and in time, 100 wave probes are installed uniformly from 'Inlet' boundary to 'Outlet' boundary. The initial wave velocities and the free-surface positions are evaluated from the stream function wave theory [10, 11]. The wave condition is identical to previous studies [1, 12] and it is tabulated in Table 3.

Same spatial discretization schemes are applied for all solvers and only the temporal discretization scheme for convection and momentum equation is different. Applied spatial discretization schemes are tabulated at Table 2. Also, note that Crank-Nicolson time-stepping method used off-centering coefficient $C_{oc} = 0.95$ due to the stability reason. In order to minimize the iterative uncertainty, 8 outer (SIMPLE) and 2 inner (PISO) loop is applied.
<table><caption>TABLE 2. SPATIAL DISCRETIZATION SCHEMES</caption><thead><tr><th>Name</th><th>Mathematical form</th><th>Spatial discretization scheme</th></tr></thead><tbody><tr><td>Gradient</td><td>∇<sub>a<sub>P</sub></sub></td><td>Gauss linear</td></tr><tr><td>Surface normal gradient</td><td><b>n</b>·(∇<sub>a</sub>)<sub>f</sub></td><td>Linear uncorrected (Structure grid)</td></tr><tr><td>VOF convection</td><td>∇·(<b>u</b>α)</td><td>Gauss Vanleer</td></tr><tr><td>VOF compression</td><td>∇·(<b>u</b><sub>r</sub>α(1 − α))</td><td>Gauss linear</td></tr><tr><td>Momentum convection</td><td>∇·(ρ<b>u</b><b>u</b>)</td><td>Gauss linear upwind for Vector</td></tr><tr><td>Momentum laplacian</td><td>∇·(μ(∇<b>u</b>))</td><td>Gauss linear uncorrected</td></tr><tr><td>Pressure laplacian</td><td>∇·(1/a<sub>P</sub>(∇<b>u</b>))</td><td>Gauss linear uncorrected</td></tr><tr><td>Under-relaxation</td><td>Velocity <b>u</b></td><td>under-relaxation of velocity to 0.7</td></tr></tbody></table>

<table><caption>TABLE 3. WAVE CONDITION</caption><thead><tr><th>Item</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Depth (D)</td><td>[m]</td><td>0.6</td></tr><tr><td>Period (T)</td><td>[s]</td><td>0.7018</td></tr><tr><td>Height (H)</td><td>[m]</td><td>0.05753</td></tr><tr><td>Wave length (λ)</td><td>[m]</td><td>0.8082</td></tr><tr><td>Wave steepness (kA)</td><td></td><td>0.24</td></tr></tbody></table>

<table><caption>TABLE 4. GRIDS FOR THE CONVERGENCE STUDY</caption><thead><tr><th>Case</th><th>λ/Δx</th><th>H/Δz</th><th>T/Δt</th><th>Δx/Δz</th></tr></thead><tbody><tr><td>Grid 1</td><td>15</td><td>3</td><td>50</td><td>2.8</td></tr><tr><td>Grid 2</td><td>25</td><td>5</td><td>100</td><td>2.8</td></tr><tr><td>Grid 3</td><td>50</td><td>10</td><td>200</td><td>2.8</td></tr><tr><td>Grid 4</td><td>100</td><td>20</td><td>400</td><td>2.8</td></tr><tr><td>Grid 5</td><td>200</td><td>40</td><td>800</td><td>2.8</td></tr></tbody></table>

To evaluate the order of convergence of each solver, a set of simulation with five different refinement level of grids is carried out. The information on the grid size and time-step (Δt) are given in Table 4. The refinement ratio between each grid is r = 2 (except between Grid 1 and Grid 2). The reference mesh is Grid 4 which have 100 cell for one wave length and this is a usual mesh resolution for the naval application. Regardless to the refinement level of grid, similar initial Courant Number (Co# ≈ 0.17) is assumed.

## Results and discussions

This section provides the results and discussions on the convergence study for three solvers (Table 1). To examine the performance of each solver, mainly three items are compared.

1. Time averaged first harmonic wave amplitudes.

2. Time averaged phase velocities.

3. Free surface profile at time 10T, 20T, 30T and 40T.

Since the test case is wave propagation with periodic boundary condition, it is difficult to compare the first harmonic phase difference. Therefore, phase velocity is used to compare the phase difference. The averaged time interval is 10 periods, thus four time averaged wave amplitude and phase velocities are compared. The most important time interval is the first one [1T − 10T]. This is because 10 periods of periodic wave propagation are enough for the practical purpose. Note that, 10 periods of periodic wave propagation is very similar to 10 wave length wave propagation with relaxation zones.

From the time averaged first harmonic amplitude, the order of convergence and its uncertainties are evaluated. This computation is performed with the open access tool developed by Eça and Hoekstra [13]. No convergence analysis has performed for
<figure><img src="image_4.png" /><figcaption>FIGURE 3. TIME AVERAGED AMPLITUDE AND PHASE VELOCITY (EULER)</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>FIGURE 4. FREE SURFACE PROFILE WITH GRID 4 AND GRID 5 (EULER)</figcaption></figure>

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

phase velocity since it showed oscillatory behavior. The value *P* and *U* inside the figures (Fig 5 to Fig 7) represent the order of convergence of wave amplitude and its uncertainty, respectively.

Fig. 3 and Fig. 4 show the results of wave propagation using Implicit Euler time-stepping method. With Grid 1 and Grid 2, the wave is fully damped after 10 periods. Therefore, the order of convergence is evaluated only for the time interval $[1T - 10T]$ and for Grid 3 to Grid 5. Higher phase velocity and large phase shift is captured during the simulation. For time interval $[1T - 10T]$ the phase velocity showed convergence, however, it is still higher than analytic phase velocity.

The wave propagation results using Crank-Nicolson temporal discretization scheme are illustrated at Fig. 5 and Fig. 6. The second-order convergence ($P = 2$) of wave amplitude is achieved for time interval $[1T - 10T]$ and $[11T - 20T]$. The order of convergence decrease to ($P = 1.5$) and ($P = 1.2$) for the time interval $[21T - 30T]$ and $[31T - 40T]$. With Grid 5, 'CN' solver loses about 1% of wave amplitude for every 10 periods, while Grid 4 loses 2.5% of wave amplitude. The phase velocity clearly shows the convergence to analytic solution. It might be useful to

note that there has a gradual increase of phase velocity respect to time.

Fig. 7, Fig. 8 give the results of wave propagation using BFD2 time integration scheme. With Backward solver, rather low accuracy is achieved with course mesh. Also, the order of convergence less then 2 is obtained even if it used second-order time integration scheme.

Table 5 gives the summary of order of convergence of wave amplitude of each solver respect to time intervals. For time interval $[1T - 10T]$, only CN solver obtained second order convergence. And CN solver get the highest order of convergence for later time interval. From this it is be possible to consider that the Crank-Nicolson time-stepping scheme is more stable than others.

Fig. 9, Fig. 10 illustrate the comparison of different solvers for the time interval $[1T - 10T]$ and $[11T - 20T]$. For course mesh (Grid 1 - Grid 3), CN solver shows better conservation of the wave amplitude. However, for the simulation set-up Grid 4 and Grid 5, the difference between CN solver and Backward solver is very small.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 5. TIME AVERAGED AMPLITUDE AND PHASE VELOCITY (CN)</figcaption></figure>

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

<figure><img src="image_11.png" /><figcaption>FIGURE 6. FREE SURFACE PROFILE WITH GIRD 4 AND GRID 5 (CN)</figcaption></figure>

<table><caption>TABLE 5. SUMMARY OF ORDER OF CONVERGENCE OF WAVE AMPLITUDE</caption><thead><tr><th></th><th>T 01-10</th><th>T 11-20</th><th>T 21-30</th><th>T 31-40</th></tr></thead><tbody><tr><th>Euler</th><td>0.7</td><td>x</td><td>x</td><td>x</td></tr><tr><th>CN</th><td>2.0</td><td>2.0</td><td>1.5</td><td>1.2</td></tr><tr><th>Backward</th><td>1.2</td><td>0.5</td><td>0.8</td><td>0.6</td></tr></tbody></table>

## Conclusion

Different temporal discretization schemes for simple naval application are presented. Implicit Euler, Crank-Nicolson and BFD2 temporal discretization scheme are applied. A new formulation for BFD2 time-stepping with explicit limiter are suggested. Three solvers are compared with wave propagation simulation which have period boundary conditions.

Implicit Euler scheme show very large damp of free surface elevation and showed this method is not appropriate for two-phase naval application. Crank-Nicolson scheme showed best propagation of wave and best order of convergence compare to other solvers. The BFD2 temporal discretization scheme showed inferior preservation of wave amplitude and phase velocity for coarse mesh. However, for the fine meshes and small time-step size, such as Grid 4 and Grid 5, BFD2 scheme showed a good wave propagation properties.

## ACKNOWLEDGMENT

This work has been performed in the framework of the Chaire Hydrodynamique et Structure Marines CENTRALE NANTES - BUREAU VERITAS.
<figure><img src="image_12.png" /><figcaption>FIGURE 7. TIME AVERAGED AMPLITUDE AND PHASE VELOCITY (BACKWARD)</figcaption></figure>

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

<figure><img src="image_14.png" /><figcaption>FIGURE 8. FREE SURFACE PROFILE WITH GIRD 4 AND GRID 5 (BACKWARD)</figcaption></figure>

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

<figure><img src="image_16.png" /><figcaption>FIGURE 9. COMPARISON OF TIME AVERAGED WAVE AMPLITUDE FOR TIME INTERVAL [1T – 10T]</figcaption></figure>

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE 10. COMPARISON OF TIME AVERAGED WAVE AMPLITUDE FOR TIME INTERVAL [11T – 20T]</figcaption></figure>

REFERENCES

[1] Kim, Y. J., Choi, Y. M., Bouscasse, B., Seng, S., and Touzé, D. L., 2018. “Influence of interface treatment in wave propagation in cfd”. *16mes Journes de l'Hydrodynamique*.

[2] Márquez Damián, S., 2013. “An extended mixture model for the simultaneous treatment of short and long scale interfaces”. PhD thesis, Universidad Nacional del Litoral, Santa Fe, Argentina.

[3] Monroy, C., and Seng, S., 2017. “Time-stepping schemes for seakeeping in openfoam”. *VII International Conference on Computational Methods in Marine Engineering*.

[4] Weller, H. G., Tabor, G., Jasak, H., and Fureby, C., 1998. “A tensorial approach to computational continuum mechanics using object-oriented techniques”. *Computers in Physics, 12*(6), pp. 620–631.

[5] Jacobsen, N. G., Fuhrman, D. R., and Fredse, J. “A wave generation toolbox for the open-source cfd library: Openfoam”. *International Journal for Numerical Methods in Fluids, 70*(9), pp. 1073–1088.

[6] Issa, R., 1986. “Solution of the implicitly discretised fluid flow equations by operator-splitting”. *Journal of Computational Physics, 62*(1), pp. 40–65.

[7] Patankar, S. V., 1980. *Numerical heat transfer and fluid flow*. Series on Computational Methods in Mechanics and Thermal Science. Hemisphere Publishing Corporation (CRC Press, Taylor & Francis Group).

[8] Henric, R., 2002. “Computational fluid dynamics of dispersed two-phase flows at high phase fractions”. PhD thesis, Imperial College of Science, Technology and Medicine.

[9] Berberović, E., van Hinsberg, N. P., Jakirlić, S., Roisman, I. V., and Tropea, C., 2009. “Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution”. *Phys. Rev. E, 79*, Mar, p. 036306.

[10] Rienecker, M. M., and Fenton, J. D., 1981. “A fourier ap-

proximation method for steady water waves”. *Journal of Fluid Mechanics, 104*, pp. 119–137.

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

[12] 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 the nonlinear potential wave solution at the boundaries”. *International Conference on Offshore Mechanics and Arctic Engineering, Volume 2: CFD and FSI*, pp. OMAE2018–78077.

[13] Eça, L., and Hoekstra, M., 2014. “A procedure for the estimation of the numerical uncertainty of cfd calculations based on grid refinement studies”. *Journal of Computational Physics, 262*, pp. 104–130.