A fully nonlinear potential flow wave modelling procedure for simulations of offshore sea states with various wave breaking scenarios

An accurate representation of a given sea state is crucial for the study of hydrodynamic loads on offshore structures. It is straightforward to check the quality of the reproduced regular waves in a numerical wave tank (NWT). However, many more parameters need to be considered to ensure the quality of irregular waves. In this paper, a fully non-linear potential flow (FNPF) wave model is used to reproduce irregular sea states with different severity of wave breaking. The numerical model solves the velocity potential from the Laplace equation and the free surface boundary conditions using a finite difference method on a σ-coordinate grid. A comprehensive procedure is introduced to ensure the quality of the reproduced fullscale sea states. The effect of wave spectrum discretisation techniques and breaking wave algorithms are compared for an optimal performance. The evaluation of the simulation results takes into account the kurtosis and the wave crest distribution in addition to the wave spectra. The vertical arrangement of the σ-coordinate grid plays an important role in representing the dispersion relation and a grid optimisation method is explained. The current study provides a working procedure that reproduces high-fidelity irregular sea states with breaking waves in an efficient FNPF NWT.


Introduction
and a working procedure is produced.
where η is the free surface elevation, φ = φ(x, η, t) stands for the velocity potential at the 149 free surface, x = (x, y) represents the coordinates at the horizontal plane and w is the vertical 150 velocity at the free surface.
where h = h(x) is the water depth measured from the still water level to the seabed.
The Laplace equation is discretised using second-order central differences on a σ-coordinate 169 system (Engsig-Karup and Bingham, 2009) and solved using a parallelised geometric multi-170 grid preconditioned conjugated gradient solver provided by Hypre (van der Vorst, 1992).

171
The gradient terms of the free-surface boundary conditions are discretised with the 5th-172 order Hamilton-Jacobi version of the weighted essentially non-oscillatory (WENO) scheme 173 (Jiang and Shu, 1996). The WENO stencil consists of three local essentially non-oscillatory 174 (ENO)-stencils and each is assigned with a smoothness indicators IS (Jiang and Shu, 1996). 175 A large IS indicates a non-smooth solution in a local stencil. The scheme is designed so that 176 the local stencil with the highest smoothness (i.e. smallest IS) is assigned with the largest 177 weight ω i and therefore contributes the most. In this way, the scheme is able to handle large 178 gradients up to shock with good accuracy. The WENO approximation for Φ is a convex The discretisation is formulated as the following: The three stencils are defined as: and The weights are written as with the regularisation parameter e = 10 −6 and the following smoothness indicators: For the temporal discretisation of the time-dependent terms in the free surface boundary 190 conditions, a third-order accurate total variation diminishing (TVD) Runge-Kutta scheme 191 (Shu and Osher, 1988) is used that consists of three Euler steps (Shu and Osher, 1988). The where x is scaled to the length of the relaxation zone. The velocity potential Φ and the 209 surface elevation η are increased to the analytical values in the wave generation zone: Similarly, the velocity potential and the surface elevation η are reduced to still water 211 values in the wave energy dissipation zone or numerical beach to eliminate wave reflection of 212 the outlet boundaries.

213
The irregular wave generation is achieved by a linear superposition of a finite number 214 of individual regular wave components with different amplitudes, frequencies and phases.

215
Assuming that each wave component is a linear wave, the first-order free surface η (1) is 216 defined as where A i is the amplitude of each wave component and θ i is the phase of each component, 218 which is defined as where k i is the wave number, ω i is the angular frequency, ε is the phase of each wave 220 component.

221
A power spectrum is used to describe the energy distribution over the frequencies. A end ω e , a finite number of discrete frequencies (N f ) are needed to sufficiently represent the  The relationship among ω s ,ω p , ω e , N s , N f and N e is summarised in Eqn. 25: Instead of evenly distributing the frequency components, the EEM method discretises the 238 power spectrum so that the energy at each frequency interval is kept constant. As a result, the After the discretisation of the power spectrum, the irregular wave free surface elevation is 249 calculated in Eqn. 26: water.

