Simulating Breaking Focused Waves in CFD: Methodology for Controlled Generation of First and Second Order

The constructive interference at a certain point in space and time of numerous wave components of varying frequencies and amplitudes results in the generation of a large focused wave. For wave components with large enough amplitudes, a wave-breaking event at the location of the focus is produced. When simulating extreme hydrodynamic conditions in numerical tanks, such a breaking wave possesses comparative advantages. It is significantly higher and steeper than any other wave within the propagating wave group, it occurs at a predefined point in space and time, and it is produced by a transient wave group with a short duration instead of, e.g., a long irregular wave sequence. Focused waves have been previously suggested to be suitable candidates for design waves in investigations of wave loading on marine structures (Tromans et al. 1991).

Ning et al. (2008) simulated focused nonbreaking waves using a higher-order boundary element model. Focused wave groups that were smaller in amplitude were generated using linear theory, while for the steepest, but still nonbreaking, waves, a Stokes wave formed by the second-order interaction of the wave components considered was used as input at the source surface at every time step. The actual focal position and time were confirmed to differ from predictions, with increasing differences for increased input wave amplitudes. To achieve the interaction of the steepest nonbreaking wave with a cylinder and compare CFD simulations with experimental results Paulsen et al. (2014) set the position of the focus behind the structure.

In experimental wave tanks, Rapp and Melville (1990), among others, used linear wave theory and Chaplin (1996) proposed an iterative process to calculate the input phases required to bring all wave components into phase at the focus location. Schmittner et al. (2009) extended Chaplin’s method to include amplitude modification as well, and more recently Fernández et al. (2014) suggested a self-correcting method employing a potential flow solver. Although they are effective for small-amplitude waves, the efficiency of these methods is reduced as the nonlinearity of the wave group increases. As a result, the focal point is shifted downstream and the quality of focus considerably deteriorates.

Nevertheless, significant advantages are entailed if the steepest focused wave, and in particular the steepest breaking wave, is accurately generated at the predetermined focusing point and time. By virtue of composite modelling, once the CFD simulations are validated against experimental results in identical conditions, CFD can be applied to test multiple scenarios with a minimum of effort. Unlike physical modelling, which takes significant amounts of labour and time, a numerical wave tank can be used to study the interaction of any structure with the same breaking wave.

This technical note presents a new methodology for the controlled generation of unidirectional focusing wave groups in a CFD model wave tank. The methodology can be used for nonbreaking waves but its application is used in particular for the generation of focused waves that break at the desired phase–focus location and time. First- and second-order input source signals are created with the proposed approach, which is described first. An example of the application is then given and the simulation results are compared with experimental measurements.

Given a desired target spectrum and focusing point and time, the empirical methodology proposed consists of the following steps:

Continuous spectra are obtained by applying a Fourier transform (FT) to each wave elevation time history. The FT and the inverse FT are defined as follows:

Since $F\left(t\right)$ is a real signal, $f\left(\omega \right)$ will be symmetric with an even real part and an odd complex part. Therefore, only positive values of $\omega$ can be used to define the inverse FT, the real part of the Eq. (3) integral from 0 to $+\infty$, multiplied by 2. Writing $f\left(\omega \right)$ in exponential form, we obtain the following:

It is seen that the absolute value of a complex FT [$\left|f(\omega )\right|$] is a function proportional to a spectral density of wave amplitude at frequency $\omega$, while the argument of $f(\omega )$ is a phase shift of the corresponding wave component. In the remainder we refer to $S\left(\omega \right)=\left|f(\omega )\right|$ as the amplitude spectrum and to $\mathit{\Phi }\left(\omega \right)=Arg(\omega )$ as the phase spectrum of the corresponding wave record.

For example, $s_n$ are spectra of fully nonlinear surface elevation signals with phase shifts πn/2, n = 0, 1, 2, 3; the absolute of $S_0$ and $S_{\mathrm{1,2},3}$ is the amplitude spectra of the second-order minus-component (subharmonics) and the nonlinear superharmonics for the first (linear), second, and third orders, respectively. We note in passing that all four surface elevation measurements should refer to the same location.

The methodology described here has two main advantages over all previous methodologies of wave focusing. The linearized part and not the fully nonlinear signal is used, and amplitude Eq. (2) and phase Eq. (3) corrections are applied at different locations in the wave flume. If amplitude corrections are conducted near the inlet, then the target spectrum is allowed to evolve and reshape as the wave group propagates over the numerical domain. A linearized input signal is the natural choice for any inlet employing linear wave theory. Since the full spectrum of a nonlinear wave group is uniquely defined by its linear part, attempting to correct a fully nonlinear spectrum requires correcting parts of the spectrum that are not produced at the inlet.

