Comparison between the Lagrangian and Eulerian Approach for Simulating Regular and Solitary Waves Propagation, Breaking and Run-Up

Abstract

The present paper places emphasis on the most widely used Computational Fluid Dynamics (CFD) approaches, namely the Eulerian and Lagrangian methods each of which is characterized by specific advantages and disadvantages. In particular, a weakly compressible smoothed particle (WCSPH) model, coupled with a sub-particle scale (SPS) approach for turbulent stresses and a new depth-integrated non-hydrostatic finite element model were employed for the simulation of regular breaking waves on a plane slope and solitary waves transformation, breaking and run-up. The validation of the numerical schemes was performed through the comparison between numerical and experimental data. The aim of this study is to compare the two modeling methods with an emphasis on their performance in the simulation of hydraulic engineering problems.

1. Introduction

Nowadays, numerical methods are becoming the main tool used by scientists and researchers to study fluid dynamics problems. Even if experimental investigations provide a lot of reliable data useful to calibrate numerical models, experimentation requires a lot of time to be planned and executed, usually at a high economical cost. Moreover, the results of experimental tests are limited to just a fixed number of points, while computational methods can provide the desired quantities in any point of the computational domain, allowing a better understanding of the physical phenomena under investigation.

The numerical models for the study of a fluid motion could be based on two different approaches, the Eulerian or the Lagrangian method. The main difference between the two methodologies is that while a Eulerian method discretizes space into a fixed grid and finds the unknown variables at the nodal points, the latter discretizes the continuum into a certain number of moving nodes. Moreover, while in a Lagrangian approach the properties of the single fluid particle are studied in accordance with the initial spatial and temporal position, the Eulerian approach represents the distribution of these properties in a domain without referring to the history of the particle’s trajectories. For example, from the experimental point of view, the PIV (Particle Image Velocimetry) technique allows reconstruction of the velocity vectors that occupy an instantaneous velocity field, based on the average particle motion in space [1,2,3]. Instead, from a Lagrangian viewpoint, PTV (Particle Tracking Velocimetry) traces the pathway of an individual particle from a sequence of images in a system. This method is better than PIV for handling non-stationary flow [4].

Eulerian approaches or grid-based methods such as the finite element method (FEM), finite difference method (FDM) and finite volume method (FVM) have been applied to a wide range of engineering applications such as breaking waves [5,6,7], eddies [8], hydraulic jumps [9,10,11,12], sediment transport [13] and current and oil drift modeling [14,15,16,17].

Despite their great success, grid-based methods suffer from some difficulties, in particular for those cases that present discontinuities or whose mesh needs to be modified during run time according to topological changes.

For these reasons, during the last few years, research has been focused on the development of Lagrangian and meshless techniques, such as the discrete element method (DEM) [18], immersed particle method (IPM) [19], smoothed particle hydrodynamics (SPH) [20] and finite volume particle method (FVPM) [21].

Lagrangian methods have been successfully applied in different fields of fluid dynamics, such as free surface flows [22], the breaking and impact of waves [23,24,25,26,27], multiphase flows [28,29], sediment transport [30,31,32], fluid-structure interaction problems [33,34,35], oscillating jets that induce breaking waves [36,37,38,39] and hydraulic jumps [40,41,42,43,44,45,46].

New research challenges in Computational Fluid Dynamics (CFD) are currently evolving; in particular, recently, research has been focused on coupling the Eulerian– Lagrangian techniques to combine the advantages of the individual models in a single model, thus increasing the accuracy, efficiency and regime of validity. For example, the SPH code is also being coupled with the discrete element Method [47] or with the distributed contact discrete element method (DCDEM) [48] for fluid–solid modeling.

The flow is described by the Navier-Stokes Equations and can be solved by using the two main classes of methods: finite element methods and meshfree methods. Both methodologies discretize a partial differential system into a set of algebraic equations that can be easily solved by computer implementation.

To simulate a turbulent flow, the RANS( Reynolds-Averaged Navier–Stokes) modeling approach has been widely employed leading to the Reynolds-averaged N-S, where the instantaneous values of velocity and pressure are defined as the sum of a mean and fluctuating component [49]. The simplest RANS models are the eddy viscosity models (based on the Boussinesq hypothesis) where the Reynolds stresses are modeled using an eddy (or turbulent) viscosity ν_t. The eddy viscosity can be determined from a turbulent time-scale (or velocity scale) and a turbulent length-scale using the turbulent kinetic energy (k) and the turbulence dissipation rate (ε).

The most commonly-applied turbulence model for its reasonable accuracy to estimate ν_t for a wide range of flows, is the two-equation k-ε model [50]. The two-equation k-ε model has been utilized in Shao [51,52] and De Padova et al. [53], showing a good comparison between the numerical and experimental turbulent quantities.

Another more recent approach is the Large Eddy Simulation (LES). The main difference between the RANS and LES approaches is that the N-S Equations are, in the former, time-averaged and, in the latter, space filtered.

Within Eulerian methods, the RANS approach is the most widely used, followed by LES due to its lower computational cost. Recently, hybrid RANS-LES methods, such as detached eddy simulation (DES), are growing in popularity [54].

The purpose of this paper is to show the versatility of a grid-based and meshless method to simulate free surface hydrodynamic problems. In Section 2, focus is put on giving the basic concepts of both numerical formulations; in Section 3, as an application, both a weakly compressible smoothed particle (WCSPH) model, coupled with a sub-particle scale (SPS) turbulence model and a non-hydrostatic discontinuous/continuous Galerkin model (non-hydrostatic D/C Galerkin) are used to study the physics of regular and solitary waves propagation, breaking and run-up.

2. Lagrangian and Eulerian Methods

2.1. SPH Formulation

SPH is a Lagrangian meshfree particle method developed by Lucy [55] and Gingold and Monaghan [56]; the reader is referred to [57,58,59,60] for a general description of SPH.

In the present paper, a WCSPH model, coupled with a sub-particle scale (SPS) approach to modeling turbulence [61] has been used. Description of the WCSPH can be found in Monaghan [62] and Gomez-Gasteira et al. [63]. WCSPH has been widely used for many applications such as wave breaking [23,25], dam break flows [64,65], wave generation and absorption [66], hydraulics jumps [42] and the interaction between jets and waves [38]. The alternative approach, developed by Lo and Shao [67,68], is defined as incompressible SPH. The reader can find a comparison between the two algorithms in [69,70,71]. The main difference between the two approaches exists in the way the pressures are calculated: the former explicitly uses an equation of state [23,62], while the latter solves a Poisson pressure equation with a source term as a function of velocity divergence or time variation density.

