Position explicit and iterative implicit consistent incompressible SPH methods for free surface flow

Accepted Manuscript

Position Explicit and Iterative Implicit Consistent Incompressible SPH Methods for Free Surface Flow Saeed Farzin , Rouhollah Fatehi , Yousef Hassanzadeh PII: DOI: Reference:

S0045-7930(18)30731-X https://doi.org/10.1016/j.compfluid.2018.10.010 CAF 4027

To appear in:

Computers and Fluids

Please cite this article as: Saeed Farzin , Rouhollah Fatehi , Yousef Hassanzadeh , Position Explicit and Iterative Implicit Consistent Incompressible SPH Methods for Free Surface Flow, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.10.010

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights

 

AC

CE

PT

ED

M

AN US



Consistent explicit and iterative implicit consistent incompressible SPH methods Improved conditions for recognition of the particles on the free surface Improvement in stability and accuracy of the pressure distribution results in dam-break problem New schemes with substituting the kernel gradient by fractions of the kernel function.

CR IP T



1

ACCEPTED MANUSCRIPT

Position Explicit and Iterative Implicit Consistent

CR IP T

Incompressible SPH Methods for Free Surface Flow

Saeed Farzina, Rouhollah Fatehib,*, and Yousef Hassanzadehc a

Department of Civil Engineering, University of Semnan, Semnan, Iran.

b

Department of Mechanical Engineering, School of Engineering, Persian Gulf University,

c

AN US

Bushehr 75169, Iran. *Email: [email protected]

Faculty of Water Engineering, Department of Civil Engineering, University of Tabriz, Tabriz,

M

Iran.

Abstract:

ED

In this paper, we propose two consistent incompressible Smoothed Particle Hydrodynamics (ISPH) methods, namely explicit ISPH (EISPH) and iterative implicit

PT

ISPH (IISPH) for simulating free surface flow as challenging complex and highly

CE

nonlinear problem. Both methods employ new consistent spatial derivatives schemes, boundary conditions free of dummy or mirror particle. In the explicit method, a

AC

projection approach is used with one pressure Poisson equation solution at each timestep. In this work, it is shown that this approach does not guarantee the incompressibility. The suggested iterative implicit approach is a remedy for this difficulty. Both methods showed satisfactory results in comparison with the available

2

ACCEPTED MANUSCRIPT

experimental and numerical solutions of the well-known 2-D dam breaking problem over dry and wet beds. However, the proposed IISPH method yields more accurate and stable results rather than EISPH even with a lower number of particles and greater time-

CR IP T

step sizes.

Keywords: Smoothed Particle Hydrodynamics, Incompressible flow, Projection

AN US

method, Consistency, Implicit method, Dam breaking.

1. Introduction

Free surface flows occur in numerous aspects in hydrodynamics, such as water sloshing

M

in tanks, breaking waves in coastal areas and dam-break problems. The study of such highly nonlinear flows with large deformations, discontinuities, and pressure

ED

oscillations is a challenging topic in fluid dynamics.

PT

Using numerical simulations of some suitable simplified models, the researchers can get a lot of significant information for variously complicated flows of two

CE

immiscible fluids. However, dealing with the moving free surfaces especially with large deformations is still a challenge for the mesh-based computational fluid dynamics

AC

methods. Eulerian solvers of Navier-Stokes equations apply a computational grid to solve the problem, and it should be coupled with a mathematical treatment of the free surface. On the other hand, Lagrangian mesh-free methods have been extended to a wide variety of problems including free surface fluid flows involving large

3

ACCEPTED MANUSCRIPT

deformations and with moving boundaries. Such new-generation methods are inherently well-suited to simulate the convection-dominated problems and discrete the governing equations without numerical diffusion. Moreover, these methods improve the flexibility of the code but bring new phenomena to be solved [1, 2].

CR IP T

The Smoothed Particle Hydrodynamics (SPH), due to the Lagrangian nature, is the most well-known particle-based method. SPH is a mesh-free, fully Lagrangian method which is originally introduced in astrophysical research by Lucy [3] and Gingold and

AN US

Monaghan [4] in 1977 and then developed to model a wide range of physical processes. The SPH method utilizes a specified set of interpolated points, named particles, to describe the physical behavior of movement and a kernel function to determine the range of interactions of particles. For each particle in the domain, in each time level, the

M

hydrodynamic properties including density, pressure, and velocity are approximated. There is no need to meshes for calculation of the spatial derivatives, description of

ED

continuous media by discrete particles, and the ability to model morphological changes

PT

are among the most important advantages of this method over other common methods. For example, for modeling multi-phase fluid flows and associated transport phenomena,

CE

the SPH possesses the unique features in comparison to the conventional mesh-based methods [5, 6].

AC

SPH simulation of the incompressible flow is more challenging than the

compressible case and is of main research focus [7]. In SPH literature, several algorithms have been proposed to represent the incompressible flow simulation and can be classified into two types:

4

ACCEPTED MANUSCRIPT

(1) Monaghan [8] in 1994 introduced the first approach, the so-called weakly compressible SPH (WCSPH) algorithm which approximates incompressibility by selecting the smallest possible artificial speed of sound which gives a low Mach number and density fluctuations within 1% [8, 9]. Using the slightly compressible artificial

CR IP T

fluid, Monaghan et al. [10] and Monaghan and Kos [11] made it possible to extend the method to free surface flows. Although this approach simulates some highly transient flows well, notably dam-break flows, pressure distribution is noisy, and the method is

AN US

highly dissipative [12].

(2) In 1999, Cummins and Rudman [13] proposed another approach for incompressible flow using the projection method of Chorin [14] which solves a pressure Poisson equation at every time-step with a source term proportional to the divergence of the

M

velocity field or the density variation. After Shao and Lo [15], this approach is called the incompressible SPH (ISPH) algorithm.

ED

So far, the comparison of the two above-mentioned methods has been done by some

PT

researchers. Lee [16] stated that the CPU time in WCSPH depends on the artificial speed of sound, but there is a linear solver in ISPH. So, the total CPU time of ISPH is

CE

less, by a factor of 2 to 8 depending on the cases. Furthermore, Lee et al. [17] simulated the 2-D lid-driven cavity flow to compare WCSPH and ISPH. Short range pressure

AC

variations produced by WCSPH were so larger than mean pressure variation. Also, the CPU time required by ISPH is shorter in all cases. Xu et al. [7] stated that in many problems, the ISPH method with projection-based pressure is more accurate and stable than the WCSPH method. Furthermore, the pressure field in ISPH is virtually noise-

5

ACCEPTED MANUSCRIPT

free. Another modified ISPH formulation was introduced by Hu and Adams [18] who used their method to simulate multi-phase incompressible flows with a large density difference. They also argued that their method gives more accurate solutions, and it is computationally more efficient than the standard WCSPH method.

CR IP T

WCSPH is a simple and practical approach. However, it may result in non-physical oscillations especially in the pressure and density fields as a result of unwanted decoupling between primitive variables. Fatehi and Manzari [19] have shown that this

AN US

difficulty is not a natural consequence of using the WCSPH method but originates from the inappropriate choice of spatial discretization scheme used for pressure term when solving the mass conservation equation [20].

There are several works on ISPH in the literature including low Reynolds number

M

flows in the channel and around bluff bodies [21], and high Reynolds number free surface flows, e.g. [22, 23]. Shao and Lo [15] used an incompressible scheme, similar to

ED

that in Koshizuka et al. [24], to simulate the free surface flow in 2-D dam-break

PT

problems. Ataie-Ashtiani et al. [25] presented a new form of the source term for the Poisson equation and also a modified Poisson equation of pressure to simulate free

CE

surface problems by an ISPH formulation. This method is used to predict 2-D dambreak and wave propagation. Khayyer et al. [1] proposed a new criterion for more

AC

