1   Introduction

With the recent advancements in High Performance Computing technology and more active share of knowledge on modeling practices, application of computational fluid dynamics (CFD) become more available and popular in many branches of engineering. In offshore Engineering, CFD-based computer simulation has been ready to replace physical model tests such as wind tunnel test for wind load (SNAME 2021; Kim et al. 2019; JIP 2019a; JIP 2020) and towing tank test for current load (Koop et al. 2020; JIP 2019b) and vortex induced motion (Jang et al. 2017; JIP 2019c). However, for the wave basin test for the global performance of the floating structures in extreme sea states, CFD simulation still has challenges to overcome (Wu et al. 2016; Baquet et al. 2019).

Scaled model test in wave basin test involves combined environmental loadings of wind, wave and current. Wind load is usually applied as prescribed force derived from the wind tunnel test or CFD simulation for the superstructure of the floaters. But wave and current are directly generated by wave maker and current flume after careful calibration process to make sure the generated wave and current parameters are within the acceptable range of intended design condition. Other than environmental modeling the wave basin test also require proper modeling of the mooring and riser system.

In many cases the physical limit of the wave basin such as basin size, water depth and capacity of wave maker current flume forces engineers to compromise similarity principles in the scaled model test. In most cases the geometric similarities of hull, and dynamic similarities for wind load, inertial and gravity force (Euler and Froude scaling) from wind, wave and current loads are kept but geometric similarity for the mooring lines and risers and dynamic similarity for the viscous force (Reynolds scaling) are sacrificed considering them less important than similitude in hull geometry and inertial and gravity forces. In deep water applications, mooring lines and risers have to be truncated to match the similitude of static force but not dynamic forces, which resulted in the usage of the model test results only for the calibration of the analytic models but not for the direct prediction of the full-scale response of the floaters.

In CFD simulation, theoretically, there is no physical limitation in modeling geometry and dynamics in full-scale hull and environment. However, one need to consider available computing resources and computing cost when designing the CFD simulation setups. The following setups are usually made to keep the CFD simulation budget withing manageable size: 1) Apply wind force the same as in the physical model test. 2) Minimize the CFD domain as small as possible without causing wave reflection from the boundary. 3) Model full layout of mooring and riser if possible.

To achieve the last two items, wave and current field away from the floater is modelled by potential theory, and mooring lines and risers are modeled by coupled line dynamics model that are commonly used in the time-domain global performance analyses in offshore engineering (Wu et al. 2016). Fully nonlinear potential-based waves are generated in a large domain that covers full water depth and footprints of the mooring and riser system, and waves in the smaller CFD domain is forced by the potential wave along the outer edge volume of the CFD domain.

Several wave models have been used to provide potential-based nonlinear waves. Spectral models such as HOS (Ducrozet et al. 2012; Canard et al. 2020) and HAWASSI (Klopman et al. 2010) models show fast and accurate simulation of nonlinear waves. However, the spectral methods show inefficiency in obtaining wave velocity field at discrete points in CFD Domain and along the mooring and risers. The other line of wave model is variational models where the vertical mode of approximation is independent of the horizontal approximation. IGN (Kim and Ertekin 2000; Webster and Zhao 2018) and TPNWT (Baquet et al. 2017) are such models. These methods are relatively slower than the spectral method and requires finer mesh than the spectral models. But the calculation of the velocity inside the CFD domain is easier and faster.

Because of the inefficiencies in either wave evolution or velocity field calculation, the wave solutions are usually stored before applied to the CFD simulations. ECN developed Grid2Grid (Choi et al. 2017) and Technip developed UDW2p protocol for this purpose (Bouscasse et al. 2021). However, for the long-time simulations such as 3-hour simulation that is required for the extreme sea states to get proper statistics of the extreme values the amount of data to be stored is significant and may delay CFD simulation because of heavy data traffic during the simulation.

