An inflow-boundary-based Navier-Stokes wave tank: verification and validation for waves propagating over flat and inclined bottoms

TL;DR

This paper introduces a novel inflow-boundary Navier-Stokes wave tank that accurately models wave generation and propagation over various bottoms, overcoming previous design challenges with high volume preservation and minimal damping.

Contribution

The paper proposes a volume-preserving inflow-boundary wave tank model using a two-phase VOF setup, addressing key issues in wave simulation such as vorticity interaction and boundary modeling.

Findings

Accurately simulates steep and irregular waves

Demonstrates excellent agreement with analytical and experimental data

Successfully models wave transformation over submerged structures

Abstract

Development of mass-source function based numerical wave tank (NWT) algorithms in the Navier-Stokes (NSE) framework is impeded by multiple design issues such as: (a) optimization of a number of variables characterizing the source region, (b) wave-vorticity interactions and (c) a mandatory requirement of modeling the domain on both sides of the wavemaker. In this paper, we circumvent these hurdles by proposing a volume-preserving inflow-boundary based Navier-Stokes wave tank. Wave generation and propagation is modeled in a two-phase PLIC-VOF set-up. Near-exact volume preservation is achieved (at arbitrarily large steepness) using kinematic stretching that is aimed towards balancing the streamwise momentum between points lying above and below the still water level. Numerical damping of steep waves is prevented by using blended third-order and limiter schemes for momentum advection. In addition, a mesh stair-stepping strategy has been adopted for modeling non-Cartesian immersed boundaries on a staggered variable arrangement. The proposed NWT model is rigorously benchmarked against various wave-propagation scenarios. These include the simulation of: (a) monochromatic waves of various steepnesses, (b) monochromatic waves superimposed with free harmonics, (c) irregular waves in deep water and (d) wave transformation occurring over a submerged trapezoidal bar. Excellent agreement with analytical, numerical and experimental data is reported with both validation as well as verification of the proposed NWT model being established.

Figures

Tables

Table 1. Table 1: Characteristics selected for monochromatic wave generation simulations: k ​ h 𝑘 ℎ kh is relative depth; H / λ 𝐻 𝜆 H/\lambda is the steepness; 𝖴𝗋 ( = H ​ λ 2 / h 3 ) annotated 𝖴𝗋 absent 𝐻 superscript 𝜆 2 superscript ℎ 3 \mathsf{Ur}(=H\lambda^{2}/h^{3}) is the Ursell number.

Case	$h\left(m\right)$	$T\left(s\right)$	$H\left(cm\right)$	$\lambda \left(m\right)$	$kh$	$H/\lambda$	$\mathrm{𝖴𝗋}$
SH[27]	$0.40$	$1.25$	$2.50$	$2.055$	$1.22$	$0.012$	$1.65$
SLH[30]	$0.40$	$2.00$	$4.00$	$3.711$	$0.68$	$0.011$	$8.61$
C[13]	$0.30$	$1.50$	$9.00$	$2.421$	$0.78$	$0.037$	$19.5$
${𝐂}_{𝐚}$[10]	$0.30$	$1.50$	$12.0$	$2.476$	$0.76$	$0.048$	$27.2$
${𝐈}_{𝐛}$	$0.75$	$0.70$	$6.00$	$0.807$	$5.84$	$0.074$	$0.09$

Table 2. Table 2: Quantification of regular wave generation in the NWT at t = 20 ​ T 𝑡 20 𝑇 t=20T . Height ( ℋ ​ ℰ ) ℋ ℰ (\mathcal{HE}) , wavelength ( 𝒲 ​ ℰ ) 𝒲 ℰ (\mathcal{WE}) and volume ( 𝒱 ​ ℰ ) 𝒱 ℰ (\mathcal{VE}) errors are shown in terms of percentages (values > 5 % absent percent 5 >5\% are underlined ; -ve values indicate under-prediction of H , λ 𝐻 𝜆 H,\lambda but also volume addition).

	SH	SLH	C	${𝐂}_{𝐚}$	${𝐈}_{𝐛}$
$\mathscr{H}{ℰ}_1$	$-3.64e+00$	$2.18e+00$	$6.06e+00$	$2.56e+00$	$-9.72e+00$
$\mathscr{H}{ℰ}_2$	$2.56e+00$	$-6.08e-01$	$-2.72e-01$	$-3.33e-02$	$-9.30e+00$
$\mathscr{H}{ℰ}_3$	$7.20e-01$	$2.55e+00$	$6.54e+00$	$6.35e+00$	$-7.96e+00$
$\mathscr{H}{ℰ}_4$	$-2.56e+00$	$-6.73e-01$	$2.24e+00$	$1.50e-01$	$-2.45e+00$
$\mathscr{H}{ℰ}_5$	$1.52e+00$	$1.30e+00$	$-1.85e+00$	$-5.33e+00$	$1.67e+00$
$\mathscr{H}{ℰ}_6$	$-2.80e-01$	–	$-1.08e+00$	$-6.24e+00$	$-3.82e+00$
$\mathscr{H}{ℰ}_7$	$6.40e-01$	–	$-6.66e-01$	$-5.38e+00$	$-4.81e+00$
$𝒲{ℰ}_1$	$-1.28e+00$	$8.03e-02$	$2.95e+00$	$4.66e+00$	$3.60e-01$
$𝒲{ℰ}_2$	$4.81e-01$	$4.51e-01$	$1.55e+00$	$3.63e+00$	$3.49e+00$
$𝒲{ℰ}_3$	$-2.11e-01$	$2.26e-01$	$-1.77e+00$	$-1.95e+00$	$1.99e+00$
$𝒲{ℰ}_4$	$-1.14e+00$	$2.90e-01$	$-1.30e+00$	$-7.17e-01$	$1.59e+00$
$𝒲{ℰ}_5$	$2.55e-01$	$-6.12e-01$	$7.97e-02$	$5.67e-01$	$6.83e-03$
$𝒲{ℰ}_6$	$-6.06e-01$	–	$5.22e-01$	$5.75e-01$	$-3.13e+00$
$𝒲{ℰ}_7$	$1.00e+00$	–	$9.21e-02$	$-3.50e-01$	$-3.11e+00$
$𝒱ℰ$	$1.07e-04$	$-5.97e-03$	$-5.41e-03$	$1.81e-02$	$3.91e-03$

Table 3. Table 3: Variables governing individual and group propagation characteristics of waves.

	${\lambda }_V\left(m\right)$	$T\left(s\right)$	$kh$	$C\left(m/s\right)$	$C_G\left(m/s\right)$	$C/C_G$
SH[27]	$2.055$	$1.25$	$1.22$	$1.644$	$1.174$	$1.400$
C[13]	$2.421$	$1.50$	$0.78$	$1.614$	$1.361$	$1.186$

Table 4. Table 4: Average wave height H ¯ ¯ 𝐻 \overline{H} calculated within: 𝖶𝖬 + 4 ​ λ ≤ x ≤ 𝕃 𝖶𝖬 4 𝜆 𝑥 𝕃 \mathsf{WM}+4\lambda\leq x\leq\mathbb{L} for the mesh configurations Δ 1 → 4 subscript Δ → 1 4 \Delta_{1\rightarrow 4} alongwith the medium and fine grid GCI’s for SH and case C waves.

	${\overline{𝑯}}_{{𝚫}_{\mathrm{𝟒}}}$	${\overline{𝑯}}_{{𝚫}_{\mathrm{𝟑}}}$	${\overline{𝑯}}_{{𝚫}_{\mathrm{𝟐}}}$	${\overline{𝑯}}_{{𝚫}_{\mathrm{𝟏}}}$	${\mathrm{GCI}}_{k+1,k+2}$	${\mathrm{GCI}}_{k,k+1}$
SH	$0.02579m$	$0.02442m$	$0.02504m$	$0.02516m$	$0.6997\%$	$0.1280\%$
C	$0.08732m$	$0.08592m$	$0.08567m$	$0.08997m$	$0.4430\%$	$0.0793\%$

Table 5. Table 5: Characteristics of component harmonics of deep-water irregular wave-train [ 12 ] .

Harmonic	$h\left(m\right)$	$H\left(cm\right)$	$T\left(s\right)$	$\delta \left(rad\right)$	$\lambda \left(m\right)$	$H/\lambda$	$\mathrm{𝖴𝗋}$
$𝓝=\mathrm{𝟏}$	$26.0$	$4.00$	$1.6$	$0$	$4.00$	$0.0100$	$3.64e-05$
$𝓝=\mathrm{𝟐}$	$26.0$	$6.00$	$2.4$	$0$	$9.00$	$0.0067$	$2.77e-04$
$𝓝=\mathrm{𝟑}$	$26.0$	$8.00$	$3.2$	$\pi /2$	$16.0$	$0.0050$	$1.17e-03$

Table 6. Table 6: Wave characteristics selected for wave transformation simulations.

Case	$h\left(m\right)$	$T\left(s\right)$	$H\left(cm\right)$	$\lambda \left(m\right)$	$kh$	$H/\lambda$	$\mathrm{𝖴𝗋}$
SH [27]	$0.40$	$1.25$	$2.50$	$2.055$	$1.22$	$0.012$	$1.65$
SLL [27]	$0.40$	$2.00$	$2.00$	$3.699$	$0.68$	$0.005$	$4.28$
SLH [30]	$0.40$	$2.00$	$4.00$	$3.711$	$0.68$	$0.011$	$8.61$

Equations

$$\vec{\nabla }\bullet \vec{V}=0$$

$$\frac{\partial \vec{V}}{\partial t}+\left(\vec{V}\bullet \vec{\nabla }\right)\vec{V}=-\frac{1}{{\rho }^{\ast }}\vec{\nabla }p+\frac{1}{{\rho }^{\ast }}\vec{\nabla }\bullet \left({\mu }^{\ast }\vec{\nabla }\vec{V}\right)+\vec{g}$$

$${\rho }^{\ast }=f{\rho }_w+(1-f){\rho }_a\mathsf{and}{\mu }^{\ast }=f{\mu }_w+(1-f){\mu }_a$$

$$\frac{\partial f}{\partial t}+\vec{\nabla }\bullet \left(\vec{V}f\right)=0$$

$${\mathbb{V}}_{\mathcal{I}}^{\star }={\mathbb{V}}_{\mathcal{I}}^n+\frac{\Delta t}{{\forall }_{\mathcal{I}}}\left\{\mathcal{A}{\mathcal{F}}_{\mathcal{I}}^n+\mathcal{B}{\mathcal{F}}_{\mathcal{I}}^{n-1}+\mathcal{C}{\mathcal{F}}_{\mathcal{I}}^{n-2}\right\}-\frac{\Delta t}{{\rho }_{I-\text{nicefrac}12}^{\ast }}{\mathcal{G}}_{p_I^n}^{\parallel }+\Delta t\left[\text{SRC}\right]$$

$${\mathcal{F}}_{\mathcal{I}}^n=[\text{ADV}{]}^n+[\text{DIF}{]}^n$$

$$\sum _{\text{face}}\frac{{\mathcal{G}}_{p_I^{\prime }}}{{\rho }_{\text{face}}}{A}_{\text{face}}^{\pm }=\frac{1}{\Delta t}\sum _{\text{face}}{\mathbb{V}}_{\mathcal{I}}^{\star }{A}_{\text{face}}^{\pm }$$

$${\mathbb{V}}_{\mathcal{I}}^{n+1}={\mathbb{V}}_{\mathcal{I}}^{\star }-\frac{\Delta t}{{\rho }_{I-\text{nicefrac}12}^{\ast }}{\mathcal{G}}_{p_I^{\prime }}^{\parallel }$$

$$\frac{2\pi T\varphi (\zeta ,\theta )}{{\lambda }^2}=$$

$$+\left({\mathcal{A}}^2{A}_{22}+{\mathcal{A}}^4{A}_{24}\right)\cosh (2k\zeta )\sin (2\theta )-\left({\mathcal{A}}^3{A}_{33}+{\mathcal{A}}^5{A}_{35}\right)\cosh (3k\zeta )\cos (3\theta )$$

$$-{\mathcal{A}}^4{A}_{44}\cosh (4k\zeta )\sin (4\theta )+{\mathcal{A}}^5{A}_{55}\cosh (5k\zeta )\cos (5\theta )\leftarrow \theta \equiv kx-\omega t$$

$${\mathbb{V}}_{\mathcal{I}}^{\star }\to {\mathbb{V}}_{\mathcal{I}}^{\star }-\exp \left(-\alpha {\left(1-\left(\frac{x-x_a}{{\ell }_d}\right)\right)}^R\right){\mathbb{V}}_{\mathcal{I}}^n$$

