NUMERICAL MODELING OF BREAKING SOLITARY WAVE RUN UP IN SURF ZONE USING INCOMPRESSIBLE SMOOTHED PARTICLE HYDRODINAMICS (ISPH)

This paper presents a numerical model for simulating wave run-up on rough sloping surfaces. Incompressible smoothed particle hydrodynamics (ISPH) has been utilized, which is capable of efficient tracking of free surface deformation in a Lagrangian coordinate system. Many of the existing models have focused on inviscid wave run-up on a smooth surface, but few numerical models and especially experimental studies have investigated the effect of beach roughness on the run up. In the present study two methods have been deployed to study the effect of roughness on wave run up. In the first method, the mass unit force, which is a coefficient of the fluid viscosity, and is dependent on the roughness of the solid boundary, has been used. In the second method, mass unit force obtained from the wall functions was utilized to enforce the friction on the particles near the boundaries. The comparisons of the numerical simulations with the analytical solutions and experimental data confirmed the capability of the model in simulating wave propagation and wave breaking. It was also concluded that the effect of roughness on wave run up depends on both the roughness itself and the beach slope. The results also indicated small roughness effect on waves running up over steep slopes.


INTRODUCTION
The problem of determining the run-up of solitary waves over beaches usually arises in the study of the coastal effects of tsunamis.Tsunamis are long water waves of small steepness.As is well known, tsunamis are caused by sudden displacements of sea water due to shallow submarine earthquakes whose epicenter is less than 40km below the seabed, landslides, submarine slides and volcanic activity.
Significant amounts of literature have been published that analyze long wave shoaling and run-up on slopes using the nonlinear shallow water equations (NSWE).The NSWE describe a thin layer of fluid of constant density in hydrostatic balance, bounded from below by bottom topography and from above by a free surface.In wave run-up modeling, one of difficulties in breaking wave simulation is the free surface tracking.The VOF method is one of the most flexible and robust approaches in this regard.
Recently, numerical methods which do not use any grid structures have been developed, such as Lagrangian or particle methods.The smoothed particle hydrodynamics (SPH) and moving particle semi-implicit (MPS) methods are two of the more robust approaches.In SPH, the state of a system is represented by a set of particles, which possess material properties and interact with each other within a range controlled by a kernel function.In the early simulations of fluid flows by the weakly compressible SPH, incompressibility was realized through an equation of state so that the fluid was assumed slightly compressible.In this case, a large sound speed has to be introduced, which could easily cause problems of sound wave reflections at the solid boundaries.In ISPH and MPS methods, the pressure is not a hydrodynamic variable to be obtained from the equation of state, but is obtained by way of solving a pressure Poisson equation derived from a semi-implicit algorithm.
Tsunamis and periodic wave run-up are complex phenomena.In spite of excellent studies carried out in this area, there are still many questions remaining to be solved.One area that needs further research is the effect of terrain roughness on tsunami run-up.Many of the existing models have focused on inviscid wave run-up on smooth surfaces, but there have been a few numerical and experimental studies which are carried out to investigate how beach roughness affects the run up.This paper presents an incompressible two-dimensional SPH model to simulate the wave run-up on smooth and rough slopes.Two methods have been implemented in the present model to account for the friction in the ISPH model.In the first method the friction force is taken equal to a constant coefficient of viscosity term, and in the second method the shear velocity is calculated from a wall function.
Early efforts for estimating wave overtopping focused on wave run-up based on numerous flume and basin tests.For instance, in the Netherlands, the height used for coastal embankments was the extreme water level with an allowance for the so-called 2% run up level (Vandermeer 2002).Vandermeer and Stam (1992) studied irregular wave run up.Their laboratory studies yielded formulations to assess various run up levels as a function of the surf similarity or breaker parameter.

GOVERNING EQUATIONS
The governing equations employed in the present model are the conservation of mass and momentum equations, which are written in the Lagrangian form as follows: where V is the velocity vector, p is pressure, t is time,  is the density of water, g is gravitational acceleration, and µ is the dynamic viscosity of water.
The momentum and continuity equations in the form of SPH written for a particle i are presented as follows (Lo and Shao 2002): In the above equations i is the reference particle and j is its neighboring particle.i m , i  , i  , i P and i V are mass, density, viscosity, pressure and velocity of particle i, respectively.ij r is the distance between particles i and j, ij W is the kernel function, and η=0.1h,where h is the smoothing length.