appropriate identification of free surface particles. In their work, an improved ISPH method is applied to the simulation of dam-break with impact problems. Enhanced wave impact is simulated by this improved method. Xu et al. [26] investigated the operation of the truncated kernel in a free surface simulation with ISPH algorithm. It

6

ACCEPTED MANUSCRIPT

was found that using incomplete kernel support; free surface numerical instability can be developed. This instability could be damped by applying artificial viscosity to particles on or neighboring the free surface. Asai et al. [27] presented an ISPH algorithm to simulate free surface flow problems. Their formulation is based on a new

CR IP T

pressure Poisson equation with a relaxation coefficient, which is estimated by a preanalysis calculation. This algorithm was tested by numerical examples of the dam-break problem. Liu et al. [28] developed an improved 2-D ISPH model to simulate the free

AN US

surface flow-structure interaction. In this model, the mirror particle treatment for solid boundaries is developed. The improvement of the model and pressure computations is validated by a benchmark test of dam-break flow. Hirschler et al. [29] used a new method for open boundary conditions to enable stable simulations using Incompressible

M

SPH based on the renormalization approach of Bonet and Lok [30]. Their proposed

boundary conditions.

ED

method is based on the mirror particle approach already introduced for solid wall

PT

Despite the valuable investigations done during the last decade, ISPH is still an emergent method in numerical simulation. There are still fundamental concepts, for

CE

instance, the treatment of incompressibility, shifting particles and physical boundary conditions that are challenging and need to be examined comprehensively. Moreover,

AC

the choice of spatial discretization is an extremely important step in the numerical solution of flow equations particularly ISPH. To have a convergent numerical method, the discretizing schemes must be consistent. This means that the truncation errors of the schemes must approach zero as the particle spacing decreases. So far, some ways for

7

ACCEPTED MANUSCRIPT

discretizing derivatives in ISPH equations are developed; however, the majority of them are not consistent [31]. In this study, employing the modified consistent first and second derivatives approximation, we have proposed consistent iterative implicit (IISPH) and explicit

CR IP T

(EISPH) projection-based ISPH methods which resolve the above-mentioned problems. Then, to ascertain the performance and capabilities of the presented, modified ISPH methods, dam-break flow problems including both dry and wet bed conditions as the

AN US

well-known benchmark test cases are simulated and validated against existing solutions.

2. SPH formulations

M

2.1. Smoothed Particle Hydrodynamics approximation

The basic principle of the SPH method is based on the concept of interpolation. For an v 

regarding the values at neighboring

ED

arbitrary field function v, the interpolated value

PT

particles v j can be approximated by the direct summation of the relevant quantities of its neighboring particles j as, v(r )    j v jW ( r  rj ; h),

CE

(1)

j

AC

where  j is the weight or the volume of neighboring particle j, r is the position vector, h is the radius of circular compact support, and W is a positive smoothing or kernel function representing a smoothed version of Dirac delta function. The kernel must be at

8

ACCEPTED MANUSCRIPT

least once differentiable, and its derivative must be continuous to avoid large fluctuations which would affect the numerical solution [32]. In this work, the following quadratic kernel function [33] is selected as the

 6 (1  q)2  W (q; h)    h 2 0

CR IP T

smoothing function for 2-D problems.

0  q  1.0 1 q

(2)

where q  r  r j h . It should be noted that the function and its first derivative are

AN US

continuous. The simplicity of this kernel function improves the performance of the computations.

M

2.2. Approximation of spatial derivatives

There are several discretization schemes for the spatial derivatives of the field variables

ED

in SPH [31]. In the following, the new techniques used for spatial derivatives are briefly

PT

described.

CE

2.2.1. The first derivative

For discretization of the first derivative in SPH, here, a new consistent scheme which

AC

has the property of first-order consistency is introduced. This property means that the scheme predicts the first derivative of linear functions exactly. The proposed numerical approximation of the first derivative

v  i

can be rewritten in vector notation as,

9

ACCEPTED MANUSCRIPT

 W ij v  i    j B 1i  (e ij )(v j v i ), rij j

(3)

where Wij  W ( ri  rj ; h) is the kernel function of particle i at the position of particle j, is the vector distance between two particles i and j and eij  rij rij is a unit

vector in the inter-particle direction from j to i. Also, 1

  B 1i     j W ij e ij e ij  ,  j  

CR IP T

rij  ri  r j

(4)

AN US

is a symmetric renormalization tensor for the first derivative. The use of renormalization tensor in SPH first derivative, first is appeared in the work of Randles and Libersky [34].

It is notable that in Eq. (3), W ij (e ij ) is used instead of the term W ij which is used

M

rij

commonly in the literate after Gingold and Monaghan [4]. According to the definition

ED

of the gradient, it is clear that here, W ij is replaced by W ij . As it is shown in Fig. 1, rij

rij

PT

these two terms behave similarly when the particles are far enough from each other [33]. However, in cases where rij is small, the derivative of the kernel function W ij

CE

rij

tends to zero for commonly used bell-shape functions such as B-splines and Wendland

AC

[35] functions. Although for the quadratic function defined in Eq. (2), it is a linear function, however, in calculating pressure gradient, this behavior reduces the repulsive force between two approaching particles with a very small distance. This, in turn, may lead to particle clustering. As it is observed in Fig. 1, the proposed term W ij strictly rij

10

ACCEPTED MANUSCRIPT

increases by reducing the particle distance rij . This property improves the performance of the proposed method in heterogeneous distributions of particles.

W Fig. 1. Comparisons of the kernel function W, its derivative W , and the ratio versus r for r

CR IP T

r

the quadratic function defined in Eq. (2).

AN US

Also, proportional to changes in Eq. (3), renormalization tensor (Eq. (4)) is defined differently from what is presented in the literature (for example [34] and [30]).

M

2.2.2. The second derivative

The stability properties of SPH equations are significantly depending on the second

ED

derivative of discretization which links to the tensile instability [16]. For the second derivative, most of SPH schemes are not consistent [31]. In this paper, we modified the

PT

consistent second derivative approximation scheme recently introduced by Fatehi and Manzari [31] by replacing the gradient of the kernel function with the term W ij (e ij ) like

CE

rij

AC

that used in the previous section for the first derivative. So, the modified scheme reads as,



2v  i  B 2i :  2 j j

v j v i e ij e ij (  e ij  v  i ), rij rij

W ij

11

(5)

ACCEPTED MANUSCRIPT

where the operator ‘‘:’’ denotes the inner product of second-order tensors and

v  i

is



computed according to Eq. (3). Also B 2 i is a renormalization tensor which is computed using the following set of equations:

CR IP T

     W B 2i :   j W ij e ij e ij e ij e ij  (  j ij e ij e ij e ij )  B 1i  (  j W ij r ij e ij e ij e ij )   I . rij j j   j

(6)

Eq. (6) is a set of three and six equations in two- and three dimensions, respectively. By 

solving this system, one will obtain the components of the symmetric tensor B 2 i . The

AN US

difference between Eq. (6) and the definition of the second renormalization tensor in [31] is due to the replacing gradient of the kernel with the term W ij (e ij ) . rij

Similar to the discussion in [31], using Taylor series expansion about the point i, it

M

can be shown that the proposed schemes (3) and (5) are consistent and converge linearly as x  0 for a constant ratio h , where x denotes the SPH particles spacing. This

ED

x

means that they can exactly predict the first and the second derivative of linear and

PT

quadratic functions, respectively. Although the mentioned proposed methods are more

CE

complex than the usual schemes. In the following, these schemes are used to discretize the spatial derivatives in the

AC

governing equations.

3. Governing equations

12

ACCEPTED MANUSCRIPT

The governing equations of viscous fluid flows are mass and momentum conservation equations. In an ISPH method, they are presented in the Lagrangian form and thus the Navier–Stokes equations read as, d  0 or  u  0, dt

CR IP T

(7)

