Modelling of Violent Water Wave Propagation and Impact by Incompressible SPH with First-Order Consistent Kernel Interpolation Scheme

: The Smoothed Particle Hydrodynamics (SPH) method has proven to have great potential in dealing with the wave–structure interactions since it can deal with the large amplitude and breaking waves and easily captures the free surface. The paper will adopt an incompressible SPH (ISPH) approach to simulate the wave propagation and impact, in which the ﬂuid pressure is solved using a pressure Poisson equation and thus more stable and accurate pressure ﬁelds can be obtained. The focus of the study is on comparing three different pressure gradient calculation models in SPH and proposing the most efﬁcient ﬁrst-order consistent kernel interpolation (C1_KI) numerical scheme for modelling violent wave impact. The improvement of the model is validated by the benchmark dam break ﬂows and laboratory wave propagation and impact experiments.


Introduction
The wave propagation and its impact on structure is a common natural phenomenon that is very important in the field of ocean and coastal engineering.During this process, the waves involve with large deformation of free surface, wave breaking, wave run-up and their strong interactions with the structures.Thus, modelling of wave system is of great interest but is also a difficult task in both physical experiments and the numerical simulations.
There are generally two classes of models for wave simulations.In the early years, the Laplace equation with fully nonlinear boundary conditions was used for simple wave problems [1][2][3], while the Navier-Stokes (N-S) equations with more realistic physical boundary conditions were solved in recent numerical models [4,5].As a result, various numerical schemes, such as the finite element, finite volume and finite difference methods have been used to solve the N-S equations to investigate the nonlinear water wave propagation and its interaction with the structures [6,7].To relieve the CPU expense and meanwhile provide reasonable solutions for the large domain problems, different forms of the Shallow Water Equations (SWEs) are also proposed in the modelling of wave system [8].
In the last two decades, the Smoothed Particle Hydrodynamics (SPH) method has emerged as a promising mesh-free Lagrangian modelling technique.Although the SPH model may use more CPU time than the grid counterpart in some cases, it has the advantage of tracking the free surface in an easy and accurate way.In this approach, the governing equations are discretized and solved by the individual particles within the computational domain.SPH was originally developed for the study of astrophysics [9,10] and then employed to study the wave propagating and overtopping over coastal structures [11].Based on the algorithm of the Moving Particle Semi-implicit (MPS) approach [12] and SPH projection method [13], a strict incompressible algorithm of the SPH model [14] was developed and then it was further improved to simulate wave breaking and post-breaking [15].The major difference between the standard weakly compressible SPH (WCSPH) [9,10] and the incompressible SPH [14] lies in that the former calculates the fluid pressure explicitly by using an equation of state, while the latter employs a strict incompressible formulation to solve the pressure implicitly by a pressure Poisson equation (PPE).Both the WCSPH and ISPH show capabilities and limitations.For WCSPH, it has been highlighted the better possibility to parallelize the numerical code for simulations in real conditions [16].The pressure field has also been improved by different form of diffusive terms in the continuity equation for wave-structure interaction problems [17,18].Furthermore, the acoustic components related to the use of the state equation can be eliminated by an appropriate filtering in the data post-processing to recover the incompressible solution [19].For the impact problems, it has been noticed that ISPH sometimes shows singularities since the pressure is inversely proportional to the adopted time step [20].On the other hand, the numerical time stepping length of ISPH can be larger than that of WCSPH, so the general computational efficiency could be higher [21,22].
In this paper, an incompressible SPH model will be employed to simulate solitary wave propagation and impact on a slope with different conditions.Although the ISPH model could predict stable and noise-free pressure field in many cases, the numerical accuracy and efficiency could be compromised in the free surface water waves.Among a variety of influence factors, the calculation of pressure gradient plays a very important role in the process.Thus, the current study presents an improved first-order derivative of pressure model for the violent wave simulations, from the first-order consistent kernel interpolation scheme (C1_KI).The proposed C1_KI ISPH model is first verified by the benchmark dam break flow and solitary wave propagation over a constant depth.Then, the wave propagation and impact on an inclined wall are investigated under different slope angles based on the self-designed laboratory experiment.

Governing Equations and Solution Algorithms
The Navier-Stokes (N-S) equations are used to describe the fluid motion.In the incompressible SPH method, the fluid density is considered to be a constant, and therefore the mass and momentum conservation equations are written in the Lagrangian form as follows: where ρ is the fluid density; u is the particle velocity; t is the time; P is the particle pressure; g is the gravitational acceleration; and ν 0 is the kinematic viscosity.A two-step projection method is used to solve the velocity and pressure field from Equations ( 1) and (2).The first step is the prediction of velocity in the time domain without considering the pressure term.The intermediate particle velocity u * and position r * are obtained by where u t and r t are the velocity and position at time t; ∆t is the time step; ∆u * is the velocity increment; and u * and r * are the intermediate velocity and position of particle.
Water 2017, 9, 400 3 of 17 The second step is the correction step in which the pressure term is added, and ∆u * * is the correction of particle velocity The following u t+∆t and r t+∆t represent the velocity and position of particle at new time step Combining Equations ( 1)-( 6), the following pressure Poisson equation (PPE) is obtained Similarly, Shao and Lo [14] proposed a projection-based incompressible approach by imposing the density invariance on each particle, leading to the following PPE: where ρ * is the particle density at intermediate time step.Due to ρ * /ρ 0 being very close to unity, the difference between the left-and right-hand sides of the denominator in Equation ( 10) can be ignored, and the combined PPE incorporating both the divergence-free and density-invariance terms is obtained as: where α is a blending coefficient and a value of 0.01 is adopted in this paper from the computational experience.In this paper, only a 2D ISPH formulation is used.

