Hydrodynamics of Regular Breaking Wave

Abstract

Turbulence and undertow currents play an important role in surf-zone mixing and transport processes; therefore, their study is fundamental for the understanding of nearshore dynamics and the related planning and management of coastal engineering activities. Pioneering studies qualitatively described the features of breakers in the outer region of the surf zone. More detailed information on the velocity field under spilling and plunging breakers can be found in experimental works, where single-point measurement techniques, such as Hot Wire Anemometry and Laser Doppler Anemometry (LDA), were used to provide maps of the flow field in a time-averaged or ensemble-averaged sense. Moreover, the advent of non-intrusive measuring techniques, such as Particle Image Velocimetry (PIV) provided accurate and detailed instantaneous spatial maps of the flow field. However, by correlating spatial gradients of the measured velocity components, the instantaneous vorticity maps could be deduced. Moreover, the difficulties of measuring velocity due to the existence of air bubbles entrained by the plunging jet have hindered many experimental studies on wave breaking encouraging the development of numerical model as useful tool to assisting in the interpretation and even the discovery of new phenomena. Therefore, the development of an WCSPH method using the RANS equations coupled with a two-equation k–ε model for turbulent stresses has been employed to study of the turbulence and vorticity distributions in in the breaking region observing that these two aspects greatly influence many coastal processes, such as undertow currents, sediment transport and action on maritime structures.

1. Introduction

Wave breaking is one of the most important process for coastal engineers since it greatly influences both the transport processes and the magnitudes of the forces on coastal structures [1, 2]. Wave breaking in the surf zone drives complicated turbulent structures and for this reason breaking is possibly the most difficult wave phenomenon to describe mathematically [3, 4, 5].

Pioneering experimental studies were carried out by [6, 7, 8, 9, 10, 11], who described the velocity field under plunging breakers in the outer region of the surf zone; more recently by [12, 13, 14, 15, 16, 17]. As observed experimentally by [18], during the pre-breaking stages, the maximum turbulence intensity appears at the core of the main vortex and decreases as the vortex moves downstream. Additional turbulence is then generated near the free surface during the breaking process.

Moreover, the difficulties of measuring velocity due to the existence of air bubbles entrained by the plunging jet have hindered many experimental studies on wave breaking encouraging the development of numerical model as useful tool to assisting in the interpretation and even the discovery of new phenomena. One of the great advantages of the numerical models is their ability to disclose the evolutions of undertow currents and turbulence quantities in the spatial and temporal domains, which are too expensive to be investigated by experiments. Therefore, the main emphasis for research is placed on the application and development of numerical methods. Furthermore, for consistent and accurate results, it is essential to calibrate the numerical models with experimental data.

The numerical models can be classified as Eulerian or Lagrangian method. In Eulerian method, the space is discretized into a grid or mesh and the unknown values are defined at the fix points, while a Lagrangian method tracks the pathway of each moving mass point. The Eulerian methods such as the finite difference methods (FDM), finite volume methods (FVM) and the finite element methods (FEM) have been widely applied in many fields of engineering because are very useful to solve differential or partial differential equations (PDEs) that govern the concerned physical phenomena [19, 20, 21, 22, 23]. Despite the great success, grid based numerical methods suffer from difficulties in some aspects such as the use of grid/mesh makes the treatment of discontinuities (e.g., wave breaking, cracking and contact/separation) difficult because the path of discontinuities may not coincide with the mesh lines.

Therefore, during the last years, research has been focused on Lagrangian techniques such as Discrete Element Method (DEM) [24], Smoothed Particle Hydrodynamics (SPH) [25], Immersed Particle Method (IPM) [26, 27]. The development and applications of the major existing Lagrangian methods have been addressed in some review articles such as [28, 29, 30]. In general, the Lagrangian methods provide accurate and stable numerical solutions for integral equations or partial differential equations (PDEs) with all kinds of possible boundary conditions using a set of arbitrarily distributed nodes or particles. During the last years, Smoothed Particle Hydrodynamics (SPH) has become a very powerful method for CFD problems governed by the Navier–Stokes equations such as fluid-dynamic problems with highly non-linear deformation [31, 32, 33, 34, 35, 36, 37]; multi-phase flows for coastal and other hydraulic applications with air-water mixture sand sediment scouring [38, 39, 40, 41, 42]; oscillating jets inducing breaking waves [43] and nonbuoyant jets in a wave environment [44, 45, 46]; fluid/structure/soil-interaction [47, 48, 49]; hydraulic jumps [50, 51, 52, 53]; multi-phase flows and oil spill [54, 55].

