Experimental parameters of the analyzed regular waves.
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.
- regular breaking waves
- shear stress
- kinetic energy
- physical modeling
- numerical modeling
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 , 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) , Smoothed Particle Hydrodynamics (SPH) , 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  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  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 . 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 at a generic point an interpolation is applied from the nodal values ai(t) through a kernel function as follows:
where ρ is the fluid density, and the summation is extended to all the N particles located inside the sphere of radius 2η centered on . The kernel function is continuous, non-zero only inside a sphere 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 ; Quinlan et al., , Randles and Libersky  and Bonet and Lok  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  showed that the Wendland kernel function  is more computationally convenient than the B-spline function, allowing better numerical convergence; Liu and Liu  showed that the quintic-spline function  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  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 can be approximated by making use of the gradient of the kernel function. For instance, the divergence can be approximated by:
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
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
where Wij is a shorthand notation for , renormalized through a procedure which enforces consistency on the first derivatives to the 1st order , 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 .
The momentum equation is then solved to yield an intermediate velocity field , which is then corrected through a smoothing procedure based on the values of the neighboring fluid particles.
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  in order to reduce the numerical noise in pressure evaluation which is present, in particular in WCSPH, owing to high frequency acoustic signals . The present method is applied only to the difference between the intermediate pressure field and the hydrostatic pressure gradient to ensure the conservation of total volume of the particle system for long time simulations .
where ki is the turbulent kinetic energy per unit mass, ε is the dissipation rate of turbulent kinetic energy, Pk is the production of turbulent kinetic energy depending on the local rate of deformation and is the eddy viscosity and . There are several empirical coefficients in the k–ε turbulent closure model. In this paper the set of constant values recommended by , i.e., = 0.09, , = 1.3, = 1.44 and =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  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 , 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 . 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 .
A sketch view of the experimental setup is shown in Figure 2.
Table 1 shows the main parameters of the examined waves listed for each experiment, such as the offshore wave height H0, the wave period T and the deepwater wavelength L0, 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 H0 = 11 cm and a period T = 2.0 s for the spilling breaker case (T1), while H0 = 6.5 cm and T = 4.0 s were used to generate a plunging breaker (T2). Table 1 shows also the Irribarren number ξ0, computed for the two tests from the following equation
|H0 (cm)||T (s)||L0 (m)||d (m)||ξ0||Breaking type|
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)|
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).
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 .
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 . 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 = 0.01, recommended by De Padova et al. , 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.
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 . 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.
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.
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).
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 . As noted by Dabiry and Gharib , 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  that the vorticity is convected due to the sharp velocity gradient of the fluid near the free surface with respect to the fluid below.
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  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  that the vorticity was convected due to the sharp velocity gradient of the fluid near the free surface with respect to the fluid below.