Calculation of Spatial Derivatives
A common approach to calculate the gradient of pressure and the divergence of velocity is through the following equations: where u ij = u i − u j is defined; ∇ i W is the gradient of SPH kernel function and a cubic spline kernel [10] is used; m is the particle mass; h is the kernel smoothing length; N is the total neighbouring particle number; and i and j indicate the reference and neighbouring particles, respectively.The viscosity term in Equation (2) adopts the following form where η 2 = 0.01h 2 is a small parameter to avoid singularity; and r ij = r i − r j is defined.The Laplacian term in Equation ( 11) is discretised by combining the SPH gradient and divergence rules to obtain where P ij = P i − P j is defined.

Free Surface and Solid Boundary Conditions
The dynamic free surface conditions require a prescribed pressure to be imposed on the surface particles, such as through P = 0.In this paper, we use three auxiliary functions combined with the ratio of particle number density to accurately identify the free surface particles.This has shown to be more robust than using either the density or the divergence rules in conventional ISPH practice.The particles on the solid boundary should satisfy the pressure boundary condition, which is represented by the following momentum balance as where n is the unit vector normal to the solid boundary.More detailed procedures to implement the free surface and solid boundary conditions can be referred to Zheng [21].

Improved First-Order Derivative Scheme
In this section three pressure calculation models are compared and an optimum is chosen to use in the practical water wave simulations.

Simplified Finite Difference Interpolation (SFDI) Scheme
In the traditional SPH calculation of pressure gradient in Equation (12), the computational results are heavily affected by the particle distributions and the shape of solid boundary.To improve this, the first-order derivative of pressure on the solid boundary was formulated by the Simplified Finite Difference Interpolation (SFDI) scheme originally proposed by Sriram and Ma [23] in their MLPG_R approach.SFDI is a second-order accurate scheme based on the Taylor series expansion, which can also be used for the inner fluid particles.In the 2D case, the pressure gradient model formulated by SFDI can be found in [24].The relevant key formulas are summarized as follows: where m = 1 and k = 2 or m = 2 and k = 1; N is the number of neighbouring particles affecting particle i; x 1 = x and x 2 = y; r j,x m is the component of position vector in x (or y) direction.

Moving Least Square (MLS) Method
MLS is a widely used interpolation scheme for meshless kernel approximation, which can get very high accuracy for the function estimation.Atluri and Shen [25] and Zheng et al. [26] have done many comparisons of improved meshless interpolation schemes with the MLS method.More details of MLS method for the function estimation and first order derivative calculation can refer to Zheng et al. [26].The accuracy of MLS scheme is very high and the calculation process is more complex than the other traditional meshless interpolation methods.Although the dimension of matrix is not large, the computational cost is demanding as it includes successive multiplication of several matrices and matrix inversions.Although MLS has been widely used for the meshless interpolation comparison [24,26,27], it is less documented in the wave propagation simulations.The key formulations of MLS are summarized as follows.
The unknown function is represented by f (x) as where N is the number of nodes that affect the function at x; and Φ j (x) is the interpolation or shape function given by Assume the basis function to be linear and define the matrix B(x) and A(x) as where W(x) and Ψ are respectively expressed by and w(x − x j ) is a weight function which can adopt different forms.
The gradient of the unknown function Equation ( 20) is estimated by The partial derivative of shape function with respect to x can be directly differentiated as where A −1 ,x is the partial derivative of A −1 , i.e., A −1 ,x = −A −1 A ,x A −1 ; and B j is the j-th column of B, for which the partial derivative is calculated as 3.3.First-Order Consistent Kernel Interpolation (C1_KI) Scheme The standard SPH practice widely uses the following symmetric summation form to calculate the first-order derivative as However, the above equation was based on the assumption of (∆ j is the particle volume), which cannot be exactly satisfied when the particles are disorderly distributed or near the solid boundary.To improve this, the present paper applies the Taylor series expansion of f (x) at point x j , then multiplies this by a kernel function W and integrates over its support domain.
Eventually it gives the values of f (x i ), f x (x i ) and f y (x i ) in the following matrix form This method is easy to be expanded to 3D and higher-order problems.As known from Equation (32), when the distributed particles are symmetric about x i , the error terms become so all of f (x i ), f x (x i ) and f y (x i ) can achieve the second-order accuracy of O(h 2 ).On the other hand, when the distributed particles are asymmetric about x i , the error terms become can reach the second-order accuracy of O(h 2 ), while f x (x i ) and f y (x i ) merely obtain the first-order accuracy of O(h 1 ).Higher-order accuracy results can always be made available by keeping more terms in the Taylor series expansions.More details on the relevant analysis can be found in [26].The pressure gradient scheme proposed in this section is named as C1_KI, because it can keep the first-order consistency of the gradient estimation.The scheme has the following attractive features.First, it avoids the kernel gradient calculation.The matrix is strictly symmetric, in which only the upper triangle elements need to be calculated and the other elements can be obtained from the symmetric relationships.Besides, it is the first time to combine the C1_KI concept with the truly incompressible SPH scheme, which is expected to inherit the merits of both and therefore provide more accurate violent wave simulations.