The motion is represented by the Navier-Stokes Equations for a weakly compressible SPH semi-discrete form:

$$\left\{\begin{array}{c}\frac{d\rho }{dt}=-\rho \nabla \mathit{v}\\ \frac{d\mathit{v}}{dt}=-\frac{1}{\rho }\nabla P+\Gamma +\mathit{f}\end{array}\right.$$

(1)

where $\rho$ is the density, v is the velocity vector, P is the pressure, Γ denotes the dissipation terms and f represents accelerations due to external forces, such as gravity. In the SPH formalism, the discrete form of the continuity equation at a point i can be written as follows [72]:

$$\frac{d{\rho }_i}{dt}={\rho }_i{\displaystyle \sum }_j\frac{m_j}{{\rho }_j}\left({\mathit{v}}_i-{\mathit{v}}_j\right)·{\nabla }_{i}W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)+D_i$$

(2)

The summations are extended to all of the particles j at a distance from i smaller than 2h (h defined as the smoothing length), i.e., lying within the circle where the adopted cubic-spline kernel function $W$ [60] is defined. The term $D_i$ represents a numerical diffusive term [73], which can be introduced in order to reduce density fluctuation. The general form of a density diffusion term is:

$$D_i=\delta h{c}_i{\displaystyle \sum }_j{\psi }_{ij}·\nabla W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)V_j$$

(3)

where $\delta$ is the Delta-SPH coefficient, which controls the magnitude of the diffusion term, c_i is the numerical speed of sound, V_j is the associated volume of the j-th particle and ${\psi }_{ij}$ is the artificial dissipation term; in the present paper, the artificial dissipation term proposed by Molteni and Colagrossi [74] was chosen.

The discrete form of the momentum equation, written in an SPH formalism, can be written as it follows [71]:

$$\frac{d{v}_i}{dt}=-{\displaystyle \sum }_j{m}_j\left(\frac{P_i+P_j}{{\rho }_i{\rho }_j}\right)\nabla W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)+⟨{\Gamma }_i⟩+\mathit{f}$$

(4)

The momentum dissipation term Γ_i is obtained coupling the viscous dissipation in the laminar regime, as approximated by Lo and Shao [68], with a sub-particle scale (SPS) model [71]. The former can be expressed with the following formula:

$${\nu }_o{\nabla }^2{\mathit{v}}_i={\displaystyle \sum }_j{m}_j\frac{4{\nu }_0({\mathit{x}}_i-{\mathit{x}}_j)·\nabla W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)}{\left({\rho }_i+{\rho }_j\right)\left({\left({\mathit{x}}_i-{\mathit{x}}_j\right)}^2+{\eta }^2\right)}\left({\mathit{v}}_i-{\mathit{v}}_j\right)$$

(5)

where 𝜐₀ is the kinematic viscosity and $\eta$ is a parameter that guarantees a non-singular operator. In the DualSPHysics solver, $\eta$ is equal to 0.001h², with $\left\{\eta \in ℝ;x_i-x_j>\eta \right\}$.

The SPS contribution to the momentum dissipation can be expressed as follows [71]:

$$\frac{1}{\rho }\nabla {\tau }_{ij}^{*}={\displaystyle \sum }_j{m}_j\left(\frac{{\tau }_i^{*}+{\tau }_j^{*}}{{\rho }_i{\rho }_j}\right){\nabla }_{i}W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)$$

(6)

where τ* is the sub-particle scale stress tensor, which can be modelled by an eddy viscosity closure as:

$$\frac{{\tau }^{ij}}{\rho }=2{\nu }_t\left(S^{ij}-\frac{1}{3}{S}^{ij}{\delta }^{ij}\right)-\frac{2}{3}{C}_l{\mathsf{\Delta }}^2{\delta }^{ij}{\left|S^{ij}\right|}^2$$

(7)

where ${\nu }_t={\left(C_s\mathsf{\Delta }\right)}^2{\left|S^{ij}\right|}^2$, $C_s$ = 0.12 is the Smagorinsky constant, $C_l$ = 0.00066, $\mathsf{\Delta }$ is the initial particle spacing and $\left|S^{ij}\right|=1/2{\left(2S^{ij}{S}^{ij}\right)}^{1/2}$ is the local strain rate. A more detailed description of the LES-SPS model using Favre averaging in a weakly compressible approach can be found in [23].

2.2. Non-Hydrostatic Discontinuous/Continuous Galerkin Model

The Eulerian grid-based method used in this paper is the non-hydrostatic D/C Galerkin, presented in Calvo et. al [7]. The model is a two-dimensional depth-integrated non-hydrostatic model with the capacity of simulating wave propagation, breaking and run-up in a finite element framework.

Eulerian grid-based depth-integrated models using the Boussinesq-type Equations (BTEs) have been traditionally utilized to simulate wave propagation, but the high order partial derivative terms included in BTEs cause numerical instabilities The shallow water models with non-hydrostatic pressure distribution have shown their capacity for accurate modeling of nonlinear and dispersive waves since their introduction by Casulli and Stelling [75] and Stansby and Zhou [76]. The vertical momentum is considered in these models by the introduction of a non-hydrostatic pressure term into the Reynolds-averaged Navier-Stokes Equations. The non-hydrostatic models for water waves are classified as either single-layer models (two-dimensional depth-integrated non-hydrostatic) or multilayer models (three-dimensional non-hydrostatic), depending on the number of layers in the vertical discretization. The capabilities of the non-hydrostatic models for water waves with single or multiple layers in the vertical direction has been researched by Stelling and Zijlema [77]; Zijlema et al. [78]; Zijlema and Stelling [79] and Zijlema et al. [80], among others.

The non-hydrostatic D/C Galerkin model is constructed using a combination of discontinuous and continuous Galerkin methods, where the depth-integrated non-hydrostatic Equations are separated into hydrostatic and non-hydrostatic parts. The hydrostatic part corresponds to the depth-integrated shallow water

Equations and is solved with a discontinuous Galerkin FEM for the simulation of discontinuous flows, wave breaking and run-up. The non-hydrostatic part results in a Poisson-type equation where the non-hydrostatic pressure is resolved by a continuous Galerkin finite element method for the modeling of wave propagation and transformation. The model utilizes linear quadrilateral finite elements for horizontal velocities, water surface elevations and non-hydrostatic pressures that permit local refinement and complex boundaries.