du 1   p  g  0 2u , dt 

(8)

where u is the velocity vector, p is the pressure, t is the time, ρ is the density, g is the

AN US

acceleration due to body forces (e.g., gravity), and 0 is the viscosity of the fluid.

To prevent numerical instabilities an artificial viscosity a  ca u x is added to the physical viscosity 0 in which

ca

is constant and chosen as 0.2. This proposed relation

M

is similar to the artificial viscosity which is used in Lax-Wendroff central scheme [36]. Of course, this is different from the scheme suggested by Monaghan [37] in the context

ED

of compressible flow. Further, the use of artificial viscosity eliminates the need for a turbulent model in this method. It is noteworthy that the commonly used the sub-

PT

particle scale (SPS) turbulence model [38] also adds an amount of diffusion to the

CE

momentum equation through eddy viscosity. However, in practice, the resolution of particles in the test cases of the SPS model (for instance [38-40]) is not high enough to

AC

be comparable with LES.

4. The pressure projection method

13

ACCEPTED MANUSCRIPT

In the present study, an improved solution for free surface flows, based on the ISPH methodology is described in the following sections. The base approach is similar to the SPH projection method proposed by Cummins and Rudman [13] for the flows without a free surface. ISPH employs a strict incompressible SPH formulation, in which the

CR IP T

pressure is not an explicit thermodynamic variable from an equation of the state, but obtained by solving a pressure Poisson equation from a hydrodynamic formulation [28]. The pressure projection method is used to enforce incompressibility by

AN US

decomposing the velocity field. The implementation of the pressure projection method in a Lagrangian frame of reference introduces mutual dependencies between velocity, discretization points and spatial stencils that are not present in a Eulerian frame of reference [41].

M

To solve these equations, here, two explicit and iterative implicit methods are

ED

applied.

PT

4.1. Explicit Incompressible SPH (EISPH) method

CE

In the incompressible SPH computation, the governing Eqs. (7) and (8) are solved using a two-step prediction and correction processes. In most works done in this field, ISPH is

AC

implemented similarly to the following: 1) The first step is an explicit integration of velocity in time-based on viscous and gravity forces in the momentum equation (Eq. (8)). For this purpose, first, an

14

ACCEPTED MANUSCRIPT

(u )

is predicted without including the pressure gradient

u   u n  t ( g  eff  2u n n ),

(9)

intermediate particle velocity term in Eq. (8).

t n , t

is the time increment, and eff is the sum

of physical and artificial viscosities. It is notable that the superscript n in

un

CR IP T

where u n is particle velocity at the time

indicates the value for the current time,

while in   n , n represents the SPH approximated derivatives (Eqs. (3) and (5)) based

AN US

on the position of the particles at the current time. Distinguishing these two points is important. In Eulerian grid-based methods, the position is independent of time. However, in Lagrangian methods, like SPH, the values of spatial derivatives are approximated according to particles’ positions. Therefore, it is necessary that the time at

M

which particles’ positions used in approximation be determined.

ED

2) Now, the pressure at time-step n+1 is calculated using the following Poisson

Therefore,



t

 u   n .

CE

2 p n 1  n 

PT

equation so that its Laplacian, neutralizes the divergence of the intermediate velocity.

( p n 1 ) ,

a system of linear equations must be

AC

To obtain the new time-level pressure

(10)

solved. It is also necessary that the Eq. (5) be rewritten so that the pressure coefficient matrix can be extracted easily and efficiently. One can write:

W ij   W ij  W  2 p  i  B 2i :  2 j  2 e ij e ij  ( k  ik e ik e ik e ik )  B 1i  e ij ( p j  p i ). rik rij j k  rij 

15

(11)

ACCEPTED MANUSCRIPT

Compared to conventional Eulerian methods, the coefficient matrix of these linear equations has more (between 10 to 30) non-zero elements in each row. The matrix is not symmetric, so the equations are solved using a Generalized Minimal Residual (GMRes)

employed to solve the Poisson pressure equation.

CR IP T

method. In this work, a sparse solver based on the well-known code BLAS [42] is

3) Next, in the correction processes, the new pressure

( p n 1 )

is used to enforce

incompressibility in the calculations. The velocity at time-step n+1,

u n 1  u   t (

1



AN US

as, p n 1  n ).

(u n 1 )

is obtained

(12)

The intermediate particle velocity (u  ) is usually not divergence-free but if the new

M

pressure exactly satisfies Eq. (10) one has  u n 1  n  0.

(13)

ED

4) Using the corrected velocity field at time-step n+1, the new positions of all fluid

t n (u  u n 1 ). 2

(14)

CE

r n 1  r n 

PT

particles are calculated as,

Finally, at the end of each round of computations, neighbor particle lists are

AC

updated.

According to the methodologies of calculating the spatial derivatives, that all of them

are calculated using the position at the current time, the above method is called Explicit

16

ACCEPTED MANUSCRIPT

ISPH (EISPH) in the present paper. It should be noted that here, the term "Explicit" is not related to the calculated pressure in Eq. (10). 4.2. Iterative Implicit Incompressible SPH (IISPH) method

CR IP T

Noting Eq. (13), it can be seen that the field variable and derivative position are calculated at two different times, while the incompressibility condition necessitates that   u n 1 n 1

to be zero. This condition is not achievable using any explicit method.

Although one can update the position during the processes (after step (1)), however, the

AN US

position of differentiation is just slightly modified.

Nevertheless, to achieve the incompressibility condition, an iterative approach is necessary. This approach has been previously pointed out by Hu and Adams [18] and

M

then Xu et al. [7]. Hu and Adams indicated that if only a divergence-free velocity is enforced, a large density variation will produce. The algorithm used in [7] is based on a

ED

density-invariant field without considering the effect of the position of differentiation. During the simulation, the particle density will be recalculated and the pressure gradient

PT

will be corrected in an iteration loop. However, in both references, the effect of changes

CE

of spatial derivatives of the pressure and velocity field is not considered in modifying the particle position.

AC

Here, to resolve this problem, a remedy called Iterative Implicit ISPH (IISPH) method is introduced and implemented as the followings. In the following equations, the first super-scripts denote the time level and the second super-scripts are iteration number.

17

ACCEPTED MANUSCRIPT

1) First, the values of position, velocity and pressure in the first iteration step k=0, are set to the values of the current time as,

r n1,0  r n , un1,0  un , pn1,0  pn .

(15) (u  )

is predicted using the values of iteration

CR IP T

2) Then, the intermediate particle velocity k as, u  un  t ( g  eff 2 un1,k  n1,k ).

(16)

3) The new pressure for iteration k+1, ( p n1,k 1 ) is calculated by solving the following

2 p n1,k 1  n1,k 

 t

AN US

system of equations as,   u  n1,k .

(17)

4) Now, the values of velocity and position of all fluid particles in the new iteration are



p n 1,k 1  n 1,k ),

t n ( u  un1,k 1 ). 2

(18)

(19)

PT

r n1,k 1  r n 

1

ED

un 1,k 1  u  t (

M

calculated as,

CE

5) If the iteration is converged, leaves the loop. For the convergence criterion in the above iteration loop, there are different choices. Here, we have chosen p n1,k 1  p n1,k  

AC

as the convergence criterion. In this condition, one may conclude that the velocities and positions of particles are converged as well. If the condition is not satisfied, by updating r k 1   r r k 1  (1   r )r k , un 1,k 1  u un 1,k 1  (1  u )un 1,k , p n 1,k 1   p p n 1,k 1  (1   p ) p n 1,k ,

18

(20)

ACCEPTED MANUSCRIPT

It returns to step 2. Here,  r , u , and  p are under-relaxation factors. To avoid the divergence of trial and error, under-relaxation is used with coefficients 0.4, 0.4 and 0.3, for the position, velocity, and pressure respectively.

particle at the new time-level is obtained as,