NUMERICAL PROCEDURE
The fractional step scheme utilized for the ISPH method consists of two steps.In the first step the velocity field is computed using momentum equation in the absence of the pressure gradient term.Incompressibility is not satisfied in this step and the fluid density ( *  ) is calculated based on the temporary positions of the particles.An intermediate particle velocity and position are obtained as follows (Lo and Shao 2002): Here ( t V ) and ( t r ) are particle velocity and position at time t respectively, ( * V ) and ( * r ) are intermediate particle velocity and position respectively and ( * V  ) is the change of particle velocity during the first step and ∆t is the time increment.
In the second correction step, the pressure term is used to enforce incompressibility by continuity equation where a Poisson equation is solved, and the new velocity values and positions of the particles are then computed.
where ( * * V  ) is the change of the particle velocity during the correction step.The incompressibility is enforced through the following equation: where *  is the intermediate particle density.By combining equations 8 and 9: Laplacian will lead to the second derivative of the kernel function that is very sensitive to particle disorder, and can cause pressure instability.Lo and Shao (2002) introduced a formulation for a stable Laplacian as follows: where . By combining equations 10 and 11 the Poisson equation, which is a linear system of equations, is obtained as follows: To solve the Poisson equation, a combination of direct and iterative methods has been utilized.In the first time-step, the pressures are computed by the direct method.This pressure is then used for the next time-step as an initial value for the iterative procedure.The solution continues by using the past time-step pressure for each time-step as an initial value for the iteration process.Since a limited number of particles around the reference particle contribute to the Poisson equation, the sparse property of the matrix of pressure coefficients has been used to decrease the computational time.Poisson equation of pressure can also be described as follows (Ataie-Ashtiani and et al. 2008): The Poisson equation (Eq.14) is a modified Poisson equation of pressure in ISPH method which uses free divergence velocity constraint instead of invariant density.( * * V  ) is computed by using equation ( 15) as follows: Finally particle velocity and position at time (t+∆t) are computed by using equations 16 and 17 as follows:

WALL BOUNDARY AND FREE SURFACE
A major challenge in SPH method is the solid boundary simulation.In this study, a line of particles are used to represent the boundary particles, and two lines of ghost particles are considered outside the solid boundary.The pressure of the ghost particles is set to that of corresponding wall particles in the normal direction of the solid wall (Ataie- Ashtiani and et al. 2008).
At the water surface the particles are considered as free surface particles, where they satisfy the following condition (Lo and Shao 2002) in which β = constant (0.97).A condition of the pressure P = 0 is given to the water-surface particles.
Since there is no particle in the outer region of free surface, the particle density decreases on this boundary, which leads to spurious pressure gradients.To overcome this problem, it is assumed that s is a surface particle with zero pressure and i is an inner fluid particle with pressure i P .In order to calculate the pressure gradients between the two particles, a mirror particle (m) with pressure i P is placed in the direct reflection position of inner particle i through the surface particle s.The gradient of the pressure between the surface particles, mirror particle (m) and inner particle (i) is expressed as follows (Lo and Shao 2002):

NUMERICAL CONVERGENCE
Because the individual fluid particles are discrete points and cannot deform as the real fluid does, the number of particles employed in the computations must be sufficiently large to lead to numerical convergence.The time step is controlled in the computations to satisfy the following condition (Lo and Shao 2002): where max v is the maximum particle velocity in the computations and 0 l is the initial particle spacing.The factor 0.1 ensures that the particle can move only a fraction (in this case 0.1) of the particle spacing per time step.In addition, the constraint of the time step due to viscous diffusion must also be satisfied by equation ( 25) as follows (Lo and Shao 2002): where α is a coefficient depending on the choice of the kernel type and particle arrangement.α is usually in the order of 0.1, which is determined from SPH numerical trials.

APPLICATION OF SURFACE ROUGHNESS IN ISPH 6.1. The first method
No slip condition sets the relative speed of fluid equal to zero on the boundary.In other words, fluid undergoes friction from the underlying boundary.In this paper an approach similar to that of Muller et al. ( 2004) has been applied.By using the viscosity equation and considering the boundary particles (Muller et al. 2004;Krištof et al. 2009): where b f is the boundary friction force, k is a constant and b n is the number of the boundary particles in the support domain, which are represented with red particles in figure 1.
Boundary friction force has been considered as a diffusive term in equation ( 26), where the variable k applies the boundary roughness effect.For the larger number of the boundary particles in the support domain, more frictional force is applied to the fluid particles.The number of the particles inside the support domain reduces away from the boundary.The boundary friction applied to the fluid particle has a direct correlation with its distance from the boundary.Due to the instability caused by the second derivatives, equation 26 is discretized as follows: The force per unit of mass is expressed as follows:

The second method
Surface roughness typically leads to an increase in turbulence production near the wall.This can result in significant increase in the wall shear stress.For an accurate prediction of near wall flows, the proper modeling of surface roughness is essential for a good agreement with experimental data.Wall roughness increases the wall shear stress and breaks up the viscous sub layer in turbulent flows.The wall boundary conditions are implemented by the following equations (Versteeg and Malalasekera 2007): where y is the distance of near wall particle to the solid surface, y is the dimensionless distance from wall and  u is the dimensionless near wall velocity.A near wall particle is taken to be laminar if  y ≤ max (11.63,  h ).The wall shear stress is assumed to be entirely viscous in origin.If  y ≥ max (11.63,  h ), the flow is turbulent and the wall function approach is used. h is the dimensionless roughness height, U is the tangential velocity on the wall and  u is the friction velocity.

 
If near wall particle is taken to be in the linear zone, then: (32) If near wall particle is taken to be in turbulent area then: where B = 5.2.The shift ∆B is a function of the dimensionless roughness height.Figure 2 shows the downward shift in the logarithmic velocity profile.

Figure 2. Downward shift of the logarithmic velocity profile
For sand grain roughness, h is the equivalent roughness height, and the downward shift can be expressed as follows (Zhang 2009): Replacing equations 31 and 34 in equation 32 leads to (Zhang 2009): By solving Equation 35,  u is obtained and w  and m f can then be calculated as follows: where a is a constant coefficient for calibration, w  is the shear stress and m f is force per unit of mass.