The governing Equations in the conservation form for depth-integrated, non-hydrostatic flow in the Cartesian coordinates system (x, y and z) are:

$$\frac{\partial \left(HU\right)}{\partial t}+\frac{\partial \left(HUU\right)}{\partial x}+\frac{\partial \left(HUV\right)}{\partial y}=-gH\frac{\partial \xi }{\partial x}-\frac{g{n}^{2}U\sqrt{U^2+V^2}}{H^{\frac{1}{3}}}-\frac{H}{2\rho }\frac{\partial q_b}{\partial x}-\frac{q_b}{2\rho }\frac{\partial \left(\xi -h\right)}{\partial x}$$

(8)

$$\frac{\partial \left(HV\right)}{\partial t}+\frac{\partial \left(HUV\right)}{\partial x}+\frac{\partial \left(HVV\right)}{\partial y}=-gH\frac{\partial \xi }{\partial y}-\frac{g{n}^{2}V\sqrt{U^2+V^2}}{H^{\frac{1}{3}}}-\frac{H}{2\rho }\frac{\partial q_b}{\partial y}-\frac{q_b}{2\rho }\frac{\partial \left(\xi -h\right)}{\partial y}$$

(9)

$$\frac{\partial H}{\partial t}+\frac{\partial \left(UH\right)}{\partial x}+\frac{\partial \left(VH\right)}{\partial y}=0$$

(10)

$$\frac{\partial W}{\partial \mathrm{t}}=\frac{q_b}{\rho H}$$

(11)

where $U$,$V$ and $W$ are the depth-averaged velocity components in the $x$, $y$ and z directions; $\rho$ is the water density; $n$ is Manning’s roughness coefficient; $g$ is the gravitational acceleration. The flow depth is $H=\xi +h$, where $\xi$ is the surface elevation measured from the still water level and $h$ is the water depth measured from the still water level. A linear distribution is assumed in the vertical direction for both the non-hydrostatic pressures and for the vertical velocities. The non-hydrostatic pressure on the free surface is assumed as zero and, in the bottom, as $q_b$. The effects of the turbulence, Coriolis force, atmospheric pressure and baroclinic pressure gradient, are traditionally ignored in the study of non-hydrostatic wave propagation. However, simulation with these effects is possible; see, e.g., Stelling and Busnelli [81]. The average vertical velocity, $W$, is $(w_{\xi }+w_b)/2$, where $w_{\xi }$ is the vertical velocity at the surface and $w_b$ is the vertical velocity at the bottom. The kinematic boundary conditions of the free surface and the bottom are:

$$w_{\xi }=\frac{\partial \xi }{\partial t}+u_{\xi }\frac{\partial \xi }{\partial x}+v_{\xi }\frac{\partial \xi }{\partial y}$$

(12)

$${\mathrm{w}}_{\mathrm{b}}=-{\mathrm{u}}_{\mathrm{b}}\frac{\partial \mathrm{h}}{\partial \mathrm{x}}-{\mathrm{v}}_{\mathrm{b}}\frac{\partial \mathrm{h}}{\partial \mathrm{y}}$$

(13)

with $u_{\xi },v_{\xi },u_b and v_b$ being the components in x and y of the velocities next to the free surface and the bottom, respectively. Due to the linearization of the vertical momentum equation, it is well recognized that the depth-integrated non-hydrostatic models are only appropriate for weakly nonlinear cases [82].

The solution begins by solving the horizontal momentum Equations (8) and (9), without the non-hydrostatic pressure terms, by means of a discontinuous Galerkin method. These horizontal momentum Equations, coupled with the mass conservation Equation (10), can be written in the conservative form of Equation (14).

$$\frac{\partial \mathit{U}}{\partial t}+\nabla .\mathit{F}\left(\mathit{U}\right)=\frac{\partial \mathit{U}}{\partial t}+\frac{\partial \mathit{E}\left(\mathit{U}\right)}{\partial x}+\frac{\partial \mathit{G}\left(\mathit{U}\right)}{\partial y}=\mathit{S}\left(\mathit{U}\right)$$

(14)

The vectors of conserved variables $\mathit{U}$, source vector $\mathit{S}$, and flux vectors $\mathit{F}\left(\mathit{U}\right)$ are given in Equations (15) to (18).

$$\mathit{U}=\left\{\begin{array}{c}{q}_x\\ q_y\\ H\end{array}\right\}$$

(15)

$$\mathit{S}=\left\{\begin{array}{c}-gH\frac{\partial \xi }{\partial x}-\frac{g{n}^{2}U\sqrt{U^2+V^2}}{H^{1/3}}\\ -gH\frac{\partial \xi }{\partial y}-\frac{g{n}^{2}V\sqrt{U^2+V^2}}{H^{1/3}}\\ 0\end{array}\right\}$$

(16)

$$\mathit{F}=\left(\mathit{E}\left(\mathit{U}\right),\mathit{G}\left(\mathit{U}\right)\right)$$

(17)

$$\mathit{E}\left(\mathit{U}\right)=\left\{\begin{array}{c}{q}_x^2/H\\ q_x{q}_y/H\\ q_x\end{array}\right\},\mathit{G}\left(\mathit{U}\right)=\left\{\begin{array}{c}{q}_x{q}_y/H\\ q_y^2/H\\ q_y\end{array}\right\}$$

(18)

In these Equations, $q_x$ and $q_y$ are discharges per unit width in the $x$ and $y$ directions and are equal to $HU$ and $HV$, respectively.

The discontinuous Galerkin formulation of the governing Equation (14) is obtained by multiplying it with a shape function $\phi$ and integrate over an element, ${\Omega }_e$. The flux term $\mathit{F}$ is integrated using Gauss theorem, resulting in Equation (19). In this equation, $\mathit{n}=\left(n_x,n_y\right)$ is the outward unit normal vector at an element boundary ${\Gamma }_e$.

$${{\displaystyle \int }}_{{\Omega }_e}^{}\phi \frac{\partial \mathit{U}}{\partial t}d\Omega +{{\displaystyle \int }}_{{\Gamma }_e}^{}\phi \mathit{F}.\mathit{n} d\Gamma -{{\displaystyle \int }}_{{\Omega }_e}^{}\mathit{F}.\nabla \phi d\Omega ={{\displaystyle \int }}_{{\Omega }_e}^{}\phi \mathit{S} d\Omega$$