$${\mathbb{V}}_a^n=\mathcal{S}{{\mathbb{V}}_a^n|}_{\text{FOU}}+(1-\mathcal{S}){{\mathbb{V}}_a^n|}_{\text{QUICK}}$$

$$\xi (t)=-{\xi }_o\cos (\omega t)$$

$$\eta (x,t)=-a\sin \left(k_{p}x-\omega t\right)-a_{2p}\cos 2\left(k_{p}x-\omega t\right)+a_{2f}\cos \left(k_{f}x-2\omega t\right)$$

$$a=\frac{{\xi }_o\tanh \left(k_{p}h\right)}{n_1}$$

$$a_{2p}=\frac{k_p{a}^2}{4}\cdot \frac{\left(2+\cosh \left(2k_{p}h\right)\right)\cosh \left(k_{p}h\right)}{{\sinh }^3\left(k_{p}h\right)}$$

$$a_{2f}=\frac{a^2\coth \left(k_{p}h\right)}{2h}\cdot \left\{\frac{3}{4{\sinh }^2\left(k_{p}h\right)}-\frac{n_1}{2}\right\}\cdot \frac{\tanh \left(k_{f}h\right)}{n_2}$$

$$U(t)=\frac{\mathrm{d}\xi (t)}{\mathrm{d}t}=\omega {\xi }_o\sin (\omega t)$$

$$a=a_{th}\cdot \underset{\text{leakage}}{\underbrace{\left\{1-C\sqrt{\frac{h}{a_{th}}}\right\}}}\cdot \underset{\text{reflections}}{\underbrace{\left(1+{ϵ}_R\cos \left(2k_{o}x+\delta \right)+{ϵ}_R\cos \delta \right)}}$$

$$\eta (x,t)=\frac{1}{2}\sum _{\mathcal{N}=1}^3{H}_{\mathcal{N}}\cos ({\omega }_{\mathcal{N}}t-k_{\mathcal{N}}x+{\delta }_{\mathcal{N}})$$

$$U=\frac{1}{2}\sum _{\mathcal{N}=1}^3{H}_{\mathcal{N}}{\omega }_{\mathcal{N}}{e}^{k_{\mathcal{N}}(y-h)}\cos ({\omega }_{\mathcal{N}}t+{\delta }_{\mathcal{N}});V=\frac{1}{2}\sum _{\mathcal{N}=1}^3{H}_{\mathcal{N}}{\omega }_{\mathcal{N}}{e}^{k_{\mathcal{N}}(y-h)}\sin ({\omega }_{\mathcal{N}}t+{\delta }_{\mathcal{N}})$$

$${\hat{\eta }}^{(n)}(x)=\frac{\omega }{2\pi }{\int }_0^{2\pi /\omega }\eta (x,t)\exp (-\mathrm{i}n\omega t)\mathrm{d}t$$

$$a^{(n)}(x)=2\sqrt{{\left[\mathrm{Re}\left({\hat{\eta }}^{(n)}(x)\right)\right]}^2+{\left[\mathrm{Im}\left({\hat{\eta }}^{(n)}(x)\right)\right]}^2}$$

$$\mathrm{Re}\left({\hat{\eta }}^{(n)}(x)\right)=$$

$$\frac{\Delta t}{3T}\left\{\eta (x,0)+2\sum _{i=1}^{N/2-1}\eta (x,t_{2i})\cos \left(n\omega t_{2i}\right)+4\sum _{i=1}^{N/2}\eta (x,t_{2i-1})\cos \left(n\omega t_{2i-1}\right)+\eta (x,T)\right\}$$

$$\mathrm{Im}\left({\hat{\eta }}^{(n)}(x)\right)=$$

$$\frac{\Delta t}{3T}\left\{2\sum _{i=1}^{N/2-1}\eta (x,t_{2i})\sin \left(n\omega t_{2i}\right)+4\sum _{i=1}^{N/2}\eta (x,t_{2i-1})\sin \left(n\omega t_{2i-1}\right)\right\}$$

Full text

**An inflow-boundary-based Navier-Stokes wave tank: verification and validation for waves propagating over flat and inclined bottoms

** Shaswat Saincher*†* and Jyotirmay Banerjee*†****Corresponding author at: Department of Mechanical Engineering, Sardar Vallabhbhai National Institute of Technology, Surat - 395007, Gujarat, India. Tel: +91 261 220 4145. Fax: +91 261 222 8394. Email: jbaner@med.svnit.ac.in

†Mechanical Engineering Department, S.V. National Institute of Technology, Surat, 395 007, India

[TABLE]

Introduction

Following the advent of computational fluid dynamics (CFD), numerical wave tanks (NWTs) have emerged as a much needed secondary standard to wave flumes and basins. Since their inception, NWT algorithms have been applied towards addressing a variety of challenging problems such as wave-structure-interaction (WSI) [1, 2], wave-breaking-induced vortex dynamics and air-sea interaction [3, 4], fluid-ship-ice interactions [5] and hydrodynamic appraisal of wave energy converters (WECs) [6]. These studies represent state-of-the-art in the field of NWT development and (a majority) have been carried out using two-phase Navier-Stokes equations (NSE); especially for WSI and coastal problems.

Albeit the level of fidelity achieved (in comparison to Fully Non-linear Potential Theory (FNPT) models [7, 1]), wave generation in the NSE framework is challenging. As a matter of fact, the free surface elevation $(\eta)$ and velocity potential $(\phi)$ do not emerge as explicit variables in the NSE. Instead, $\eta$ and $\phi$ need to be implicitly linked to the continuity and momentum equations through free-surface modeling techniques and wave theories. The integration of analytic expressions for $\eta,\phi$ (from wave theory) into the NSE has to be realized by means of “numerical wavemakers” (abbreviated WM from this point onward).

Over the past two decades, several numerical wavemakers have been proposed such as inflow-boundary-based [8, 9, 10], mass-source-based [11, 12, 13], relaxation-zone-based [14], internal-inlet-based [15], moving-boundary-based [16] and momentum-source-based [17] techniques. Of these, the mass-source function and inflow-boundary techniques are (probably) the first and most extensively used NSE-based WMs.

In case of the mass-source WM [11], surface waves are generated through periodic ingestion/ejection of mass from a group of computational cells termed as the “source region” (cf. Figure 1).

Within the source region, zero-divergence condition of the velocity field $(\vec{V})$ is deliberately violated using: $\vec{\nabla}\bullet\vec{V}=s(t)$ where $s(t)\left(\propto\eta(t)\right)$ is strength of the source region. The mass-source WM is advantageous in that specification of $\eta(t)$ is sufficient for prescribing $s(t)$. Further, given that $\int\limits_{T}\eta(t)\cong 0$ for any (arbitrarily non-linear) wave theory, the source WM does not alter water-phase mass over a wave period $(T)$. Having said that, the mass-source WM is nonetheless faced with many design issues:

• •

the mass-source WM has multiple design variables associated, namely, height $(H_{s})$, width $(W_{s})$ and placement $(h_{s})$ (cf. Figure 1); this complicates geometric optimization of the source region [13].

• •

if a wide range of wave heights is to be simulated, the mass-source WM is only suitable in intermediate water $(1.0\leq kh\leq 2.5)$ [13] where $U\propto y$. With increasing relative depth $(kh>2.5)$, $U\propto{\mathrm{e}}^{y}$ which forces one to continually increase $h_{s}$ (cf. Figure 1) to prevent damping of source-injected momentum. Unfortunately,

• –

there is a limit to increasing $h_{s}$ since source height $(H_{s})$ must remain finite.

• –

$h_{s}\uparrow\,;H_{s}\downarrow$ induces wave distortion, especially when $H/\lambda>0.03$ [13].

This limitation of the source-function WM in deep water was resolved by Perić and Abdel-Maksoud [12] (for $kh\approx 40$) by proposing a modified source design which involved decreasing $h_{s}$ and over-designing the source-strength: $s(t)\rightarrow K\cdot s(t)\ni K\gg 1$ (cf. Figure 1). It was demonstrated that said design strategy is effective in eliminating wave distortion during steep wave generation in deep water [12, 13].

• •

owing to jet-like flows emanating from the source region [12], wave-vorticity interactions (WVI) become particularly intense in $kh<1$ which induce loss of source-injected momentum (to viscous effects), wave distortion and (eventual) height reduction [11, 13]. The extent of WVI induced by the source WM in $kh<1$ appears to be governed by the order of streamwise momentum induced under wave crests $(U_{y=h})$ which is in turn a function of the shallow-water non-linearity $(H/h)$.

• •

aforementioned height damping and WVI in $kh<1$ cannot be eliminated through Perić and Abdel-Maksoud’s [12] design strategy. Owing to negligible damping of wave momentum in near-shallow water $\left(\frac{\partial U}{\partial y}\sim 0\leftarrow kh<1\right)$, strong WVI would occur regardless of where the source is placed along the water column. In such a scenario, even if $s(t)$ were over-designed to compensate for height reduction, it would in fact increase WVI and aggravate height damping. In [13], it has been quantitatively demonstrated by the authors that WVI in $kh<1$ can be arrested by stretching the source region to occupy the entire water depth ($H_{s}\cong h$; cf. Figure 1). Then, wave damping (resulting from an increase in source area) could be arrested by mildly over-designing the source strength: $s(t)\rightarrow K\cdot s(t)\ni K>1$.

• •

an obvious shortcoming of the mass-source WM is the requirement of modeling the domain on both sides of the wavemaker (cf. Figure 1) which necessitates wave absorption at both ends of the NWT. It is a general practice that placement of the source region is offset from the center towards the western end of the NWT so as to facilitate a longer “wave simulation region” [13]. Hence, the portion of the NWT westward of the WM serves no research purpose and is a computational liability.

Given the above-mentioned shortcomings, various empirical design strategies have been proposed for the mass-source WM [11, 12, 13]. Despite these advances, the fact that WM design requires significant alterations with changing relative depth $(kh)$ and steepness $(H/\lambda)$ proves detrimental to the robustness of a proposed NWT algorithm. Aim of the present work is to demonstrate that said limitations of the mass-source technique could be overcome using an inflow-boundary-based WM formulation.

In case of the inflow-boundary technique, analytic expressions for velocities $U(t),V(t)$ and elevation $\eta(t)$ are directly specified at a vertical boundary of the NWT (cf. Figure 2). Thus, it is sufficient to model the domain at one end of the WM. Further, unlike the mass-source WM, water volume is influxed (under crests)/effluxed (under troughs) rather than ejected/ingested through the boundary. Hence there is no vortex formation in the WM near-field thereby precluding WVI. Most importantly, given that the specification of $U(t),V(t)$ and $\eta(t)$ is dictated by wave theory [8, 9], the baseline version of the inflow technique (cf. Figure 2) has virtually zero design variables requiring optimization.

However, the baseline inflow WM is susceptible to inducing wave setup [18, 10]. Said susceptibility is attributable to the fact that, with increasing non-linearity of the waves, there emerges a streamwise momentum discrepancy $|U_{CL}|\gg|U_{TL}|$ (subscripts $CL$ and $TL$ represent crest and trough elevations respectively) such that water volume influxed under crests is considerably greater than volume effluxed under troughs. This could be mathematically stated as: $\int\limits_{T}\int\limits_{h+\eta(t)}U(y,t)\,\mathrm{d}y\,\mathrm{d}t=\mathcal{V}_{+}$ where $\mathcal{V}_{+}$ is the water volume added through the inflow boundary over a wave period; $\mathcal{V}_{+}>0$ holds for any (arbitrarily non-linear) wave theory [19, 20]. Over multiple wave generation cycles, the added volume becomes sufficiently large to induce wave setup in the NWT.