Recently Kim et al. (2022) introduced the higher-level spectral (HLS) method that adopts the superior features of each line of ware models. Pseudo-spectral method is used to formulate the evolution of the ware elevation and surface potential following HOS formulations. The governing equation in the Fluid domain is solved by variational method where the vertical profile of the velocity field is given by small numbers of interpolation functions. As a result, the HLS model solves the wave evolution equation as effective as the HOS model and velocity field in an arbitrary location can be efficiently calculated as the variational methods. As a result, the HLS model can run parallelly with the CFD simulation while providing wave solution in runtime, to make the CFD-based numerical wave tank simulation faster and more economic.

In this paper, the HLS model originally developed for deep water in Kim et al. (2022) is extended to the finite depth of water for the application to the offshore wind turbine development. Also considered is the interaction with the uniform current. Verification of the HLS model and an example of coupling with CFD-based numerical wave tank is presented.

The simulation of CFD-based numerical wave coupled with HLS model shows number of improvements in efficiency and accuracy of the numerical simulation. To name a few 1) No needs for additional storage of wave data; 2) Minimum number of interpolation function for deep and finite-depth of water that leads to efficient calculation of wave kinematics during CFD simulation; 3) Robust wave simulation for the extreme breaking waves

2   Mathematical formulation

We consider the free-surface flow of an inviscid and incompressible fluid in a uniform water depth, h. The coordinate system is chosen such that the z-axis directs against the gravity and the Oxy-plane is the still-water level. We consider long-crested wave where wave propagates in positive x direction. The location of the free surface is denoted by $z=\zeta(x, t) $.

The velocity potential, $ \phi(x, z, t)$, and stream function, $ \psi(x, z, t)$, are introduced to describe the fluid motion. The fluid velocity is defined by

where U_c is current velocity in-line with the wave propagation direction.

The velocity potential and stream function are governed by the Laplace equation in the fluid domain and is subject to appropriate boundary conditions on the boundaries.

The governing equation in the fluid domain is therefore given by

The normal velocity of the fluid vanishes on the sea floor:

In deep water, the velocity of the fluid vanishes on the sea floor:

On the free surface, the kinematic and dynamic conditions are applied. The kinematic boundary conditions can be written in the following three different but equivalent forms:

Dynamic boundary condition has the following equivalent forms:

where $ \hat{\phi} \equiv \phi(x, \zeta, t), \hat{\psi} \equiv \psi(x, \zeta, t)$ and g is the gravity. $n_z=1 / \sqrt{1+\zeta_x^2} $ is the z-component of the normal vector on the free surface. The Equations (5) and (6) are given in terms of the free-surface values of the velocity potential and stream function. In these new expressions all mathematical operations are made on the horizontal axis and time only, to make the numerical operations simple and fast.

The above governing equations and boundary conditions are homogeneous without any external forcing. We add two forcing mechanisms into the free-surface conditions:

1) Mass and momentum source terms to force the known approximate solution for a progressive wave, $ \zeta_w$ and $ \hat{\phi}_w$ in and up- and down-stream part of the free surface. A damping function μ(x) that varies from zero at the inner domain to a constant value at the outer area near the ends of the domain

2) Turbulent viscosity $v_B(x, t) $ due to wave breaking, which is triggered when the wavefront slope is greater than a threshold value. See Baquet et al. (2017) for detail description.

Also introduced is the linearization parameter ϵ (x) that varies from zero at the outer domain to unity at the inner domain. The free surface conditions with all these correction terms are given by

The fluid domain for the governing equations is also modified as

Equations (7) and (8) provides the time evolution of ζ and $ \hat{\phi}$ but not on $ \hat{\psi}$. The stream function on the free surface, $ \hat{\psi}$, can be obtained after solving the Laplace equation given in Eq. (9) for the stream function in the whole domain, ψ, with the boundary condition given by