(19)

In the discontinuous Galerkin method, the variable vector $\mathit{U}$ is approximated over a quadrilateral element by:

$$\mathit{U}\left(x, y\right)\approx {\displaystyle \sum }_{\mathit{j}=\mathbf{1}}^{\mathbf{4}}{U}_j{\phi }_j\left(x, y\right)$$

(20)

where $U_j$ are nodal values of the variables and ${\phi }_j\left(x, y\right)$ are the bilinear approximation functions of the solution variables or shape functions. Since in a discontinuous Galerkin context the discontinuities of variables at the element boundaries are allowed, the intercell flux is taken as being dependent on the values $\mathit{U}$ in each of the two adjacent elements. Thus, the normal flux $\mathit{F}.\mathit{n}$ is not exclusively defined and is substituted by a numerical flux $\tilde{\mathit{F}}\left({\mathit{U}}_{L,}{\mathit{U}}_R\right)$ where ${\mathit{U}}_L$ and ${\mathit{U}}_R$ are the variables at the left (internal) and right (external) sides of the element boundary in the counterclockwise direction, respectively. Consequently, the second integral in Equation (19) is expressed as:

$${{\displaystyle \int }}_{{\mathit{\Gamma }}_{\mathit{e}}}^{}\phi \mathit{F}.\mathit{n} d\Gamma ={{\displaystyle \int }}_{{\mathit{\Gamma }}_{\mathit{e}}}^{}\phi \tilde{\mathit{F}}d$$

(21)

In this analysis, $\tilde{\mathit{F}}$ is chosen to be the common Harten-Lax-van Leer (HLL) numerical flux [83]. The first step of the solution process ends with the intermediate estimation of the conserved variables $q_x$ and $q_y$ from Equation (17): ${\tilde{q}}_x^{n+1} and {\tilde{q}}_y^{n+1}$.

In the second step of the solution, a Poisson equation is constructed that is then solved by a continuous Galerkin method to obtain the non-hydrostatic pressures $q_b^{n+1}$. A time approximation of the vertical moment Equation in (11) is:

$$w_{\xi }^{n+1}=w_{\xi }^n-w_b^{n+1}+w_b^n+\frac{2\mathsf{\Delta }t{q}_b^{n+1}}{\rho H^n}$$

(22)

The vertical velocity at the bottom is valued from the kinematic boundary condition (16), estimated as:

$$w_b^{n+1}=-U^n\frac{\partial h}{\partial x}-V^n\frac{\partial h}{\partial y}$$

(23)

$$w_b^n=-U^{n-1}\frac{\partial h}{\partial x}-V^{n-1}\frac{\partial h}{\partial y}$$

(24)

From the remaining part of the horizontal momentum, Equations (8) and (9) including the non-hydrostatic pressure terms:

$$\frac{\partial \left(HU\right)}{\partial t}=-\frac{H}{2\rho }\frac{\partial q_b}{\partial x}-\frac{q_b}{2\rho }\frac{\partial \left(\xi -h\right)}{\partial x}$$

(25)

$$\frac{\partial \left(HV\right)}{\partial t}=-\frac{H}{2\rho }\frac{\partial q_b}{\partial y}-\frac{q_b}{2\rho }\frac{\partial \left(\xi -h\right)}{\partial y}$$

(26)

The final horizontal velocities, influenced by the non-hydrostatic pressure, can be time approximated as

$$U_c^{n+1}=\tilde{U}{ }^{n+1}-\frac{\mathsf{\Delta }t}{2\rho }\frac{\partial q_b^{n+1}}{\partial x}-\frac{\mathsf{\Delta }t}{2\rho }\frac{q_b^{n+1}}{H_c^n}\frac{\partial }{\partial x}\left({\xi }_c^n-h\right)$$

(27)

$$V_c^{n+1}=\tilde{V}{ }^{n+1}-\frac{\mathsf{\Delta }t}{2\rho }\frac{\partial q_b^{n+1}}{\partial y}-\frac{\mathsf{\Delta }t}{2\rho }\frac{q_b^{n+1}}{H_c^n}\frac{\partial }{\partial y}\left({\xi }_c^n-h\right)$$

(28)

In the continuous Galerkin solution of the non-hydrostatic pressures, the continuous horizontal velocities $U_c^{n+1}$ and $V_c^{n+1}$, continuous water level ${\xi }_c^n$ and flow depth $H_c^n$, vertical velocities $w_{\xi }$ and $w_b,$ and non-hydrostatic pressures $q_b^{n+1}$ are approximated with nodal variables, which have the same value in all the adjoining elements. The continuous intermediate velocities $\tilde{U}{ }^{n+1}$ and $\tilde{V}{ }^{n+1}$ are obtained from the intermediate evaluation of the discontinuous discharges ${\tilde{q}}_x^{n+1} and {\tilde{q}}_y^{n+1}$ integrating over the entire domain $\Omega$.

$${{\displaystyle \int }}_{\Omega }^{}\phi \tilde{U}{ }^{n+1}d\Omega ={{\displaystyle \int }}_{\Omega }^{}\phi {\tilde{U}}_{dc}^{n+1}d\Omega$$

(29)

$${{\displaystyle \int }}_{\Omega }^{}\phi \tilde{V}{ }^{n+1}d\Omega ={{\displaystyle \int }}_{\Omega }^{}\phi {\tilde{V}}_{dc}^{n+1}d\Omega$$

(30)

where the discontinuous intermediate velocities ${\tilde{U}}_{dc}^{n+1}$ and ${\tilde{V}}_{dc}^{n+1}$ are approximated in an element using the values of $\frac{{\tilde{q}}_x^{n+1}}{H^n}$ and $\frac{{\tilde{q}}_y^{n+1}}{H^n}$, respectively, evaluated at the nodes. Similarly, the continuous water level elevations ${\xi }_c^n$ are obtained from:

$${{\displaystyle \int }}_{\Omega }^{}\phi {\xi }_c^{n}d\Omega ={{\displaystyle \int }}_{\Omega }^{}\phi {\xi }^{n}d\Omega$$

(31)

where ${\xi }^n=H^n-h$. Finally, the continuous flow depth, $H_c^n$, is $H_c^n={\xi }_c^n+h$. To obtain a correct solution between the velocity field and the non-hydrostatic pressures, the continuity equation is applied directly to the water column