To the best of the authors’ knowledge, there haven’t been diligent efforts towards improving volume preservation characteristics of inflow-boundary-based WMs in the literature. In this paper, we propose to improve volume preservation characteristics of the baseline inflow WM using a novel Kinematic Stretching Technique (abbreviated KST from this point onward). Broadly speaking, a KST alters $U(y,t)$ (predicted by wave theory) by replacing $y$ with a “stretched” coordinate $\zeta(\propto\eta(t))$ which effectively stretches/compresses the velocity profile at points (considered along the wave) lying above/below the SWL. The central motivation underlying the development of KSTs is towards improving Airy theory-based estimations of hydrodynamic loads induced by irregular waves on offshore structures [21, 22]. In this regard, several attempts have been made to improve the estimation of $U(y,t)$ above the SWL using KSTs such as: (a) Wheeler stretching [21], (b) Vertical stretching [18], (c) Extrapolation stretching [9], (d) Effective node elevation method [22] and (e) Effective water depth method [22]. It can be verified that each of the aforementioned KSTs introduces: $U_{CL}^{KST}<U_{CL}^{Airy}\,\,;\,\,\left|U_{TL}^{KST}\right|\geq\left|U_{TL}^{Airy}\right|$. These modifications would invariably result in a reduction in ${\mathcal{V}}_{+}$. This indicates that a KST could be employed to design inflow-boundary based WMs with superior volume conservation characteristics. However, an existing KST-based inflow-boundary WM would prove severely limited during generation of strongly non-linear waves which necessitate higher order wave theories (abbreviated HoT from this point onward) for accurate kinematic description. Given that $U_{CL}^{Airy}<U_{CL}^{HoT}$, it can be anticipated that $U_{CL}^{KST}\ll U_{CL}^{HoT}$. Hence, applying existing KSTs to strongly non-linear wave generation would under-predict $U_{CL}$ thereby defeating the very purpose of a HoT and (possibly) leading to height reduction. Another concern is regarding the extent upto which $\mathcal{V}_{+}$ could be minimized for strongly non-linear waves through existing KSTs.

From the above discussion, a need emerges to develop a KST-based inflow boundary WM that: (a) is applicable to non-linear waves governed by HoTs, (b) does not alter $U^{HoT}_{CL}$ and (c) preserves water volume for any given set of wave characteristics $h,H,T$. In the present work, we achieve this by proposing the “modified inflow technique” (cf. Figure 2) as a robust wave generation methodology for NSE-based NWTs [10]. In case of modified inflow, $U(y,t)\,\forall\,y\geq h$ is predicted using a HoT (namely Stokes V) whilst $U(y,t)\,\forall\,y<h$ is predicted using a modified form of Wheeler stretching (KST). Hence, during wave generation, only the troughs are kinematically strengthened leaving the crest velocities (predicted by HoT) unaltered. It is demonstrated that the modified inflow WM can be used to achieve $\mathcal{V}_{+}\cong 0$ for any (arbitrarily non-linear) target wave design. When compared against the mass-source WM, the modified inflow technique is particularly advantageous in that the latter involves only one WM design variable necessitating parametric optimization compared to four in case of the former. Robustness of the proposed NWT model is further improved through the inclusion of blended schemes for momentum advection to prevent numerical damping of steep waves $(H/\lambda>0.03)$. The proposed NWT algorithm is benchmarked against various monochromatic and polychromatic wave generation as well as wave-transformation scenarios. It is demonstrated that the simulations show excellent agreement with analytic predictions as well as existing numerical findings and experimental measurements reported in the literature.

Rest of the paper is structured as follows: the mathematical model of the wave tank is described in section 2, validation against monochromatic wave generation is presented in section 3; uncertainty quantification of monochromatic wave simulations is reported in section 4; benchmarking for polychromatic wave generation over flat bottoms is presented in section 5 whilst benchmarking against wave-transformation over an inclined bottom is reported in section 6. Salient aspects of the work are summarized in section 7.

Numerical wave tank

The present work is aimed towards improving robustness of our existing NWT algorithm [13] through a reduction in number of WM design variables which is accomplished by replacing the mass-source generator with the modified-inflow technique [10]. Mathematical model of the NWT is detailed in the following subsections.

Governing equations and solution methodology

Wave hydrodynamics in the NWT has been modeled using the two-phase Navier-Stokes equations following a “one-fluid” approach [23]. Given that air and water are individually incompressible (within the range of velocities induced by wave motion), the NSE are retained in a “non-conservative” form to avoid generation of unrealistically large velocities at the interface [23]. The NSE would thus be comprised of the continuity and momentum equations which are represented as,

[TABLE]

Here, $U$ and $V$ denote the streamwise $(x)$ and vertical $(y)$ components of velocity $\vec{V}$, $p$ is the pressure, $t$ is time, $\vec{\mathrm{d}A}$ is the area of surface surrounding the control volume $\mathrm{d}\forall$, $\rho^{*}$ and $\mu^{*}$ are the mixture density and viscosity respectively and $g$ is the acceleration due to gravity. The mixture properties $\rho^{*}$ and $\mu^{*}$ are in turn determined using the volume of fluid (VOF) method [24],

[TABLE]

where subscripts $w$ and $a$ denote water and air phases respectively and $f$ is the volume fraction. Transport of $f$ is governed by the pure advection equation,

[TABLE]

The transport Section 2.1 and Eq. 3 have been discretized on a finite-volume staggered grid [23] with the solution advanced explicitly in time. At the beginning of a new $(n+1)$ time level, $\rho^{*}$ and $\mu^{*}$ are determined from $f-$field of the previous $(n)$ time level which is followed by solution of Eq. 3 using the $n^{th}$ level $\vec{V}-$field. Eq. 3 is solved geometrically using Youngs PLIC-VOF technique [24] which is based on the recurrence of three steps; interface reconstruction, interface advection and material redistribution:

• •

in the first step, the air-water interface is reconstructed from the $n^{th}$ level $f-$field which is used to evaluate the interface normal $\vec{\mathfrak{n}}$ on a $3\times 3$ cell stencil $\left(\vec{\mathfrak{n}}=\vec{\nabla}f\right)$. Once $\vec{\mathfrak{n}}$ is obtained, volume conservation is invoked to estimate the placement of the interface within a cell. Here, the interfacial location is directly obtained by means of pre-computed analytical solutions of a set of sixteen possible orientation cases.

• •

in the second step, the reconstructed fluid region is advected across the velocity field $\vec{V}$. Advection is performed in the Eulerian framework through geometric calculation of fluid fluxes across the cell faces. Again, the fluxes are directly obtained from a set of thirty-two pre-calculated analytical solutions. Eq. 3 is advanced in an operator-split manner to eliminate “double-fluxing errors” and to allow a maximum Courant number of $C_{max}\leq 1$ for $f-$transport. Second-order accuracy in fluid advection is retained by alternating the splitting direction every time step.

• •

as a final step, a conservative redistribution algorithm is run after each split to eliminate overshoots $(f>1)$ and undershoots $(f<0)$ in the volume fraction field.

Thus, barring the (iterative) redistribution algorithm, the casewise selection structure ensures extremely fast interface reconstruction and advection computations in the NWT algorithm [13]. The PLIC-VOF solution yields the $f^{n+1}$ field which is followed by solution of the NSE using the velocity projection method through a predictor-corrector approach. The following discrete form is adopted for Section 2.1 during the predictor $(\star)$ step;

[TABLE]

where $\mathbb{V}$ represents any velocity component, $I$ denotes a pressure cell, $\mathcal{I}$ denotes a momentum cell, $\forall_{\mathcal{I}}$ is volume of a momentum cell, $\Delta t$ is time-step size, $\mathcal{G}^{\parallel}_{p^{n}_{I}}$ is the gradient of pressure at $\mathcal{I}$ along $\mathbb{V}$ and $\left[\textsf{SRC}\right]=0\leftarrow\mathbb{V}\equiv U$ whilst $\left[\textsf{SRC}\right]=-g\leftarrow\mathbb{V}\equiv V$. The function $\mathscr{F}$ in Eq. 4 further expands as,

[TABLE]

Here, $[\textsf{ADV}]^{n}\equiv-\sum\limits_{\textsf{face}}\mathbb{V}^{n}_{a}\mathbb{V}^{n}_{\mathfrak{N}}A_{\textsf{face}}^{\pm}$ is momentum advection in discrete form where, $A_{\textsf{face}}^{\pm}$ is signed face area, $\mathbb{V}_{a}$ is the “advected quantity” determined using a blend of first and third order advection schemes [25] (see section 2.4 for details) whilst $\mathbb{V}_{\mathfrak{N}}$ is the “advecting agency” determined using linear interpolation. Further, momentum diffusion in discrete form is given by: $[\textsf{DIF}]^{n}\equiv\dfrac{1}{\rho^{*}_{I-\nicefrac{{1}}{{2}}}}\left\{\sum\limits_{\textsf{face}}\mu^{*}\mathcal{G}^{\parallel}_{\mathbb{V}^{n}_{\mathcal{I}}}A_{\textsf{face}}^{\pm}+\sum\limits_{\textsf{face}}\mu^{*}\mathcal{G}^{\perp}_{\mathbb{V}^{n}_{\mathcal{I}}}A_{\textsf{face}}^{\pm}\right\}$ where $\mathcal{G}^{\parallel}_{\mathbb{V}^{n}_{\mathcal{I}}}$ is the gradient of $\mathbb{V}^{n}_{\mathcal{I}}$ parallel to itself whilst $\mathcal{G}^{\perp}_{\mathbb{V}^{n}_{\mathcal{I}}}$ is the gradient of $\mathbb{V}^{n}_{\mathcal{I}}$ normal to itself. Because momentum cells are back-staggered by half a cell-size, $\mu^{*}$ is directly available at pressure cells whilst evaluating parallel gradients $(\mathcal{G}^{\parallel})$, however, the same has to be obtained through bi-linear interpolation for perpendicular gradients $(\mathcal{G}^{\perp})$. It should be further noted that $\mathcal{G}$ has been determined using the second-order central difference scheme in all cases. In the present work, $\mathbb{V}^{\star}_{\mathcal{I}}$ in Eq. 4 is estimated using the forward Euler method [26]: $\left(\begin{smallmatrix}\mathscr{A}\\ \mathscr{B}\\ \mathscr{C}\end{smallmatrix}\right)\equiv\left(\begin{smallmatrix}1\\ 0\\ 0\end{smallmatrix}\right)$ for monochromatic wave-propagation scenarios (section 3) whilst the same is estimated using the third-order Adams-Bashforth method [26]: $\left(\begin{smallmatrix}\mathscr{A}\\ \mathscr{B}\\ \mathscr{C}\end{smallmatrix}\right)\equiv\left(\begin{smallmatrix}23/12\\ -16/12\\ 5/12\end{smallmatrix}\right)$ for polychromatic wave-propagation (section 5 and section 6).

The predictor $(\star)$ step is followed by a corrector/velocity projection step which involves solution of the following equation of pressure correction;

[TABLE]

where $p^{\prime}_{I}$ is the correction in pressure necessary to make $\mathbb{V}^{\star}_{\mathcal{I}}$ solenoidal subject to a condition that the root-mean-square value of the imbalance between both sides of Eq. 6 is less than $1e-06$. Said imbalance is iteratively reduced using the Gauss-Seidel method. Once the $p^{\prime}_{I}$ field is obtained, the $n+1^{th}$ level velocity field is computed using,

[TABLE]

which, owing to a staggered grid, is directly obtained at $\mathcal{I}$ without any interpolation; this ensures a tight coupling between velocity and pressure [23]. Domain and mesh design strategies adopted for the NWT are presented in the next subsection.

Domain designing and meshing strategy

Computational model of the proposed inflow-boundary based NWT is shown in Figure 3. Given the fact that the NWT would be benchmarked against a variety of wave propagation scenarios, the domain and mesh configurations have been depicted here in a generalized manner with problem specific designs illustrated throughout section 3-section 6. Referring to Figure 3, $\mathbb{L}(\propto\lambda)$ is the length of the wave simulation region, $\ell_{d}(\propto\lambda)$ is length of the east sponge layer (ESL) and $\mathbb{H}(\sim h+2H)$ is height of the domain. For wave-transformation simulations, the ESL has been replaced by a $1:25$ dissipative beach whose numerical design is described in section 6.1. In all cases, the wave simulation region is uniformly divided into $nxm$ cells. To achieve reflection-free wave absorption, the ESL is discretized using a successively stretched grid of $nxr(\sim 50-100)$ cells. For the vertical mesh, the region $0\leq y\leq h$ is divided into $nyd$ cells whilst $h\leq y\leq\mathbb{H}$ is divided into $nyu$ cells. As evident from Figure 3, two different vertical meshing strategies have been adopted.

For monochromatic waves, the vertical mesh is successively stretched directly from $y=h$ [13]. For polychromatic waves, the vertical mesh is uniformly clustered within $\textsf{TL}\leq y\leq\textsf{CL}$ and the cells are successively stretched away from this region111for polychromatic waves, we consider TL and CL as minimum and maximum displacements (respectively) of the free-surface about $y=h$ over $0\leq x\leq\mathbb{L}$. towards the boundaries (cf. Figure 3). Lastly, variables $nx_{\lambda}$ (cells per wavelength) and $ny_{H}$ (cells per wave-height) have been defined to analogize the mesh design process with steepness $(H/\lambda)$ of the incident waves as per guidelines proposed in [13]. The wave generator and absorber designs are described in the next subsection.