Convergence and Accuracy of Different Pressure Gradient Models
In this section, several investigations are made into the convergence and accuracy behaviours of the three pressure gradient models, under the uniform and random particle distributions.The error of the numerical results is quantified as where E_a x is the defined mean error; f x (x j ) is the numerical value of the gradient component; f x (x j ) is the analytical solution; and M is the number of sampling points in the inner fluid domain or near the boundary area.Here, the near boundary area is defined as the inner region at a distance of 2DX (twice of the particle spacing) from the computational boundary and the remaining area inside is the inner fluid area.The computational domain of the test is chosen as a square with side length being of 1 m as shown in Figure 1a-d.The test function uses a polynomial function of second order expressed by f (x, y) = exp (2x+3y) .The calculation nodes are irregularly distributed in the domain by using the quasi-random number generator.The sample node distributions are illustrated in Figure 1, for the particle number of 400, 1600, 6400 and 25,600, respectively, corresponding to the particle size DX of 0.05 m, 0.025 m, 0.0125 m and 0.00625 m.
Water 2017, 9, 400 7 of 17 quasi-random number generator.The sample node distributions are illustrated in Figure 1, for the particle number of 400, 1600, 6400 and 25,600, respectively, corresponding to the particle size DX of 0.05 m, 0.025 m, 0.0125 m and 0.00625 m.When the particles are uniformly distributed, the convergence rate of the inner domain and near boundary area is shown in Figure 2a,b, respectively.It is shown from Figure 2a that _ x E a can reach the second-order accuracy for all three pressure gradient methods in the inner fluid domain.Among these SFDI can get the highest accuracy, while the errors of C1_KI and MLS are almost the same.As shown in Figure 2b near the boundary area, MLS can achieve the best results with minimum error while C1_KI performs similarly to SFDI to achieve first-order accuracy.On the other hand, Figure 3a suggests that when the particles are randomly distributed, C1_KI can achieve the most promising results in the inner fluid domain.In the area near the boundary as shown in Figure 3b, both C1_KI and MLS obtain similar results which are better than those computed by SFDI.When the particles are uniformly distributed, the convergence rate of the inner domain and near boundary area is shown in Figure 2a,b, respectively.It is shown from Figure 2a that E_a x can reach the second-order accuracy for all three pressure gradient methods in the inner fluid domain.Among these SFDI can get the highest accuracy, while the errors of C1_KI and MLS are almost the same.As shown in Figure 2b near the boundary area, MLS can achieve the best results with minimum error while C1_KI performs similarly to SFDI to achieve first-order accuracy.On the other hand, Figure 3a suggests that when the particles are randomly distributed, C1_KI can achieve the most promising results in the inner fluid domain.In the area near the boundary as shown in Figure 3b, both C1_KI and MLS obtain similar results which are better than those computed by SFDI.
Water 2017, 9, 400 8 of 17 Convergence and error analysis of first-order derivative calculation under uniform particle distributions in the: (a) inner fluid area; and (b) near boundary area.Generally speaking, Figure 3 implies that all of the three pressure gradient models can only achieve the first-order accuracy when the particles are irregularly distributed.In practical water wave modelling, the situation of particle randomness should increase from the wave propagation to wave-structure interaction cases because there are frequent exchanges of location within neighbour particles and surface-to-inner particles.Computationally, MLS is the most expensive one since it includes quite a few inverse calculations and matrix multiplications [24,26].Overall, C1_KI scheme performs quite satisfactory among the three in view of computational accuracy under all tested conditions, thus it will be used to study the wave propagation and impact in the next model applications.

Dam Break Wave Impact on a Vertical Wall
The dam break flow has always been a benchmark violent free surface flow to validate the numerical models.In this example, a rectangular column of water is initially confined.The width of water column is a and the height is H .At the beginning of the computation, the dam is instantaneously collapsed.A schematic setup of the dam break flow domain is given in Figure 4, where L is the distance between the two vertical walls.There are two pressure sensors 1 p and 2 p located on the left and right walls, respectively, with a distance of 1 h and 2 h from the horizontal Water 2017, 9, 400 8 of 17 Convergence and error analysis of first-order derivative calculation under uniform particle distributions in the: (a) inner fluid area; and (b) near boundary area.Generally speaking, Figure 3 implies that all of the three pressure gradient models can only achieve the first-order accuracy when the particles are irregularly distributed.In practical water wave modelling, the situation of particle randomness should increase from the wave propagation to wave-structure interaction cases because there are frequent exchanges of location within neighbour particles and surface-to-inner particles.Computationally, MLS is the most expensive one since it includes quite a few inverse calculations and matrix multiplications [24,26].Overall, C1_KI scheme performs quite satisfactory among the three in view of computational accuracy under all tested conditions, thus it will be used to study the wave propagation and impact in the next model applications.

Dam Break Wave Impact on a Vertical Wall
The dam break flow has always been a benchmark violent free surface flow to validate the numerical models.In this example, a rectangular column of water is initially confined.The width of water column is a and the height is H .At the beginning of the computation, the dam is instantaneously collapsed.A schematic setup of the dam break flow domain is given in Figure 4, where L is the distance between the two vertical walls.There are two pressure sensors 1 p and 2 p located on the left and right walls, respectively, with a distance of 1 h and 2 h from the horizontal Generally speaking, Figure 3 implies that all of the three pressure gradient models can only achieve the first-order accuracy when the particles are irregularly distributed.In practical water wave modelling, the situation of particle randomness should increase from the wave propagation to wave-structure interaction cases because there are frequent exchanges of location within neighbour particles and surface-to-inner particles.Computationally, MLS is the most expensive one since it includes quite a few inverse calculations and matrix multiplications [24,26].Overall, C1_KI scheme performs quite satisfactory among the three in view of computational accuracy under all tested conditions, thus it will be used to study the wave propagation and impact in the next model applications.