$$\left(\frac{\partial U_c^{n+1}}{\partial x}+\frac{\partial V_c^{n+1}}{\partial y}\right)+\left(\frac{w_{\xi }^{n+1}-w_b^{n+1}}{H_c^n}\right)=0$$

(32)

Substituting Equations (22), (23) and (25) in Equation (28), applying Gauss’s theorem in some terms, and neglecting others, a Poisson equation is established from which the hydrostatic pressure in the bottom, $q_b^{n+1}$, is obtained.

A continuous Galerkin finite element model to solve the Poisson equation on an element ${\Omega }_e$ is

$$\begin{gathered} {{\displaystyle \int }}_{{\Omega }_e}\phi \frac{2\Delta t{q}_b^{n+1}}{\rho {\left(H_c^n\right)}^2}d{\Omega }_{}+{{\displaystyle \int }}_{{\Omega }_e}\frac{\partial \phi }{\partial x}\frac{\Delta t}{2\rho }\frac{\partial q_b^{n+1}}{\partial x}d{\Omega }_{}+{{\displaystyle \int }}_{{\Omega }_e}\frac{\partial \phi }{\partial y}\frac{\Delta t}{2\rho }\frac{\partial q_b^{n+1}}{\partial y}d{\Omega }_{} \\ -{{\displaystyle \int }}_{{\Omega }_e}\phi \frac{\Delta t}{2\rho H_c^n}\frac{\partial q_b^{n+1}}{\partial x}\frac{\partial }{\partial x}\left({\xi }_c^n-h\right)d{\Omega }_{}-{{\displaystyle \int }}_{{\Omega }_e}\phi \frac{\Delta t}{2\rho H_c^n}\frac{\partial q_b^{n+1}}{\partial y}\frac{\partial }{\partial y}\left({\xi }_c^n-h\right)d{\Omega }_{} \\ = -{{\displaystyle \int }}_{{\Omega }_e}\phi \begin{array}{c}\left(\frac{\partial \tilde{U}{ }^{n+1}}{\partial x}+\right.\end{array} \frac{\partial \tilde{V}{ }^{n+1}}{\partial x} \left.+\frac{\left(w_{\xi }^n+w_b^n-2w_b^{n+1}\right)}{H_c^n}\right) d{\Omega }_{} \end{gathered}$$

(33)

After all of the variables in Equation (29) have been approximated, as in Equation (20), and all of the elements of the domain have been assembled, the resulting linear system is solved to obtain the non-hydrostatic pressures, $q_b^{n+1}$.

In the last step of the solution, once the non-hydrostatic pressures $q_b^{n+1}$ are known, the Equations (27) are solved on an element to obtain the final solutions for the discharges $q_x^{n+1}$ and $q_y^{n+1}$. A Galerkin finite element model to solve Equation (24) on an element ${\Omega }_e$ is:

$${{\displaystyle \int }}_{{\Omega }_e}^{}\phi q_x^{n+1}d\Omega ={{\displaystyle \int }}_{{\Omega }_e}^{}\phi {\tilde{q}}_x^{n+1}d\Omega -{{\displaystyle \int }}_{{\Omega }_e}^{}\phi \frac{\mathsf{\Delta }t{H}^n}{2\rho }\frac{\partial q_b^{n+1 }}{\partial x}d\Omega -{{\displaystyle \int }}_{{\Omega }_e}^{}\phi \frac{\mathsf{\Delta }t{q}_b^{n+1 }}{2\rho }\frac{\partial \left({\xi }^n-h\right)}{\partial x}d\Omega$$

(34)

$${{\displaystyle \int }}_{{\Omega }_e}^{}\phi q_y^{n+1}d\Omega ={{\displaystyle \int }}_{{\Omega }_e}^{}\phi {\tilde{q}}_y^{n+1}d\Omega -{{\displaystyle \int }}_{{\Omega }_e}^{}\phi \frac{\mathsf{\Delta }t{H}^n}{2\rho }\frac{\partial q_b^{n+1 }}{\partial y}d\Omega -{{\displaystyle \int }}_{{\Omega }_e}^{}\phi \frac{\mathsf{\Delta }t{q}_b^{n+1 }}{2\rho }\frac{\partial \left({\xi }^n-h\right)}{\partial y}d\Omega$$

(35)

Finally, with the discharges $q_x^n$, $q_y^n$, $q_x^{n+1}$, $q_y^{n+1}$ and the flow depth $H^n$ known, the unknown flow depth $H^{n+1}$ is determined from the continuity Equation (13) using the discontinuous Galerkin formulation in Equation (19). Details of the entire solution process can be found in Calvo et. al [7]. After each of the three solution steps, a special slope limiter for quadrilateral elements is applied. Details of the developed slope limiter for quadrilateral elements can also be found in Calvo et. al. [7].

3. Applications

A total number of four cases of breaking waves have been simulated in the present study, named T1, T2, T3 and T4, respectively. Among these, there is one case of spilling/plunging breaker (T1), one case of plunging breaker (T2), one case of a solitary wave propagation over a fringing reef (T3) and one case of a solitary wave run-up on a beach of constant slope (T4).

The validity of the numerical schemes described above, has been checked against the experimental observations. In particular, the first two data sets (T1 and T2) were based on experiments by De Serio and Mossa [84]; the third data set was chosen out of the experiments by Roeber [85] and Roeber et al. [86]; the last data set (T4) was chosen out of the experiments by Titov and Synolakis [87].

3.1. Experimental Set Up

3.1.1. Regular Breaking Waves on a Plane Beach

The experiments by De Serio and Mossa [84] were carried out in a wave flume 45 m long and 1 m wide, located at the Department of Civil, Environmental, Land, Building Engineering and Chemistry of the Polytechnic University of Bari, Italy (Figure 1). The flume has an initial flat bottom of 12 m extending from the wave paddle to section 73, while from section 73 up to the shoreline the bottom has a 1/20 slope, obtained with a wooden panel. The wavemaker is a piston-type, whose paddles produce the desired wave by moving the water mass according to a proper input signal. The wave elevations were measured with a resistance probe, while the instantaneous Eulerian velocities were acquired by a backscatter, two-component, four beam Laser Doppler Anemometer (LDA) system and a signal processor (Figure 1).

Table 1. Experimental parameters of the analysed regular waves.