Modified-inflow WM and wave absorption

As detailed in section 1, the baseline-inflow WM adds a net volume of water $\mathcal{V}_{+}$ into the NWT domain; for strongly non-linear waves, the added volume induces setup which hampers the fidelity of the simulation222it is obvious that the magnitude of the induced wave setup increases with decreasing domain length.. The issue of volume addition is addressed here by proposing a WM based on a KST similar to Wheeler stretching [21] but within the framework of a higher-order wave theory (namely Stokes V). At fifth order, the following velocity potential $\phi(\zeta,\theta)$ is prescribed at the inflow boundary,

[TABLE]

where, $k$ is the circular wavenumber, ${\mathcal{A}}$ is a topological parameter, $\omega$ is the circular frequency and the coefficients $A_{11}$, $A_{13}$ etc. are weights assigned to component harmonics of the Stokes V wave. Here, the potential $\phi(y,\theta)$ in Section 2.3 has been modified by replacing the vertical coordinate $y$ with a stretched coordinate $\zeta$ such that $\zeta\equiv y\cdot\frac{h+\wp}{h+\eta(t)}$ where $\wp>0\,\,\forall\,\,y<h\,;\,\wp=\eta(t)\,\,\forall\,\,y\geq h$. We term the proposed method as the modified inflow technique. It is noteworthy that $\wp=\eta(t)$ yields baseline inflow [8] (cf. Figure 4 (a)) whilst $\wp=0$ corresponds to Wheeler’s method [21] which kinematically over-designs the troughs but also under-designs the crests. In the case of modified inflow $(\wp>0)$, kinematic predictions of Stokes V theory are preserved above the SWL with kinematic over-design only applied below the SWL (cf. Figure 4(b)). Further, prescription of $\wp$ is flexible and depends on the wave design in question; the added volume $\mathcal{V}_{+}$ can hence be directly controlled through $\wp$ [10]. Most importantly, the proposed modified inflow WM is characterized by only one design variable: $\wp$. In section 3.1.2 it is demonstrated that the modified-inflow technique outperforms both baseline-inflow as well as Wheeler stretching in terms of the fidelity of waves being generated.

In addition to the wave generator, accuracy of NWT simulations is also contingent on reflection-free absorption of incident wave energy at the “far end” opposite to the WM. In the present work, sponge layers are used for efficiently absorbing waves propagating over horizontal bottoms. For waves propagating over inclined bottoms, physical beaches have been modeled as partially submerged/emergent boundaries for wave dissipation [10]. This is done so that wave-transformation simulations (reported in section 6) bear a close resemblance to the experimental conditions of Beji and Battjes [27] against whose measurements the former are validated. To avoid confusion, the numerical model of the dissipative beach is contextually described in section 6.

Incident waves are absorbed in the ESL through introduction of spatially-varying sink terms [13] in $\mathbb{V}^{\star}_{\mathcal{I}}$ over a purposefully coarsened streamwise mesh (cf. Figure 3),

[TABLE]

where $x_{a}(\equiv\mathbb{L})$ marks the beginning of the ESL, $x(>x_{a})$ is measured along the direction of wave propagation, $\ell_{d}$ is the length of the ESL whilst $\alpha$ and $R$ are sponge layer strength and rate of absorption respectively. It should be noted that numerical optimization of ESL design parameters is not attempted here. Instead, $\ell_{d}\geq 4.0\lambda$; $\alpha=10$; $R=0.4-1.0$ are directly adapted from our previous works [13, 10] for which excellent wave absorption performance was demonstrated over a wide range of steepness $(0.004\leq H/\lambda\leq 0.052)$ and relative depths $(0.76\leq kh\leq 6.2)$. The need for and implementation of blended advection schemes in the present NWT framework is elucidated in the next subsection.

Addressing numerical damping of waves

In our previous work [13], we have extensively discussed the mechanics of steep wave generation $(H/\lambda>0.03)$ in context to the mass-source WM being applied to both deep $(kh\approx 6.2)$ and near-shallow $(kh\approx 0.8)$ water conditions. This was motivated by an observation that in the NSE framework, keeping $h$ and $T$ fixed, the generated waves become increasingly susceptible to damping as $H\uparrow$ and also by an observation that, barring few studies (besides our own) [12, 28], this issue (even though encountered in multiple studies) has received limited attention in the literature. Whilst said susceptibility to damping (especially the “rate of height damping” along the length of the NWT) is closely linked to type of numerical WM selected for simulation [13, 28], it is now realized that the former is also controlled by the strategy adopted for numerically modeling wave propagation (which is, for the most part, inviscid) in a viscous framework such as the NSE. The act of numerically approximating the NSE itself introduces several errors which cause the numerical model to exhibit certain traits that are actually absent from the physical system [29]. In case of ocean waves, one such (numerical) artifact that might be mistaken for a physical trait is spurious/numerical diffusion $(\Gamma^{num})$. We justify this assertion through a pilot assessment of the magnitude of momentum advection $[\mathsf{ADV}]$ and diffusion $[\mathsf{DIF}]$ terms in the discrete NSE (cf. Eq. 4) in comparison to gravity $[\mathsf{SRC}]$ which is the dominant restoring force [19].

The evaluation, shown in Figure 5, is done for two different waves: case A [13] and case $\mathbf{C_{a}}$ [10]; the latter is twelve times higher in comparison to the former. The results indicate a shift in balance between $[\mathsf{ADV}]$ and $[\mathsf{DIF}]$ terms such that the problem transforms from a (conventional) advection-diffusion scenario at low steepness to an advection-dominant scenario at large steepness. Hence, the selection of the numerical scheme used for approximating “advected momentum” $\mathbb{V}_{a}^{n}$ (cf. Eq. 5), being inconsequential for case A, is crucial for case $\mathbf{C_{a}}$. This is intimately linked to the fact that $\Gamma^{num}\propto\rho U\Delta x$ considering a FOU-based evaluation of $\mathbb{V}_{a}^{n}$ [29]. Then, for ocean waves it can be proved that $\Gamma^{num}\propto\rho HT^{-1}\mathrm{e}^{H\lambda^{-1}}\Delta x$ using Airy theory333in fact, streamwise velocity at the crest level can be approximated as: $U_{CL}^{Airy}\cong K_{1}HT^{-1}\mathrm{e}^{K_{2}H\lambda^{-1}}$ where $K_{1}\neq K_{2}>\pi\leftarrow kh<\pi$ and $K_{1}=K_{2}=\pi\leftarrow kh\geq\pi$.; numerical diffusion increases rapidly with wave height. In fact, for a FOU-based evaluation of $\mathbb{V}_{a}^{n}$, case $\mathbf{C_{a}}$ would require a twelve times finer mesh compared to case A to restrict $\Gamma^{num}$ within the same order of magnitude; this is computationally unjustifiable.

The trouble with nullifying $\Gamma^{num}$ in two-phase NSE-based NWT simulations is that FOU (which is inherently bounded and non-oscillatory [29]) cannot simply be replaced by “pure” higher order schemes (such as CD, SOU or QUICK); the latter class of algorithms induce dispersive oscillations in $\mathbb{V}_{a}^{n}$, especially at the interface [23]. In fact, in the authors’ experience, even TVD schemes in pure form tend to produce spurious momentum near the interface (leading to wave distortion). In the present work, we circumvent the issue of spurious momentum generation by implementing a blend of higher-order (QUICK) and “limiter” (FOU) schemes to estimate advected momentum [25],

[TABLE]

where $0.5<\mathscr{S}\leq 1$. Said blending strategy has been observed to be extremely effective at restricting numerical damping of steep waves without inducing topological distortion. This is demonstrated in the next section on monochromatic wave generation.

Monochromatic waves: validation

For regular wave generation, both sinusoidal $(\mathsf{Ur}<10)$ as well as trochoidal $(\mathsf{Ur}\gtrapprox 20)$ wave designs have been considered from literature. Three important capabilities of the proposed NWT are tested: (a) (topological) quality of generated waves, (b) minimization of setup [10] and (c) minimization of height error $\mathcal{HE}$ [13] (especially for $H/\lambda>0.03$). The characteristics considered for monochromatic wave generation simulations are listed in Table 1. These include sinusoidal high frequency (SH) waves, sinusoidal low frequency “high” (SLH) waves, steep waves C and $\mathbf{C_{a}}$ and an extreme-steepness design $\mathbf{I_{b}}$. In terms of the Ursell number, case $\mathbf{C_{a}}$ represents the limit of application of Stokes theory [20]. NWT simulations of said wave designs are presented next.

Topology of generated waves: qualitative assessment

In the present section, the monochromatic wave-generation performance of the proposed NWT model is qualitatively assessed against Stokes V theory in terms of the topology of the generated waves by means of spatial and temporal free-surface elevation profiles.

Small steepness waves $(H/\lambda\approx 0.01)$

It can be appreciated that SH and SLH are “borderline” low steepness waves ($H/\lambda\approx 0.01$; cf. Table 1). In this situation, NSE-based wave generation can be treated as a combined advection-diffusion problem (cf. section 2.4).

Thus, pure FOU is retained ($\mathscr{S}=1$; cf. Eq. 10) for momentum advection in both cases. Further, some simulation parameters (common to both wave designs) have been directly selected based on guidelines established in our previous work [13]: $\mathbb{L}=20.0\,m$, $\mathbb{H}=0.44\,m$, $\ell_{d}\sim 5\lambda$; $nx_{\lambda}\approx 84$, $ny_{H}\approx 18$; $R=0.4$. The pressure field $(p)$ in the water phase has been initialized following the hydrostatic law $p=\rho_{w}g(h-y)\,\,\forall\,\,y\leq h$. At this juncture, it is worth recalling that a refined temporal resolution is necessary for NSE-based wave generation at $H/\lambda\lessapprox 0.01$ [13]. It is further anticipated that $\wp_{\textbf{SH}}\neq\wp_{\textbf{SLH}}\leftarrow\mathsf{Ur}_{\textbf{SH}}\ll\mathsf{Ur}_{\textbf{SLH}}$ (cf. Table 1). Hence, optimum values of $\Delta t$ and $\wp$ have been determined through parametric selection which is depicted in Figure 6. It is evident that the topology of (comparatively shorter) SH waves is strongly governed by $\Delta t$ whilst that of (comparatively non-linear) SLH waves is largely governed by $\wp$. Nonetheless, $\Delta t=T/5000$ is observed to yield sufficiently accurate wave topology within $+5\lambda$ from the wavemaker whilst a near-exact agreement with Stokes V theory is achieved for $\Delta t=T/10000$ (cf. $\eta(x)$ profiles in Figure 6). Considering the placement of the toe of the structure from the WM, $\Delta t=T/5000$ is deemed sufficient for the wave-transformation simulations reported later in section 6. Moreover, $\wp_{\textbf{SH}}=0.05$ and $\wp_{\textbf{SLH}}=0.18$ are observed to effectively nullify wave-setup such that $\mathcal{VE}<0.01\%$ is achieved in both cases; interestingly, increasing $\wp$ beyond these values is observed to induce wave set-down through volume subtraction (cf. Figure 6). The small-steepness wave simulations reported here demonstrate that criteria chosen for selecting spatio-temporal resolution based on source-function WM simulations [13] are equally applicable to the inflow-based NWT. Thus, the selection of $nx_{\lambda}$, $ny_{H}$ and $\Delta t$ in a NWT is contingent on the mathematical framework (FNPT, NSE etc.) used for modeling the wave propagation and is rather independent of the wavemaker design.

Large steepness waves $(0.03\leq H/\lambda\leq 0.05)$

The ability of the proposed modified inflow WM to accurately generate steep, trochoidal waves is assessed next by means of two designs C and $\mathbf{C_{a}}$ which are described in Table 1. Nearly identical NWT setups have been considered in both cases [10]: $\mathbb{L}=19.0\,m$, $\mathbb{H}=0.6\,m$, $\ell_{d}\approx 4\lambda$; $nx_{\lambda}\approx 170$, $ny_{H}\approx 6$; $R=1$. Given the steep nature of the target wave designs, time is non-uniformly advanced using the forward Euler method: $T/\Delta t\geq 750\ni C_{max}\leq 0.25$. Further, $\left(\begin{matrix}\wp\\ \mathscr{S}\end{matrix}\right)\equiv\left(\begin{matrix}0.265\\ 1\end{matrix}\right)$ and $\left(\begin{matrix}\wp\\ \mathscr{S}\end{matrix}\right)\equiv\left(\begin{matrix}0.335\\ 0.5\end{matrix}\right)$ are parametrically chosen for cases C and $\mathbf{C_{a}}$ respectively with the pressure initialized using the hydrostatic law.