r n1  r n1,k 1 , un1  un1,k 1 , pn1  pn1,k 1.

CR IP T

6) Finally, after convergence, the values of position, velocity, and pressure of each

(21)

AN US

4.3. Numerical treatment of initial and boundary conditions 4.3.1. Initial conditions

For all the cases considered in this work, the initial velocity and pressure of the fluid

M

particles are considered zero in the computational domain at time t=0. The pressure field tends quickly to the hydrostatic profile after a few time-steps.

ED

Fluid particles were initially created in the form of a fixed rectangular grid with

PT

equidistant particle spacing. For boundary particles, the particle spacing is the same as the initial particle spacing between fluid particles. For each particle, a quadratic kernel

CE

is used with the same smoothing length h  3x where x is the initial particle spacing. Since in this method, the cut-off length of the kernel function is equal to the smoothing

AC

length, choosing h  3x gives rather a small number (i.e., 10 to 30) of neighboring particles. This leads to lower cost of computation and better performance.

4.3.2. Wall boundary conditions

19

ACCEPTED MANUSCRIPT

In SPH context, there are several methods in the literature for modeling solid wall boundaries such as the curtain breakwater and sea bottom. This condition is modeled by fixed solid or dummy particles or by mirror particles method that mirrors the properties of inner fluid particles. However, the solid boundaries require special care and the no-

CR IP T

slip condition must be imposed. A suitable solid boundary condition approach should lead to a fairly smooth and accurate pressure field around the solid surface.

In this paper, a recent no-slip solid boundary treatment method [20, 43] is used. In

AN US

the case of the no-slip condition and fixed walls, as the fluid particles on the wall have the same velocity as the wall, their acceleration is zero. Thus, neglecting the normal component of the force due to viscous stresses on the wall, one can write

p .n   g .n ,

(22)

M

where n is the unit vector normal to the wall. For the boundary particle i, ni can be

ED

calculated from the summation of the kernel functions multiplied by unit vectors in the inter-particle direction as,

 

j

 jW ij e ij

.

(23)

PT

ni 

 jW ij e ij j

CE

The SPH particles arrangement and examples of normal vectors n for particles on

AC

the wall boundary are shown in Fig. 2.

Fig. 2. A schematic of the SPH particles arrangement near a wall boundary and free surfaces

Discretizing the pressure gradient term in Eq. (22) for particle i, leads to

20

ACCEPTED MANUSCRIPT

  W ij   j j  B 1i  r (eij )   ni ( p j  pi )   g  ni . ij  

(24)

It must be noted that the implementation of the proposed solid boundary treatment will be possible only if the kernel truncation problem in the neighborhood of the solid

CR IP T

boundaries is solved. In the current work, the use of consistent renormalized first and second derivatives schemes helps to avoid this problem. The good performance of this algorithm has been shown in [20, 43] and there is no need for dummy or mirror

AN US

particles.

4.3.3. Free surfaces condition

In the current study, free surfaces condition is simplified to

(25)

M

τ  n  0.

That is

ED

(u   u T )  n  0,

(26)

PT

which is treated as a condition on the velocity field. Also, a condition for the pressure field is necessary. In the absence of the surface

CE

tension effect, the pressures of both phases are equal on the interface. Here, by

AC

neglecting the effect of air pressure, the condition is the zero pressure for free surface particles

p  p0  0.

(27)

In this work, a particle is identified to be on the free surface if satisfies the four following equations:

21

ACCEPTED MANUSCRIPT

W

ij



0.80

i

j

W

(28 a)

,

e  0.25,

(28 b)

ij ij

j

S   14,

(28 c)

CR IP T

B2 

B 1  0,

where S

(28 d)













AN US

B2



 ( B 2 xx  B 2 yy  B 2 xy ) in which B 2 xx , B 2 yy and B 2 xy are elements of the

second renormalization tensor. It should be noted that the constants in the right-hand side of equations (28 a) to (28 d) may depend on the size of the smoothing radius h.

M

The above four conditions are interpreted as follows: a) In the case of a particle located within the fluid,

W is approximately equal to ij

1 i

.

ED

j

There are no dummy or mirror particles for the treatment of the free surface boundary.

PT

So, for the particles located on the free surface, because of having less neighboring

CE

particles, the summation of neighbors kernel function is less than the above amount. b) In the case of a particle is located within the fluid, neighboring particles are normally

AC

distributed uniformly so that the value of

 W j

e

ij ij

is almost zero. Relatively high

j

values for this expression mean that the particle is near the free surface.

22

ACCEPTED MANUSCRIPT





c) For a particle within the fluid, tensors B 1 and B 2 are generally close to unit tensor. But for particles near the boundary or on the surface, the components of these two 

tensor increases. Condition (28 c) is a criterion which indicates that the tensor B 2 is

CR IP T

large enough. d) For those particles on the wall where there is no fluid particle nearby and also for individual or group of detached particles, it is possible that none of the previous

4.4. Determination of time-step size

AN US

conditions is satisfied. Condition (28 d) is employed to consider these particles.

M

Since explicit time integration is performed using the ISPH projection method, the sizes of time- steps must be limited by stability constraints to get accurate results. For stable

ED

time integration, three time-step size criteria must be satisfied, including the Courant– Friedrichs–Lewy condition, the viscous-diffusion condition and a limitation due to the

PT

body force. In this work, due to relatively large Reynolds numbers, only the first

CE

condition must be checked.

x t   ( ), U max

AC

(29)

in which 0    1.0 is a constant, x is the initial distance between two neighboring SPH particles, and U max is the maximum velocity in the flow. However, because of the nonlinearity of the governing equations, the case   1.0 may not reachable, and the upper

23

ACCEPTED MANUSCRIPT

limit depends on the solution method. Nevertheless, for the cases in this paper, it is between 0.05 and 0.25.

CR IP T

4.5. Shifting SPH particles In the implementation of the pressure projection algorithm as a Lagrangian method, moving the particles according to the fluid velocity field may introduce some particle splitting and stability problems. Especially in the present method which uses less

AN US

neighboring particles for calculation of spatial derivatives [20]. A remedy is using a modified velocity closer to the average velocity of the neighboring particles and another is to shift particles after each time-step. In this study, a slightly different approach from what is used in [20] and [7] is employed. This approach leads to obtaining an

M

approximately uniform distribution of particles in the computational fluid domain.

ED

Furthermore, application of this technique helps to overcome spurious oscillations and instabilities. The shifting distance should be large enough to prevent instability and

PT

small enough not to cause inaccuracy due to the Taylor series correction [7]. For particles on the free surface, this algorithm is not suitable. Thus a modified method

CE

which includes only surface particles in the computation of the shifting distance is used.

AC

This approach also has a tuning parameter  [7] which in this study is between 0.01 and 0.1.

The procedures of the current explicit and iterative implicit projection methods are

summarized in Tables 1 and 2.

24

ACCEPTED MANUSCRIPT

Table 1 Summary of the proposed explicit algorithm (EISPH).

Initialize particles’ velocity, position, pressure, and the time: u0, r0, P0 at time t0; for each time-step n do

// The prediction step: for each particle i do Calculate the intermediate particle velocity

(u  ) by Eq.

(9);

AN US

end for

CR IP T

Find the neighboring particles;

Obtain the pressure by solving Poisson Eq. (10); // The correction step: for each particle i do

(u n 1 ) by Eq.

(12);

M

Calculate the new particle velocities

ED

Calculate the new position of particles (r n 1 ) by Eq. (14); Shift the position of particles;

PT

Interpolate the velocities and pressures at the new position; end for

AC

CE

end for

Table 2 Summary of the proposed iterative implicit algorithm (IISPH).

Initialize particles’ velocity, position, pressure, and the time: u0, r0, P0 at time t0;

for each time-step n do

25

ACCEPTED MANUSCRIPT

Initialize the values of position, velocity and pressure at the first iteration step k=0 for each iteration k do Find the neighboring particles;

