1 INCOMPRESSIBLE SMOOTHED PARTICLE HYDRODYNAMICS ( ISPH ) MODELLING OF BREAKWATER OVERTOPPING

This paper describes an investigation into using incompressible smoothed particle hydrodynamics (ISPH) to simulate the overtopping of a coastal structure such as a breakwater. The paper presents the ISPH formulation that employs the multiple boundary tangent method and the latest developments such as particle shifting that produce noise-free pressure fields. The numerical model is compared with experimental overtopping data for a solitary wave and a crest-focussed wave group approaching a trapezoidal breakwater. The ISPH model is shown to produce close agreement for the free-surface evolution for both types of wave and generates overtopping volumes in satisfactory agreement with experimental data. Closer agreement with experimental data is obtained for ISPH compared to more popular weakly compressible SPH for the same resolution or particle size. Future work identifies conducting a convergence study and using more sophisticated boundary treatments.


INTRODUCTION
At present, engineers possess mainly empirical means to predict the overtopping of coastal protection structures.The modelling and simulation of such highly nonlinear and potentially violent free-surface motion is an extremely challenging task with numerous difficulties.In this work, we examine the behaviour of violent wave overtopping processes using truly incompressible smoothed particle hydrodynamics (SPH).
SPH describes a fluid by replacing its continuum properties with locally smoothed quantities at discrete Lagrangian locations.Thus, the domain can be multiply-connected with no special treatment of the free surface, making it ideal for examining complicated flow situations.
SPH has become increasingly popular in recent years as a novel technique to model the violent hydrodynamics in wave breaking, etc. (Rogers and Darlymple 2006).However, the vast majority of SPH schemes are based on the weakly compressible SPH (WCSPH) formulation where the density is allowed to vary slightly so that compressibility effects are kept within 1%.This is generally acceptable for engineering computations (Rogers et al. 2010), but can lead to severe problems when predicting impact pressures while the propagation of waves can exhibit significant decay and dissipation.
In Manchester, we have been developing a truly incompressible SPH (ISPH) scheme where a pressure Poisson equation is solved to predict pressure fields.A novel shifting algorithm has effectively eliminated unphysical diffusion at the free-surface leading to accurate noise-free pressure fields for a range of test cases including periodic waves, impulsively moving plates, and slam problems (Lind et al. 2012, Skillen et al. 2013).For the simulation of overtopping, previous work using the weakly compressible approach (Stansby et al. 2008) produced close agreement with experimental data for solitary waves.The application of ISPH to overtopping is the next important step in the development of this method.
In this paper, first we introduce the SPH discretisation.This includes a short description of the methodology for enforcing strict incompressibility within SPH.The next section presents the numerical results for two cases of a solitary wave and a crest-focussed wave group overtopping a breakwater.Comparisons with experimental data for overtopping the breakwater are presented.

SPH Methodology
Smoothed Particle Hydrodynamics (SPH) is based on the approximate representation of continuous interpolations or integrals by a discrete particle representation.The value of a flow quantity, ϕ, at a position vector r is approximated by a local summation where V j is the volume of the jth particle located at r j with scalar quantity ϕ j , and   j r r   is the weighting function called the smoothing kernel.The kernel is specified prior to simulation by an analytical expression making its computation straightforward.In all results presented, a fifth-order spline kernel has been used (Lind et al. 2012).The kernel has a characteristic smoothing length, h, which defines the region of influence, and is preset to be 1.3x for the duration of all simulations where x is the initial particle spacing.
In SPH function derivatives can be expressed as another summation by simply using a derivative of the smoothing kernel, i.e.
. This avoids complicated expressions for calculating derivatives since an analytical expression is known for the specified kernel.Expression (2) also ensures that there is zero divergence for a uniformly constant field.

Governing Equations
Herein, we solve the incompressible Navier-Stokes equations in Lagrangian form such that the divergence of the velocity field is zero where u is the velocity.The conservation of momentum is therefore written as: where P is pressure, ρ is density, ν is viscosity, t is time and g is gravity.Incompressibility is enforced using a pressure Poisson equation.Using a fractional step method, the pressure gradient term is discretized according to: where W represents the corrected kernel gradient as used by Lind et al. (2012).The viscous and gravity terms are discretised as