Results of the NWT simulations are depicted in Figure 7. In Figure 7(a,b), local variation of free surface elevation $(\eta(t))$ is compared for baseline $(\wp=\eta(t);\zeta\equiv y)$, Wheeler stretched $(\wp=0;\zeta\equiv y)$ and modified $(\wp>0\,\,\forall\,\,y<h;\wp=\eta(t)\,\,\forall\,\,y\geq h)$ inflow-based WMs. Although Wheeler stretching over-designs $\left|U_{TL}\right|$, the technique only manages to correct wave-setup by a small amount. In contrast, the proposed “modified inflow” WM clearly outperforms the conventional inflow-boundary based WMs as setup induced in $\eta(t)$ due to volume addition is convincingly nullified in both cases. In addition, volume preservation characteristics of the three inflow-boundary WMs have been quantified in Figure 7(c). Quantification is based on monitoring the percentage change in primary phase (water) volume $(\mathcal{VE})$ occurring within the NWT [13] during the last five wave periods of the simulation. A net volume addition over a wave period $(\mathcal{V}_{+})$ is clearly observable in the baseline and Wheeler formulations whilst modified inflow WM is seen to exactly preserve volume over a wave period. Figure 7(c) also establishes that $\mathcal{V}_{+}\propto\mathsf{Ur}$ (cf. Table 1) which explains the stronger setup induced for case $\mathbf{C_{a}}$ (cf. Figure 7(b)). It naturally follows that “optimum” $\wp$ required for (exactly) balancing $\mathcal{V}_{+}$ would also increase with Ursell number $(\wp\propto\mathsf{Ur})$ which is established from Figure 7.

The ability of the proposed scheme-blending strategy (cf. Eq. 10) to prevent numerical damping of steep waves is exhibited in Figure 7(d) for case $\mathbf{C_{a}}\,(\mathsf{Ur}>27)$. Given that $H/\lambda\approx 0.05$ for case $\mathbf{C_{a}}$, pure FOU based momentum advection $(\mathscr{S}=1)$ leads to considerable height damping $(\mathcal{HE}\approx 15\%)$ as the waves propagate towards the ESL. In such a situation, the results substantiate the computational reasonableness of using a blended scheme $(\mathscr{S}=0.5)$ in lieu of mesh refinement $(nx_{\lambda}\approx 170)$ to control numerical damping. However, it is worth observing from Figure 7(d) that large amounts of QUICK in the blend $(\mathscr{S}=0.2)$ introduces dispersion errors which induce wave distortion; $\mathscr{S}\geq 0.5$ is thus retained for all validation cases considered in section 5-section 6. Here, FOU plays the important role of limiting spurious momentum generated due to a “density discontinuity” existing at the air-water interface [23].

Extreme steepness waves $(H/\lambda>0.05)$

The final wave design considered for benchmarking the proposed NWT model against regular wave generation is the extreme-steepness case $\mathbf{I_{b}}$ in deep water (cf. Table 1). The following NWT setup has been considered for case $\mathbf{I_{b}}$ simulations: $\mathbb{L}=8.0\,m$, $\mathbb{H}=1.0\,m$, $\ell_{d}\approx 6\lambda$; $nx_{\lambda}\approx 177$, $ny_{H}\approx 66$; $R=0.4$. Time is non-uniformly advanced using the forward Euler method: $T/\Delta t\geq 700\ni C_{max}\leq 0.25$. Further, $\left(\begin{matrix}\wp\\ \mathscr{S}\end{matrix}\right)\equiv\left(\begin{matrix}0.032\\ 0.5\end{matrix}\right)$ are parametrically chosen with the pressure initialized using the hydrostatic law. Extreme-steepness wave-generation capabilities of various inflow-boundary-based WMs are depicted in Figure 8. A thorough appreciation of the versatility of the modified-inflow WM can be gained by assessing the $\eta(t)$ profiles reported in Figure 8(a).

Of the three inflow-boundary WMs considered, modified-inflow is observed to attain the closest agreement with Stokes V predictions. The baseline-inflow and Wheeler-stretching-based simulations are observed to suffer from both height damping as well as wave-setup. Interestingly, the height-damping observed in the case of Wheeler stretching $(\wp=0)$ is chiefly attributable to the fact that the technique under-designs streamwise momentum at points lying above the SWL. Hence, a preservation of Stokes V momentum prediction above the SWL combined with momentum over-design below the SWL helps arrest numerical damping of waves (even for $H/\lambda\approx 0.075$) in the case of the modified-inflow-based NWT.

Owing to $\mathsf{Ur}\sim 0.1$, the topology of case $\mathbf{I_{b}}$ is predominantly sinusoidal; this results in a comparatively lower net volume addition at $t=20T$ even in the case of the baseline inflow WM (cf. Figure 8(b)). This substantiates the hypothesis previously put forth in section 3.1.2 that $\mathcal{V}_{+}\propto\mathsf{Ur}$. Quite interestingly, one may observe from Figure 8(b) that the topological symmetry of the wave design gets reflected in the corresponding $\mathcal{VE}(t)$ record over a wave-generation cycle. Owing to extreme steepness, the generated waves necessitate an “optimally blended” advection scheme which is evident from Figure 8(c).

Topology of generated waves: quantitative assessment

Regular wave generation performance of the modified-inflow WM is quantified in the present section.

The quantification is done in terms of percentage errors in height and wavelength of individual wave packets [13] as well as the net percentage change in water volume (over twenty wave-generation cycles) for all five wave designs considered in section 3.1.1-section 3.1.3. It is noticeable that the proposed WM is capable of simulating a target wave-topology with sufficient accuracy. This statement is especially true for the wavelength $(\lambda)$ as $\mathcal{WE}<5\%$ in all cases. This in turn demonstrates that the kinematic over-design introduced below the SWL (cf. Figure 4) does not hamper the ability of the WM to correctly capture amplitude dispersion. Achieving the target wave-height $(H)$ proved to be more challenging, especially for $H/\lambda\gtrapprox 0.05$ (cases $\mathbf{C_{a}}$ and $\mathbf{I_{b}}$). Nonetheless, $\mathcal{HE}\approx 5\%$ for $x>{\mathsf{WM}}+4\lambda$ is certainly acceptable considering the large steepness of the designs considered. It is also evident from Table 2 that the modified-inflow WM exhibits excellent volume preservation characteristics: $\mathcal{VE}<0.02\%$ is sustained for all five designs even after twenty cycles of wave-generation.

In a novel attempt, the propagation characteristics of the generated waves are rigorously benchmarked against Stokes V predictions in the next section.

Graphical verification of the phase and group velocities

In the present section, the phase $(C)$ and group velocities $(C_{G})$ of the waves generated by the modified-inflow WM are validated. Such an assessment is considerable more rigorous in comparison to topological validation (reported in section 3.1 and section 3.2) because the $C,C_{G}$ validation mandates a simultaneous agreement between the wave profile, it’s propagation speed $(C)$ and amplitude dispersion. In the present work, $C$ and $C_{G}$ of regular SH and case C waves have been graphically evaluated by means of “wave-waterfall diagrams” [31] and validated against Stokes V theory. The waterfall diagrams are essentially $x-t$ plots of $\eta(x)$ profiles with shoreward distance (NWT length) represented along the abscissa and discrete time instances (being integer multiples of the wave period $T$) of said $\eta(x)$ profiles represented along the ordinate. Then, the inverse slope of a line passing through the same wave phase across adjacent wave periods represents the phase speed of the waves.

Further, if the waves propagate in deep $(kh\geq\pi\vdash C=2C_{G})$ or shallow water $(kh\leq\pi/10\vdash C=C_{G})$, the slope of the $C_{G}$ line would be respectively twice or equal to the slope of the $C$ line. Hence, the deep and shallow-water limits of wave theory present certain simplifications [31] to the analysis in that the same waterfall plot could be used for graphically representing both $C$ as well as $C_{G}$. In the present evaluation, however, the peculiar choice of wave design(s) precludes citation to such simplifying assumptions of relative depth (cf. Table 3). It is evident from Table 3 that $C$ isn’t an integer multiple of $C_{G}$ for either wave design. Thus, it is not possible to graphically verify $C_{G}$ on a waterfall diagram of a regular wave train in either case. This limitation is evident from Figure 9. Whilst an excellent graphical validation is achieved for $C$ in both cases, the line representing $C_{G}^{-1}$ proves rather immaterial in either case. Given that $C_{G}$ is the propagation speed of the envelope comprising a group of waves [19], it becomes necessary to generate a wave envelope that evolves in a spatio-temporal sense along the shoreward direction.

Once formed, shoreward evolution of such an envelope could then be traced by means of waterfall diagrams similar to those of Figure 9. In order to generate “groups” of regular SH and case C waves, the inflow WM was run on an zero$\rightarrow$up-ramping$\rightarrow$steady$\rightarrow$down-ramping$\rightarrow$zero cycle with a periodicity of $10T$ for twenty wave periods (cf. Figure 10). The up-ramping and down-ramping of the WM strength was achieved using appropriate cosine functions of the form: $\frac{1}{2}\left\{1\pm\cos\left\{\frac{\pi(t-\mathscr{X}T)}{3T}\right\}\right\}$, where $\mathscr{X}$ is an integer constant that was incremented by $10$ between consecutive cycles. Further, given that the beginning as well as end of a “group generation cycle” would invariably involve small steepness wave generation, $\Delta t$ was conservatively set to $T/10000$ in both cases to preserve numerical stability of the small waves [13]. In order to prevent height damping of case C waves, a $50-50$ FOU-QUICK blend $(\mathscr{S}=0.5)$ was adopted for momentum advection. The regular SH and case C wave groups generated in the NWT (corresponding to the WM signals reported in Figure 10) are depicted by means of wave-waterfalls in Figure 11. It can be appreciated that the wave groups are numerically stable and correctly propagate with the group speed $C_{G}$.

Propagation with $C_{G}$ is most strikingly manifested in that the $C_{G}^{-1}$ line consistently follows the highest wave in the center of the group at all times. Thus, the $C_{G}^{-1}$ line intersects the envelope of the wave group at constant phase across consecutive wave periods. However, it should be noted that waves at the leading edge of the group are in a state of perpetual replacement owing to amplitude dispersion. Considering that this replacement occurs with the phase speed, the $C^{-1}$ line correctly intersects the leading edge of the group at constant phase across consecutive wave periods (cf. Figure 11). Thus, the waves generated using the modified-inflow technique corroborate Stokes V predictions in both topological as well as a kinematic sense.

Formation of evanescent modes

The proposed inflow-boundary-based NWT model would be eventually applied to wave transformation occurring over a submerged shoal. Thus, it is important to estimate the minimum separation necessary between the toe of the structure and the WM such that the waves incident on the structure are free from topological distortions induced during generation.

At the WM, the wave train gets distorted due to the emergence of evanescent modes which are essentially “newborn standing waves”. The evanescent waves mark the region across which the velocity field transitions from “wavemaker-induced” to “wave-induced”. A formal mathematical treatment aimed towards analyzing the spatial evolution of evanescent components induced by the modified-inflow technique is well beyond the scope of the present paper. Instead, the formation of evanescent modes at the inflow boundary is assessed graphically through generation of envelopes representing wave propagation during the final five periods $(15T-20T)$ in the NWT. The wave envelopes thus obtained have been plotted, for successive wave periods, in Figure 12 for both regular SH as well as case C waves. In context to Figure 12, the effect of standing evanescent waves is manifested in the wave topology (thus envelope) getting displaced from the SWL [31, 32]. It is evident that the proposed modified-inflow WM does exhibit evanescent waves (which are highly prominent in the case C envelopes $(kh\approx 0.8)$). The existence of evanescent modes in case of the modified-inflow technique is attributed to the fact that, unlike the works of Maguire [31] and Keaney [32], the explicit goal of designing the WM was volume preservation and not elimination of evanescent waves. Nonetheless, given that the evanescent waves dampen exponentially as one moves away from the WM [31], their effect is observed to be significant only within a distance of one wavelength from the wavemaker. This is nonetheless significant when one considers the fact that $\lambda\approx 8h$ for case C waves. Thus, it is advisable to place the structure atleast $+2\lambda$ from the WM, especially if $kh\leq 1$ for the incident waves. Steeper incident waves in $kh<1$ (such as cases C or $\mathbf{C_{a}}$) might necessitate a minimum distance of $+4\lambda$ between the inflow-boundary and the (toe of the) structure.