Dam Break Wave Impact on a Vertical Wall
The dam break flow has always been a benchmark violent free surface flow to validate the numerical models.In this example, a rectangular column of water is initially confined.The width of Water 2017, 9, 400 9 of 17 water column is a and the height is H.At the beginning of the computation, the dam is instantaneously collapsed.A schematic setup of the dam break flow domain is given in Figure 4, where L is the distance between the two vertical walls.There are two pressure sensors p 1 and p 2 located on the left and right walls, respectively, with a distance of h 1 and h 2 from the horizontal bed.In the following analysis, all the variables are non-dimensionalised by using the dam width a and gravitational acceleration g, such as t = t g/2a or H = H/(2a).
Water 2017, 9, 400 9 of 17 bed.In the following analysis, all the variables are non-dimensionalised by using the dam width a and gravitational acceleration g , such as In the present simulations it is assumed a = 0.5 m, / 2.0 H a = and / 4 L a = .The total particle number is 60 × 120 with a particle size of 0.00833 m.The non-dimensionalised time step is given by / (2 ) 0.003 t t g a Δ = Δ =  . To validate the proposed C1_KI SPH computations, Figure 5a,b gives the comparisons of dam break flow leading edge and water column height, respectively, with the experimental results from [28].A good agreement is found in both cases.are 1600, 6400, 19600, respectively, for the study of the model convergence.Correspondingly, the particle sizes of three cases are 0.025 m, 0.0125 m and 0.00714 m.The computational results of pressure time history on sensor point 1 p are shown in Figure 6a-c for the different dam height to width ratios / H a .The vertical height of sensor point 1 p is from the bed. Figure 6 shows that the peak pressure value changes with the height of water column, which is 0.4, 0.8 and 1.6, respectively, for / H a = 1, 2 and 4. Thus the relationship between maximum pressure and water column height seems to follow a linear correlation.Besides, C1_KI SPH computations demonstrate good convergence behaviours and different particle numbers In the present simulations it is assumed a = 0.5 m, H/a = 2.0 and L/a = 4.The total particle number is 60 × 120 with a particle size of 0.00833 m.The non-dimensionalised time step is given by ∆ t = ∆t g/(2a) = 0.003.To validate the proposed C1_KI SPH computations, Figure 5a,b gives the comparisons of dam break flow leading edge and water column height, respectively, with the experimental results from [28].A good agreement is found in both cases.from the bed. Figure 6 shows that the peak pressure value changes with the height of water column, which is 0.4, 0.8 and 1.6, respectively, for / H a = 1, 2 and 4. Thus the relationship between maximum pressure and water column height seems to follow a linear correlation.Besides, C1_KI SPH computations demonstrate good convergence behaviours and different particle numbers To carry out more robust model validations, three tests are run by changing the height of dam with H = a, H = 2a and H = 4a, and the simulations continue to the stage when the dam break wave front hits on the right vertical wall.For analysis purpose, the non-dimensional pressure is expressed as p = p/2ρga.The total particle numbers N used in case of H = a are 400, 1600, 4900, in H = 2a are 800, 3200, 9800, and in H = 4a are 1600, 6400, 19600, respectively, for the study of the model convergence.Correspondingly, the particle sizes of three cases are 0.025 m, 0.0125 m and 0.00714 m.
The computational results of pressure time history on sensor point p 1 are shown in Figure 6a-c for the different dam height to width ratios H/a.The vertical height of sensor point p 1 is h 1 = 0.1a from the bed. Figure 6 shows that the peak pressure value changes with the height of water column, which is 0.4, 0.8 and 1.6, respectively, for H/a = 1, 2 and 4. Thus the relationship between maximum pressure and water column height seems to follow a linear correlation.Besides, C1_KI SPH computations demonstrate good convergence behaviours and different particle numbers lead to almost identical results.Although the numerical results demonstrate some kinds of oscillation in the time history, they tend to become smoother with the refinement of particle size.
Water 2017, 9, 400 10 of 17 lead to almost identical results.Although the numerical results demonstrate some kinds of oscillation in the time history, they tend to become smoother with the refinement of particle size.; and (c) To quantify the accuracy of pressure predictions made by C1_KI SPH, the comparisons with experimental data of Zhou et al. [29] are also made.In this case, the water column dimension is a = 2 m and 0.5 , and the distance between the two vertical walls is 5.367 L H = .On the right wall, there is a pressure sensor point 2 P with a height of 2 0.133 h a = from the bed.From the C1_KI SPH simulations, the particle snapshots with pressure distributions are shown in Figure 7a-f at different time instants after the dam break.Besides, Figure 8 gives the time history of dam break flow impact pressures on the right wall computed by using different particle numbers.Compared with the experimental pressure data [29], not only a good agreement has been found, but also the pressure history of C1_KI SPH computations becomes more reasonable with a higher particle resolution.To quantify the accuracy of pressure predictions made by C1_KI SPH, the comparisons with experimental data of Zhou et al. [29] are also made.In this case, the water column dimension is a = 2 m and H = 0.5a, and the distance between the two vertical walls is L = 5.367H.On the right wall, there is a pressure sensor point P 2 with a height of h 2 = 0.133a from the bed.From the C1_KI SPH simulations, the particle snapshots with pressure distributions are shown in Figure 7a-f at different time instants after the dam break.Besides, Figure 8 gives the time history of dam break flow impact pressures on the right wall computed by using different particle numbers.Compared with the experimental pressure data [29], not only a good agreement has been found, but also the pressure history of C1_KI SPH computations becomes more reasonable with a higher particle resolution.
Water 2017, 9, 400 10 of 17 lead to almost identical results.Although the numerical results demonstrate some kinds of oscillation in the time history, they tend to become smoother with the refinement of particle size.To quantify the accuracy of pressure predictions made by C1_KI SPH, the comparisons with experimental data of Zhou et al. [29] are also made.In this case, the water column dimension is a = 2 m and 0.5 , and the distance between the two vertical walls is 5.367 L H = .On the right wall, there is a pressure sensor point 2 P with a height of 2 0.133 h a = from the bed.From the C1_KI SPH simulations, the particle snapshots with pressure distributions are shown in Figure 7a-f at different time instants after the dam break.Besides, Figure 8 gives the time history of dam break flow impact pressures on the right wall computed by using different particle numbers.Compared with the experimental pressure data [29], not only a good agreement has been found, but also the pressure history of C1_KI SPH computations becomes more reasonable with a higher particle resolution.To show more clearly the improvement of proposed model, Figure 9 gives the pressure errors between the experimental data [29] and ISPH results computed by C1_KI, MLS and SFDI models using different particle numbers N. Again it fully demonstrates that C1_ KI scheme can achieve the lowest error and have the fastest convergence rate.

