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

## *Hydrodynamics of Regular Breaking Wave DOI: http://dx.doi.org/10.5772/intechopen.94449*

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:

$$a(\bar{\mathbf{x}},t) \approx < a(\bar{\mathbf{x}},t)> = \sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} a(\mathbf{x}\_j, t) \mathcal{W}(\mathbf{x}\_j - \overline{\mathbf{x}}, \eta) \tag{1}$$

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*. The *kernel* 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:

$$\nabla \cdot \boldsymbol{a}(\bar{\mathbf{x}},t) \approx \boldsymbol{\varsigma} \cdot \nabla \cdot \boldsymbol{a}(\bar{\mathbf{x}},t) \\ > \sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} \left[ \boldsymbol{a}(\bar{\mathbf{x}},t) - \boldsymbol{a}(\boldsymbol{x}\_j,t) \right] \nabla \boldsymbol{W} \left( \mathbf{x}\_j - \bar{\mathbf{x}}, \eta \right) \tag{2}$$

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

$$\begin{cases} \frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{v} = 0 \\\\ \frac{D\mathbf{v}}{Dt} = \frac{1}{\rho} \nabla p + \frac{1}{\rho} \nabla \bullet T + \mathbf{g} \\\\ p - p\_0 = c^2 (\rho - \rho\_0); T = \mu\_T \mathbf{S} \end{cases} \tag{3}$$

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, *μ<sup>T</sup>* 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

$$\begin{cases} \left< \frac{\mathbf{D}\rho\_i}{Dt} \right> = \sum\_j m\_j (\boldsymbol{v}\_i - \boldsymbol{v}\_j) \nabla W\_{\vec{\eta}} \\\\ \left< \frac{\mathbf{D}\mathbf{v}\_i}{Dt} \right> = -\sum\_j m\_j \left( \frac{p\_i}{\rho\_i^2} + \frac{p\_j}{\rho\_j^2} \right) \nabla W\_{\vec{\eta}} + \sum\_j \frac{m\_j}{\rho\_j} (T\_i - T\_j) \cdot \nabla W\_{\vec{\eta}} + \mathbf{g} \\\\ p\_i - p\_0 = c\_i^2 (\rho\_i - \rho\_0); T\_i = \mu\_{T\_i} \mathbf{S}\_i \end{cases} \tag{4}$$

where *Wij* is a shorthand notation for *<sup>W</sup> <sup>x</sup><sup>i</sup>* � *<sup>x</sup> <sup>j</sup>*, *<sup>η</sup>* � �, 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 *v*^*,* which is then corrected through a smoothing procedure based on the values of the neighboring fluid particles.

$$\boldsymbol{\nu}\_{i} = (\mathbf{1} - \boldsymbol{\phi}\_{v})\boldsymbol{\hat{\nu}}\_{i} + \boldsymbol{\phi}\_{v}\frac{\sum\_{j}\frac{m\_{j}}{\rho\_{j}}\boldsymbol{\hat{\nu}}\_{j}\boldsymbol{\mathcal{W}}\_{ij}}{\sum\_{j}\frac{m\_{j}}{\rho\_{j}}\boldsymbol{\mathcal{W}}\_{ij}}\tag{5}$$

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 *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 *<sup>μ</sup><sup>T</sup>* <sup>¼</sup> *<sup>c</sup><sup>μ</sup> <sup>k</sup>*<sup>2</sup> *<sup>ε</sup>* in Eq. (4) was evaluated through a *k-ε* model by [74]:

$$\begin{cases} \frac{Dk\_i}{Dt} = P\_{k\_i} + \frac{1}{\sigma\_k} \sum\_j m\_j \frac{\nu\_{T\_i} + \nu\_{T\_j}}{\rho\_i + \rho\_j} \frac{k\_i - k\_j}{r\_{\vec{\eta}}^2 + 0.01h^2} r\_{\vec{\eta}} \cdot \nabla W\_{\vec{\eta}} - \varepsilon\_i\\ \frac{D\varepsilon\_i}{Dt} = \frac{1}{\sigma\_\varepsilon} \sum\_j m\_j \frac{\nu\_{T\_i} + \nu\_{T\_j}}{\rho\_i + \rho\_j} \frac{\varepsilon\_i - \varepsilon\_j}{r\_{\vec{\eta}}^2 + 0.01h^2} r\_{\vec{\eta}} \cdot \nabla W\_{\vec{\eta}} + + C\_{\varepsilon\_1} \frac{\varepsilon\_i}{k\_i} P\_{k\_i} + C\_{\varepsilon\_2} \frac{\varepsilon\_i}{k\_i} \sum\_j \frac{m\_j}{\rho\_j} \varepsilon\_j W\_{\vec{\eta}} \end{cases} (6)$$

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 *ν<sup>T</sup>* is the eddy viscosity and *rij* ¼ *x<sup>i</sup>* � *x <sup>j</sup>*. 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<sup>μ</sup>* = 0.09, *σ<sup>k</sup>* ¼ 1, *σε* = 1.3, *C<sup>ε</sup>*<sup>1</sup> = 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" 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.