The Eqs. (5) and (6) governs the time evolution of the two canonical variables ζ and $\hat{\phi} $. The right-hand side of the two equations consist of terms that are either given by known forcing terms or spatial derivatives of the canonical variables, except the free-surface value of the stream function, $ \hat{\psi}(x, t)$, that can be obtained from the solution of the boundary value problem (BVP) governed by the Laplace equation for the stream function $ \psi(x, z, t)$ and the boundary condition Eq. (10). Note that the formulation given in this section is valid regardless of water depth. The effect of water depth is realized through the solution of the BVP given in different water depth. Next section describes how to solve the BVP numerically.

3   Numerical method

Bai and Kim (1995) solved the problem using a finite-element method based on Hamilton's principle. The vertical structures of the stream function and velocity potential are interpolated by analytic functions and finite-element shape functions are used for the horizontal interpolation of the solutions. Similar approaches are used for the IGN and TPNWT models. One of the drawbacks of this method is numerical error in handling the first-order derivative terms in Equations (7) through (10). It is well known that finite-element and finite-difference methods has dispersion error for the advection terms, which are the first-order derivative terms. The most accurate numerical method to handle this advection term is the spectral method, which provides exact derivative terms of any order for the numerical solutions interpolated by Fourier series. The high-order spectral (HOS) method, which is a specialized Pseudo-Spectral Method for the water-wave problem, is a good example of this approach. In HOS, the velocity eigenfunctions that satisfy the Laplace equation and all boundary conditions except the free-surface boundary conditions. The boundary conditions on the free-surface are approximated by the boundary conditions on the mean water level based on Taylor expansion of the velocity potential at the mean-water level. For small to moderate wave steepness, the HOS provides efficient and accurate solution with minimum dispersion error. However, as the wave become steep more Taylor expansion terms are necessary and sometimes solutions are not convergent.

HLS code implements a hybrid numerical method that combines the HOS and FEM. The spatial distribution of the surface variables $\zeta, \hat{\phi} \text { and } \hat{\psi} $ are all represented by the Fourier series as in HOS for accuracy. The governing Equation (9) and boundary condition (10) of the stream function in the whole domain, ψ, are solved by the FEM for the robustness.

In this paper we introduce two different ways of solving the BVP for deep and finite-depth of water to minimize the number of interpolation functions in vertical direction. In deep water we use exponential function, and polynomials in finite depth of water.

3.1   Pseudo-spectral method for canonical variables

The Fourier series representation of the surface elevation, for example, is given by

It is assumed that the surface elevation is periodic in space with the period L. The discrete Fourier coefficients of the surface variables are given by

where the $\zeta_n(t) $ are the values of $ \zeta(x, t)$ at the discrete points in between x = 0 and x = L, or

We denote the discrete Fourier transformation (DFT) operator with the N discrete points as F_N, or

and its inverse $F_N^{-1} $ as

The free-surface conditions (5) and (6) involve operations of x-derivative and quadratic and cubic nonlinear terms to evaluate time derivatives of the canonical variables. The x-derivatives can be easily evaluated in the Fourier space. On the other hand, the nonlinear terms are easier to evaluate in physical space. But the evaluation of the nonlinear terms requires de-aliasing operation for the inverse Fourier transformation before performing multiplications to avoid aliasing error (Orszag 1971; Patterson Jr and Orszag 1971).

The DFT of x-derivatives are given by

where the diagonal matrix D_N is given by

For the quadratic and cubic nonlinear terms, the multiplication operations are performed in physical space after de-aliasing. Since we have cubic nonlinear term, 2N de-aliasing is used for the quadratic and cubic terms:

The inverse Fourier transformation with 2N de-aliasing, $F_{2 N}^{-1}\left[\left\{\tilde{f}_k\right\}_{k=0, 1, \cdots, N-1}\right] $, is defined by

where $ \tilde{f}_k^{\prime}$ is the Fourier coefficients obtained by zero-padding of the original Fourier coefficients $ \tilde{f}_k$, i.e.,

Then the two canonical equations in Eqs. (5) and (6) are given in the discrete Fourier space as

3.2   Finite-Element Method for Stream Function

Since the fluid domain varies because of the surface elevation, $ \zeta(x, t)(\epsilon(x) \zeta(x, t)$ to be precise), it is more convenient to adopt a transformed coordinate system $ (x, \gamma, t)$ such that