Monochromatic waves: uncertainty quantification

The proposed PLIC-VOF NWT algorithm has been successfully validated against Stokes V theory for the topological quality as well as kinematics of propagation of the generated waves. Whilst a “validation exercise” ensures that the NWT algorithm correctly models the hydrodynamics of wave propagation, the same doesn’t guarantee whether the simulations achieve a sufficiently high order of error convergence. Evaluation of the order of convergence of an algorithm falls within the purview of “model verification” which can be undertaken by means of a Grid-Convergence-Index (GCI) assessment. More specifically, a validation exercise ensures that the set of PDEs adopted correctly models the physics of the problem whilst a verification exercise ensures that the numerical solution (of the set of PDEs) approaches the exact solution as the spatio-temporal resolution is “infinitely refined”. In the present work, the GCI methodology [33] has been implemented as a code verification exercise for the generation of regular SH and case C waves (cf. Figure 13). The mean wave height $\overline{H}$ of the waves, four wavelengths away from the wavemaker, has been selected as the parameter for GCI-based verification. The following procedure has been adopted for carrying out the GCI-based verification of the NWT algorithm:

four sets of meshes $(\Delta_{1\rightarrow 4})$ have been selected for both regular SH and case C waves such that the cell size is uniform along both shoreward and depthward directions. The number of cells considered within the wave propagation region for the GCI-assessment are reported in Figure 13 where, $\Delta_{4}$ and $\Delta_{1}$ represent the coarsest and finest meshes respectively. The grid has been successively refined by a fixed, integer rate $\mathscr{R}=2$. Further, the GCI methodology is only applicable to uniform meshes and hence the mesh structure presented here is vastly different from that adopted in the validation studies reported in section 3.1.1 and section 3.1.2. Regular SH and case C waves have been simulated on each mesh configuration for twenty wave periods; the corresponding $\eta(x)$ profiles are reported in Figure 13. 2. 2.

at $t=20T$, the average numerical wave height $\overline{H}$ has been evaluated within the region: $\mathsf{WM}+4\lambda\leq x\leq\mathbb{L}$. The $\overline{H}$ values thus obtained are reported in Table 4. 3. 3.

of the meshes $\Delta_{1\rightarrow 4}$, only three are necessary for the GCI analysis [33]. Considering the fact that the order of convergence of the model has to be a real number, the values of $\overline{H}$ selected should be such that they: (a) pertain to three consecutive meshes and (b) do not contain inflexions. Based on this reasoning, $\overline{H}$ values pertaining to the mesh sets $\Delta_{1\rightarrow 3}$ and $\Delta_{2\rightarrow 4}$ have been selected for the GCI analysis of SH and case C waves respectively. 4. 4.

the formal order of convergence of the method is evaluated using the expression [33]: \mathfrak{p}={\ln\left(\frac{\overline{H}_{k+2}-\overline{H}_{k+1}}{\overline{H}_{k+1}-\overline{H}_{k}}\right)}\Big{/}{\ln(\mathscr{R})} where $k=1$ for SH waves and $k=2$ for case C waves. With reference to the $\overline{H}$ values reported in Table 4, one obtains $\mathfrak{p}=2.444$ and $\mathfrak{p}=2.485$ for SH and case C waves respectively. It is thus established that regular wave generation within the NWT is atleast second-order accurate. 5. 5.

the average wave height obtainable on an infinitely refined mesh $\overline{H}_{0}$ is predicted using Richardson extrapolation [33]: $\overline{H}_{0}=\overline{H}_{k}+\frac{\overline{H}_{k}-\overline{H}_{k+1}}{\mathscr{R}^{\mathfrak{p}}-1}$. The extrapolation (depicted graphically in Figure 13(a,b) by means of $\overline{H}-\Delta x$ plots) reveals $\overline{H}_{0}=0.0252\,m$ for SH waves and $\overline{H}_{0}=0.0856\,m$ for case C waves. 6. 6.

the GCI has been evaluated for the medium and fine grids [33]: \mathrm{GCI}_{k+1,k+2}=F_{s}\cdot\left|\frac{\overline{H}_{k+1}-\overline{H}_{k+2}}{\overline{H}_{k+1}}\right|\Big{/}{\left(\mathscr{R}^{\mathfrak{p}}-1\right)}\times 100\% ; \mathrm{GCI}_{k,k+1}=F_{s}\cdot\left|\frac{\overline{H}_{k}-\overline{H}_{k+1}}{\overline{H}_{k}}\right|\Big{/}{\left(\mathscr{R}^{\mathfrak{p}}-1\right)}\times 100\% using a “moderately conservative” safety factor [33]: $F_{s}=1.25$. The medium and fine grid GCI values are reported in Table 4 for both wave designs considered. 7. 7.

the medium and fine grid GCI’s have finally been used to verify that the mesh set $\Delta_{k\rightarrow k+2}$ lies within the asymptotic range of convergence [33]: $\frac{\mathrm{GCI}_{k+1,k+2}}{\mathscr{R}^{\mathfrak{p}}\cdot\mathrm{GCI}_{k,k+1}}\cong 1$. It is found that asymptotic convergence of the set $\Delta_{k\rightarrow k+2}$ indeed holds with $\frac{\mathrm{GCI}_{2,3}}{\mathscr{R}^{\mathfrak{p}}\cdot\mathrm{GCI}_{1,2}}=1.0046$ and $\frac{\mathrm{GCI}_{3,4}}{\mathscr{R}^{\mathfrak{p}}\cdot\mathrm{GCI}_{2,3}}=0.9979$ for SH and case C waves respectively.

This concludes the GCI-based verification of the proposed NWT model. It is demonstrated that $\mathfrak{p}=2.444$ ($+0.732\%$ error in $\overline{H}$) and $\mathfrak{p}=2.485$ ($-4.867\%$ error in $\overline{H}$) for regular SH and case C waves respectively. Being atleast second-order accurate, the generation of regular waves using the modified-inflow technique is indeed of high fidelity. The NWT is evaluated for polychromatic wave generation in the next section.

Polychromatic wave generation over flat bottoms

In the present section, the proposed PLIC-VOF NWT is benchmarked against analytical, experimental and numerical studies reported in the literature for polychromatic wave trains propagating over even (flat) bottoms. Two scenarios have been considered: (a) simulation of free-wave generation during piston-type WM motion in near-shallow water [34, 35] and (b) simulation of a short irregular wave train in deep water $(kh>10)$ [12].

Regular wave train superimposed with free harmonics

The first case under consideration is free-wave generation occurring when a piston-type wave paddle executes sinusoidal motion in near-shallow water $(kh<1)$. This scenario is challenging to accurately simulate due to the added task of capturing free harmonics that are super-imposed on carrier waves. Piston-type motion of a wave paddle in an experimental wave tank can be mathematically prescribed as,

[TABLE]

where, $\omega$ is the angular frequency and $\xi_{o}$ is the amplitude. Then, the free-surface profile $\eta(x,t)$ of waves generated in the flume would be a super-position of (a) the fundamental harmonic, (b) a bound Stokes II harmonic and (c) a free harmonic. Using wavemaker theory (WMT), Madsen [34] developed an analytic (inviscid) solution to this problem that yielded the three harmonic amplitudes. Following Madsen’s solution, $\eta(x,t)$ would be governed by the expression,

[TABLE]

where $a$, $a_{2p}$ and $a_{2f}$ are amplitudes of the primary, second bound and second free harmonics respectively. The harmonic amplitudes are obtained using [34],

[TABLE]

where, $n_{1}\equiv\frac{1}{2}\left\{1+\frac{2k_{p}h}{\sinh\left(2k_{p}h\right)}\right\}$, $n_{2}\equiv\frac{1}{2}\left\{1+\frac{2k_{f}h}{\sinh\left(2k_{f}h\right)}\right\}$ and $k_{p}$, $k_{f}$ are wavenumbers of the primary and second free harmonics respectively governed by the dispersion relation: $(n\omega)^{2}=gk\tanh(kh)$. Solving the dispersion relation for the fundamental yields $k_{p}h\cong 0.47$. Then, the task at hand is a NSE-based simulation of the piston-type WM motion and comparison of the resultant $\eta(x,t)$ profile against Equation 12. The inflow velocity for piston-type WM motion is obtained from,

[TABLE]

It is evident from Equation 16 that $\frac{\partial U}{\partial y}=0$ and $\int\limits_{T}U(t)\,\mathrm{d}t=0$ (for piston-type motion) and, as a consequence, $\mathcal{V}_{+}=0$ from the inflow boundary; baseline inflow is thus retained for wave generation. The following NWT setup (cf. Figure 3) has been considered: $\mathbb{L}=40.0\,m$, $\mathbb{H}=0.6\,m$, $\ell_{d}\sim 5.0\lambda_{p}$; $nx_{\lambda_{p}}\approx 180$, $ny_{H}\approx 10\,(H\equiv 2a)$; $T/\Delta t=5000$, $R=1$; $\wp\equiv\eta(t)$ (baseline inflow) and $\mathscr{S}=1$.

The validation study itself is split into two sub-problems: the first part deals with comparison of $\eta(x)$ profiles with the simulations of [35] whilst the second part involves validating local $\eta(t)$ signals with those measured experimentally by [34]. It is worth noting at this juncture that waves generated in a experimental flume are slightly smaller than those predicted by wavemaker theory owing to leakage of water past the wave paddle [34]. Thus, one is forced to consider separate validation studies because experimental wave generation entails a “loss” in wave-height whilst numerical wave generation is “lossless” with regard to no water leaking past the wave paddle.

For the first sub-problem, the WM parameters are selected as $\xi_{o}=6.1\,cm$ and $T=2.75\,s$ [35]; the resulting $U(t)$ (cf. Equation 16) is directly input to the inflow boundary. Dong and Huang [35] have reported their simulation results in the form of $\eta(x)$ profiles at $t=21.8\,s$. For the sake of validation, their $\eta(x)$ plot was spatially shifted against the present simulation results such that the first troughs matched in both cases. The resulting (superimposed) $\eta(x)$ profiles are shown in Figure 14(a-c); an excellent agreement is observed between the present simulations and those reported in [35]. However, when the aforementioned simulation framework was directly applied to the second sub-problem, it was observed that the numerical waves were $\approx 16\%$ larger than experimental measurements reported in [34]. This happens because for a given $\xi_{o}$, wave generation in the NWT would be “lossless” whilst that in an actual flume would be “lossy”. Hence, “$a$” appearing in Madsen’s theory actually denotes different amplitudes; $a$ in Equation 12 is the actual amplitude of the fundamental whilst Equation 13 yields a theoretical carrier amplitude $a_{th}$ corresponding to lossless generation. Thus, in order to validate Madsen’s experiments, slightly smaller waves need to be generated in the NWT. Madsen proposed the following modified equation for determining actual amplitude $a$ of the first harmonic [34];

[TABLE]

where $a_{th}$ is obtained from Equation 13 and $C(=0.04)$, $\epsilon_{R}(=0.057)$, $\delta(=0.96\,rad)$ are experimentally determined variables [34]. The actual amplitude $a$ is then re-substituted on the LHS of Equation 13 to yield $\xi_{*}$ as a “diminished amplitude” of the wave paddle. It is noteworthy that the present calculations yield $\xi_{*}=5.55\,cm$ as opposed to $\xi_{o}=6.1\,cm$ used by Dong and Huang [35]. Comparison of presently obtained $\eta(t)$ signals with analytic and experimental profiles of Madsen [34] is shown at the bottom in Figure 14. A decent match is obtained at both stations with the PLIC-VOF NWT signals following the analytic curve (Equation 12) more closely than Madsen’s experimental signals. Validation of the NWT against the generation of a deep-water irregular wave-train over a flat bottom is presented in the next section.

Irregular wave train in deep water

The proposed NWT model is validated against a polychromatic deep-water wave propagation scenario conceived by Perić and Abdel-Maksoud [12] which was simulated using their novel “over-designed source” (OVD) technique. Analytically, the free surface elevation of the irregular wave train at various locations in the NWT is given by,

[TABLE]