At least two wave profile time-histories with a constant phase shift of 180° corresponding to crest- and trough-focused wave groups are required to be isolated even from odd harmonics in the measured surface elevation. This technique, known as spectral decomposition, is used for isolating harmonic components corresponding to the Stokes expansion orders, see, e.g., Orszaghova et al. (2014). With this approach, however, the linearized part is not separated from third and higher order terms (odd harmonics), and therefore cannot be used for calculating new input. Combining waves with four phase shifts, namely, crest, trough, positive, and negative slope-focused waves, improves the isolation of the linearized part, and although, e.g., fifth- and higher-order terms are still present, their amplitude is negligible, see, e.g., Buldakov et al. (2017). It is however noted that if breaking occurs before the location of phase correction then neither the spectral decomposition nor the proposed methodology can be effectively applied.

All experiments used for the validation of the numerical simulations are conducted in a 20 m × 1.2 m × 1 m wave flume with a water depth of 0.5 m. The flume is equipped with two wavemakers located at each end of the facility; when waves are generated from one wavemaker the other acts as an active absorber. An Edinburg Design Limited wavemaker resembles but is not a typical piston type wavemaker and therefore differs from a wavemaker of a numerical flume, later to be introduced. The main differences relate to the gaps between the physical wavemaker, the flume’s bed and the side walls, the absence of water behind the numerical wavemaker, and the presence of a wedge-shaped back part in the physical wavemaker, which moves in a direction opposite to the front part and prevents the generation of waves upstream.

The wave flume has been replicated with the numerical model olaFlow, a new development from of the well-known ihFoam model (Higuera et al. 2015). Both technologies are based on OpenFOAM; thus, they solve the Navier–Stokes equations for two incompressible phases (water and air, by means of the volume of fluid technique) using a finite volume discretization. One of the advantages of the new model is that it enables the simulation of fully controllable piston and flap-type wavemakers.

Mesh covers a portion of the experimental flume. Its dimensions are 12.5 m × 0.65 m, enough for the wave to propagate without significant effects from the boundaries. The base mesh resolution is 1 cm × 1 cm for the whole length and up to 0.4 m depth. Above that, the area of interest where the wave propagates and evolves has double that resolution, 2.5 mm × 2.5 mm, to obtain a more detailed flow description. The Courant number is also chosen to be 0.10, which from experience is judged to be low enough to very accurately simulate wave celerities. The mesh totals 475,000 cells. Due to the very small displacements of the piston during the initial instants (on a mm scale), the numerical simulation time is shortened to 22 s, in which the focusing event takes place at t = 20 s. The simulations are performed in a dual Xeon workstation (2.5 GHz); each takes slightly less than 72 h using 10 parallel processes. Depending on the machine, up to four cases, sufficient for the full set of simulations in each iteration, can be run simultaneously.

Surface elevation in the physical and numerical flume is measured with seven resistance-type wave probes and seven numerical probes located 4 m, 5.7 m, 6.9 m, 7.7 m, 8.2 m, 8.45 m, and 8.7 m from the initial position of the wavemaker. The first and last probes are used for amplitude Eq. (5) and phase Eq. (6) corrections, respectively; in the remainder, the probe used for amplitude correction is referred to as the amplitude matching point (AMP), while the probe used for phase corrections is referred as the focus point (FP).

The proposed methodology is used to generate a focused wave breaking 8.7 m from the wavemaker, using a Gaussian target spectrum with peak frequency of f_peak = 0.9 Hz. The depth of the water is set to 0.5 m. Considering a range of target spectra with the same peak frequency, Buldakov et al. (2017) present experimental results, which demonstrate a correlation between the target spectrum and the onset and type of breaking for focusing wave groups; for example, nondimensional elevation time histories for limited nonbreaking and breaking waves created with the Gaussian and JONSWAP target spectra are presented in the Appendix. The former spectrum is, however, chosen for simulations as the shape (wave tails on each side of the main crest) of the focused wave for the latter spectrum (JONSWAP) would require a longer and thereby more computationally expensive numerical domain; a longer domain would be necessary to avoid the coexistence of the focused wave with long wave reflections at the focal point and time.

The focus time is set to 64 s and the repeat period to 128 s entailing a df = 1/128 Hz, where df is the difference between each distinct wave frequency. For this df, the target spectrum used has 256 wave frequencies. The linearized part of the spectrum measured in the first wave probe is matched to the target 4 m away from the wavemaker; thus, this is the location of the AMP. The FP is set 8.7 m from the wavemaker. For the presentation of all results, the focus location is selected as the origin of the coordinate system; therefore, the x-coordinate of the wavemaker is 8.7 m and the wave group propagates in the negative direction toward 0 m. Accordingly, the time reference is also shifted so that the focussing event takes place at t = 0 s.

Fig. 1 presents the numerical surface elevation at AMP and FP and the corresponding amplitude spectrum and phases. For this first step, linear wave theory is used to calculate the phases at the inlet and the motion of the paddle is computed with the target spectrum as input. Four wave groups are generated in total by adding a constant phase shift to the crest-focused case [Figs. 1(a and b)]. The decomposed spectrum for the fully nonlinear signal measured at AMP and FP is shown in Figs. 1(c and d). Near the wavemaker, the linearized part of the measured spectrum (thick line) is smaller than the target (black crosses), especially around the peak frequency, while the wave is clearly out of focus at FP [Figs. 1(b and d)].