where $ \theta(x, t)=\zeta(x, t)+h$ is the thickness of the water layer.

In the transformed domain we can then separate the depth-wise and horizontal variations of the stream function as

where $ \left\{f_m(\gamma)\right\}_{m=1, \cdots, M}$ is the set of orthogonal interpolation functions in the vertical direction that satisfy

where δ_mn is the Kronecker's delta. The eigenfunctions can be obtained from the general non-orthogonal interpolation functions by solving the matrix eigenvalue problem

The velocity fields are given by

in deep water and

in finite water depth.

Then the Laplace equation for the stream function $\psi(x, z, t) $ with the Neumann-type boundary condition (10) can be reduced to the following system of differential equations for $\psi_1(x, t), \psi_2(x, t), \cdots, \psi_M(x, t) $.

Deep water:

Finite depth of water:

where

The Eqs. (33) and (34) can be solved for the stream function coefficients $\psi_m(x, t) $ for given $ \zeta \text { and } \frac{\partial \hat{\phi}}{\partial x}$ that are btained from the time integration of the canonical Equations (18) and (19). The boundary-value problem is solved by the finite-element method with 2N nodes. The discrete values of $ \zeta \text { and } \frac{\partial \hat{\phi}}{\partial x}$ at the 2N nodes are de-aliased inverse DFT:

From the solution of the BVP the free-surface value of the stream function, $\hat{\psi}(x, t) $, is obtained as

such that all terms in the right-hand side of the canonical Equations (22) and (23) are now defined for the time integration to the next step.

4   Numerical results

The HLS deep water model has been verified for the extreme wave in the Gulf of Mexico (Kim et al. 2022), satisfying all qualification criteria proposed from the 'Reproducible Offshore CFD JIP' (Fouques et al. 2021). The JIP the CFD modeling practices of various topics are collected and consolidated by the work group for each topic. The consolidated CFD modeling practices are verified their reproducibility by at least 3 verifiers from each work group. For the numerical wave tank, the qualification criteria of the numerically generated waves are developed as a part of the CFD modeling practices. In this paper, application of the HLS model has been made for an offshore wind turbine site in the Dong-Hae Sea located at the east cost of the Republic of Korea. The extreme wave of 50-yr return period is used for the verification of the HLS model and CFD wave simulation using HLS input wave. The water depth and wave condition for the tested wave is given in Table 1.

For HLS verification 100 realizations of irregular waves are simulated, and the following statistical quantities are compared with available qualification criteria:

• Spectral density S(f)

• Significant wave height H_S

• Spectral peak period T_P

• Skewness λ₃ and kurtosis λ₄ of the surface elevation

• Probability distribution of crest height

• Probability distribution of crest rise-velocity

The verified HLS waves are coupled with the CFD-based numerical wave tank to verify the wave elevation and velocity profile generated in the numerical wave tank.

4.1   HLS model setup

The following parameters are used for the HLS simulations for the two irregular wave cases.

Finite-depth wave model is used considering half wavelength of the peak period wave is longer than the water depth. The vertical interpolation of the stream function is made by the odd-order polynomial functions following (Kim and Ertekin 2000):

The order of polynomial, N_p, is determined based on the convergence of wave statistics and velocity profile, which will be addressed later in this paper. The input wave signal is the linear combination of linear sinusoidal waves of 1 000 components with randomly chosen phases, and amplitudes initially generated from the wave energies within the partitions of JONSWAP wave energy spectrum with equal period spacing. The amplitudes of the input wave components are iteratively corrected to fit the wave spectrum of the simulated nonlinear wave to the target JONSWAP wave spectrum, following the procedure described in Baquet et al. (2017).