Test	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

shows the main parameters of the two examined waves, in particular, the offshore wave height H₀, the wave period T, the offshore wavelength L₀, the Irribarren number ξ₀, which has been estimated in section 76, located where the bottom is flat with water depth d equal to 0.70 m. Based on the Irribarren number ξ₀, the two regular tested waves were characterized by a spilling/plunging and plunging breakers, respectively.

The sketch of Figure 2 shows the different sections named 76, 55, 49, 48, 47, 45, which have been used to make the comparisons between the simulated and measured free surface elevations.

3.1.2. Solitary Waves

The experiment by Roeber [85] and Roeber et al. [86] was performed in a wave channel 48.8 m long and 2.16 m wide, equipped with a piston-type wavemaker for solitary wave generation. The test involves a steep solitary wave with a height A = 0.5 m and a water depth of h = 1 m, resulting in A/h = 0.5. The wave was transformed over a fringing reef with a 1:5 slope onto a dry reef flat. Further details on this experimental test can be found in Roeber [85] and Roeber et al. [86].

The experimental study by Titov and Synolakis [87], reproduced a solitary wave run-up with wave height A/h = 0.3 (where A is the solitary wave height and h is the still water depth) on a plane beach with a slope of 1:19.85. The solitary wave had a height A = 1 m with a water depth of 3.33 m, resulting in A/h = 0.3.

3.2. Numerical Parameters

In this section, the parameters used for all four tests (T1-T4) in both numerical approaches have been described.

3.2.1. SPH Parameters

The SPH results discussed here were obtained using the hardware accelerated DualSPHysics code [88]. Since the first release back in 2011, DualSPHysics has been shown to be robust and accurate for simulating free surface flows yet requiring high computational cost.

Therefore, recently, the high-performance computing of SPH has mainly focused on Graphical Processing Units (GPUs) [89], which are superior in terms of price and energy consumption compared to traditional CPUs. On DualSPHysics, the C++ programming language is used to code the SPH formulation for CPU execution, while GPU executions are based on NVIDIA’s CUDA architecture [72]. The specific features of the open-source code are detailed in [72,89].

For the first two tests (T1–T2), following the experimental setting [84], the numerical wave tank was 22.5 m long and 0.97 m high, while the initial water depth was equal to 0.70 m. In De Padova et al., [27] the convergence analysis for T1 showed that only d/Δx ≥ 35 (Δx≤ 0.02 m) guaranteed the independence of the SPH result from the resolution and yielded a result in accordance with the experiments. However, in this study, simulations were performed by choosing a higher resolution than shown by De Padova et al. [25,27], in order to show the efficient and reliable use of DualSPHysics code executed on a GPU architecture also with a large number of particles. In particular, the 2D flow (T1 and T2) was simulated by discretizing the computational domain through a particle distribution with initial particle distance Δx =Δz varying from to 0.02 to 0.003 m. Figure 3 shows that d/Δx ≥ 35 (Δx≤ 0.02 m) guaranteed the independence of the SPH result from the resolution. This again confirms the results by [27]. Following, the results and performance analysis for T1 and T2 were carried out considering an initial particle distance Δx = Δz = 0.01 m.

In all four tests T1-T4, the solid boundary conditions needed for idealizing the seabed, are discretized by a set of boundary particles that differ from the fluid particles. When a fluid particle approaches the boundary, the density, and hence the pressure exerted by the boundary particles, increase. This generates the necessary repulsive force on the water particles.

Instead, the offshore boundary condition is treated as a dynamic boundary condition [90], i.e., the numerical wave paddle is composed of wall particles whose density is computed with the continuity equation and pressure obtained from the equation of state [72]. For both of the regular waves’ tests, the first order numerical wavemaker is located at a distance of 0.5 m offshore on the side of section 76.

For the third test (T3), following the experimental setting by (Roeber [85] and Roeber et al. [86]), the numerical flume was 45 m long, with an initial water depth of 1 m. The fringing reef, whose slope was 1/5, was 5.1 m long and its toe was at a distance of 17 m from the numerical wavemaker. The latter generated the solitary wave according to the Rayleigh solitary wave generation theory; in this case, an initial particle distance of Δx = Δz = 0.02 m with a total number of particles equal to 50,242.

In the last test (T4), following the experimental setting by Titov and Synolakis [87], the numerical wave tank was 154.83 m long, with an initial water depth of 3.33 m. The toe of the inclined beach was at 34.83 m from the piston, which was of the same type as that described for the third case. The initial particle distance was Δx = Δz = 0.025 m, resulting in a total number of 369,885 particles.

For all of the cases, the initial conditions are described by setting the velocities of all the particles of the computational domain u = 0 at time t = 0. Moreover, the pressure distribution is assumed to be hydrostatic at the beginning of the computation. In the SPH computations, the free surface is easily and accurately tracked by the particles. To simplify, a pressure equal to zero is given to each of the surface particles.

According to De Padova et al. [25], the ratio of the smoothing length h to the initial particle spacing Δx was maintained to a constant value of h/Δx = 1.5 for all four simulations (T1-T4).

Table 2. Main characteristics of the SPH simulations.

Test	Time Simulation (s)	Δx(m)	h/Δx	Nparticles (−)
SPH_T1	100	0.01	1.5	62,892
SPH_T2	100	0.01	1.5	62.892
SPH_T3	13	0.02	1.5	50,242
SPH_T4	45	0.025	1.5	369,885

shows the main characteristics of the SPH simulations.

3.2.2. Eulerian Model Parameters

Among the most important parameters in the simulation of wave propagation, breaking and run-up using the non-hydrostatic D/C Galerkin model are the breaking waves parameters. Wave overturning and small-scale processes related to wave breaking cannot be reproduced in depth-averaged models, consequently some kind of closure is required. Usually, this closure involves two steps. The first is a trigger mechanism to locate in space and time the limits of breaking. The second is a mechanism that generates the dissipation of energy in the model.

The dissipation of energy is accomplished by shifting locally to the depth-integrated hydrostatic shallow water Equations, modeling breaking wave fronts as shocks. This procedure was firstly presented by Smit et al. [91] as the hydrostatic front approximation (HFA) method for the modeling of wave breaking.