CR IP T

// The prediction step: for each particle i do

Calculate the intermediate particle velocity

AN US

end for

(u  )

by Eq. (16);

Obtain the pressure by solving Poisson Eq. (17); // The correction step: for each particle i do

M

Calculate new particle velocities

(un1,k 1 ) by Eq.

ED

Calculate new position of particles end for

PT

if the iteration is converged, then Leave the loop;

CE

else

AC

Update the iteration values

end if end for for each particle i do

26

(18);

(r n1,k 1 ) by Eq.

(19);

ACCEPTED MANUSCRIPT

Obtain the values of position, velocity and pressure at the new time-level by Eq. (21); Shift the position of particles; Interpolate the velocities and pressures at the new position;

CR IP T

end for end for

AN US

5. Results and discussions

In the present work, dam-break flow problem, as a well-known benchmark test case for two-dimensional incompressible flow problems are solved to demonstrate the performance and capabilities of the proposed ISPH methods. The main cause of the

M

failure of a dam is due to the flood caused by the rupture of a dam and sudden emptying

ED

of the reservoir which is basically different from natural floods for the following reasons: production of very high peak discharge values in a very short time period,

PT

occurrence of shock waves, creation of important dynamic effects and transport of a

CE

significant volume of solid materials. So, before the occurrence of this calamitous event, preparation of a flood map with maximum water levels and the water front arrivals for

AC

each prototype is necessary [44]. First of all, we show the effectiveness of the schemes presented in equations (3) to (6) by solving a simple test case. 5.1. Comparison of kernel functions

27

ACCEPTED MANUSCRIPT

Here, we solve the fully developed flow in a straight 2D channel. This very simple benchmark usually is used to prevent the difficulties such as treatment of the free surface, shifting particles, oscillating pressure field, to make effects on the solution. So, only the discretization schemes and the kernel function are influential. The test case consists of two constant parallel solid walls

CR IP T

making a 2D conduit. The fluid flows in the horizontal direction by exerting a constant body force. After a while, the flow reaches its steady state with a velocity profile equals to (

) in which

and

( )

are dimensionless position and horizontal velocity,

respectively.

AN US

To solve the problem numerically, a channel of 2h length with periodic boundary conditions in the horizontal direction is considered. Here,

is the smoothing radius and cut-off length

of the kernel functions as well. Also, just 20 particles are placed in the height of the channel with uniform spacing

to facilitate the evaluation and comparison of the accuracy.

M

Six commonly used functions are considered for the kernel. Two of them are splines with 3rdand 5th- order polynomials. The super Gaussian function is a combination of a polynomial and

ED

the Gaussian function which was introduced by Monaghan [37]. The other functions are simple

PT

polynomials which are also used in SPH methods. Table 3 compares the results of the method using different kernel functions. To compare the

CE

results, all the numerical parameter is set identical. Then, the error at each particle (

) is

calculated by subtracting numerical results of the dimensionless horizontal velocity from its

AC

exact value and L∞-norm is Max(| |) and

∑| | . |

L1-norm is ∑ |

Since the error are small, the

difference cannot recognized visually. So, they are reported quantitatively. The results in table 3 show that the quartic function of Lucy [3], cubic spline by Monaghan and Lattanzio [45], and the quadratic function of equation (2) performed better in term of accuracy. The good accuracy of the quadratic function of Johnson and Beissel [33] is rather unexpected. It

28

ACCEPTED MANUSCRIPT

is usually considered as a case with lack of precision, because as it was shown in Fig. 1, its first derivative is not smooth at r=h. However, since in this paper, we use the alternative form of Eq. (3) instead if the commonly used kernel gradient, this low-order function leads to results as

CR IP T

accurate as those of cubic spline and quartic function.

Table 3 CPU time and errors of the results of the presented SPH method for the fully developed flow in a channel using various kernel functions.

Relative error norms of velocity

CPU Time Kernel Function (s)

Quartic Lucy [3]

127

Cubic spline [45]

131

Quintic spline [9]

129

AN US

134

1.1E-03

1.1E-03

6.2E-04

1.3E-03

7.4E-04

6.7E-03

3.4E-03

125

1.7E-03

9.4E-04

125

3.0E-03

1.7E-03

PT

Quintic Wendland [35]

ED

Quadratic [33]

L∞

2.2E-03

M

Super Gaussian [37]

L1

In Table 3, the CPU times of the method are also compared. It can be observed that the

CE

computational cost is almost equivalent. The reason is that evaluating the kernel function has a small share in the computations in comparison with those of renormalization tensors. However,

AC

as expected, the splines and the super Gaussian have slightly more computational costs. Based on this investigation, the quadratic function of equation (2) is used as the kernel in the following numerical simulations.

29

ACCEPTED MANUSCRIPT

5.2. Simulation of the dam breaking on a dry bed In this section, to ascertain the capability of the proposed schemes in dealing with large deformation of the free surface, they are applied to simulate dam-break problems

CR IP T

over the dry and wet bed in two spatial dimensions and then test cases are used to validate the codes. The initial conditions are plotted in Fig. 2.

5.2.1. Martin and Moyce [46] experiment

AN US

We first used dam-break experiments have been carried out by Martin and Moyce [46] where the initial reservoir height in Fig. 3 is H=L=5.7cm, and the length and height of

M

the channel are d=5.366H and D=3.0H.

ED

Fig. 3. Schematic sketch of the initial set-up for two-dimensional dam-break on the dry bed

PT

Many characteristics of dam-break flow depend strongly on whether the front propagates on dry channel bed or over an ambient fluid layer [47]. The position of the

CE

surge front toe is named xfront. The solutions of this problem reported in the literature are including experimental data [46], analytical Ritter solution by the assumption of

AC

shallow-water conditions [48] and numerical solutions by Colagrossi and Landrini [49] and Ferrari et al. [50]. Our EISPH and IISPH numerical results of xfront is obtained using 3898 and 9320 fluid particles and plotted versus time in Fig. 5, compared with the experimental data

30

ACCEPTED MANUSCRIPT

[46] and the numerical solution by Ferrari et al. [50]. Here, xfront is computed by the fluid particle with the maximum x position, and t=0 is the time instant for the dam break. It is worth mentioning that Ferrari et al. [50] simulated this problem assuming inviscid flow in three dimensions using 2M SPH particles.

CR IP T

Here, the constant of shifting SPH particles and the maximum Courant number are considered as   0.1 and   0.1 , respectively. Also, the physical viscosity of the fluid in

AN US

all cases in this paper is set to 0.001 Pa.s.

Fig. 4. Comparison of time history of the Surge front position from the literature and the

M

proposed methods EISPH (left) and IISPH (right)

The experimental data in Fig. 4 is faced with a slower water front advancement. It is

ED

most likely due to ignoring the roughness effects on the solid boundaries of the tank in the numerical methods. This alters the propagation velocity and triggers the

PT

development of turbulence near the water front [51]. Another reason could be the

CE

wetting of the surface, where surface tension damps the collapsing dynamics. It can be seen in Fig. 4 that all the numerical results of the proposed ISPH methods

AC

are in reasonable agreement even with the relatively small number of particles. Over time, the roughness of the walls will produce more difference between experimental and simulated propagation velocities.

31

ACCEPTED MANUSCRIPT

5.2.2. Zhou et al. [52] experiment To more examine the proposed methods, we used another test case reproduces the experiments detailed by Zhou et al. [52], where the initial reservoir height H in Fig. 3 is equal to 60 cm, L=2H, and the length and height of the channel are d=5.366H and

CR IP T

D=3.0H. Our explicit and iterative implicit numerical solutions have been computed for water using 5155 and 11187 fluid particles. Figs. 5 and 6 show our EISPH and IISPH results for particles’ pressure and the evolution of the free surface at different times

AN US