The periodicity of the input wave signal is enforced by multiplying a tapering function τ(x) to the linear input wave solution. One third of the computational domain, 500 m at both ends of the domain, are used for the transition area for the linear-nonlinear transition, wave forcing and tapering of the wave solution. The remaining 2 000 m of the domain is the fully nonlinear domain. Total 13 wave gauges separated by 100 m spacing are put in the fully nonlinear domain starting from x =900 m. The wave signal at gauge #7 (WG7, hereafter) located at x =1 500 m is used for the calibration of wave input wave signal. The transition functions and the locations of the wave gauges are depicted in Figure 1.

Figure  1  Nonlinearity ϵ(x), damping μ(x) and tapering τ(x) functions and 13 wave gauge locations in HLS simulations

Download: Full-Size Img

For each irregular wave condition 100 realizations are made for 12 000 s that includes 1 200 s transition time and 10 800 s to obtain wave statistics of the 3-h wave to evaluate the qualification of the HLS solution. The averaged wave spectrum from all realizations at WG7 is used for the wave calibration.

4.2   CFD Model Setup

The size of CFD domain is determined considering typical size of the offshore wind turbine hull dimensions. A two-dimensional CFD domain is made from the sectional area normal to the wave propagation. The following parameters are used for the CFD domain, mesh size and time step:

Figure 2 shows the CFD domain and grid system that are used for the CFD simulation presented in this paper. The grid system and time step used in the simulation has been verified for convergence in Bouscasse et al. (2021).

Figure  2  CFD domain and grid for wave simulation

Download: Full-Size Img

CFD simulation is made by a commercial CFD software Star-CCM+ from Siemens, running in the Lonestar 6 supercomputer in Texas Advanced Computing Center (TACC). For the air-water interface, volume of fluid (VOF) method is used to model air and water phases. Both air and water are considered as incompressible fluid. Detached eddy simulation (DES) model is used for the turbulence model.

On the bottom and lateral boundaries, inlet velocity is specified from the HLS solution. On the top boundary, pressure outlet condition prescribing zero-gauge pressure is used. In addition to the inlet boundary conditions, the coupling between HLS wave and CFD wave solutions are ensured by the Euler-Overlay Method, where the differences between the two solutions are damped out by the spatially-varying damping function, μ(x, z), which varies from zero to μ_max in the overlay zone where the damping is applied. Figure 3 shows the contour plot of the damping function that is used in this study. Typical value of μ_max for irregular wave simulation is 2π/T_p.

Figure  3  Damping function for the HLS-CFD coupling

Download: Full-Size Img

4.3   Summary of HLS wave simulation results

Figure 4 shows the wave spectra from 100 realizations at the calibration wave gauge WG7. The wave spectra at each frequency are averaged to compare with the target JONSWAP wave spectrum. The recommended qualification criterion for the wave spectrum from the reproducible offshore CFD JIP (Fouques et al. 2021) is "Spectral density S(f) at calibration wave gauge should be within 5% of the target spectrum for all points within the following frequency range $ f \in\left[\frac{3}{4 T_P}, \frac{3}{2 T_P}\right]$". The frequency covers from 0.05 Hz to 0.10 Hz. The criterion is well satisfied in the frequency range wider than the required range.

Figure  4  Wave energy spectrum at the calibration wave gauge (Case 1)

Download: Full-Size Img

Figure 5 shows the ensembled probability distribution of the crest height at WG7. The crest height of all wave events from the 100 realizations are collected and sorted to obtain the probability of exceedance (POE) of the crest height. Forristall (2000) distribution from the 2^nd-order theory (Huang and Guo 2017) and Huang's distribution from physical and numerical wave tank data are also compared. The ±5% range from Huang's distribution (Huang and Zhang 2018; Bouscasse et al. 2021), which is suggested as a qualification criterion, is also shown. HLS result shows very good agreement with Huang's distribution.

Figure  5  Probability distribution of wave crest height at the calibration wave gauge, ensembled from 100 realizations

Download: Full-Size Img