The present chapter is organized as follows. First, an WCSPH method is developed using the RANS equations and a two-equation k–ε model is formulated using the particle approach. Then the numerical model is employed to reproduce breaking in spilling and plunging waves in a sloped wave channel. The experimental data by [14] are used to check the model results. This reveals the importance of experimental data in these studies. The present chapter is aimed to describe some recent results obtained within the frame of numerical and experimental analyses of wave breaking. The new insight is the investigation of the ability of WCSPH with a k–ε turbulence closure model to disclose the turbulence dispersion and the temporal and spatial evolutions of turbulence quantities in different types of breakers.

2. Mathematical formulation

A Lagrangian numerical model is developed to solve free surface turbulent flows. The flow field is governed by the Reynolds Averaged Navier–Stokes (RANS) equations and the k–ε turbulence equations [56]. These equations are solved by the WCSPH method in which an artificial compressibility is introduced to solve explicitly in time the equations of motion of an incompressible fluid.

Using the SPH approach, the fluid flow domain is initially represented by a finite number of particles. These particles can be viewed as moving numerical nodes, which move according to the governing equations and boundary conditions. Each discrete point is associated to an elementary fluid volume (or particle) i, which has position xi and constant mass mi.

To find the value of $\mathit{a}\left(\overline{\mathit{x}},t\right)$ at a generic point $\overline{\mathit{x}}$ an interpolation is applied from the nodal values ai(t) through a kernel function $W=\left(\mathit{x}-\overline{\mathit{x}},\eta \right)$ as follows:

$$\mathit{a}\left(\overline{\mathit{x}},t\right)\approx <\mathit{a}\left(\overline{\mathit{x}},t\right)>=\sum _{j=1}^N\frac{m_j}{{\rho }_j}\mathit{a}\left({\mathit{x}}_j,t\right)W\left({\mathit{x}}_j-\overline{\mathit{x}},\eta \right)$$E1

where ρ is the fluid density, and the summation is extended to all the N particles located inside the sphere of radius 2η centered on $\overline{\mathit{x}}$. The kernel function is continuous, non-zero only inside a sphere $\mathit{x}-\overline{\mathit{x}}<2\eta$ and tends to the Dirac delta function when η (defined as the smoothing length) tends to zero. There are different available kernel functions and the kernel operations can be inaccurate for cases where the particle distribution is non-uniform or the support for the interpolations is incomplete [55]; Quinlan et al., [57], Randles and Libersky [58] and Bonet and Lok [59] introduced a Kernel correction which ensure at least first-order consistency; however the corrected kernel is non-symmetric which leads to non-conservative interpolations.

Dehnen and Aly [60] showed that the Wendland kernel function [61] is more computationally convenient than the B-spline function, allowing better numerical convergence; Liu and Liu [62] showed that the quintic-spline function [63] is more effective in interpolating the second-order derivatives. The SPH computations discussed in the present paper were based on the cubic-spline kernel function proposed by [64] that is more effective in the simulation of several different hydraulic flows [65, 66].

The advantage of SPH approach is that differential operator applied to $\mathit{a}\left(\overline{\mathit{x}},t\right)$ can be approximated by making use of the gradient of the kernel function. For instance, the divergence $\nabla ·\mathit{a}\left(\overline{\mathit{x}},t\right)$ can be approximated by:

$$\nabla ·\mathit{a}\left(\overline{\mathit{x}},t\right)\approx <\nabla ·\mathit{a}\left(\overline{\mathit{x}},t\right)>=\sum _{j=1}^N\frac{m_j}{{\rho }_j}\left[\mathit{a}\left(\overline{\mathit{x}},t\right)-\mathit{a}({\mathit{x}}_j,t)\right]\nabla W\left({\mathit{x}}_j-\overline{\mathit{x}},\eta \right)$$E2