HYDRODYNAMIC TESTS 7.1. Solitary wave propagation
Solitary wave propagation on a constant water depth is simulated in order to test the accuracy of the numerical scheme for the ISPH model.The analytical solution for the wave profile can be derived from the Boussinesq equation as follows (Lo and Shao 2002): )) ( 4 where η is the water surface elevation, H the wave amplitude, h the water depth and the solitary wave celerity.The horizontal velocity under the wave profile is given A solitary wave with a wave height of H = 0.04cm is considered.The water depth h is 0.18m and the initial particle spacing 0 l is 1cm. 11245particles are involved in the computations.Figure 3 shows the comparisons of simulation of a solitary wave with analytical solution.The results show that the wave maintains its form during the wave propagation, and the agreements with the analytical solution are convincing.

Regular wave propagation
The incompressible SPH method has also been simulated for regular wave propagation.The analytical solution for the wave profile can be derived from the equation ( 40) as follows: where k is the wave number, L is the wave length, w is the wave frequency and T is the wave period.The kinematic condition must be satisfied on the wave maker.If S is the horizontal displacement of a piston-type wave maker, then the wave maker position may be described as follows (Dean and Dalrymple 1991) The range of wave maker motion (S) can be computed from the following equation (Dean and Dalrymple 1991): To absorb the wave reflection from the wall at the end of the flume, an exponential damping zone is placed over a distance of at least a wave length.In the damping zone, the velocity of fluid particles will be damped according to the formulation as follows (Xu 2010): where V is the fluid velocity, and: δx is the damping zone length which is equal to 1.5 in the present model, and 0 x is the starting point of damping zone.A regular wave with the wave height of H = 0.04cm is considered.The water depth h is 0.24m, the initial particle spacing 0 l is 1cm, the wave period is 0.9s and the wave length is 1.108m. 13917particles are involved in the computations.Figure 4 shows the regular wave propagation, and the performance of the wave absorber and damping zone.The numerical simulations are in good agreement with the analytical solution.

WAVE BREAKING
The laboratory breaking solitary wave experiments provided appealing tests for demonstrating the capability of the ISPH model to simulate a breaking solitary wave run up on a mild, 1:20 (Synolakis 1987), and a 1:15 slope (Grilli et al. 1997).

Comparison with the laboratory measurements on a slope of 1:20
In the experiments the still water depth was d = 21.168cm, the slope of the beach 1:20 and the incident wave height H=5.88cm (Synolakis 1987).The initial particle spacing was considered to be 0 l =1.176cm, and 5457 particles were used in the simulations.In figure 5, the numerical model simulations are compared with the experimental measurements of Synolakis (1987).In figure 5-a wave height increases, and figure 5-b shows the final stage of the solitary wave shoaling.The wave form becomes highly asymmetric due to decrease of water depth, and figure 5-c shows the wave run up.In figure 5-b, the amplitude increases about 2.22cm higher than the incident wave height.Further steeping of the wave front causes the wave to start breaking.

Comparison with the laboratory measurements on a slope of 1:15
In this experiment the still water depth was d = 30.36cm, the slope of the beach was 1:15 and the incident wave height was H = 13.8cm(Grilli et al. 1997).The initial particle spacing was considered to be 0 l =1.38cm and 9948 particles were used in the simulations.Figures 6-a to 6-d show the comparisons of the numerical simulations of a solitary wave breaking on a 1:15 slope with the experimental measurements.In figure 6-b, process of wave shoaling has been finished.Wave front takes a vertical form and breaking process starts.In figures 6-c and 6-d the wave breaking is completed, and in figure 6-d, the differences between the numerical and experimental values are apparent.

NUMERICAL MODELING OF WAVE RUN-UP ON A ROUGH SURFACE
In the present study, simulation of wave run up is carried out by two different methods for investigating the influence of boundary roughness.The run-up simulations were conducted for both smooth and rough (n = 0.02, n = 0.024, n = 0.033; n is the Manning coefficient) beaches with two different slopes (β), namely, β = 5 and β = 10 degrees.Numerically simulated values were compared with the laboratory results of Teng and Feng (2000).The water depth ranged from 0.127m to 0.203m.
In the first method, f k in equation ( 28) is the boundary friction constant and is related to the surface roughness.The numerical results for wave run-up on both smooth and rough surfaces with β = 5 degrees are presented in figure 7. The values of f k are calibrated for a specified value of a/h, and then the run up is computed for the other values of a/h using the computed constant roughness.In figure 7, a is the wave amplitude, h is the water depth and R is wave run-up.The values of f k are calibrated for different Manning coefficients as presented in Table 1.The Manning coefficient for a smooth surface is 0.009.In the second method, the a constant in equation 37 is calibrated for n = 0.02.a is equal to 0.00316 and the numerically simulated results for wave run-up for both smooth and rough slopes for β = 10 degrees are presented in figure 8. Table 2 shows the values of RMSE for the first and second methods, where RSME is defined as follows: where exp F and num F are values for the experimental and numerical run up.

CONCLUSIONS
ISPH method has been utilized to simulate wave propagation and run up.The model is capable of simulating the propagation of solitary and regular waves.It is also capable of simulating the whole process of wave propagation, shoaling, breaking and run-up on beach slopes.
Two methods were used to evaluate the run up.In the first method the approach was based on computing the frictional force.In second method, the shear velocity was calculated from the wall function.
It is concluded that run up to water depth ratio increases with the increase of wave amplitude to water depth ratio.The results showed that the effect of terrain roughness on wave run-up depends on both the roughness itself and the beach slope.These results showed that for wave run up over rough slopes, roughness can cause a reduction in run up height compared with wave run up on smooth slopes and run up height increases with increasing slope.The roughness effect on wave for a slope of 10 degrees is smaller than that for a slope of 5 degrees.Generally, when the beach slope becomes milder, the effect of roughness increases.Both methods showed good agreements with the experimental data, with closer agreements for the second method for smooth surfaces, and surfaces with n = 0.02 and 0.024, and also for the first method with n = 0.033 for a slope of 5°.For the slope of 10°, the first method showed better agreements.

Figure 1 .
Figure 1.Boundary particles in the support domain

Figure 5 .
Figure 5. Comparisons between the simulated results and experimental data for water surface profile of a solitary wave breaking on a 1:20 slope at a) t = 1.7s, b) t = 2.4s, c) t = 3.7s.

Figure 6 .
Figure 6.Comparisons of the computed and measured surface profile of solitary wave breaking on a 1:15 slope at a) t = 1.95s, b) t = 2.05s, c) t = 2.1s, d) t = 2.15s.

Figure 7 .Figure 8 .
Figure 7.Comparison of numerical results and experimental data of run-up using the first and second methods on a surface with β = 5 degrees: a) smooth surface, b) surface with n = 0.02, c) surface with n = 0.024, d) surface with n = 0.033.

Figure 9
Figure9shows wave run up for both smooth and rough surfaces with a Manning coefficient equal to n = 0.02 for both methods.

Figure 9 .
Comparison for run up over smooth and rough surfaces on slopes of: a) 5 degrees, b) 10 degrees, c) 15 degrees.