Solitary Wave Propagation over a Constant Depth
Recently the correlation in between the energy conservation properties and applied pressure gradient models has been extensively explored [30].The solitary wave propagation over a constant water depth is considered in this section.The interaction between tsunami wave and coastal structure is the topic related to practical coastal and ocean engineering problems and therefore attracts increasing attentions.In most cases, the solitary wave is used to represent certain characteristics of the tsunami wave [31,32].Long-distance wave propagation is still a huge challenge to SPH models since the wave form cannot be well reserved due to the particle disorders and numerical dissipations.To verify the robustness of the proposed C1_KI SPH scheme, the computed solitary wave profiles are compared with the analytical solutions derived from the Boussinesq equation referred to Zheng et al. [33].
Here, consider a solitary wave with the wave amplitude a = 0.05 m, water depth d = 0.25 m and total length of the tank L = 30.0m.The computational time step is kept constant as 0.001 s.The To show more clearly the improvement of proposed model, Figure 9 gives the pressure errors between the experimental data [29] and ISPH results computed by C1_KI, MLS and SFDI models using different particle numbers N. Again it fully demonstrates that C1_ KI scheme can achieve the lowest error and have the fastest convergence rate.To show more clearly the improvement of proposed model, Figure 9 gives the pressure errors between the experimental data [29] and ISPH results computed by C1_KI, MLS and SFDI models using different particle numbers N. Again it fully demonstrates that C1_ KI scheme can achieve the lowest error and have the fastest convergence rate.

Solitary Wave Propagation over a Constant Depth
Recently the correlation in between the energy conservation properties and applied pressure gradient models has been extensively explored [30].The solitary wave propagation over a constant water depth is considered in this section.The interaction between tsunami wave and coastal structure is the topic related to practical coastal and ocean engineering problems and therefore attracts increasing attentions.In most cases, the solitary wave is used to represent certain characteristics of the tsunami wave [31,32].Long-distance wave propagation is still a huge challenge to SPH models since the wave form cannot be well reserved due to the particle disorders and numerical dissipations.To verify the robustness of the proposed C1_KI SPH scheme, the computed solitary wave profiles are compared with the analytical solutions derived from the Boussinesq equation referred to Zheng et al. [33].
Here, consider a solitary wave with the wave amplitude a = 0.05 m, water depth d = 0.25 m and total length of the tank L = 30.0m.The computational time step is kept constant as 0.001 s.The

Solitary Wave Propagation over a Constant Depth
Recently the correlation in between the energy conservation properties and applied pressure gradient models has been extensively explored [30].The solitary wave propagation over a constant water depth is considered in this section.The interaction between tsunami wave and coastal structure is the topic related to practical coastal and ocean engineering problems and therefore attracts increasing attentions.In most cases, the solitary wave is used to represent certain characteristics of the tsunami wave [31,32].Long-distance wave propagation is still a huge challenge to SPH models since the wave form cannot be well reserved due to the particle disorders and numerical dissipations.To verify the robustness of the proposed C1_KI SPH scheme, the computed solitary wave profiles are compared with the analytical solutions derived from the Boussinesq equation referred to Zheng et al. [33].
Here, consider a solitary wave with the wave amplitude a = 0.05 m, water depth d = 0.25 m and total length of the tank L = 30.0m.The computational time step is kept constant as 0.001 s.The numerical solitary wave is generated by a piston-type wave maker according to the theory given in [34], in which the motion of wave maker is defined in a dimensionless form by with the dimensionless wave celerity c = √ 1 + a. Figure 10 shows the computed free surface profiles with the analytical solutions for different particle numbers at different times of the wave propagation.It shows that the numerical wave surfaces approach to the analytical ones with the increasing particle number in vertical direction Ny.Meanwhile, the convergence of model computations is evidenced by the close agreement among the three numerical results.Especially it is promising to note that both the wave height and wave shape are well maintained during the wave propagation even if the wave has travelled nearly 30 m over a water depth of 0.25 m.This implies that the dampening of wave height over long-distance travel could be attributed to the influence of pressure gradient calculation schemes.Moreover, in addition to the possible numerical dissipations it has also been shown in the literature that some functions to reproduce the solitary waves lead to a progressive decay in the wave amplitude along the channel as well as the presence of trailing waves due to instantaneous truncation in the motion law of the wave maker.Finally, to show the model convergence more clearly, Figure 11 gives the errors of wave surface elevation between C1_KI SPH results computed with different particle numbers and the analytical results at time t = 15.32 s.
Water 2017, 9, 400 12 of 17 numerical solitary wave is generated by a piston-type wave maker according to the theory given in [34], in which the motion of wave maker is defined in a dimensionless form by wave celerity 1 c a = + .Figure 10 shows the computed free surface profiles with the analytical solutions for different particle numbers at different times of the wave propagation.It shows that the numerical wave surfaces approach to the analytical ones with the increasing particle number in vertical direction Ny .Meanwhile, the convergence of model computations is evidenced by the close agreement among the three numerical results.Especially it is promising to note that both the wave height and wave shape are well maintained during the wave propagation even if the wave has travelled nearly 30 m over a water depth of 0.25 m.This implies that the dampening of wave height over long-distance travel could be attributed to the influence of pressure gradient calculation schemes.Moreover, in addition to the possible numerical dissipations it has also been shown in the literature that some functions to reproduce the solitary waves lead to a progressive decay in the wave amplitude along the channel as well as the presence of trailing waves due to instantaneous truncation in the motion law of the wave maker.Finally, to show the model convergence more clearly, Figure 11 gives the errors of wave surface elevation between C1_KI SPH results computed with different particle numbers and the analytical results at time t = 15.32 s.Water 2017, 9, 400 12 of 17 numerical solitary wave is generated by a piston-type wave maker according to the theory given in [34], in which the motion of wave maker is defined in a dimensionless form by wave celerity 1 c a = + .Figure 10 shows the computed free surface profiles with the analytical solutions for different particle numbers at different times of the wave propagation.It shows that the numerical wave surfaces approach to the analytical ones with the increasing particle number in vertical direction Ny .Meanwhile, the convergence of model computations is evidenced by the close agreement among the three numerical results.Especially it is promising to note that both the wave height and wave shape are well maintained during the wave propagation even if the wave has travelled nearly 30 m over a water depth of 0.25 m.This implies that the dampening of wave height over long-distance travel could be attributed to the influence of pressure gradient calculation schemes.Moreover, in addition to the possible numerical dissipations it has also been shown in the literature that some functions to reproduce the solitary waves lead to a progressive decay in the wave amplitude along the channel as well as the presence of trailing waves due to instantaneous truncation in the motion law of the wave maker.Finally, to show the model convergence more clearly, Figure 11 gives the errors of wave surface elevation between C1_KI SPH results computed with different particle numbers and the analytical results at time t = 15.32 s.