(t=0.7, 1.2, 1.4, 1.6, 1.8 and 1.9 s).

Fig. 5. Dam-break flow and impacts against the vertical wall by EISPH (left) and IISPH (right) methods at times t=0.7, 1.2, 1.4, 1.6, 1.8 and 1.9 s from top to bottom using 5155 particles

M

Fig. 6. Dam-break flow and impacts against the vertical wall by EISPH (left) and IISPH (right)

ED

methods at times t=0.7, 1.2, 1.4, 1.6, 1.8 and 1.9 s from top to bottom using 11187 particles

In Fig. 5 it can be seen that the pressure distributions for both EISPH and IISPH

PT

methods are smooth, oscillation-free and in acceptable ranges. Nevertheless, in EISPH

CE

method, some pressure overshoots are observable at for example t=1.2, 1.6 and 1.9 s. However, the EISPH results illustrate a significant improvement in solution compared

AC

to the classical WCSPH results. In contrast, the results of IISPH are away from this defect. Furthermore, this separation of particles after hitting the wall which is observed in EISPH results at times after 1.2 s (see Fig. 6), is not present in the case of IISPH.

32

ACCEPTED MANUSCRIPT

On the other hand, the comparison of Fig. 5 with Fig. 6 shows that increasing the number of particles to about double, the results are somewhat improved . Also, In the first row of Fig. 6, which corresponds to the time 0.7 s, the wave fronts is ahead of the ones at the same time in Fig. 6. Moreover, comparing Figs. 5 and 6 at the time of 1.9 s,

CR IP T

one can observe that wave breaking after the impact in EISPH method results is visibly improved.

In the following, the time history of free surface elevation at locations x=2.228 m

AN US

and x=2.725 m from the left vertical wall is investigated. The obtained results of the numerical simulation using the proposed EISPH and ISPH methods have been compared with the experimental data [52], WCSPH solution of Ferrari et al. [50], and the commercial code Fluent by Ferrari et al. [50] using the Volume of Fluid method in

M

Figs. 7 and 8 for x=2.228 m and in Figs. 9 and 10 for x=2.725 m.

PT

ED

Fig. 7. Time evolution of the free surface height at location x=2.228 m using EISPH method

CE

Fig. 8. Time evolution of the free surface height at location x=2.228 m using IISPH method

AC

Fig. 9. Time evolution of the free surface height at location x=2.725 m using EISPH method

Fig. 10. Time evolution of the free surface height at location x=2.725 m using IISPH method

33

ACCEPTED MANUSCRIPT

As it is shown in Figs. 7 to 10, the obtained numerical free surface heights of the proposed ISPH methods follow the experimental measurements well in comparison with other illustrated solutions. The differences between numerical and experimental results may be due to details of the initial condition in the experiments and the friction

CR IP T

on the bottom.

Although the iterative implicit method leads to more computational costs than the explicit one, it is more flexible and suitable for demonstrating the characteristics of

twice of that of the explicit method.

AN US

these violent flows. Also, it is more stable and can be used with time-step sizes up to

Next, for more accurate assessment of the capability and efficiency of the proposed methods, the value of pressure on the vertical wall during the impact is studied. In Fig.

M

11, the experimental data [52], numerical simulation by Colagrossi and Landrini [49] and Ferrari et al. [50], as the examples of the well-known solutions in the case, are

ED

presented. The number of particles used for modeling is 4900 and 39072 particles,

PT

respectively for free surface and two-phase flow methods in [49], and 2M particles in [50]. In the experiment, a circular shaped gauge of 9 cm diameter has been used, located

CE

on the vertical wall with the center 0.16 cm above the bed [49]. The experimental sensor cannot determine the pressure evolution exactly at a point on the wall, but the impact

AC

pressure measurement is averaged over time and location. Instead, it can be directly calculated numerically at any point of the wall. In impact problems, as the spatial and temporal pressure gradients are high, any of the approaches to average the pressure provides significantly different results [50]. However, it is observable in Fig. 11 that

34

ACCEPTED MANUSCRIPT

over time, the impact of the water front against the vertical wall produces a jump in the pressure field and a second pressure peak has been induced by the backward plunging water front [49, 50]. It seems that the differences between numerical and experimental data, particularly

CR IP T

at the initial times of the motion, are related to details of the initial experimental conditions and neglecting the roughness effects on the bottom. Also, it should be noted that as it was shown in Figs. 5 and 6, the pressure results are smooth enough, and the

AN US

fluctuations in the results in Fig. 11 are of the physics itself. This may also be concluded from the experimental values in this figure.

In Fig. 12 the obtained results of averaged pressure at the sensor’s position on the wall using EISPH are shown compared to the experimental data [52] and in Fig. 13 the

M

related results of IISPH is compared with the same experimental data and numerical results by Ferrari et al. [50]. Here, the averaged pressure is calculated similar to the

ED

experiment by averaging the pressure of the particles located in the interval of 115 to

PT

205 mm from the bottom of the right wall. Considering Fig. 12, in our EISPH approach, increasing the number of particles will dramatically reduce pressure fluctuations. Also,

CE

for the case of 11187 particles, the peak values of the pressure with the experiments are in good agreement.

AC

Another improving capability of the proposed IISPH method is presented in Fig. 13.

As expected, employing IISPH method leads to much less unphysical pressure fluctuations in the case of 5155 particles compared to the EISPH method of Fig. 12. Also, the first and second pressure peak values agree well with the experimental

35

ACCEPTED MANUSCRIPT

modeling. Furthermore, in Fig. 13, it is visible that increasing the number of computational particles causes more accurate and smoother results for pressure.

methods and experimental data

CR IP T

Fig. 11. Comparison of pressure evolution on the vertical wall; numerical results of the SPH

AN US

Fig. 12. Pressure evolution on the vertical wall using the EISPH method

Fig. 13. Pressure evolution on the vertical wall using the IISPH methods

M

It should be noted that besides its better accuracy, the implicit method has more

ED

computational costs than the explicit method. In our numerical experiments, typically, 10 to 15 iterations are needed in each time-step of the implicit method. Although each

PT

iteration approximately has a cost of a time-step in the explicit method, however, the

CE

overall computational cost of the implicit method is less than 10 times of the explicit method because here, largest time-step sizes (up to twice of those of the explicit one)

AC

are applicable.

5.3. Simulation of the dam breaking on a wet bed

36

ACCEPTED MANUSCRIPT

Dam-break on the wet bed, because of its distinctive, interesting features, is a powerful numerical benchmark test to assess the performance of numerical methods in free surface flows simulations [53, 54]. Here, a dam-break experiment on a wet bed considered by Janosi et al. [47] is simulated to ascertain the validity of the proposed

CR IP T

methods. Fig. 14 shows the schematic sketch of the initial conditions, where the initial upstream and downstream of reservoir height are du=0.15 m and dd= 0.038 m, respectively and upstream and downstream reservoir length are lu=0.38 m and ld= 3.55

AN US

m, respectively.

Fig. 14. Schematic sketch of the initial set-up for two-dimensional dam-break on the wet bed

The aforementioned experiment has been simulated using both EISPH and IISPH

M

methods on 11905 particles. Here, In EISPH method, the constant of shifting SPH

ED

particles and the maximum Courant number are chosen   0.1 and   0.05 , respectively. Also, in IISPH case, these parameters are set to be   0.08 and   0.12 , respectively.

PT

The obtained results of the free surface evolution by the proposed methods comparing to experimental results [47] are illustrated in Fig. 15 for various times.

CE

Comparison of the results shows that both methods reproduce the experimental surfaces

AC

mostly accurate and particularly over the time, the accuracy increases. Also, in this test case, like the earlier ones, the results of IISPH method are significantly more accurate than those of EISPH. It should be stated that the difference detected between our methods and experimental data at time 0.28 s are most likely due to relatively low removal velocity

37

ACCEPTED MANUSCRIPT