The linearized amplitude spectrum at AMP and the phases of the linearized spectrum at FP are then used with the target as inputs to Eqs. (2) and (3). Figs. 2(a and b) present the new, corrected input amplitudes and phases, which are substituted into the analytical solution of Sharma and Dean (1981) to calculate the second order elevation at the inlet and then the motion of the numerical paddle. To the second order, the amplitude of the highest crests in the group increases slightly and the motion of the numerical paddle is adjusted to create the required water level suppression [Figs. 2(c and d)]; the first-order paddle motion, the gray line in Fig. 2(d), can of course be used. It is, however, noted that for second-order generation, the computational time does not increase per se, as the only change derives from the different velocities of the wavemakers, to which the model will react, changing the time step to fulfil Courant number restrictions.

The amplitudes and phases at AMP and FP produced with the corrected paddle displacement are presented in Fig. 3. The measured amplitude spectrum approaches the target [Fig. 3(a)], and the quality of focus improves [Fig. 3(b)]. Unlike the results of the first simulation, the wave now occurs at 0 s, and the main crest at FP is higher. An additional correction is required for the amplitudes and phases to fully converge to the target and zero.

The fully nonlinear surface elevation for the focused breaking wave generated to the first and second order is validated against experimental measurements in Fig. 4. The overall comparison between numerical and physical results is excellent. The highest crests of the experimental waves have slightly increased amplitudes, but at focus the agreement is nearly perfect. Furthermore, the crest of the focused wave was observed to plunge in the experiments, which also agrees with the plunging crest reproduced in the simulations; see Fig. 7.

The effects of second-order generation for numerical waves reflect mainly on the elevation of the waves surrounding the main crests, which are seen to be lower than those created with first-order input. The result of not generating second order would be that the long wave bounded to the wave group would not be well represented; thus, the spurious free waves would be generated, as stated in Sand (1982). In this sense, the numerical model olaFlow already demonstrated capabilities to accurately reproduce second order waves, both for static and dynamic boundary wave generation (Higuera et al. 2013, 2015).

In Fig. 5, the elevation of the linearized and the nonlinear parts is reconstructed from the results of spectral decomposition. Numerical and physical measurements are presented and the measured elevation of the linearized part near the wavemaker (AMP) and at focus (FP) is nearly identical for all three cases. Second-order sum and higher-order components are also seen to be in good agreement in principle. Slightly deeper troughs and higher adjacent crests are observed for experimental measurements, while the use of a second order accurate signal leads to focused higher order waves with troughs symmetric around the main crest; the differences in the amplitude of the numerical crests with and without second order correction were negligible (<3%). The strongest discrepancy reported refers to the second-order difference/long wave components, where a parasitic crest is present in the experimental and numerical results. Nonetheless, the higher-order motion of the numerical paddle produces a deeper suppression of the water level under the wave group (AMP) and the focused wave at FP and prevents the generation of the parasitic crest. For studies of extreme overtopping events using focused wave groups, long error waves have been shown to yield erroneously enhanced run-up distances and overtopping volumes (Orszaghova et al. 2014).

An empirical methodology is described for the generation of breaking focused waves in computational fluid dynamics simulations. The application of the methodology is illustrated for a numerical wave flume created using the CFD model olaFlow. Numerical waves are created with a piston-type wavemaker, and the results are validated against experimental measurements. The overall comparison of the fully nonlinear surface elevation signals measured in the physical and numerical flume is excellent. The signals are decomposed into their linear and higher-order components and a high degree of agreement between experiments and computations is reported.

The proposed methodology is shown to result in the controlled generation of a numerical wave breaking at the predefined location of focus. There is a need to produce wave groups with four phase shifts to apply spectral decomposition, and to iterate puts the experimenter at a disadvantage. For this work, the results converged to the target within three iterations. Practical experience with other target spectra in experimental facilities shows that a fourth iteration is seldom required (Buldakov et al. 2017).

Control is provided over the amplitude spectrum of the focusing wave group, which in the example presented is accurately reproduced 4 m from the wavemaker. The successful use of different probes for amplitude (AMP) and phase (FP) iterations is a key advantage of the suggested methodology, as is the use of a paddle-driving signal with second-order corrections. The former allows for parametric studies, where a wave group is accurately created with a target spectrum and evolves to focus in the wave tank. For second-order generation, computational time does not increase per se and parasitic waves do not pollute the measurements. Once the corrected input is calculated the CFD model can be applied to test multiple scenarios, of, e.g., the breaking wave interacting with different structures, with minimum effort.

Figs. 6 and 7 present additional experimental results and illustrate the type of breaking crest in the simulations, respectively.

The authors are thankful to EPSRC for supporting this project within the Supergen Marine Technology Challenge (Grant EP/J010316/1).