257
Deepwater wave steepness-induced breaking is initialised with a steepness criterion: where β is the threshold of wave slope at the wave front.

268
When wave breaking is detected, the free surface boundary conditions Eqn. 2 and Eqn. 3 then 269 become: where ν b is the artificial turbulence viscosity.
where α is the stretching factor and i and N z stand for the index of the grid point and 291 the total number of cells in the vertical direction.

292
The stretching approach reduces the computational cost as well as the numerical errors 293 when the same number of vertical cells are used. However, the dispersion relation and phase 294 velocity are sensitive to the stretching arrangement. A methodology for an optimal vertical 295 grid distribution is adopted to ensure an accurate representation of the wave propagation.

296
Since each wave component is a linear wave, the Airy wave theory is used to explain the 297 method. According to the Airy wave theory for infinite water depth, the amplitude of the 298 velocity potential only depends on the vertical coordinate z: where C = A g/k, A is the wave amplitude and k is the wave number. 300 expansion of the normalised function Φ(z) at the free surface ζ up to 3rd order can be expressed 302 as the following: 303 e kz = e kζ + e kζ k∆z + 1 2 e kζ k 2 ∆z 2 + 1 6 e kζ k 3 ∆z 3 + O(∆z 4 ).

304
The truncation error E is then defined as the difference between the original analytical 305 expression and its Taylor expansion up to the 3rd order: It is seen that the truncation error of a finite difference scheme depends on ∆z = z − ζ.

307
Alternatively, a more general form of Eqn. 34 can be expressed for a finite difference scheme 308 of order O as: assumption, based on which the irregular wave theory is derived. In this manuscript, the 327 low-frequency end (ω s ) is chosen so that it truncates off only 0.05% of the total wave energy.

328
For the high-frequency end (ω e ), it is chosen either at a frequency that truncates 99% of the 329 wave energy or at the high-frequency limit that the waves fulfil the linear wave theory. The 330 linear wave theory requirement in deep water is: Where k is the wave number, ω e is the angular frequency, g is gravitational acceleration 332 and A is the wave amplitude. Here, A = H s /2. From eqn. 37, the upper limit of the high-

333
frequency ω e is determined. In the tested cases, this frequency limit corresponds to 2.5ω p to 334 3.5ω p depending on the spectrum and sea state. train. The horizontal cell size is decided as the follows: where L e is the wavelength corresponding to the high frequency limit ω e . Such resolution 342 is used for all the following simulations.
Where C s is phase velocity corresponding to the longest wave with the low frequency limit 350 (ω s ). As a result, the time step is determined by using the phase velocity of the longest wave 351 (lowest frequency ω s ) and the cell size that is determined from the shortest wave (highest 352 frequency ω e ).   Table. 1. Here, γ is the shape parameter of the JONSWAP spectra, d stands 366 for still water level, 0.5H s K p is the characteristic wave steepness and ω is the frequency range 367 used in the simulations.   Table. 2.   With the parallel computation algorithm, the numerical model is able to utilise the com-  As can be seen, the first two wave gauges G1 and G2 tend to show much smaller kur-    is concentrated near the peak frequency, such as the JMB case, the EEM method is seen 550 to produce the spectrum shape well. When a spectrum is more evenly distributed over the 551 frequency, such as the JNB case, the EEM method has a tendency to show higher energy 552 concentration near the peak frequency than theory. The viscous damping algorithm is seen 553 to be more stable in preserving the kurtosis and wave spectrum, while the filtering algorithm 554 produces more varying results. The combined use of the PEM method and the viscous damp-555 ing breaking algorithm seems to give a more universal solution considering kurtosis, spectra 556 shape and wave crest distribution. With increasing severity of wave breaking, more energy 557 loss is seen at the high frequency part of the spectra, and the wave crest distribution has a 558 trend to exceed the upper limit.

559
In conclusion, the methodology described in the manuscript fills the gap of a guideline