of the gate. As reported in [55], the velocity of the gate in the experiment was 1.5ms-1. This means that the removal of the gate takes about 0.1 s which is comparable to the time of the simulation.

CR IP T

Fig. 15. Profiles of dam breaking flow over a wet bed by the proposed EISPH and IISPH methods and experimental data at times t= 0.28, 0.34, 0.46 and 0.59 s from left to right

AN US

6. Conclusions

In the present work, we proposed explicit (EISPH) and iterative implicit (IISPH) consistent ISPH methods for solving free surface flow problems. Some of the methods'



M

properties are listed below;

In both methods, consistent discretization schemes are employed. Therefore, for a

ED

problem in which the methods are stable in time, it is expected that increasing the



PT

number of computational particles leads to a convergent solution. The proposed EISPH method, apart from the discretization of spatial derivatives, is

CE

similar to the prevalent ISPH methods in other features, especially, the use of a pressure Poisson equation to satisfy the incompressibility condition. Here, it is

AC

shown that in such methods the time at which the positions of particles are estimated in the process of approximating the divergence of the new time-level velocity is the old time-level. So, satisfying .u n 1  n  0 does not guarantee the incompressibility. Although, this difference reduces with the time-step size, however, very small time-

38

ACCEPTED MANUSCRIPT

step sizes, in turn, may destroy an important advantage of the ISPH. One of the main advantages which can make ISPH method more appropriate than WCSPH is that in ISPH, the artificial speed of sound is not used and thus a larger time-step size can be considered. To resolve the above difficulty, here, an iterative implicit approach has been

CR IP T



suggested. In this approach, in each time step, an iteration loop is executed. Thus, at the end of each time-step, both velocity and the positions used in divergence



AN US

estimation are at the new time-level so, one can expect .u n 1  n 1  0 .

It is clear that this iterative implicit method has more computational cost than the explicit one. However, it provides the significant improvement in the stability of the solutions and pressure distribution results.

M

To display the capabilities and performance of our improved ISPH methods, the well-

ED

known dam-break flow problems including both dry and wet bed conditions were simulated and the suggested methods were validated against analytical, experimental

PT

and existed numerical solutions. It was shown that both presented methods using the modified discretization schemes can handle the violating free surface flow problems

CE

accurately. The results of IISPH method, even with the relatively small number of

AC

particles, are smooth, especially for pressure distribution. In comparison with EISPH, the proposed IISPH method yields more accurate and stable results even with the lower number of particles and greater time-step sizes.

References

39

ACCEPTED MANUSCRIPT

[1] Khayyer, A., Gotoh, H., & Shao, S. (2009). Enhanced predictions of wave impact pressure by improved incompressible SPH methods. Applied Ocean Research, 31(2), 111-131. [2] Lobovsky, L., Vimmr J. (2007). Smoothed Particle Hydrodynamics and finite volume modelling of incompressible fluid flow. Mathematics and Computer in Simulation, 76: 124-

CR IP T

131. [3] Lucy, L. B. (1977). A numerical approach to the testing of the fission hypothesis. The astronomical journal, 82, 1013-1024.

[4] Gingold, R. A., & Monaghan, J. J. (1977). Smoothed particle hydrodynamics: theory and

AN US

application to non-spherical stars. Monthly notices of the royal astronomical society, 181(3), 375-389.

[5] Shadloo, M., Oger, G., & Le Touzé, D. (2016). Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and

M

challenges. Computers & Fluids, 136, 11-34.

[6] Fatehi, R., Shadloo, M. S., & Manzari, M. T. (2014). Numerical investigation of two-phase

ED

secondary Kelvin–Helmholtz instability. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 228(11), 1913-1924.

PT

[7] Xu, R., Stansby, P., & Laurence, D. (2009). Accuracy and stability in incompressible SPH

CE

(ISPH) based on the projection method and a new approach. Journal of Computational Physics, 228(18), 6703-6725.

AC

[8] Monaghan, J. J. (1994). Simulating free surface flows with SPH. Journal of Computational Physics, 110(2), 399-406.

[9] Morris, J. P., Fox, P. J., & Zhu, Y. (1997). Modeling low Reynolds number incompressible flows using SPH. Journal of Computational Physics, 136(1), 214-226.

40

ACCEPTED MANUSCRIPT

[10] Monaghan, JJ., Thompson, MC., and Hourigan, K. (1994). Simulation of free surface flows with SPH. Proceedings of ASME Symposium on Computational Methods for Fluid Dynamics, Lake Tahoe, 19-23. [11] Monaghan, J., & Kos, A. (1999). Solitary waves on a Cretan beach. Journal of waterway,

CR IP T

port, coastal, and ocean engineering, 125(3), 145-155. [12] Lind, S., Xu, R., Stansby, P., & Rogers, B. D. (2012). Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. Journal of Computational

AN US

Physics, 231(4), 1499-1523.

[13] Cummins, S. J., & Rudman, M. (1999). An SPH projection method. Journal of Computational Physics, 152(2), 584-607.

[14] Chorin, A. J. (1968). Numerical solution of the Navier-Stokes equations. Mathematics of

M

computation, 22(104), 745-762.

[15] Shao, S., & Lo, E. Y. (2003). Incompressible SPH method for simulating Newtonian and

ED

non-Newtonian flows with a free surface. Advances in water resources, 26(7), 787-800. [16] Lee, E.-S. (2007). Truly incompressible approach for computing incompressible flow in

PT

SPH and comparisons with the traditional weakly compressible approach: University

CE

of Manchester.

[17] Lee, E.-S., Moulinec, C., Xu, R., Violeau, D., Laurence, D., & Stansby, P. (2008).

AC

Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. Journal of Computational Physics, 227(18), 8417-8436.

[18] Hu, X., & Adams, N. A. (2007). An incompressible multi-phase SPH method. Journal of Computational Physics, 227(1), 264-278.

41

ACCEPTED MANUSCRIPT

[19] Fatehi, R., & Manzari, M. T. (2011). A remedy for numerical oscillations in weakly compressible smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids, 67(9), 1100-1114. [20] Fatehi, R., & Manzari, M. (2012). A consistent and fast weakly compressible smoothed

Numerical Methods in Fluids, 68(7), 905-921.

CR IP T

particle hydrodynamics with a new wall boundary condition. International Journal for

[21] Shadloo, M. S., Zainali, A., Sadek, S. H., & Yildiz, M. (2011). Improved incompressible smoothed particle hydrodynamics method for simulating flow around bluff bodies.

AN US

Computer methods in applied mechanics and engineering, 200(9-12), 1008-1020.

[22] Rafiee, A., Pistani, F., & Thiagarajan, K. (2011). Study of liquid sloshing: numerical and experimental approach. Computational Mechanics, 47(1), 65-75. [23] López, Y. R., Roose, D., & Morfa, C. R. (2013). Dynamic particle refinement in SPH:

Mechanics, 51(5), 731-741.

M

application to free surface flow and non-cohesive soil simulations. Computational

ED

[24] Koshizuka, S., Nobe, A., & Oka, Y. (1998). Numerical analysis of breaking waves using the moving particle semi-implicit method. International Journal for Numerical Methods in

PT

Fluids, 26(7), 751-769.

CE

[25] Ataie-Ashtiani, B., Shobeyri, G., & Farhadi, L. (2008). Modified incompressible SPH method for simulating free surface problems. Fluid Dynamics Research, 40(9), 637-

AC

661.

[26] Xu, R., & Stansby, P. (2010). The influence of the truncated kernel to free-surface predictions with ISPH and a new solution. Paper presented at the Proceedings of the 5th International SPHERIC workshop.

42

ACCEPTED MANUSCRIPT

[27] Asai, M., Aly, A. M., Sonoda, Y., & Sakai, Y. (2012). A stabilized incompressible SPH method by relaxing the density invariance condition. Journal of Applied Mathematics, 2012. [28] Liu, X., Xu, H., Shao, S., & Lin, P. (2013). An improved incompressible SPH model for