Solitary Wave Impacting on Vertical and Inclined Walls
In order to further show the effectiveness of improved C1_KI SPH technique, an investigation is made on the numerical results of a solitary wave propagation and impact on a solid wall with different Water 2017, 9, 400 13 of 17 inclination angles.The physical experiment was carried out by Zheng et al. [35] in a 3-D wave flume with piston wave maker in Harbin Engineering University (HEU).The schematic diagram of the wave tank is shown in Figure 12.The wave tank length is L = 10 m and the water depth is d = 0.25 m.The solitary wave height is a = 0.15 m and the wave nonlinearity is ε = a/d = 0.6.As shown in Figure 12, two pressure sensor points are located on the right wall at a distance of 0.05 m and 0.15 m from the tank bottom to monitor the pressure time history, i.e., h 1 = 0.05 m and h 2 = 0.10 m.Besides, two wave elevation gauges are located at section of Set1 and Set2, which are 2.0 m away from the left and right boundaries, respectively.In the C1_KI SPH computations, the initial particle spacing is chosen to be 0.01 m and the time step is kept constant as 0.001 s.

Solitary Wave Impacting on Vertical and Inclined Walls
In order to further show the effectiveness of improved C1_KI SPH technique, an investigation is made on the numerical results of a solitary wave propagation and impact on a solid wall with different inclination angles.The physical experiment was carried out by Zheng et al. [35]    In order to show the pressure time histories more clearly, only the dynamic parts of the pressures are provided, while the total pressure should be obtained by adding to the part of static pressure.From the results in both figures, it shows that C1_KI SPH computations can achieve good results as compared with the experimental data [35].Since there is a reflected wave after the solitary wave impacts on the slope, it shows double peaks in the wave elevation time history at Set2 in Figure 13.Due to the wave running up and down, there are some discrepancies in the numerical results after the solitary wave reflects from the slope.The pressure time histories at point P1 in Figure 14 are not as good as those at point P2, as there are some small oscillations in the peak domain.This is due to that some small water drops hit upon the water surface and make the pressure of inner particles change rapidly.The pressure amplitude at P1 and P2 are almost the same.It should be noted that the double peak pressure patterns in Figure 14 have also been reported in a latest ISPH work [36].
For generality, Figures 15a,b and 16a,b give the comparisons of wave surface elevations at Set1 and Set2, and time histories of wave impact pressures at P1 and P2, respectively, when the slope angle is α = 120°.As there exist the overturning and re-entering waves, the amplitude of reflected waves becomes smaller than the case of Figure 13.Due to being difficult to simulate the air-pocket and small bubbles in present model, the amplitude and phase of waves demonstrate some differences between the experimental data [35] and numerical ISPH results.With the increase of slope inclinations, i.e., from α = 90° to α = 120°, the amplitude of reflected wave elevations and wave impact pressures becomes smaller, and the duration of wave impact becomes longer, following the comparisons between Figures 13 and 15, or between Figures 14 and 16.Due to the 3D effect of wave propagation and small vibrations from the piston wave maker, there is a phenomenon of trailing waves appearing in Figures 13a and 15a.
Finally, Figure 17a,b gives a snapshot of the computed wave profile with the experimental photo at time t = 6 s, including the pressure contour distributions in the fluid domain.It shows again that the wave elevation profile obtained by C1_KI SPH can achieve a good agreement with the laboratory photograph, and the pressure distribution of wave field is quite stable and noise-free.In order to show the pressure time histories more clearly, only the dynamic parts of the pressures are provided, while the total pressure should be obtained by adding to the part of static pressure.From the results in both figures, it shows that C1_KI SPH computations can achieve good results as compared with the experimental data [35].Since there is a reflected wave after the solitary wave impacts on the slope, it shows double peaks in the wave elevation time history at Set2 in Figure 13.Due to the wave running up and down, there are some discrepancies in the numerical results after the solitary wave reflects from the slope.The pressure time histories at point P1 in Figure 14 are not as good as those at point P2, as there are some small oscillations in the peak domain.This is due to that some small water drops hit upon the water surface and make the pressure of inner particles change rapidly.The pressure amplitude at P1 and P2 are almost the same.It should be noted that the double peak pressure patterns in Figure 14 have also been reported in a latest ISPH work [36].
Water 2017, 9, 400 14 of 17   For generality, Figure 15a,b and Figure 16a,b give the comparisons of wave surface elevations at Set1 and Set2, and time histories of wave impact pressures at P1 and P2, respectively, when the slope angle is α = 120 • .As there exist the overturning and re-entering waves, the amplitude of reflected waves becomes smaller than the case of Figure 13.Due to being difficult to simulate the air-pocket and small bubbles in present model, the amplitude and phase of waves demonstrate some differences between the experimental data [35] and numerical ISPH results.With the increase of slope inclinations, i.e., from α = 90 • to α = 120 • , the amplitude of reflected wave elevations and wave impact pressures becomes smaller, and the duration of wave impact becomes longer, following the comparisons between Figures 13 and 15, or between Figures 14 and 16.Due to the 3D effect of wave propagation and small vibrations from the piston wave maker, there is a phenomenon of trailing waves appearing in Figures 13a and 15a.Finally, Figure 17a,b gives a snapshot of the computed wave profile with the experimental photo at time t = 6 s, including the pressure contour distributions in the fluid domain.It shows again that the wave elevation profile obtained by C1_KI SPH can achieve a good agreement with the laboratory photograph, and the pressure distribution of wave field is quite stable and noise-free.