Pressure Projection Method
The solution of the pressure gradient term follows the approach first suggested by Cummins and Rudman (1999) for SPH.To solve Eq. ( 3), the solution first advects the particles to an intermediate location on viscosity and body forces is then computed This produces a matrix equation which is solved using a Bi-CGSTAB method.(van der Vorst 1992).This enables the projection of * i u onto the divergence-free space giving the velocity at time n+1: The particles' positions are then updated using a centred time integration scheme

Particle Shifting
The pressure projection approach of Cummins and Rudman (1999) has been shown to be unstable for irregular particle distributions (Xu et al. 2009).Lind et al. ( 2012) and Skillen et al. (2013) have shown that noise-free pressure fields for free-surface flow can be computed by shifting particles a small distance δr i according to the gradient of the local concentration of the particles C: where D is a local diffusion coefficient.The diffusion coefficient is a numerical parameter evaluated according to the maximum timestep: where ∆t max is provided by a CFL condition: where U max is the maximum particle velocity.In order to prevent the particles being shifted either beyond or too close to each other, the shifting distance, δr i , is restricted to a maximum value of 0.2h.
A measure of the particle concentration can be computed from the sum of the kernel function: So the gradient of the concentration is given by where f ab is a tensile correction term (Lind et al. 2012).The hydrodynamic variables are finally corrected by a Taylor series approximation.

Boundary Treatment
The original SPH method does not include the presence of boundaries naturally since the particle approximation in Eq. ( 1) is only a volume integral and does not account for the missing kernel support in the vicinity of boundaries.There have been many different treatments for boundaries proposed in SPH (see Ferrand et al. 2012 for a comparison of some popular methods).In this work, we use the multiple boundary tangent (MBT) method of Yildiz et al. (2009).In the MBT approach, a fluid particle will interact with a wall particle, w i .A tangent at w i is computed and then a localized mirror image of the fluid particles is generated using that tangent line as shown in Fig. 1 for different cases.For overlapping mirror images, the mass of each particle is adjusted so as to maintain a physical correct volume.

Solitary Wave Overtopping
In Fig. 3, we can see a comparison for the free surface at different time instants with (a) a Boussinesq-type model, a weakly compressible SPH model (Stansby et al. 2008) and (b) ISPH for solitary wave.The solitary wave has a height of 0.1 m and is generated in the SPH simulations using the Goring (1978) wavemaker.The agreement between the computed ISPH results and other numerical predictions results is close, but there are differences which will be the focus of future investigation.Fig. 8 shows a comparison of the measured free surface at each wave gauge and the free surface predicted by ISPH using the same resolution of 0.01 m (note that the vertical scales shown for the experimental and numerical results are not the same).The agreement is reasonable except when the focused wave group reaches WG5 which is located in the breaking region where the measurements will have included the aeration, where there is some oscillation evident in the ISPH results.Finally, Fig. 9 shows the time history of the overtopping volume predicted by ISPH where there are two clear overtopping events in agreement with experimental observations and the final overtopping volume predicted by ISPH is 9.8 litres/m compared to 12 litres/m measured in the experiments.The equivalent WCSPH scheme produced no overtopping whatsoever for this resolution.

CONCLUSIONS
This paper has presented the application of an incompressible SPH numerical scheme to the overtopping of breakwaters using the latest formulation improvements for ISPH including particle shifting.The paper has presented investigations into the modelling of different wave types including solitary waves and crest focussed wave groups looking at the hydrodynamics of the overtopping process in the vicinity of the obstacle with an improved ISPH formulation.The results for only one resolution (or particle size) have been presented but the results compare favourably with more established WCSPH models, and in the case of focussed wave group overtopping, the ISPH predicted overtopping in satisfactory agreement with the experimental data where the WCSPH scheme produced no overtopping whatsoever.In future work, more recent developments in solid wall boundary treatments will be introduced which will improve the robustness of the scheme.

Figure 1
Figure 1 Multiple boundary tangent (MBT) method (Yildiz et al. 2009) for generating localised mirror images for different cases

Figure 2
Figure 2 Geometry for UKCRF experiments and ISPH simulations

Figure 3
Figure 3 Surface profiles for a solitary wave at t = 8s, t = 9s, t = 10s

Figure 4 Figure 6
Figure 4 Solitary wave free surface time histories at toe of beach

Figure 7
Figure 7 wave gauge location for focussed wave group overtopping

Figure 8
Figure 8 Focussed wave group free-surface profiles at wave gauges WG1-5 (note: vertical axes are not on same scale, bottom axes represent time in seconds)

Figure 9
Figure 9 Overtopping volume history for focussed wave groups