Evolution of wave spectrum and wave crest height distribution along the wave gauges are shown in Figure 6 and Figure 7, respectively. The wave spectrum from HLS simulation is indistinguishable from the JONSWAP spectrum except at WG13, where slight shift of peak towards lower frequency is observed. Wave crest distribution is within 5% error bound except for WG3. But wave crest distribution for the probability of exceedance higher than 5×10^-4 of the Huang's distribution is within 5% error bound at all wave gauges.

Figure  6  Averaged wave spectra at different wave gauges (WG)

Download: Full-Size Img

Figure  7  Ensembled probability distribution of wave crest height at different wave gauge (WG)

Download: Full-Size Img

Figure 9 the variation of H_s, T_p, maximum wave height, maximum crest height, skewness and kurtosis along the wave gauges. Mean values of all statistical values are within 3% difference from the target or reference values. The kurtosis values show relatively higher deviation from the reference value 3.1 from the 2^nd-order theory for the narrow-banded wave. The calculated values are around 3.15, which are consistent with the other fully nonlinear wave models shown in (Fouques et al. 2021).

Figure  8  Wave Statistics of the selected wave realization for CFD simulation

Download: Full-Size Img

Figure  9  Variation of wave statistics along wave gauges

Download: Full-Size Img

In summary, the statistical values of the wave numerically simulated by HLS model all satisfy the qualification criteria proposed by the Reproducible Offshore CFD JIP.

The HLS results presented above is from the HLS model with N_p = 4. HLS simulations with higher values of N_p are also performed to confirm the convergence of the results. But N_p higher than 4 showed almost identical statistical results. Figure 10 shows the probability distribution of the wave crest height from N_p = 4 and N_p =6. The wave crest distributions from the two simulations are almost identical.

Figure  10  Wave crest height comparison between HLS simulations with N_p = 4 and N_p= 6

Download: Full-Size Img

4.4   Selection of Wave Event for CFD Simulation

CFD simulation has been made for a selected wave event from the 100 realizations. Specifically, a P90 (probability of non-exceedance = 0.9) extreme value in a 3-h sea state is selected. In terms of probability of exceedance shown in Figure 10, extreme event of probability level of 10^-4 is selected for the CFD simulation. The selected event, indicated by an arrow in Figure 10, belongs to the 6^th realization among the 100 realizations simulated. Figure 10 shows the wave statistics of the realization.

As shown in Figure 10, the wave statistics at the calibration wave gauge, WG7, showed converged results with N_p = 4. However, the wave elevation time history and velocity profile show that higher number of vertical interpolation function is needed. Figure 11 shows the time histories of wave elevation near the P90 extreme wave crest height event from HLS and CFD simulations. N_p = 6 provides converged wave crest height and wave elevation for both HLS and CFD simulations. Figure 12 shows the velocity profiles from N_p = 4, 6 and 8. The HLS simulation shows converged velocity profile at N_p = 6. CFD simulation shows minor difference between the velocity profiles from N_p = 6 and N_p = 8. The difference is mainly near free surface, where the velocity is sensitive to the air-water mixing and difficult to converge. For the air gap calculation, simulation with N_p = 6 seems enough.

Figure  11  Wave elevation near the P90 extreme wave crest height event

Download: Full-Size Img

Figure  12  Velocity profile at the P90 extreme wave crest height event

Download: Full-Size Img

The selected wave event and HLS set up can be used for the CFD simulation with a floating wind turbine for the air gap calculation. The CFD simulation will be published in the near future.

5   Conclusions

A new numerical wave model, HLS, is developed to simulate fully nonlinear deep-water wave. The developed HLS model has been verified following the calibration and qualification procedure proposed by 'Reproducible Offshore CFD JIP'.

The HLS model has 1) deep and shallow water wave modeling; 2) wave-current interaction; 3) coupling with CFD as a run-time library to provide realistic wave kinematics for the design of offshore wind turbines, which are main operating in shallow-to-intermediate depth of water. The HLS model has been used to select input wave condition for an air gap calculation for the Penta Semi, a 15 MW floating wind turbine developed by KRISO. It has been shown that HLS model can efficiently simulate and identify extreme wave event, to provide input wave for CFD simulation.