For the trigger mechanism of wave breaking, two different criteria are used. The first criterion is used for solitary waves and establish that wave breaking occurs when $\frac{\mathsf{\xi }}{\mathrm{h}}>0.8$ [92,93,94]. This criterion is referred to as the local criterion in Bacigaluppi et al. [95]. All of the elements containing a node with the wave breaking condition $\frac{\mathsf{\xi }}{\mathrm{h}}>0.8$ are shifted to the hydrostatic mode by making ${\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }=0, \frac{\partial {\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }}{\partial \mathrm{x}}=0$, $\frac{\partial {\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }}{\partial \mathrm{y}}=0$ and ${\mathrm{w}}_{\mathsf{\xi }}$ = ${\mathrm{w}}_{\mathrm{b}}$ = $0$ in the element.

The second criterion is used for regular waves and was first introduced by Kazolea et al. (2014) [96] as the surface variation criterion. A node is flagged if $\frac{\frac{\partial \mathsf{\xi }}{\partial \mathrm{t}}}{\sqrt{\mathrm{g}\left|\mathrm{h}\right|}}>Ɣ$, with $Ɣ$ ∈ [0.3, 0.65] depending on the type of breaker. Flagged nodes are clustered to form a breaking region. All of the elements containing a flagged node are shifted in the element to hydrostatic mode by making ${\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }=0, \frac{\partial {\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }}{\partial \mathrm{x}}=0$, $\frac{\partial {\mathrm{q}}_{\mathrm{b}}^{\mathrm{n}+1 }}{\partial \mathrm{y}}=0$ and ${\mathrm{w}}_{\mathsf{\xi }}$ = ${\mathrm{w}}_{\mathrm{b}}$ = $0$. More details regarding the implementation of this detection criterion can be found in Bacigaluppi et al. (2019) [95].

The boundaries conditions in all tests are: at wall boundaries the gradients of water elevations and the velocities normal to the boundary are set to zero; at open boundaries, the water elevations are determined by a Sommerfeld radiation condition; at wave inflow boundaries, the water elevations and velocities are calculated with the corresponding inflow wave analytical formula; no boundary conditions are required for the non-hydrostatic pressure.

In the simulation of tests T1 and T2, the experiments of De Serio and Mossa [84] were replicated with a numerical domain of 32 m long and 0.24 m wide and the initial water depth was 0.70 m. Square elements with sizes Δx = Δy = 0.04 m were considered, resulting in 5607 nodes. The surface variation criterion was used as the trigger mechanism of wave breaking. The parameter Ɣ was chosen to replicate the waves before and after the breaking zone (consequently trying to match the total energy dissipation of wave breaking with the HFA method).

In test T3, we analyzed the tests of Roeber [85] and Roeber et al. [86] with a numerical domain of 45 m long and 0.2 m wide and a water depth of 1 m. Square elements with ∆x = ∆y = 0.05 m were used, leading to 4505 nodes. A Manning coefficient of n = 0.01 was established to represent the finished concrete surface of the reef model. The local criterion was utilized as the trigger mechanism of wave breaking.

Finally, in test T4, the experimental setting of Titov and Synolakis [87] was reproduced in a numerical domain of 220 m long and 2.4 m wide and a water depth of 3.33 m. Square elements with sizes Δx = Δy = 0.4 m were considered, resulting in 3857 nodes. A Manning coefficient n = 0.01 was used to define the bed surface roughness. The local criterion was used to activate wave breaking.

Table 3. Main characteristics of the non-hydrostatic D/C Galerkin simulations.

TEST	Time Simulation (s)	Δx, Δy (m)	Wave Breaking Criterion	Nnodes(−)
D/C Galerkin_T1	60	0.04	$\frac{\frac{\partial \xi }{\partial t}}{\sqrt{g\left|h\right|}}>0.2$	5607
D/C Galerkin_T2	60	0.04	$\frac{\frac{\partial \xi }{\partial t}}{\sqrt{g\left|h\right|}}>0.3$	5607
D/C Galerkin_T3	25	0.05	$\frac{\xi }{h}>0.8$	4505
D/C Galerkin_T4	40	0.4	$\frac{\xi }{h}>0.8$	3857

shows the main characteristics of the non-hydrostatic D/C Galerkin simulations.

4. Results and Performance Analysis

The efficiency and performance of both numerical approaches have been analysed in this section. The numerical results have been compared with laboratory experiments by De Serio and Mossa [84] in T1 and T2, by Roeber [85] and Roeber et al. [86] in T3 and Titov and Synolakis [87] in T4.

4.1. Spilling/Plunging and Plunging Breaking Waves on a Plane Beach

Both numerical approaches reproduce the breaking waves quantitatively well in both tests (T1 and T2). In fact, in the horizontal direction of the surf zone, the following zones can be distinguished [97]. Initial wave deformation occurs in what has been termed the shoaling zone (section 76–49), where wave profile is characterized by a rapid change in shape. Subsequently, the wave reaches the breaking point in the outer zone (section 48–47) originating an overturning jet whose strength depends on the type of breaker. In the inner surf zone (section 45), the wave undergoes a gradual transformation into a turbulent bore, until reaching the swash zone. Figure 4a,b show the snapshot of the free surface for both the cases of a spilling/plunging (T1) and a plunging breaking wave (T2).

As an example, in Figure 5a,b, both laboratory and numerical wave surface elevations, at vertical sections 49, 48 and 45, have been plotted for T1 and T2.

Table 4. Statistical analysis of the comparison of the free surface for T1 and T2.

	IW
Sect.	49	48	45
SPH_T1	0.92	0.90	0.89
D/C Galerkin_T1	0.82	0.78	0.72
SPH_T2	0.95	0.89	0.90
D/C Galerkin_T2	0.91	0.85	0.75

shows the statistical analysis of the comparison between the two numerical models with the experimental results in terms of the index $I_W$ proposed by Wilmott [98].

$$I_W=1-\frac{{{\displaystyle \sum }}_{k=1}^N{\left(X_{C_k}-X_{m_k}\right)}^2}{{{\displaystyle \sum }}_{k=1}^N{\left[\left|X_{C_k}-\overline{X_m}\right|+\left|X_{m_k}-\overline{X_m}\right|\right]}^2}$$

(36)

where X_c and X_m are the modeled and measured values, respectively, while the bar denotes the average of the modeled and measured values. It takes into account the relative error among experimental and output values, and it will exhibit values closer to one for higher levels of accordance.