The characteristics of component harmonics (governed by the deep-water dispersion relation: $\omega=\sqrt{gk}$) of the wave train are listed in Table 5. The following expressions for velocities from Airy theory [19] are used as input to the WM $(x=0)$;

[TABLE]

where a deep-water approximation $(\tanh{k_{\mathscr{N}}h}\rightarrow 1)$ has been invoked. The NWT configuration and a local-clustering-based vertical meshing strategy adopted for the simulation have been reported in Figure 3. The following setup is adopted for simulation: $\mathbb{L}=60.0\,m$, $\mathbb{H}=28.0\,m$, $\ell_{d}=10.0\lambda_{1}$; $nx_{\lambda_{1}}=88$, $ny_{H_{1}}=10$; ${T_{1}}/\Delta t=1000$ (following [12]), $R=1$.

Evidently, the length of the domain has been selected based on the second harmonic $(\lambda_{2})$ to reduce computation time. Further, considering the fact that the component harmonics are low steepness and (very) low Ursell number waves, $\wp\equiv\eta(t)$ (baseline inflow) and $\mathscr{S}=1$ have been retained. The wave train is simulated upto a physical time of $t=42T_{1}$. It is also worth mentioning that the WM strength is gradually increased over a duration of $t=2T_{1}$ using a cosine ramping function to maintain numerical stability during the initial stages of the simulation [13].

Results of irregular wave generation have been validated against Airy theory in terms of both spatial $(\eta(x))$ and temporal $(\eta(t))$ wave profiles. In Figure 15, $\eta(x)$ profiles are depicted at intervals of $t=16.8\,s$ whilst $\eta(t)$ profiles have been recorded at four locations uniformly-spaced every $15.0\,m$ from the WM. The spatio-temporal intervals are selected such that neither of the three component periods $(T_{1},T_{2},T_{3})$ nor component wavelengths $(\lambda_{1},\lambda_{2},\lambda_{3})$ be multiples of $16.8$ or $15.0$ respectively. Then, owing to frequency dispersion, the $\eta(x)$ and $\eta(t)$ profiles of interest would themselves vary amongst the chosen space-time intervals; this makes validation exceptionally challenging. Nonetheless, the simulations show excellent agreement with Airy theory in all cases (cf. Figure 15). It is also worth mentioning that, even with baseline inflow, the net percentage change in water volume at $t=67.2\,s$ was minimal: $\mathcal{VE}\cong-1.38e-02\%$ which is jointly attributable to $H/\lambda\leq 0.01$ and $\mathsf{Ur}\lll 1$ for individual harmonics (cf. Table 5).

Thus, the simulations reported in section 5 demonstrate that the proposed NWT algorithm is capable of generating high-fidelity polychromatic waves over flat bottoms.

Wave transformation over an inclined bottom

The proposed NWT algorithm is finally benchmarked against scenarios that involve a natural transition of monochromatic waves into a polychromatic wave train (comprised of both free as well as bound harmonics) owing to bathymetric variations. To this effect, the problem of wave transformation over a submerged trapezoidal bar has been considered [27]. The problem represents an effective coastal protection strategy in which the structural design is aimed towards transforming incident long waves to transmitted short waves thereby arresting beach erosion. In simulating such scenarios, the challenge lies in that the structure (which generally has a non-Cartesian geometry) needs to be modeled in the NWT and appropriate momentum and pressure boundary conditions need to be prescribed at the structure. Further, the wave topology itself undergoes dramatic changes over the structure. Achieving a simultaneous phase-agreement amongst vastly different wave topologies coexisting in the NWT during wave transformation is challenging.

Description of the problem, corresponding NWT setup and strategy adopted for treating non-Cartesian submerged/emergent boundaries are presented in section 6.1 and section 6.2. This is followed by systematic validation of the NWT simulations against weak as well as strong wave transformation in section 6.3 through section 6.6.

Problem description and NWT setup

A representative sketch of Beji and Battjes’ [27] experimental setup is shown in Figure 16. A servo-controlled piston-type wave paddle was employed for generating both regular as well as JONSWAP-spectrum based wave trains. The waves pass over a submerged trapezoidal bar with a $1:20$ upslope and a $1:10$ downslope. A $2\,m$ long stretch of water over the bar crest acts as a non-dispersive medium for the waves which eventually get dissipated over a beach with a gentle $1:25$ slope. In the experiments, topological changes occurring in the wave train were recorded by an array of seven wave gauges (see Figure 16): WG1 $(x=5.7\,m)$, WG2 $(x=10.5\,m)$, WG3 $(x=12.5\,m)$, WG4 $(x=13.5\,m)$, WG5 $(x=14.5\,m)$, WG6 $(x=15.7\,m)$ and WG7 $(x=17.3\,m)$. A constant water depth of $h=0.4\,m$ was maintained during the experiments [27]. For a rigorous validation of the NWT, three different incident wave designs have been considered which are reported in Table 6. These include: (a) short, sinusoidal high frequency (SH) waves [36] and (b) comparatively longer, sinusoidal low frequency (SL) waves [36]. We further classify the latter category as: (a) SL “low” $\equiv$ SLL waves and (b) SL “high” $\equiv$ SLH waves. The SLH waves were introduced by Huang and Dong [30] as a “steeper version” of SLL waves and were thus not a part of the original experimental paradigm of Beji and Battjes [27]. Further, it is clear that $\mathsf{Ur}_{SH}<\mathsf{Ur}_{SLL}<\mathsf{Ur}_{SLH}$; the relative amplitudes $\left(a^{(n)}/H\right)$ of higher order bound harmonics also increase in that order [37].

Thus, $\mathsf{Ur}$ of the upstream wave train governs the “intensity” of wave transformation over the submerged obstacle [27, 30]. The upslope/downslope of the bar and upslope of the beach represent non-Cartesian immersed boundaries in the NWT. A mesh stair-stepping strategy has been adopted for immersed boundary treatment [10] which is presented in the next subsection.

Numerical treatment of immersed boundaries

Immersed boundaries have been modeled following an “obstacle approach” which is illustrated in Figure 17. Immersed boundary treatment involves the following steps:

• •

pressure/volume-fraction $(p,f)$ $C\forall$s falling “inside” the bar/beach boundary $(\Omega)$ are flagged; the flagged cells are skipped during $U,V,p,f$ calculations;

• •

elimination of flagged cells from the calculation leads to a characteristic “stair-stepped” approximation [30] $\Omega^{\sqcap}$ to $\Omega$ (cf. Figure 17);

• •

backward staggering of $U,V$ $C\forall$s with respect to $p,f$ $C\forall$s ensures an exact placement of momentum cell-centers along $\Omega^{\sqcap}$ (cf. Figure 17);

• •

no-slip boundary conditions are locally imposed at $U,V$ cell-centers lying on $\Omega^{\sqcap}$;

• •

the flagged $p,f$ $C\forall$s adjacent to $\Omega^{\sqcap}$ are employed as “ghost cells” (see cell $i,j$ in Figure 17) for imposing zero gradient conditions in $p^{\prime},f$ over $\Omega^{\sqcap}$ ($\mathcal{G}^{\perp}_{p^{\prime}_{I},f_{I}}=0$).

It should be noted that no local modifications in cell sizes were attempted for improving the approximation $\Omega^{\sqcap}$ and hence $\Omega^{\sqcap}\rightarrow\Omega$ would result only upon mesh refinement. The proposed methodology facilitates a simplified treatment of non-Cartesian geometries in the NWT even when the placement of solution variables is staggered. Results of wave transformation simulations are discussed in the sections that follow.

Weak transformation: SH waves $(T=1.25\,s;H=2.5\,cm)$

The first scenario considered is that of SH wave transformation over the trapezoidal bar. As waves transform over the submerged obstacle, $H/\lambda\propto x$ owing to generation of free and bound harmonics. Small, moderate and large $\mathsf{Ur}$ waves would thus co-exist along the NWT. In such a situation, an optimal spatial resolution $(nx_{\lambda})$ cannot be decided solely based on the regular wave generation guidelines provided in [13]; a mesh dependence analysis becomes necessary. In order to assess the mesh dependence of the solution, normalized free surface elevation $(\eta(t)/(0.5H))$ profiles of SH waves, recorded at WG2, WG4 and WG6, have been compared amongst five spatial resolutions: $nx_{\lambda}\cong 26,53,105,211,422$ during the interval $t\in[15T:16T]$. The following NWT framework has been considered (cf. Figure 16): $\mathbb{L}=29.95\,m$, $\mathbb{H}=0.44\,m$; $ny=100$, $ny_{H}\approx 12$; $T/\Delta t=5000$, $\wp=0.5H$; $\mathscr{S}=0.8$ is selected to prevent height damping of high $\mathsf{Ur}$ waves propagating over the bar crest. Results of the parametric analysis are shown in Figure 18. As anticipated, $H/\lambda\propto x$ means that “early” mesh independence is achieved at WG2 but not at WG4 where the waves are steeper. Interestingly, the $\eta(t)$ profiles in Figure 18 indicate that $\left.\frac{H}{\lambda}\right|_{\mathrm{WG2}}\cong\left.\frac{H}{\lambda}\right|_{\mathrm{WG6}}$ yet solution dependence on grid size is stronger at WG6. This finding establishes that mesh selection for wave transformation cannot be solely based upon guidelines extrapolated from regular wave generation simulations. From Figure 18, the $3072\times 100$ mesh is selected for further simulations and validation.

Spatial development of the SH wave train is presented at the top in Figure 19. It is apparent from the $\eta(x)$ profiles that the spatial evolution of SH waves is slow. This is because the “envelope” of SH waves evolves with a group velocity $(C_{G})$ which is less than the phase velocity $(C)$ of the individual wave packets; $C/C_{G}=1.4$ (cf. Table 3). In the initial stages of development, the train is essentially sinusoidal. For $t\geq 9T$, the SH wave train begins to shoal over the bar upslope which is characterized by $\lambda\downarrow$ and $H\uparrow$ and there is a visible onset of non-linearity. The wave topology resolved at $t=15T$ over the bar crest reaffirms the observation of [27] that SH waves shoal to more closely resemble higher-order Stokes waves. Subsequently, the waves “de-shoal” [27] over the lee-side of the bar with negligible super-harmonic generation. Nonetheless, the de-shoaling is accompanied by topological distortion $(t=20T)$ which is inturn attributed to a superposition of the transmitted carrier wave with free waves (that have gained amplitude owing to energy transfer from bound harmonics [38]). However, in the absence of harmonic generation, the lee-side wave train is essentially in intermediate water $(kh>1)$ and the free waves gradually dissipate. The SH wave train tends to regain it’s incident character after the bar downslope; wave transformation is hence termed “weak” in this case.

For the sake of validation, normalized $\eta(t)$ profiles recorded in the NWT are compared against experimental measurements of [27] during $t\in[15T:25T]$; this is reported at the bottom in Figure 19. The simulations demonstrate good agreement with experiments. The proposed NWT model is benchmarked against strong wave transformation in the following subsections.

Strong transformation: SLL waves $(T=2.0\,s;H=2.0\,cm)$

The SL waves are approximately $1.8\times$ longer than the SH waves (cf. Table 6). Long wave propagation over a submerged obstacle is characterized by extensive short wave generation on the lee side [38, 27]; wave transformation is thus “strong” in the case of both SLL as well as SLH waves. Both $H/\lambda\propto x$ and $kh\propto x$ hold during strong transformation owing to extensive harmonic generation on the lee side of the obstacle. Thus, the SLL waves have also been subjected to a mesh dependence analysis using the same NWT setup and range of mesh sizes as considered in the SH case. The normalized free surface elevation $(\eta(t)/(0.5H))$ profiles of SLL waves, recorded at WG2, WG4 and WG6, have been compared for different mesh sizes during the interval $t\in[10T:11T]$ and reported at the bottom in Figure 20. It should be noted that the SLL waves propagate in near-shallow water with $\approx 40\%$ larger group celerity $(C_{G}\approx 1.62\,m/s)$ in comparison to the SH waves $(C_{G}\approx 1.17\,m/s)$. Thus, local $\eta(t)$ signals reach a “fully developed state” much earlier (in terms of $t/T$) in the SLL case. It is evident from Figure 20 that the wave topology becomes mesh independent beyond $1536\times 100$ at WG2 and WG4 but $\eta(t)$ continues to evolve beyond $3072\times 100$ in the downslope region at WG6. Stronger mesh dependence observed at WG6 is largely attributable to: (a) larger steepness and (b) polychromatic character of waves transmitted to the lee side of the bar [39]. Nonetheless, in an attempt to maintain consistency with the SH setup, $3072\times 100$ is selected as the independent mesh for SLL wave simulations.