For further details on the different methods for SPH approximations of all the vector operators, the reader can see [67, 68]. In a Lagrangian frame, the Reynolds-averaged Navier–Stokes (RANS) equations and the k–ε turbulence equations take the following form

$$\left\{\begin{array}{l}\frac{\mathit{D\rho }}{\mathit{Dt}}+\rho \nabla ·\mathit{v}=0\\ \frac{D\mathit{v}}{\mathit{Dt}}=\frac{1}{\rho }\nabla p+\frac{1}{\rho }\nabla \bullet T+\mathit{g}\\ p-p_0=c^2\left(\rho -{\rho }_0\right);T={\mu }_{T}S\end{array}\right.$$E3

where v = (u, v) is the velocity vector, p is pressure, g is the gravity acceleration vector, T is the turbulent shear stress tensor, c is the speed of sound in the weakly compressible fluid, μ_T is the dynamic eddy viscosity, S is rate-of-strain tensor and the subscript 0 denotes a reference state for pressure computation.

The RANS equations (3) in the SPH semi-discrete form become

$$\left\{\begin{array}{l}⟨\frac{D{\rho }_i}{\mathit{Dt}}⟩=\sum _j{m}_j\left(v_i-v_j\right)\nabla W_{\mathit{ij}}\\ ⟨\frac{D{\mathit{v}}_{\mathit{i}}}{\mathit{Dt}}⟩=-\sum _j{m}_j\left(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}\right)\nabla W_{\mathit{ij}}+\sum _j\frac{m_j}{{\rho }_j}\left(T_i-T_j\right)\bullet \nabla W_{\mathit{ij}}+\mathit{g}\\ p_i-p_0=c_i^2\left({\rho }_i-{\rho }_0\right);T_i={\mu }_{T_i}{S}_i\end{array}\right.$$E4

where W_ij is a shorthand notation for $W\left({\mathit{x}}_i-{\mathit{x}}_j,\eta \right)$, renormalized through a procedure which enforces consistency on the first derivatives to the 1st order [69], leading to a 2nd order accurate discretization scheme in space. The semi-discretized system (4) is then integrated in time by a 2nd order two-step XSPH explicit algorithm [70].

The momentum equation is then solved to yield an intermediate velocity field $\widehat{\mathit{v}}$, which is then corrected through a smoothing procedure based on the values of the neighboring fluid particles.

$${\mathit{v}}_i=\left(1-{\varphi }_v\right){\widehat{\mathit{v}}}_i+{\varphi }_v\frac{{\sum }_j\frac{m_j}{{\rho }_j}{\widehat{\mathit{v}}}_j{W}_{\mathit{ij}}}{{\sum }_j\frac{m_j}{{\rho }_j}{W}_{\mathit{ij}}}$$E5

using a velocity smoothing coefficient φ_v. The corrected velocity value is then used to update the particle position and to solve the continuity equation. The new density values are finally used to compute pressure, according to the equation of state.

A pressure smoothing procedure is also applied to the difference between the local and the hydrostatic pressure values [71] in order to reduce the numerical noise in pressure evaluation which is present, in particular in WCSPH, owing to high frequency acoustic signals [72]. The present method is applied only to the difference between the intermediate pressure field $\widehat{p}$ and the hydrostatic pressure gradient to ensure the conservation of total volume of the particle system for long time simulations [73].

The eddy viscosity coefficient ${\mu }_T=c_{\mu }\frac{k^2}{\epsilon }$ in Eq. (4) was evaluated through a k-ε model by [74]:

$$\left\{\begin{array}{l}\frac{D{k}_i}{\mathit{Dt}}=P_{k_i}+\frac{1}{{\sigma }_k}\underset{j}{\sum m_j}\frac{{\nu }_{T_i}+{\nu }_{T_j}}{{\rho }_i+{\rho }_j}\frac{k_i-k_j}{r_{\mathit{ij}}^2+0.01h^2}{\mathit{r}}_{\mathit{ij}}·\nabla W_{\mathit{ij}}-{\epsilon }_i\\ \frac{D{\epsilon }_i}{\mathit{Dt}}=\frac{1}{{\sigma }_{\epsilon }}\sum _j{m}_j\frac{{\nu }_{T_i}+{\nu }_{T_j}}{{\rho }_i+{\rho }_j}\frac{{\epsilon }_i-{\epsilon }_j}{r_{\mathit{ij}}^2+0.01h^2}{\mathit{r}}_{\mathit{ij}}·\nabla W_{\mathit{ij}}++C_{{\epsilon }_1}\frac{{\epsilon }_i}{k_i}{P}_{k_i}+C_{{\epsilon }_2}\frac{{\epsilon }_i}{k_i}\sum _j\frac{m_j}{{\rho }_j}{\epsilon }_j{W}_{\mathit{ij}}\end{array}\right.$$E6

where k_i is the turbulent kinetic energy per unit mass, ε is the dissipation rate of turbulent kinetic energy, P_k is the production of turbulent kinetic energy depending on the local rate of deformation and ${\nu }_T$ is the eddy viscosity and ${\mathit{r}}_{\mathit{ij}}={\mathit{x}}_i-{\mathit{x}}_j$. There are several empirical coefficients in the k–ε turbulent closure model. In this paper the set of constant values recommended by [74], i.e., $c_{\mu }$ = 0.09, ${\sigma }_k=1$, ${\sigma }_{\epsilon }$ = 1.3, $C_{{\epsilon }_1}$ = 1.44 and $C_{{\epsilon }_2}$=1.92, is adopted.

According to Eq. (6) both the production term and dissipation term for ε become singular when k approaches zero. Furthermore, no turbulence energy can be produced if there is no turbulent kinetic energy initially. Thus, it is necessary to “seed” a small amount of k in both the initial condition and inflow boundary condition. In this paper the initial seeding of turbulent kinetic energy recommended by [75] is adopted.

3. Validation and application

3.1 Experimental set up

The results obtained from the numerical model outlined in the previous section have been validated against extensive experimental data [14], and then used to obtain further insight in the physics of the flow here analyzed.

The detailed experimental setup has been given in De Serio and Mossa [14]. Here only some important parameters are summarized.

Experiment was carried out in the wave flume 45 m long and 1 m wide of the Department of Civil, Environmental, Land, Building Engineering and Chemistry (DICATECh) of the Polytechnic University of Bari (Italy). A beach with constant slope of 1/20 is connected to a region with constant water depth of h = 0.7 m. The wave generating system is a piston-type one, with paddles producing the desired wave by providing a translation of the water mass, according to the proper input signal. The instantaneous Eulerian velocities were acquired by a backscatter, two-component, four beam Laser Doppler Anemometer (LDA) system and a Dantec LDA signal processor (58 N40 FVA Enhanced) based on the covariance technique. The wave elevations were measured with a resistance probe placed in the transversal section of the channel crossing the laser measuring volume.

Figure 1a–f show the different parts of the complex experimental apparatus, which comprises the LDA system, the resistance wave gauge system and the wavemaker system. Further details about the experimental tests can be found in [14].

Figure 1.

Experimental apparatus: (a) LDA probe; (b) DANTEC FVA signal processor and process computer; (c) laser coherent Innova and Dantec 2D fiber flow optics; (d) process computer with a AD/DA board for the wavemaker control; (e) a part of the wave channel; (f) the wavemaker.

A sketch view of the experimental setup is shown in Figure 2.

Figure 2.

Sketch view of experimental setup.

Table 1 shows the main parameters of the examined waves listed for each experiment, such as the offshore wave height H₀, the wave period T and the deepwater wavelength L₀, estimated in section 76, where the bottom is flat and the mean water depth h is equal to 0.70 m. In the experiments, the regular wave had a height H₀ = 11 cm and a period T = 2.0 s for the spilling breaker case (T1), while H₀ = 6.5 cm and T = 4.0 s were used to generate a plunging breaker (T2). Table 1 shows also the Irribarren number ξ₀, computed for the two tests from the following equation

	H0 (cm)	T (s)	L0 (m)	d (m)	ξ0	Breaking type
T1	11	2	4.62	0.70	0.37	Spilling/plunging
T2	6.5	4	10.12	0.70	0.74	Plunging

Table 1.

Experimental parameters of the analyzed regular waves.

$${\xi }_0=\frac{\mathit{tan}\beta }{\sqrt{\frac{H_0}{L_0}}}$$E7

in which β is the bottom slope angle.

Water surface elevations and velocities were measured at six different sections along the longitudinal axis of symmetry of the wave channel named 76, 55, 49, 48, 47, 46 and 45 (see Table 2). Specifically, for all two tests, section 48 was in the pre-breaking region, section 47 was where the incipient breaking occurred, while in sections 46 and 45, the wave re-arranged into a bore.

Investigated section	Distance from paddles (m)	d (cm)
Section 76	10.56	70.0
Section 55	19.80	31.0
Section 49	22.44	16.5
Section 48	22.88	14.0
Section 47	23.32	11.3
Section 45	24.20	8.5

Table 2.

Location of measurement sections.

3.2 Numerical model setup and validation

The WCSPH method coupled with a k–ε turbulence model has been employed to reproduce the above experiment. The computational domain has been reduced to be 20.0 m long so as to save computing expenses (Figure 3).

Figure 3.

Sketch of the computational domain wave channel with location of the seven investigated sections, used to calibrate the numerical model.

The adopted offshore boundary condition guarantees a regular development of the wave train before the sloping section of the channel and, therefore, does not influence the quality of the numerical solution, as shown by [32].

For both the two tests, the offshore boundary condition has been treated as dynamic boundary condition modeled by a numerical wave paddle also composed of ghost particles whose motion has been forced to obtain the frequency and amplitude of the wave paddle needed to generate the desired sinusoidal wave [76]. The initial water depth was set equal to 0.70 m. In the present simulations, the initial particle spacing Σ = 0.022 m, the value of η/Σ = 1.5 and $\phi$ = 0.01, recommended by De Padova et al. [32], have been adopted.

The instantaneous SPH particle distribution and velocity magnitude snapshots of the breaking wave are shown in Figures 4a–e and 5a–e, respectively, for the spilling and plunging breakers. These results show that the general features of wave breaking, collapsing and a turbulent bore propagating have been well captured by the SPH computations.

Figure 4.

Instantaneous SPH velocity field in the SPH simulation of spilling wave (T1): (a) before; (b)–(c) during and (d)–(e) after breaking.

Figure 5.

Instantaneous SPH velocity field in the SPH simulation of plunging wave (T2): (a) before; (b)–(c) during and (d)–(e) after breaking.

In order to further verify the accuracy of the SPH model the time series of wave elevations, horizontal and vertical velocities at the investigated sections (Figure 3) have been compared with the experimental data of De Serio and Mossa [14]. As an example, in Figure 6a–c both laboratory and numerical wave surface elevations, and velocities at vertical sections 48 and 45 are plotted for T1, referring to the point located at 1 cm from the bottom. The agreement between the calibrated numerical results and the laboratory measurements is fairly good.

Figure 6.

Instantaneous computed and measured (a) wave elevations, (b) horizontal and (c) vertical velocities in section 48 and section 45 for T1.

4. Results and discussion

One of the great advantages of the numerical models is their ability to show the evolutions of vorticity and turbulence quantities in the spatial and temporal domains, which are too expensive to be investigated by experiments. Using the SPH computational results, the turbulent kinetic energy distributions are shown in Figures 7a–e and 8a–e, respectively, for the spilling (T1) and plunging (T2) waves. For both breakers, the turbulence quantity has the largest values near the free surface and decreases into the water column. However, the results highlight that there exist fundamental differences in the dynamics of turbulence between the spilling and plunging breakers, which can be related to the processes of wave breaking production.

Figure 7.

Instantaneous turbulence intensity distributions in the SPH simulation of spilling wave (T1): (a) before; (b)–(c) during and (d)–(e) after breaking.

Figure 8.

Instantaneous turbulence intensity distributions in the SPH simulation of plunging wave (T2): (a) before; (b)–(c) during and (d)–(e) after breaking.

For the spilling wave (T1), higher turbulence levels are mainly concentrated in the breaking wave front and the highest turbulence level appears in the roller region (Figure 7d). After the breaking, as the wave propagates forward, the turbulence kinetic energy decreases (Figure 7e). Instead, the turbulence levels increase rapidly after the wave breaking for the plunging case (T2) as shown in Figure 8c–e. The maximum turbulence level is generated as the plunging jet touches down on the wave trough (Figure 8d) in sections 46–45 (Figure 2); After the breaking, the roller continues to spread downwards and therefore high turbulence levels are generated beneath the free surface after breaking (Figure 8e).

Using the SPH computational results, the vorticity maps are shown in Figures 9a–e and 10a–e, respectively, for the spilling and plunging waves. Vorticity is defined as

Figure 9.

Instantaneous values of ω distributions in the SPH simulation of spilling wave (T1): (a) before; (b)–(c) during and (d)–(e) after breaking.

Figure 10.

Instantaneous values of ω distributions in the SPH simulation of spilling wave (T2): (a) before; (b)–(c) during and (d)–(e) after breaking.

$$\omega =\left(\frac{\partial u}{\partial z}\right)-\left(\frac{\partial v}{\partial x}\right)$$E8

and is computed using instantaneous values of the horizontal and vertical velocity.

As noted by several authors [77, 78], for both breakers (T1 and T2), when the breaking begins, positive vorticity occupies the whole region of the surface roller and spreads out over the whole water column. However, the vorticity levels increase rapidly after the wave breaking for the plunging case (T2) due to the strong impingement of the jet on the forward trough, inducing a propagation of the positive vorticity towards the bottom (Figure 10c–e).

Moreover, the results highlight that there exist differences in the dynamics of vorticity between the spilling and plunging breakers. In fact, only during spilling formation (T1), small structures of negative vorticity are generated, instead when the plunging breaker (T2) occurs the fluid is relatively free of negative vorticity regions.

Figures 11 and 12 show the comparison between the instantaneous map of vorticity and of the surface parallel convective acceleration for the spilling and plunging waves (T1 and T2) when the breaking begins at time step of Figures 9b and 10b, respectively. The surface parallel convective acceleration here has been computed following [79]. As noted by Dabiry and Gharib [80], for both breakers (T1 and T2), a flow deceleration (Figures 11b and 12b) occurs in the same location where peaks of positive vorticity appear (Figures 11a and 12a). Therefore, the present results confirm the findings by Dabiri and Gharib [80] that the vorticity is convected due to the sharp velocity gradient of the fluid near the free surface with respect to the fluid below.

Figure 11.

SPH simulation of spilling wave (T1): (a) Vorticity map and (b) surface-parallel convective acceleration map.

Figure 12.

SPH simulation of plunging wave (T2): (a) Vorticity map and (b) surface-parallel convective acceleration map.

5. Conclusions

In the present chapter a WCSPH method has been developed using the RANS equations and a two-equation k–ε model has been formulated using the particle approach. Then the numerical model has been employed to reproduce breaking in spilling and plunging waves in a sloped wave channel. The experimental data by [14] have been used to check the model results. Finally, we have fully exploited the advantages of numerical modeling to disclose fundamental differences between the different types of breakers by investigating the temporal and spatial evolutions of turbulence quantities and vorticity field.

For the spilling wave (T1), during the pre-breaking and breaking stages, the maximum turbulence intensity has been generated near the free surface and decreases as the vortex moves downstream. Instead, the turbulence levels increased rapidly after the wave breaking for the plunging case (T2). In fact, the maximum turbulence level was generated as the plunging jet touches down on the wave trough and after the breaking, the roller continued to spread downwards and therefore high turbulence levels were generated beneath the free surface after breaking.

For both breakers (T1 and T2), analyzing the instantaneous vorticity distributions, when the breaking begins, positive vorticity has occupied the whole region of the surface roller and has spread out over the whole water column. However, only during spilling formation (T1), small structures of negative vorticity have been generated, instead when the plunging breaker (T2) occurs the fluid was relatively free of negative vorticity regions.

Furthermore, comparing the instantaneous map of vorticity and of the surface parallel convective acceleration for the spilling and plunging waves (T1 and T2), the present results confirmed the findings by Dabiri and Gharib [80] that the vorticity was convected due to the sharp velocity gradient of the fluid near the free surface with respect to the fluid below.