Hydrodynamics of Regular Breaking Wave

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 timeaveraged 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.


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.
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.

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 a x, t ðÞ at a generic point x an interpolation is applied from the nodal values ai(t) through a kernel function W ¼ x À x, η ðÞ 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 x.Thekernel function is continuous, non-zero only inside a sphere x À x < 2η 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 a x, t ðÞ can be approximated by making use of the gradient of the kernel function. For instance, the divergence ∇ Á a x, t ðÞ 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 Reynoldsaveraged 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 W ij is a shorthand notation for W x i À x j , η ÀÁ , 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 fieldv, 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 [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 fieldp 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 μ T ¼ c μ k 2 ε in Eq. (4) was evaluated through a k-ε model by [74]: 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 ν T is the eddy viscosity and r ij ¼ x i À 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 μ = 0.09, σ k ¼ 1, σ ε = 1.3, C ε 1 = 1.44 and C ε 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" asmall 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.

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, twocomponent, 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].
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 H 0 , the wave period T and the deepwater wavelength L 0 , 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 0 = 11 cm and a period T = 2.0 s for the spilling breaker case (T1), while H 0 =6.5cm 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 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 prebreaking region, section 47 was where the incipient breaking occurred, while in sections 46 and 45, the wave re-arranged into a bore.

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 [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 φ = 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.
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.

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  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 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

11
Hydrodynamics of Regular Breaking Wave DOI: http://dx.doi.org /10.5772/intechopen.94449 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), aflowdeceleration (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.

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.