The spatial development of the SLL wave train is depicted at the top in Figure 21. The incident waves initially shoal over the windward/weather face of the bar and become progressively asymmetrical $(t>4T)$. This indicates an energy-gain by higher bound harmonics.

As the waves travel over the bar crest (which is a non-dispersive medium) triplet resonance occurs [40] and a prominent free-wave appears behind the steepened crest $(t=6T)$. As the waves propagate over the leeward face of the bar, the train “de-shoals” [40] and the primary wave breaks up into several smaller amplitude waves. Through Fourier decomposition, Huang and Dong [30] discovered that during de-shoaling, energy transfer is largely directed from the higher bound harmonics to free waves owing to a dramatic reduction in topological non-linearity in deep(er) water. This amplitude-gain by free waves on the lee side of the obstacle is clearly observed in Figure 21 for $t\geq 9T$. Unlike the bound harmonics, the free waves propagate relative to the carrier waves and are governed by the dispersion relation $(n\omega)^{2}=gk_{n}\tanh(k_{n}h)$ [38]. The leeward side is thus characterized by extensive relative motion amongst the waves.

The SLL simulations have also been validated against experiments [27, 36] in terms of normalized $\eta(t)$ profiles measured at six wave gauge locations (WG2-WG7); the results are reported at the bottom in Figure 21. The overall agreement between the proposed NWT model and experiments is observed to be good. As a matter of fact, long-time $(>5T)$ validation against wave transformation experiments is seldom attempted in the literature ([41] being an exception). Thus, the consistency retained between the present NWT simulations and experiments [27] over nine wave periods is commendable and is, in fact, superior to some of the recent simulations reported in the literature [42, 39, 43].

Strong transformation: SLH waves $(T=2.0\,s;H=4.0\,cm)$

Sinusoidal low frequency high (SLH) wave propagation over the trapezoidal bar is one of several cases considered by Huang and Dong [30] whilst numerically simulating wave deformation and vortex generation over submerged breakwaters. The SLH waves are twice as high as the SLL waves (cf. Table 6). Hence, stronger wave transformation and greater topological non-linearity (over the bar crest) are expected in the SLH case. In fact, there is some evidence that the SLH waves are susceptible to breaking during transformation. This is because the SLH waves are equivalent in steepness $(H/\lambda={0.011})$ to another wave design $(T=2.5\,s,H=5.2\,cm,h=0.4\,m\,\vdash{\mathsf{Ur}}\approx 18,H/\lambda={0.011})$ considered by Beji and Battjes [40] in their experiments, that yielded plunging breakers over the bar. Unfortunately, there is a lack of evidence whether breaking occurred in Huang and Dong’s [30] simulations or not; this seems highly unlikely given a (relatively) coarse spatial resolution of $nx_{\lambda}=93$ adopted by them in the NWT. In this respect, we demonstrate here that breaking of SLH waves in the NWT is suppressed by numerical damping of wave height over the bar crest and breaking can (rather) be “triggered” by increasing the blending parameter $\mathscr{S}$ (cf. Equation 10). Existence of such a relationship is expected; accurate simulation of the extent of wave-steepening over the bar crest necessitates higher-order treatment of momentum advection (cf. section 2.4).

Given the non-availability of $\eta(x)$ or $\eta(t)$ signals for SLH wave transformation in [30], it would be contributive to first illustrate spatio-temporal development of the SLH wave train against SLL waves in a comparative framework. Given much higher incident waves, the domain height is increased to $\mathbb{H}=0.75\,m$ which involves lengthening the propagation region to $\mathbb{L}=37.7\,m$ (which is the actual length of the experimental wave tank used in [27]). Rather than resorting to mesh dependence assessment, we directly select: $nx=1992\vdash nx_{\lambda}\approx 196$ and $ny=100\vdash ny_{H}\approx 13$. Further, $\left(\begin{matrix}\wp\\ \mathscr{S}\end{matrix}\right)\equiv\left(\begin{matrix}0.18\\ 1\end{matrix}\right)$ and $\left(\begin{matrix}0.18\\ 0.6\end{matrix}\right)$ have been considered for non-breaking and breaking SLH simulations respectively; the value of $\wp$ is adopted from the parametric selection procedure reported in Figure 6. Time is uniformly advanced using the forward Euler method with $\Delta t=T/5000$ retained from SH and SLL simulations. The above simulation setup helped ensure that $C_{max}\leq 0.25$ even when waves broke (as weak plungers [44]) over the bar crest.

Spatial topologies $(\eta(x))$ of both breaking and non-breaking SLH wave trains at $t=15T$ are shown in Figure 22 and compared with SLL waves generated using the same NWT setup (albeit with $\wp=0.5H$). It is demonstrated that, during SLH wave transformation, much stronger topological non-linearity is induced over the bar crest. Zoomed-in views of wave profiles over the bar crest (see bottom in Figure 22) indicate that, unlike SH waves, shoaled SL trains resemble shallow-water cnoidal waves [27]. Interestingly, doubling the incident wave height yields very similar $\eta(x)$ distribution on the lee side of the obstacle; albeit with higher waves. However, an examination of normalized $\eta(t)$ signals in both cases (cf. Figure 23) reveals that $\eta(t)/H$ actually decreases on the lee side of the bar in the SLH case. This finding is substantiated by a general observation made in [38, 30] and recently in [45] that $\eta(t)/H\downarrow$ on the lee side as (incident) $H\uparrow$ on the weather side of a submerged obstacle.

Following an assessment of vorticity fields (presented later in section 6.6) it has been hypothesized that $\eta(t)/H\downarrow$ on the lee side stems from rotational dissipation of packet energy towards inducing pervasive vortex generation in air, especially over the bar crest. Furthermore, Figure 22 and Figure 23 jointly support the claim that breaking can be triggered in the SLH wave train by $\downarrow{\mathscr{S}}$. It is obvious that wave breaking induces considerable energy dissipation and height damping past the breaking point ($x\gtrapprox 12\,m$; cf. Figure 22). Moreover, as evidenced from Figure 22 and Figure 23, breaking induces extensive short wave generation over the leeward-side of the bar [40, 45]. Short wave formation is manifested in packet energy transfer from $n\omega\,\,(n=1-3)$ to higher harmonic frequencies $n\geq 4$ [45] and since $C_{n\geq 4}\ll C_{1}$, the higher frequency waves almost appear “frozen in space” when visualized against the carrier wave train.

Vorticity dynamics and harmonic decomposition

One of the key advantages of the proposed NWT model is the two-phase NSE philosophy adopted for modeling the waves. As evident from our previous work [13], the NSE paradigm facilitates interrogation of the solution in both air and water phases, thus providing greater insight into momentum dynamics. In context to submerged bars, the design goal is to induce short wave behavior on the lee side, thereby preventing beach erosion [45]. In addition, vorticity dynamics induced over the weather and lee faces as well as over the bar crest is of great interest from a design point of view since persistent vortical activity (in water) would induce hydrodynamic scour in the long-term. Hence, for the strong wave transformation scenarios, it is contributive to investigate the vorticity dynamics in air as well as the water phases. The spatial variation of amplitudes $a^{(n)}\,\,;\,\,n=1-3$ of the first three harmonics is also presented as a rigorous validation exercise.

As the incident waves encounter the obstacle, energy transfer to higher harmonics occurs from the fundamental; the higher harmonics consist of both bound as well as free waves [30, 38]. However, separation of free and bound components for $n>2$ is not our concern here; a single wave gauge based Discrete Fourier Transform (DFT) is hence sufficient. For DFT, the free surface elevation $\eta(x,t)$ has been recorded in each cell falling within $x\in[4\,m:22.5\,m]$ every $\Delta t=0.0004\,s$ during the time interval $t\in[14T:15T]$. We proceed by defining the Fourier transform $\hat{\eta}^{(n)}(x)$ of the free-surface elevation,

[TABLE]

where $\mathrm{i}=\sqrt{-1}$ and $n=1,2,3,\dotsc,$. The integral in Equation 20 is evaluated using the composite Simpson’s rule. Then, $a^{(n)}(x)$ could be directly obtained using [19],

[TABLE]

where the real and imaginary parts of the complex elevation $\hat{\eta}^{(n)}(x)$ are given by,

[TABLE]

where $N\equiv T/\Delta t$ is the number of samples and $t\equiv N\Delta t$ is the physical time.

The spatial variation of harmonic amplitudes $a^{(n)}$ obtained using Equation 21 is reported in Figure 24. Evidently, an acceptable validation is obtained against the simulations reported in the literature [46, 30]. Energy transfer from the fundamental $(n=1)$ to higher harmonics $(n=2,3)$ over the bar crest is evident from the plots. Interestingly, a reduction of $\left|\eta(t)/H\right|$ on the lee side for SLH waves is also reflected in the amplitude decomposition. During breaking, the greatest energy loss occurs from $n=2$ followed by the fundamental $(n=1)$. It is reaffirmed from Figure 24 that wave-breaking does not induce any upstream influence in the NWT.

The contours of vorticity $\vec{\Omega}$ (superimposed with velocity vectors and the free-surface contour) alongwith $a^{(n)}(x)$ decomposition are shown in Figure 24 for SLL, SLH and breaking SLH waves. The velocity vectors clearly indicate a region of intermediate/deep water $(kh>1)$ on the lee side of the bar in all cases (which is desired). Unlike rectangular dikes, gradual shoaling over the bar upslope precludes flow separation on the lee side; this exhibits reduced susceptibility of the lee side to hydrodynamic scour [30]. The $\vec{\Omega}$ contours indicate that the SLH waves act to expand the region of interaction between air and water by inducing strong regions of rotation, especially over the bar crest. It is hypothesized that the observed vortical activity would drain packet energy to viscous effects thereby causing relatively smaller waves to be transmitted to the lee side of the bar $(\uparrow H\vdash\downarrow\left|\eta(t)/H\right|)$ [38]. The region of interaction around the air-water interface gets further expanded as SLH waves break over the bar crest with consecutive breaking events leading to the formation of a surface current (cf. bottom in Figure 24).

To the best of the authors’ knowledge, the vorticity dynamics of breaking/non-breaking wave transformation has never been reported in the literature and is therefore one of the novel aspects of the present work. At this juncture, the authors would also like to assert the validity of the vortex dynamics reported in Figure 24 from the standpoint of the CFD paradigm used for simulation. Based on the instantaneous Navier-Stokes equations, the NWT solves for the instantaneous velocity field $(\vec{V})$ induced by the waves and not the mean ($\overline{V}$; unlike RANS) or filtered ($\tilde{V}$; unlike LES) fields. It is obvious that both $\overline{V}$ and $\tilde{V}$ are extracted from $\vec{V}$ with the residual terms being accounted for using a turbulence closure model. Hence, it is argued that the present NWT model automatically accounts for the generation of “two-phase turbulence” at the wave surface and within the air and water phases during breaking. Thus, there has been a deliberate effort to “not model” the effects of turbulence but rather to allow the flow solver to capture turbulent and vortical scales to the extent that is allowed by the mesh. Whilst this approach may strike one as being overtly simplistic, the same in reality preserves the flow simulation from turbulence-modeling-induced artifacts such as over-production of turbulent viscosity [47] and unrealistic dampening of fluctuations caused as a consequence.

Summary and outlook

We propose a robust numerical wave tank (NWT) algorithm (and accompanying design paradigm) for simulating ocean wave propagation and wave-structure interaction scenarios in a two-phase, Navier-Stokes framework. Robustness of the proposed algorithm is primarily manifested in:

• •

a novel kinematic-stretching based modified inflow boundary wavemaker which is volume-preserving and is characterized by a single design variable,

• •

a strategy for blending low and higher-order momentum advection schemes that is effective in restricting numerical damping of steep waves without inducing additional wave distortion at the air-water interface and,

• •

a simplified methodology for treating non-Cartesian submerged/emergent boundaries on staggered grids which ensures sustenance of a tight coupling between computed pressure and velocity fields.

The proposed NWT model has been successfully benchmarked against several problems including generation of monochromatic waves of varying steepness, simulation of free-wave generation during piston-type wavemaker motion, generation of a deep-water irregular wave train and (weak and strong) wave transformation (including weak plunging of waves) occurring over a submerged trapezoidal bar. Very good agreement with analytical, numerical and experimental studies reported in the literature is obtained. The proposed NWT algorithm has already been parallelized using MPI [48] and is currently being applied toward simulating the hydrodynamics of OWC devices.