Conclusions and Discussions
An improved ISPH method is proposed to simulate violent water wave propagation and impact upon coastal structures such as vertical and inclined walls.After the comparisons of different first-order derivative calculation models, such as SFDI, MLS and C1_KI, C1_KI SPH has found to perform the best considering both the inner fluid domain and the near boundary region when the particles are randomly distributed.In practical model applications, the results of dam break flow computed by C1_KI SPH show good agreement with the experimental data.Then the robustness of C1_KI SPH is further verified by the solitary wave propagation over a long-distance with almost no wave height dampening, which implies that the accurate pressure calculation could be the key issues in modelling of wave system.Finally, the model is applied to the solitary wave impact on a slope with different orientation angles.After compared with the experimental data for two different slopes angles, i.e., 90° and 120°, the computational results of C1_KI SPH demonstrate its great potential in predicting the wave surface elevations and impact pressures in water wave system.
The main contributions of the paper lie in the following aspects: Firstly, it proposes the improved numerical scheme of C1_KI, which effectively avoids the process of kernel gradient calculation.The relevant matrix is strictly in symmetric distribution, so only the upper triangle elements need to be calculated while others can be obtained from the symmetric relationship.Secondly, a comprehensive analysis of SFDI, MLS and C1_KI benchmarks indicated that the proposed one achieves the most accuracy regardless of whether the particles are regularly or irregularly distributed in the inner fluid area.Considering the solid boundary effect, the accuracy of

Conclusions and Discussions
An improved ISPH method is proposed to simulate violent water wave propagation and impact upon coastal structures such as vertical and inclined walls.After the comparisons of different first-order derivative calculation models, such as SFDI, MLS and C1_KI, C1_KI SPH has found to perform the best considering both the inner fluid domain and the near boundary region when the particles are randomly distributed.In practical model applications, the results of dam break flow computed by C1_KI SPH show good agreement with the experimental data.Then the robustness of C1_KI SPH is further verified by the solitary wave propagation over a long-distance with almost no wave height dampening, which implies that the accurate pressure calculation could be the key issues in modelling of wave system.Finally, the model is applied to the solitary wave impact on a slope with different orientation angles.After compared with the experimental data for two different slopes angles, i.e., 90° and 120°, the computational results of C1_KI SPH demonstrate its great potential in predicting the wave surface elevations and impact pressures in water wave system.
The main contributions of the paper lie in the following aspects: Firstly, it proposes the improved numerical scheme of C1_KI, which effectively avoids the process of kernel gradient calculation.The relevant matrix is strictly in symmetric distribution, so only the upper triangle elements need to be calculated while others can be obtained from the symmetric relationship.Secondly, a comprehensive analysis of SFDI, MLS and C1_KI benchmarks indicated that the proposed one achieves the most accuracy regardless of whether the particles are regularly or irregularly distributed in the inner fluid area.Considering the solid boundary effect, the accuracy of

Conclusions and Discussions
An improved ISPH method is proposed to simulate violent water wave propagation and impact upon coastal structures such as vertical and inclined walls.After the comparisons of different first-order derivative calculation models, such as SFDI, MLS and C1_KI, C1_KI SPH has found to perform the best considering both the inner fluid domain and the near boundary region when the particles are randomly distributed.In practical model applications, the results of dam break flow computed by C1_KI SPH show good agreement with the experimental data.Then the robustness of C1_KI SPH is further verified by the solitary wave propagation over a long-distance with almost no wave height dampening, which implies that the accurate pressure calculation could be the key issues in modelling of wave system.Finally, the model is applied to the solitary wave impact on a slope with different orientation angles.After compared with the experimental data for two different slopes angles, i.e., 90 • and 120 • , the computational results of C1_KI SPH demonstrate its great potential in predicting the wave surface elevations and impact pressures in water wave system.
The main contributions of the paper lie in the following aspects: Firstly, it proposes the improved numerical scheme of C1_KI, which effectively avoids the process of kernel gradient calculation.The relevant matrix is strictly in symmetric distribution, so only the upper triangle elements need to be calculated while others can be obtained from the symmetric relationship.Secondly, a comprehensive analysis of SFDI, MLS and C1_KI benchmarks indicated that the proposed one achieves the most accuracy regardless of whether the particles are regularly or irregularly distributed in the inner fluid area.Considering the solid boundary effect, the accuracy of C1_KI is comparable with that of the existing approaches, but its numerical scheme is more straightforward in terms of the formulation and is more efficient in terms of CPU time costs.Thirdly, it is the first time that C1_KI scheme and truly incompressible SPH are combined for the water wave propagation and impact simulations, which has been evidenced by the high accuracy results and low dissipations of the wave peaks as well as stable pressure distributions.

Figure 3 .
Figure 3. Convergence and error analysis of first-order derivative calculation under disordered particle distributions in the: (a) inner fluid area; and (b) near boundary area.

Figure 2 .
Figure 2. Convergence and error analysis of first-order derivative calculation under uniform particle distributions in the: (a) inner fluid area; and (b) near boundary area.

Figure 3 .
Figure 3. Convergence and error analysis of first-order derivative calculation under disordered particle distributions in the: (a) inner fluid area; and (b) near boundary area.

Figure 3 .
Figure 3. Convergence and error analysis of first-order derivative calculation under disordered particle distributions in the: (a) inner fluid area; and (b) near boundary area.

Figure 4 .
Figure 4. Schematic setup of dam break flow domain.

Figure 5 .
Figure 5. Comparisons between C1_KI SPH computations and experimental results on dam break flow: (a) leading edge; and (b) water column height.

Figure 4 .
Figure 4. Schematic setup of dam break flow domain.
the following analysis, all the variables are non-dimensionalised by using the dam width a and gravitational acceleration g , such as /

Figure 4 .Figure 5 .
Figure 4. Schematic setup of dam break flow domain.In the present simulations it is assumed a = 0.5 m, / 2.0 H a = and / 4 L a = .The total particle number is 60 × 120 with a particle size of 0.00833 m.The non-dimensionalised time step is given by / (2 ) 0.003 t t g a Δ = Δ =  .To validate the proposed C1_KI SPH computations, Figure5a,b gives the comparisons of dam break flow leading edge and water column height, respectively, with the experimental results from[28].A good agreement is found in both cases.