Both codes have been shown to be robust, efficient and reliable and a good agreement has been observed for all sections for the wave elevations. However, in the shallow waters (sections 49–48), the computed Galerkin method results appear to slightly underpredict the wave crest and slightly overestimate the experimental values in the trough. After the breaking wave (section 45), a slight overestimate of the crest data is observed for the Galerkin method in T1 and T2. The statistical parameter in Table 4 confirms these conclusions.

Furthermore, skewness [99], a measure of crest-trough shape, has been computed and shown in Figure 6a,b for both tests (T1–T2). The trend of wave skewness increases as the wave shoals and breaks (sections 48–49) and decreases near the shoreline (sections 47–45). Indeed, for both experimental and computed data, the trend generally becomes more pronounced as the water depth decreases, with higher crests and shallower troughs evident. However, the SPH model shows better results than the Galerkin model in terms of the wave shape. The reason for the underestimation of the breaking waves in the Galerkin model is that the HFA method requires high vertical resolution for capturing the hydrostatic shock of the series of regular breaking waves; high vertical resolution that is obviously not present in depth-averaged models. Smit et al. [91] indicated that a multilayer version of a non-hydrostatic shallow water model would require at least 20 layers to obtain a correct solution of wave height using the HFA method. Similar underestimations are observed in other works using the HFA method in depth-averaged models for the simulation of a breaking wave series [79,95,100,101].

4.2. Solitary Waves Propagation

Both numerical approaches reproduce the principal processes well in both two tests (T3 and T4); in fact, in test T3, as shown by Roeber (2010) [85] and Roeber et al. (2010) [86], the wave evolves from subcritical flow to supercritical flow over the reef edge at x = 22 m, collapses at x = 23 m, and, subsequently, rushes over the initially dry reef flat (Figure 7).

Additionally, in T4, both of the numerical approaches reproduce the wave breaking and run-up processes well. In fact, both of the numerical results confirm the experimental findings by Titov and Synolakis [87]; in fact, the breaking of waves occurs between t(g/h)½ = 20 and 25 and the maximum run-up was observed at t(g/h)½ = 45 (Figure 8).

4.3. Performance Analysis

Both codes have been shown to be robust, efficient and reliable in simulating the analysed phenomena. However, in terms of computing cost, the Eulerian method was faster than the Lagrangian method. In fact, tracking sufficient number of particles and analyzing them in the Lagrangian approach required more time than the Eulerian iterations.

With the large numbers of neighboring particle interactions and the time step restricted by the Courant condition for weakly compressible flow, SPH is computationally demanding.

However, in this work the power computing of GPUs was used to accelerate DualSPHysics when compared with the performance achieved using a classic version of SPHysics code [25].

This important acceleration of the code using the new GPU technology can be observed in Figure 9 where the runtime of the hardware accelerated DualSPHysics code is shown by comparing its performance against the classic version of SPHysics code using a single core. In particular, in De Padova et al. [25], the open-source code SPHysics [23], was executed up to t=25s from the beginning for test T1 and T2 and it cost about 40 h of CPU time using a personal computer (CPU 2.66 GHz and RAM 4.0 GMB PC). An initial particle distance of Δx = Δz 0.022 m was chosen for these two tests, resulting in a total number of 30,000 particles.

In this study, with an initial particle distance of Δx = Δz = 0.01 m and a total number of 62,892 particles, the hardware accelerated smoothed particle hydrodynamics code DualSPHysics was executed up to t = 100 s from the beginning and cost about 12 h of CPU time using a personal computer (CPU 2.70 GHz and RAM 384 GB PC, equipped with a NVIDIA Quadro K2000 GPU). Furthermore, with an initial particle distance of Δx = Δz 0.003 m and a total number of 389,000 particles, the hardware accelerated smoothed particle hydrodynamics code DualSPHysics was executed up to t = 100 s from the beginning and cost about 225 h of CPU time using a personal computer (CPU 2.70 GHz and RAM 384 GB PC, equipped with a NVIDIA Quadro K2000 GPU).

Instead, the discontinuous Galerkin FEM was executed up to t = 100 s from the beginning and cost about 9.75 h of CPU time using a personal computer (CPU 2.66 GHz and RAM 4.0 GMB PC). Square elements with sizes Δx = Δy = 0.04 m were chosen for these two tests, resulting in a total number of 5607 nodes.

5. Conclusions

The most widely used CFD approaches, namely the Eulerian and Lagrangian methods, are characterized by specific advantages and disadvantages. The most relevant advantages of Lagrangian approaches over Eulerian methods are the exact conservation of mass and momentum and the meshless properties. Moreover, the Lagrangian method gives detailed information of individual particles that can be crucial in many applications. In terms of computing cost, the Eulerian method was faster than the Lagrangian method. In fact, tracking a sufficient number of particles and analyzing them in the Lagrangian approach required more time than the Eulerian iterations.

New research challenges on Computational Fluid Dynamics (CFD) are constantly on-going; in particular, recently, research has focused on coupling the Eulerian–Lagrangian techniques to combine the advantages of the individual models in a single model, thus increasing the accuracy, efficiency and regime of validity.

The present study dealt with the comparison of the Lagrangian and Eulerian approaches in the simulation of regular breaking waves on a plane slope and solitary waves transformation, breaking and run-up. The performance of the numerical models was evaluated in accordance with experimental data.

Both Eulerian and Lagrangian approaches were capable of simulating the analysed phenomena. However, the Lagrangian method showed more accuracy to reproduce the wave shape, both in deep waters and shoaling zones right up to the breaking region. Eulerian models cannot reproduce wave overturning exactly (e.g., only one water surface elevation value is permitted). Additionally, when using the HFA method in the case of regular waves, a high vertical resolution multilayer model is required to capture the hydrostatic shock of a breaking wave [91].

For the tests T3–T4, the numerical results demonstrated that the Eulerian and Lagrangian methods have similar accuracy on predicting the solitary wave transformation, breaking and run-up processes.

In this study, the hardware accelerated DualSPHysics code on a GPU was used; the new DualSPHysics code is a user-friendly platform designed to encourage researchers to use the SPH technique to investigate novel and real CFD problems that are beyond the scope of classical models. In recent years, DualSPHysics code has been shown to be accurate for simulating coastal engineering problems, with a continuous improvement in efficiency thanks to the exploitation of hardware such as GPUs for scientific computing.

In this study, this important acceleration of the code using the new GPU technology has been observed by comparing its performance against the classic version of SPHysics code using a single core.