CR IP T

simulation of wave–structure interaction. Computers & Fluids, 71, 113-123. [29] Hirschler, M., Kunz, P., Huber, M., Hahn, F., & Nieken, U. (2016). Open boundary conditions for ISPH and their application to micro-flow. Journal of Computational Physics, 307, 614-633.

AN US

[30] Bonet, J., & Lok, T.-S. (1999). Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Computer methods in applied mechanics and engineering, 180(1-2), 97-115.

[31] Fatehi, R., & Manzari, M. (2011). Error estimation in smoothed particle hydrodynamics

M

and a new scheme for second derivatives. computers & Mathematics with Applications, 61(2), 482-498.

ED

[32] Gingold, R., & Monaghan, J. (1982). Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics, 46(3), 429-453.

PT

[33] Johnson, G. R., & Beissel, S. R. (1996). Normalized smoothing functions for SPH impact

CE

computations. International Journal for Numerical Methods in Engineering, 39(16), 2725-2741.

AC

[34] Randles, P., & Libersky, L. (1996). Smoothed particle hydrodynamics: some recent improvements and applications. Computer methods in applied mechanics and engineering, 139(1-4), 375-408.

43

ACCEPTED MANUSCRIPT

[35] Wendland, H. (1995). Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in computational Mathematics, 4(1), 389396. [36] Lax, P., & Wendroff, B. (1960). Systems of conservation laws. Communications on Pure

CR IP T

and Applied mathematics, 13(2), 217-237. [37] Monaghan, J. J. (1992). Smoothed particle hydrodynamics. Annual review of astronomy and astrophysics, 30(1), 543-574.

[38] Edmond YM Lo, Shao S (2002) Simulation of near-shore solitary wave mechanics by an

AN US

incompressible SPH method. Applied Ocean Research, 24: 275-286.

[39] Shao, S. (2009). Incompressible SPH simulation of water entry of a free‐falling object. International Journal for Numerical Methods in Fluids, 59(1), 91-115. [40] Gotoh, H., Shao, S., & Memita, T. (2004). SPH-LES model for numerical investigation of

M

wave interaction with partially immersed breakwater. Coastal Engineering Journal, 46(01), 39-63.

ED

[41] Bøckmann, A., Shipilova, O., & Skeie, G. (2012). Incompressible SPH for free surface flows. Computers & Fluids, 67, 138-151.

PT

[42] Blackford, L. S., Petitet, A., Pozo, R., Remington, K., Whaley, R. C., Demmel, J., . . .

CE

Henry, G. (2002). An updated set of basic linear algebra subprograms (BLAS). ACM Transactions on Mathematical Software, 28(2), 135-151.

AC

[43] Hashemi, M., Fatehi, R., & Manzari, M. (2012). A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows. International Journal of Non-Linear Mechanics, 47(6), 626-638.

[44] Hassanzadeh, Y. (1997). Rapidly varied unsteady flow in a small scale dry bed model. International Journal of Engineering, 10(1), 1.

44

ACCEPTED MANUSCRIPT

[45] Monaghan, J. J., & Lattanzio, J. C. (1985). A refined particle method for astrophysical problems. Astronomy and astrophysics, 149, 135-143. [46] Martin, J., Moyce, W., Price, A., & Thornhill, C. (1952). Part V. An experimental study of the collapse of fluid columns on a rigid horizontal plane, in a medium of lower, but

CR IP T

comparable, density. Phil. Trans. R. Soc. Lond. A, 244(882), 325-334. [47] Jánosi, I. M., Jan, D., Szabó, K. G., & Tél, T. (2004). Turbulent drag reduction in dambreak flows. Experiments in Fluids, 37(2), 219-229.

[48] Ritter, A. (1892). Die Fortpflanzung der Wasserwellen. Vereine Deutcher Ingenieure

AN US

Zeitswchrift, 36, 947-954

[49] Colagrossi, A., & Landrini, M. (2003). Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics, 191(2), 448-475. [50] Ferrari, A., Dumbser, M., Toro, E. F., & Armanini, A. (2009). A new 3D parallel SPH

M

scheme for free surface flows. Computers & Fluids, 38(6), 1203-1217. [51] Colagrossi, A. (2005). A meshless Lagrangian method for free-surface and interface flows

ED

with fragmentation. These, Universita di Roma. [52] Zhou, Z., De Kat, J., & Buchner, B. (1999). A nonlinear 3-D approach to simulate green

PT

water dynamics on deck. Paper presented at the Proc. 7th International Symposium on

CE

Numerical Ship Hydrodynamics, Report. [53] Stansby, P., Chegini, A., & Barnes, T. (1998). The initial stages of dam-break flow.

AC

Journal of Fluid Mechanics, 374, 407-424.

[54] Khayyer, A., & Gotoh, H. (2010). On particle-based simulation of a dam-break over a wet bed. Journal of Hydraulic Research, 48(2), 238-249.

[55] Crespo, A. J. C. (2008). Application of the smoothed particle hydrodynamics model SPHysics to free surface hydrodynamics. Ph. D. thesis, University of De Vigo.

45

ACCEPTED MANUSCRIPT

5 Kernel function W 4 Kernel derivative

CR IP T

3

The ratio

2

AN US

1 0 0

0.2

0.4

r

0.6

0.8

1

M

W Fig. 1. Comparisons of the kernel function W, its derivative W , and the ratio versus r for r

AC

CE

PT

ED

the quadratic function defined in Eq. (2).

46

r

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

FIG. 2. A schematic of the SPH particles arrangement near a wall boundary and free surfaces

FIG. 3. Schematic sketch of the initial set-up for two-dimensional dam break on dry bed

47

CR IP T

ACCEPTED MANUSCRIPT

FIG. 4. Comparison of time history of the Surge front position from the literature and the proposed

AC

CE

PT

ED

M

AN US

methods EISPH (left) and IISPH (right)

48

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

49

CR IP T

ACCEPTED MANUSCRIPT

FIG. 5. Dam break flow and impacts against the vertical wall by EISPH (left) and IISPH (right) methods

AC

CE

PT

ED

M

AN US

at times t=0.7, 1.2, 1.4, 1.6, 1.8 and 1.9 s from top to bottom using 5155 particles

50

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

51

CR IP T

ACCEPTED MANUSCRIPT

FIG. 6. Dam break flow and impacts against the vertical wall by EISPH (left) and IISPH (right) methods

CE

PT

ED

M

AN US

at times t=0.7, 1.2, 1.4, 1.6, 1.8 and 1.9 s from top to bottom using 11187 particles

AC

FIG. 7. Time evolution of the free surface height at location x=2.228 m using EISPH method

52

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

FIG. 8. Time evolution of the free surface height at location x=2.228 m using IISPH method

53

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

FIG. 9. Time evolution of the free surface height at location x=2.725 m using EISPH method

AC

CE

PT

ED

FIG. 10. Time evolution of the free surface height at location x=2.725 m using IISPH method

54

AN US

CR IP T

ACCEPTED MANUSCRIPT

FIG. 11. Comparison of pressure evolution on the vertical wall; numerical results of the SPH methods

AC

CE

PT

ED

M

and experimental data

55

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

FIG. 12. Pressure evolution on the vertical wall using the EISPH method

56

ACCEPTED MANUSCRIPT

AN US

CR IP T

FIG. 13. Pressure evolution on the vertical wall using the IISPH methods

AC

CE

PT

ED

M

FIG. 14. Schematic sketch of the initial set-up for two-dimensional dam break on wet bed

57

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 15. Profiles of dam breaking flow over a wet bed by the proposed EISPH and IISPH

AC

CE

PT

ED

M

methods and experimental data at times t= 0.28, 0.34, 0.46 and 0.59 s from left to right

58