Figure 5 .
Figure 5. Comparisons between C1_KI SPH computations and experimental results on dam break flow: (a) leading edge; and (b) water column height.

Figure 6 .
Figure 6.Time history of computed pressures for different water column height and particle number: (a) H a = ; (b) 2 H a =

Figure 6 .
Figure 6.Time history of computed pressures for different water column height and particle number: (a) H a = ; (b) 2 H a = ; and (c) 4 H a = .

Figure 9 .
Figure 9. Comparisons of pressure errors for C1_KI, MLS and SFDI models with different particle numbers.

Figure 8 .
Figure 8. Pressure time history on right wall computed by C1_KI SPH with different particle numbers, compared with experimental data [29].

Figure 8 .
Figure 8. Pressure time history on right wall computed by C1_KI SPH with different particle numbers, compared with experimental data [29].

Figure 9 .
Figure 9. Comparisons of pressure errors for C1_KI, MLS and SFDI models with different particle numbers.

Figure 9 .
Figure 9. Comparisons of pressure errors for C1_KI, MLS and SFDI models with different particle numbers.

Figure 10 .Figure 11 .
Figure 10.Comparisons of wave surface profile computed by C1_KI SPH using different particle resolutions with analytical solutions at different wave propagation stages.

Figure 10 .
Figure 10.Comparisons of wave surface profile computed by C1_KI SPH using different particle resolutions with analytical solutions at different wave propagation stages.

Figure 10 .Figure 11 .
Figure 10.Comparisons of wave surface profile computed by C1_KI SPH using different particle resolutions with analytical solutions at different wave propagation stages.

Figure 11 .
Figure 11.Convergence study of wave surface elevations computed by C1_KI SPH with different particle resolutions.

6 . 2 h
in a 3-D wave flume with piston wave maker in Harbin Engineering University (HEU).The schematic diagram of the wave tank is shown in Figure 12.The wave tank length is L = 10 m and the water depth is d = 0.25 m.The solitary wave height is a = 0.15 m and the wave nonlinearity is a d As shown in Figure 12, two pressure sensor points are located on the right wall at a distance of 0.05 m and 0.15 m from the tank bottom to monitor the pressure time history, i.e., 1 h = 0.05 m and = 0.10 m.Besides, two wave elevation gauges are located at section of Set1 and Set2, which are2.0 m away from the left and right boundaries, respectively.In the C1_KI SPH computations, the initial particle spacing is chosen to be 0.01 m and the time step is kept constant as 0.001 s.

Figure 12 .
Figure 12.Schematic view of solitary wave tank and measurement locations.

Figure
Figure 13a,b gives the comparison of wave surface elevations at Set1 and Set2, and Figure 14a,b gives the comparison of pressure time histories at P1 and P2, for the slope inclination angle α = 90°.In order to show the pressure time histories more clearly, only the dynamic parts of the pressures are provided, while the total pressure should be obtained by adding to the part of static pressure.From the results in both figures, it shows that C1_KI SPH computations can achieve good results as compared with the experimental data[35].Since there is a reflected wave after the solitary wave impacts on the slope, it shows double peaks in the wave elevation time history at Set2 in Figure13.Due to the wave running up and down, there are some discrepancies in the numerical results after the solitary wave reflects from the slope.The pressure time histories at point P1 in Figure14are not as good as those at point P2, as there are some small oscillations in the peak domain.This is due to that some small water drops hit upon the water surface and make the pressure of inner particles change rapidly.The pressure amplitude at P1 and P2 are almost the same.It should be noted that the double peak pressure patterns in Figure14have also been reported in a latest ISPH work[36].For generality, Figures15a,b and 16a,b give the comparisons of wave surface elevations at Set1 and Set2, and time histories of wave impact pressures at P1 and P2, respectively, when the slope angle is α = 120°.As there exist the overturning and re-entering waves, the amplitude of reflected waves becomes smaller than the case of Figure13.Due to being difficult to simulate the air-pocket and small bubbles in present model, the amplitude and phase of waves demonstrate some differences between the experimental data[35] and numerical ISPH results.With the increase of slope inclinations, i.e., from α = 90° to α = 120°, the amplitude of reflected wave elevations and wave impact pressures becomes smaller, and the duration of wave impact becomes longer, following the comparisons between Figures13 and 15, or between Figures14 and 16.Due to the 3D effect of wave propagation and small vibrations from the piston wave maker, there is a phenomenon of trailing waves appearing in Figures13a and 15a.Finally, Figure17a,b gives a snapshot of the computed wave profile with the experimental photo at time t = 6 s, including the pressure contour distributions in the fluid domain.It shows again that the wave elevation profile obtained by C1_KI SPH can achieve a good agreement with the laboratory photograph, and the pressure distribution of wave field is quite stable and noise-free.

Figure 12 .
Figure 12.Schematic view of solitary wave tank and measurement locations.

Figure
Figure 13a,b gives the comparison of wave surface elevations at Set1 and Set2, and Figure 14a,b gives the comparison of pressure time histories at P1 and P2, for the slope inclination angle α = 90• .In order to show the pressure time histories more clearly, only the dynamic parts of the pressures are provided, while the total pressure should be obtained by adding to the part of static pressure.From the results in both figures, it shows that C1_KI SPH computations can achieve good results as compared with the experimental data[35].Since there is a reflected wave after the solitary wave impacts on the slope, it shows double peaks in the wave elevation time history at Set2 in Figure13.Due to the wave running up and down, there are some discrepancies in the numerical results after the solitary wave reflects from the slope.The pressure time histories at point P1 in Figure14are not as good as those at point P2, as there are some small oscillations in the peak domain.This is due to that some small water drops hit upon the water surface and make the pressure of inner particles change rapidly.The pressure amplitude at P1 and P2 are almost the same.It should be noted that the double peak pressure patterns in Figure14have also been reported in a latest ISPH work[36].

Figure 17 .
Figure 17.Wave surface profile during its run-up at t = 6 s, compared between: (a) C1_KI SPH computations; and (b) Laboratory photograph.