A Robust volume conservative divergence-free ISPH framework for free-surface flow problems

Advances in Water Resources 96 (2016) 423–437

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

A Robust volume conservative divergence-free ISPH framework for free-surface flow problems Gourabananda Pahar a, Anirban Dhar b,∗ a b

Department of Civil Engineering, Indian Institute of Technology Kharagpur, Kharagpur WB 721302, India Department of Civil Engineering, Indian Institute of Technology Kharagpur, Kharagpur WB 721302, India

a r t i c l e

i n f o

Article history: Received 28 January 2016 Revised 17 August 2016 Accepted 22 August 2016 Available online 23 August 2016 Keywords: Flow through porous media Incompressible smoothed particle hydrodynamics Density current Volume conservation

a b s t r a c t This study presents a Volume Conservative approach for resolving volume conservation issue in divergence-free incompressible Smoothed Particle Hydrodynamics (ISPH). Irregular free surface deformation may introduce error in volume computation, which has a cascading effect over time. Proposed correction decreases this numerical compressibility to a minimal value. The correction is obtained directly by solving Navier-Stokes momentum equation. Consequently, the framework does not require any parametric study for mixed source/sink term or iterative solution of pressure Poisson equations. The correction is implemented on four different types of flow: (a) pressurized flow in a closed box, (b) dambreak flow, (c) flow through porous block, (d) lock-exchange flow of immiscible fluids (both free-surface and pressurized flow). All four scenarios are shown to have minimal error compared to pure divergence-free ISPH. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Free surface flow problems are common in natural systems. Underlying physics of the problem becomes difficult to model in presence of free surface boundary. Shallow-water equations use depthaveraged approach to compute hydrodynamic variables. The fluctuation in pressure and vertical velocity cannot be captured precisely by Eulerian frameworks (Pahar and Dhar, 2014) with shallow-water equations. In recent years, mesh-free Smoothed Particle Hydrodynamics (SPH) has become popular for free-surface flow problems. However, SPH suffers some drawbacks in handling problems with incompressibility. Pure incompressibility or divergence-free continuity equation does not consider any density change. Numerical compressibility may appear for pure incompressible SPH model. The error in representative volume increases steadily over time. This may result in unrealistic value. A Volume Conservative ISPH (VC-ISPH) algorithm is proposed to overcome the compressibility problem. SPH was derived by Gingold and Monaghan (1977), and Lucy (1977) for modelling of astrophysical problems. Later, Monaghan (1994) introduced SPH in free surface flow simulation. Pressure field obtained by standard weakly compressible SPH may contain spurious variation. Various correction and renormalization algo-

∗

Corresponding author. Fax: +91 3222 282254. E-mail addresses: [email protected] (G. Pahar), [email protected] (A. Dhar). http://dx.doi.org/10.1016/j.advwatres.2016.08.010 0309-1708/© 2016 Elsevier Ltd. All rights reserved.

[email protected],

rithms need to be applied to get smoothed pressure and density field (Gomez-Gesteira et al., 2010). Incompressible SPH can produce comparatively better result with larger time step (Lee et al., 2008). Incompressibility in SPH is mainly achieved by using a semi-implicit pressure projection scheme (Hosseini et al., 2007). Pressure is derived implicitly by solving Pressure Poisson Equation (PPE). The source/sink term of PPE differs with choice of scheme. Cummins and Rudman (1999) proposed that the source/sink term should contain divergence of intermediate velocity. Their algorithm of solenoidal velocity field provides smooth pressure profile. However, particles tend to proceed along streamlines in pure divergence-free Lagrangian method (Shadloo et al., 2012). An artificial displacement term needs to be added at the end of each time step for removing tensile instability (Xu et al., 2009). Shao and Lo (2003) replaced the source/sink term with intermediate density based on particle position. Density invariance schemes are very effective to conserve proper volume of the fluid considering the sensitivity of density in SPH. The scheme also does not require any artificial displacement to maintain stability. However, purely divergence-free velocity field may not be obtained (Gui et al., 2015). On the other hand, divergence free schemes consider density to be constant at all time levels. Kernel compactness algorithm is not directly utilized unlike density-invariant schemes. Consequently, small errors in density computation may be induced. As no compactness checking algorithm is utilized, the errors may accumulate over time to obtain an infeasible value. Thus, total volume of fluid may not be conserved by the scheme. Various researchers tried to create hybrid schemes to balance both diver-

424

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

gence free and density invariance criteria. Hybrid schemes mainly solves PPE with combined source term containing weighted average of intermediate velocity and density. These schemes are mainly different variants of density invariance scheme. Asai et al. (2012) derived the weight factors based on hydrostatic pressure test. Principle of energy dissipation and height-depth ratio of the flow system were used by Gui et al. (2014) for calculating weight factor. A constant value of 0.5 was suggested by Koh et al. (2013) for sloshing problem. In Moving particle Semi-implicit method, an improved higher order source/sink term was proposed by Khayyer and Gotoh (2011). Hu and Adams (2007) solved multiple poisson equations in a single time step until a suitable convergence criteria is achieved. Their method was later extended (Hu and Adams, 2009) by solving two consecutive PPEs containing divergence-free and density invariant source/sink terms in halftime steps. Very few studies (Gui et al., 2015) have concentrated on quantification of errors in terms of density and velocity-divergence for free-surface flow simulation. However, purely divergence-free volume-conservative approach including multiphase flow through porous media is not reported in the literature. This study presents a volume conservative divergence free incompressible SPH. A corrective displacement is added at the end of each time step. The correction is derived based on total incompressibilty of continuity equation. Divergence-free incompressible SPH algorithm is adopted as base algorithm due to its capability of conserving solenoidal criteria. The corrective framework can be applied with divergence-free ISPH without any parametric study. It also does not need iterative solutions of multiple PPE within a single time-step. Current model is validated for a lid-driven cavity flow with Reynolds number 10 0 0. The proposed model is applied on three two-dimensional (2D) representative numerical examples. Change of volume is measured against basic ISPH algorithm (Xu et al., 2009) for each case. 2. ISPH Model for multiphase flow and flow through porous media 2.1. Governing equations

main with radius 2h. Every particle will interact with the neighboring particles present in the support domain. The amount of interaction is determined by kernel/weight (w) function. Value of kernel decreases with distance from the centre particle and outside support domain it is zero.

f (r ) =

∇ .u = 0

(1a)

Du η = − ∇ P + ηg + ν∇ 2 u + FR Dt ρ

(1b)

u = Darcy velocity vector, P = pressure, ρ = density of the fluid, g = acceleration due to gravity, ν = kinematic viscosity, η = effective porosity inside porous media, FR = resistive force (R/ρ ) due to presence of porous media and divergence of residual stress tensor. For pure free-surface multiphase flow through non-porous region, η becomes unity. FR will represent ρ1 ∇ .τ¯ , where τ¯ is SubParticle-Scale (SPS) tensor. Detailed derivation of the governing equations for both pure fluid and porous regions are presented in Appendix.





f (r )δ (|r − r | )dr



(2)

Dirac delta function is substituted with a distributed kernel.

f (r ) =







f (r )w(|r − r |, h )dr



(3)

Integral formulation can be expressed as summation of discrete quantities in support domain.

f (r ) =



f b w ( |r − rb |, h )vb

(4)

b

such that



w ( |r − rb |, h )vb = 1

(5)

b

where vb = volume associated with particle b. 2.3. Operators Instead of using basic kernel gradient operators, velocity and pressure operators are computed using simple algebraic calculations. Velocity divergence of a particle can be written as,

∇ u|a =

 mb b

ρb

(ub − ua ).∇ wab

(6)

To conserve linear momentum, pressure gradient can be estimated by,

   ∇P Pb Pa mb + ∇ wab |a = ρ ρb2 ρa2 b

(7)

where mb = mass of particle b and ρb = density of particle b. Adami et al. (2012) has given a modified pressure gradient for simulating multiple fluids with different density.

1 For divergence-free incompressible SPH, governing equation needs to be represented in terms of material derivative. Different forms of integrated flow models can be modelled in ISPH, provided the momentum equation contains material derivative of velocity. Cheng et al. (2012) derived an integrated equation for free-surface flow interaction with porous media for Eulerian framework. Generalized Lagrangian form of governing equations for multiphase flow including flow through porous media can be represented as,



ρ

∇ P |a =

 1  2 va + v2b P˜ab ∇ wab ma

(8)

b

where P˜ab = density weighted interparticle averaged pressure. Hu and Adams (2006) proposed that the interparticle average pressure can be expressed as,

P˜ab =

ρa Pb + ρb Pa ρa + ρb

(9)

Near free-surface, particles do not have full kernel support. The kernel derivative can be corrected by using algorithm proposed by Bonet and Lok (1999). Kernel second order derivatives are computed using a combination of first-order kernel derivatives (Morris et al., 1997).

1

ρ

∇ .(∇ P )|a =

 b

2

vb (Pa − Pb )rab .∇ wab ρb |rab |2 + χ

(10)

To avoid zero in denominator, a very small coefficient χ is added(0.0 0 01h2 ). Quintic Wendland (1995) is adopted for computing non-linear flow phenomena.

W (r, h ) =



7 r 1+ h 4 π h2



1−

r 2h

4

,

∀r ∈ [0, 2h]

(11)

2.4. Resistive force 2.2. Smoothed particle hydrodynamics(SPH) In SPH, fluid domain is considered to be constituted of a scattered set of nodes/particles. Each particle has a circular support do-

Resistive force term (FR ) denotes an extra linear and non-linear drag forces inside porous media. In pure fluid domain, the term represents only ρ1 ∇ .τ¯ , where τ¯ is Sub-Particle-Scale (SPS) tensor.

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

425

2.4.2. Turbulence modeling Present study utilizes an explicit Smagorinsky (1963) model. Large scale flow formulation is used for sub-particle scale formulation. Turbulence stress tensor (τ¯ ) can be calculated as,

1

ρ

2 3

τ¯i j = 2νt Si j − κδi j

(13)

where νt = turbulence Eddy Viscosity, δi j = kronecker delta function, Si j = strain rate.

1 Si j = 2



∂ ui ∂ u j + ∂ x j ∂ xi



Turbulence Eddy Viscosity (ν t ) can be evaluated as,

νt = Cs2 δ l 2 |S|

(14)

where δ l = characteristic length of small eddies, Cs = a constant. δ l √ can be computed as xz (Cheng et al., 2012), x, z being initial particle distance along x and z axes. Inside porous media, only volume of pores are taken into consideration. Initial particle spacing is taken as characteristic length of eddies for incompressible fluids. Local strain rate (|S|) can be computed as,

|S | = ( 2Si j Si j )0.5

Fig. 1. Initial Condition for Lid driven Cavity (Scenario-I).

(15)

Turbulent kinetic energy (κ ) can be evaluated as, 2.4.1. Drag forces In porous domain, drag force inside porous media (R/ρ ) is consisted of Darcy (linear) and Forchheimer (non-linear) drag forces (Ghimire, 2009). These parameters are derived based on empirical and/or experimental approach (Gent, 1995; Liu et al., 1999).

FR = −

νη2 us K

+

Fch η3 us |us | √ K

 +

1

ρ

∇ .τ¯

κ=

Cν 2 2 δ l |S | C

(16)

The constants are given by Pan et al. (2012) as Cν = 0.08, C = 1.0, Cs = 0.15. Turbulent kinetic energy can be incorporated into pressure term while solving momentum equation (Shao and Gotoh, 2005). Cs = 0.15 is adopted for all the simulations mentioned in this work.

(12) 2.5. Solution procedure

where us = seepage velocity, K = intrinsic permeability, Fch = Forchheimer’s inertia factor. For computation of intrinsic permeability, the analytical formula suggested by Shao (2010) is used. Resistive force term can vary depending upon choice of porous media and formulation of equations.

Density Invariance schemes in SPH maintains divergence-free velocity criteria by changing intermediate density. Divergence free schemes consider the density of the fluid to be constant at all time levels. This study utilizes a purely divergence-free formulation (Cummins and Rudman, 1999) with full pressure projection

Fig. 2. Velocity distribution in pure and VC-ISPH for lid-driven cavity (Scenario-I).

426

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Fig. 3. Comparison of instantaneous velocity distribution in horizontal and vertical centre-line with Ghia et al. (1982) (Scenario-I).

maintain proper particle arrangement. Resistive force term FR may have different values based on choice of porous media and nature of formulation (Akbari, 2014; Aly and Asai, 2015). However, different choice of FR can be incorporated into the explicit intermediate velocity step. The steps are described below. Data: Particle positions rn , Effective Porosity η, velocities un at time-step n Result: Updated rn+1 , un+1 at time-step n + 1 while t
du ∗ dt u∗ ←

← ηg + ν∇ 2 un + FR

un + t

du ∗ dt





 ∗

Solution of P from ∇ . ρ1 ∇ P n+1 = 1t ∇ . uη un+1 ← u∗ − ηρt ∇ P n+1



t n+1 + un rn+1 ← rn +  2η u



rn ← rn+1 t ← t + t n←n+1 end

Fig. 4. Comparison of kernel and velocity divergence error in pure and VC-ISPH for dambreak-flow (Scenario-I).

In purely divergence-free algorithm, particles march along streamlines to form local flow structure (Xu et al., 2009) leading to particle clumping. This can be corrected by adding a small artificial displacement with particle position at the end of each time level (Shadloo et al., 2012). 2 ra = β umax ra,o b

by keeping constant density at all time levels. First particles are moved to intermediate position using known time level velocity. Intermediate velocity of the fluid particles is computed without enforcing incompressibility on the intermediate particle position. Then the intermediate velocity field is projected into a divergence free and a curl free component. PPE is solved using biconjugate stability method with Jacobi preconditioner. Obtained pressure gradient is added with intermediate velocity to calculate divergencefree velocity at unknown time level. Particles are advected to new position using a second-order time-marching scheme. A small corrective displacement is added at the end of each time step to

ra

|rab |3

(17)

where umax = maximum flow velocity, ra,o = average interparticle spacing, β = a problem dependent constant. Artificial displacement is much smaller than actual displacement.

rab = ra − rb ra,o = b

|rab | Na

where Na = number of particles present in the neighborhood of particle a. Finally, the actual position of the particles is computed

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

427

Fig. 5. Comparison of instantaneous pressure distribution in horizontal and vertical centre-line with Botella and Peyret (1998) (Scenario-I).

Density invariance scheme uses Eqn. 19 to maintain constant density. However, true divergence free nature of the flow may not be achieved. The mass and momentum conservation equations for correction can be written as,

Fig. 6. Initial Condition for dambreak flow (Scenario-II).

as,

rn+1 ← rn+1 + rn+1

(18)

The correction ensures proper particle distribution for the system. However, numerical compressibility may not be removed using the algorithm. Asai et al. (2012) pointed out that the total volume of the fluid may not remain constant for the aforementioned scheme. 2.6. Volume correction General continuity equation for fluid flow can be written as,

1 dρ + ρ dt

∇  .u= 0

(19)

Pure Divergence-free

Divergence-free algorithm solves only the pure velocity divergence-free portion of the continuity equation. Consequently, it may induce error in terms of particle number density. The error will have a cumulative effect with time and may undergo a variation of large unrealistic value, which is not at all compatible with incompressibility of fluid. The corrective displacement (Eqn. 17) obtained explicitly by Shadloo et al. (2012); Xu et al. (2009) does not ensure proper particle number density. Here, the corrective displacement is computed by solving a second PPE. Thus, the resultant numerical compressibility may be kept within a suitable limit by considering a corrective velocity (uc ) in the continuity equation.

0 1 dρ > + ∇ .uc + ∇  . u=0 ρ dt

(20 )

1 dρ + ∇ .uc = 0 ρ dt

(21a)

1 Du c = − ∇ P∗ Dt ρ

(21b)

where uc = corrective velocity, P ∗ = corrective pressure required to minimize numerical compressibility. For a single fluid system, density based on particle position (ρ ∗ ) within support domain can be determined as,

ρ∗ =



w ( |r − rb |, h )mb

(22)

b

It has been assumed that the known particle position satisfies only divergence free condition. Ideally, for incompressible fluids, kernel value should always be equal to unity except nearby free surface and on dummy boundary particles. Instead of using density in source/sink term of the PPE, kernel value is used as an estimate of particle number density. Particle number density at known time level can be determined through kernel value (K∗ ) of each point.

K∗ =



w ( |r − rb |, h )vb

(23)

b

Taking divergence of momentum equation, and putting continuity equation of unknown time level into Eqn. 21, corrective PPE can be obtained as,

∇.

1 1 − K∗ ∗ ∇ P = ρ∗  t2

(24)

Corrective velocity can be determined by taking gradient of corrective pressure.

uc = −

t ∗ ∇P ρ∗

rnc +1 = uc t

(25) (26)

428

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Fig. 7. Comparison of kernel error in pure and VC-ISPH at t∗ = 1.50(a), t∗ = 3.00(b), t∗ = 5.70(c), t∗ = 6.45(d) of dam-break flow (Scenario-II).

Fig. 8. Comparison of Pressure Computed at right wall (Scenario-II).

rnc +1 is the corrective displacement required for maintaining density invariance. The corrective displacement is added with computed particle position at unknown time level as per the algorithm presented in Solution Procedure (2.5).

r

n+1

←r

n+1

+

rnc +1

(27)

Hydrodynamic variables are modified based on the obtained corrective displacement using Taylor’s series. The variables used in corrective algorithm (Eqn. 21) have no direct relation with specific hydrodynamic variables (velocity, pressure). The error as well as density depend only on particle position.

ghost particles do not change their position with time. Instead their property is modified by the information collected by interpolation nodes (Pahar and Dhar, 2016). Interpolation nodes are exact reflection of fixed ghost particles with respect to boundary particle layer. Direction of velocity is changed based on the nature of boundary condition (same for free-slip and opposite for no-slip). Divergence of position vector should be equal to 2 for full kernel support. Nearby free-surface, particles have truncated support domain. Thus ∇ .r have much smaller value on those particles (Lee et al., 2008). On the free surface particles, reference pressure condition can be directly imposed during solution of PPE. However, it may not conserve divergence-free condition on the free-surface particles. Bockmann et al. (2012) proposed an implicit consideration of free-surface particles. PPE can be discretized in the following manner for free surface particles.

1

ρ

∇ .(∇ P )|a =

 b

2

vb ((C + 1 )Pa − Pb )rab .∇ wab ρb |rab |2 + χ

(28)

where C = a constant. Eq. 28 combines the effect of divergencefree PPE and Dirichlet boundary condition.

 b

2CPa

vb rab .∇ wab =0 ρb |rab |2 + χ

(29)

2.7. Boundary condition No-slip and free-slip boundaries can be simulated by using either mirror (Morris et al., 1997) or dummy (Koshizuka et al., 1998) boundary particles. In this study, fixed ghost particles suggested by Federico et al. (2012) are used to simulate slip boundaries. Fixed

Implicit treatment of free surface helps to minimize error in velocity-divergence at the free-surface particles resulting from directly putting reference pressure as Dirichlet boundary condition. This implicit treatment is not suitable for other forms of Incompressible SPH .i.e., density invariant ISPH or their variants.

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

429

3. Results and discussion The developed VC-ISPH is applied for four different 2D cases of incompressible flow. The framework is validated for a lid driven cavity with Reynolds number 10 0 0. Second scenario considers a dambreak flow in dry bed. A water column passing vertically through a rectangular porous media is simulated in the third scenario. Fourth scenario assesses the effect of proposed algorithm on two cases of lock-exchange flows with two immiscible fluids. All scenarios are also simulated using pure divergence-free ISPH formulation (Xu et al., 2009). The change of volume is compared between divergence-free ISPH and proposed volume conservative algorithm. Mean value of kernel error (Ekernel ) is used as error criterion. Ekernel can be defined as,

Ekernel =

N  1   Ki  N

(31)

i=1

Fig. 9. Comparison of kernel and velocity divergence error in pure and VC-ISPH for dambreak-flow (Scenario-II).

where Ki = Error in kernel calculation (1 − Ki ) at particle i, Ki = kernel value for i-th particle and N = total number of fluid particles. Boundary and dummy particles are not included in error calculation. Mean error in velocity-divergence (Ediv ) is calculated by,

Ediv =

N  1   ∇ .ui  N

(32)

i=1

Both errors in kernel and velocity-divergence are compared for pure ISPH and proposed VC-ISPH models. Both ISPH models follow same aforementioned algorithm. Only corrective displacement at the end of each time step differs from each other [ISPH: (Eqns. 17, 18) and VC-ISPH: (Eqns. 26, 27)]. A value of 0.02 for β is assumed in pure ISPH model for particle shifting. Smoothing length h for each case is taken as 1.60x, x being the initial interparticle distance for a staggered distribution. At the no-slip boundary, 4 layers of fixed ghost particle is provided to maintain Neumann boundary condition. 3.1. Scenario-I

Fig. 10. Initial condition for flow through rectangular porous block (Scenario-III).

2.8. Time step Courant-Friedrichs-Lewy condition needs to be maintained for choosing a proper time step in ISPH.



t = min

σ1

δl umax

, σ2

δl2 ν



(30)

where umax = maximum velocity of the particles, σ 1 and σ 2 are problem-dependent constants( ≤ 0.5). The present study considers both as 0.1. In porous media, actual seepage velocity should be taken for computation instead of Darcy velocity.

Current model is applied to a benchmark numerical case of Lid driven cavity with Reynolds number 10 0 0. Both pure divergence free ISPH model and current model are utilized for the test (Fig. 1). Top lid of the square box of side L moves with a constant velocity Ulid . The size of the cavity (L) and lid velocity Ulid is set as unity. Viscosity is calculated accordingly. Initially, the fluid velocity inside the domain is set as zero. Top surface acts as a free-slip boundary. Other sides of the box are taken as no-slip boundary condition. The resolution for the simulation is taken as L/125. The steady velocity profile (U m/s) obtained from pure ISPH and VC-ISPH is shown in Fig. 2. It is evident that the variation in velocity field is negligible. The velocity profile at the centre-line is compared with benchmark data of Ghia et al. (1982) for validation in Fig. 3. The present model and pure divergence free ISPH produce identical result for velocity at centre-line. This validates the current numerical framework. The absolute errors in kernel (Ekernel ) and velocity divergence (Ediv ) obtained over time are compared with pure divergence free ISPH model in Fig. 4. Pure divergence free ISPH produces small error in kernel computation inside the computational domain. However, current framework decreases that error considerably. Improvement in velocity divergence can also be observed in the current model. The pressure profiles across horizontal and vertical centrelines of the cavity are compared with Botella and Peyret (1998) in Fig. 5. Their solution was obtained on a 513 × 513 Chebyshev-Gauss-Lobatto

430

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Fig. 11. Comparison of kernel error in pure and VC-ISPH for ∇ .r = 1.60 (a,b,c) and ∇ .r = 1.70 (d,e,f) at 0.10s, 0.50s, 5.00s (Scenario-III).

Fig. 12. Comparison of kernel and velocity divergence error in pure ISPH (a,c) and VC-ISPH (b,d) (Scenario-III).

grid. In their collocation method, velocity in both direction is calculated with polynomial of degree 128. The computed pressure is defined on a polynomial of two degree less than that of velocity. Both numerically obtained pressure profile show almost identical nature. However, minor deviations can be observed near cav-

ity wall. Similar results are also reported by Khorasanizade and Sousa (2014). Khorasanizade and Sousa used a resolution of L/130 for comparing their ISPH solution with Botella and Peyret’s (1998) result.

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

431

Fig. 13. Initial Condition for lock-exchange flow (a) without and (b) with free-surface at top (Scenario-IV).

Fig. 14. Comparison of fluid front positions for pure and VC-ISPH (Scenario-IVa).

Fig. 15. Position of heavier and lighter fluids for pure ISPH (a, c, e, g) and VC-ISPH (b, d, f, h) at different time instances (Scenario-IVa).

3.2. Scenario-II A dambreak flow in dry bed (Adami et al., 2012) is modelled with pure divergence-free and VC-ISPH. Initially water column of width 2H m and height H m is placed in a box of length 5.366H m (Fig. 6). Pressure is measured at the right wall of the box at a height 0.19H from the bottom. Initially, hydrostatic pressure distribution is assumed. A vertical barrier of negligible thickness is placed at 2H m distance from the left wall. The barrier is lifted instantaneously and water is allowed to move. Value of H is set as 1 m for the current simulations. Initial interparticle spacing for the dambreak flow is taken as 0.02m for an equidistant particle distribution. Total 7508 fluid particles including dummy and boundary particles have been considered to simulate flow for both pure divergence-free and current model. Free surface criteria is set as ∇ .r = 1.60 (Bockmann et al., 2012). Variation of kernel error in fluid particles with time are shown in Fig. 7. Time is nondimensionalized by t(g/H)1/2 . The pressure computed by the simulations at the right wall is compared with experimental data in Fig. 8. Pressure variation across time obtained from both simulations are very similar. Numerical variation in computed pressure is com-

paratively less in case of VC-ISPH. As pure ISPH does not use any packing algorithm to check particle number density, free surface may get identified at the interior along with actual free-surface. VC-ISPH maintains proper particle arrangement to avoid such situations. The highest peaks obtained in both numerical models have a time lag compared to experimental data. The delay in pressure occurs due to single phase flow simulation, where air cushion effect is not taken into consideration. Similar discrepancy in pressure computation have been previously reported in ISPH (Gui et al., 2015) and WCSPH (Adami et al., 2012; Marrone et al., 2011). Marrone et al. (2011)’s simulation can capture the highest pressure peak. However, it fails to produce smaller pressure peaks. Moreover, their model involves introduction of artificial diffusive terms into continuity equation to remove spurious pressure noises. Pure ISPH results show that the error in kernel has a cascading effect on fluid particles over time. However, deviation from unit kernel value is not observed for VC-ISPH. The mean error in kernel (Ekernel ) and velocity divergence (Ediv ) for pure divergence free ISPH and the developed VC-ISPH are plotted against time in Fig. 9. Error in kernel is reduced to a minimal level with the proposed VC-ISPH. The proposed framework also preserves the divergencefree nature of the flow. Particle number density in case of pure di-

432

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Fig. 16. Kernel Errors in pure (a, c, e, g) and VC-ISPH (b, d, f, h) (Scenario-IVa).

Fig. 17. Comparison of kernel and velocity divergence error in pure ISPH and VCISPH (Scenario-IVa).

vergence free ISPH varies depending on the value of free-surface criterion. Same case of dambreak flow was tested by Gui et al. (2015) for validating their improved ISPH model. The error in absolute density for divergence free ISPH was in order of 0.1 − 0.25. However, for current simulation maximum mean error is in order of 0.12. Similar variation of divergence error can also be found in the literature. Generally the present model performs better compared to pure ISPH for conserving divergence-free criteria. 3.3. Scenario-III Scenario-III deals with a water column passing through a porous block(Fig. 10). A numerical model similar to Akbari (2014) is solved using proposed method and compared with di-

vergence free formulation. For a continuum approach, position of porous media needs to incorporated into particle volume depending on its position. Inside porous media, fluid particles move further away from each other compared to pure fluid region as volume of solid matrix need to be accumulated. Thus inside porous media, the particles represent fluid as well as solid volume. Akbari (2014) changed apparent density of the fluid in porous media to incorporate volume change. However, pure divergence free ISPH considers constant density throughout the time period. Thus, volume of the particle is increased directly by introducing a factor ( η1 ) of porosity into it. Pressure gradient is derived using linear momentum conservation (Eq. 7). Initially, a 0.20 m water column is assumed to be kept at a height 0.60 m from datum (Fig. 10). A 0.20 m rectangular block of porous media is kept just below the water column. The water column and the porous block is separated by a horizontal gate. The gate is moved suddenly to allow the water column to pass through the porous block. Effective porosity of the block is taken as 0.5. Intrinsic permeability in both direction is assumed as 0.01m2 . The error in kernel value computed by both models at different time intervals are shown in the Fig. 11. The free surface criteria is set to be 1.60 (Bockmann et al., 2012) and 1.70 for divergence of position vector (∇ .r). Initial interparticle distance is taken as 0.0045m for staggered particle distribution.Total 5761 particles are used to simulate the scenario. There should be no change in water volume after passing through the porous media. However, pure divergence-free ISPH fails to preserve volume conservation property. There is also sudden increase of particle number density nearby interface. In VC-ISPH, particle number density remains very close to unity. (Akbari, 2014) solved the problem using density invariance scheme. Thus, particle number density was kept constant at all time levels for their model. The error in kernel calculation and velocity divergence is plotted with time in Fig. 12 for both cases. Spurious fluctuations in velocity divergence reduces in case of volume conservative ISPH compared to pure ISPH. Interior particles near interface may get wrongly identified as free-surface particles for pure divergence-free ISPH. The particles change their

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

433

Fig. 18. Comparison of density profile in pure and VC-ISPH at t∗ = 1.30(a), 3.30(b), 5.40(c), 6.30(d) of lock exchange flow (Scenario-IVb).

arrangement frequently while passing through porous and nonporous phase. In porous phase, they need to represent volume of solid matrix along with actual fluid volume. After passing through porous to non-porous phase, their velocity suddenly increases due to absence of any resistive force. The particles need to have proper compaction after transition to free-surface particle to interior particle (Fig. 11c). For pure divergence-free ISPH, the particles stop moving close to each other as soon as free surface threshold is reached. Thus, two different criteria for free-surface identification produce different results at t = 10.00s. In case of VC-ISPH, the particle compaction checking algorithm ensures proper particle distribution. For pure ISPH, as particle distribution is very close to free-surface criterion, sudden movement can trigger identification of free-surface particles. Consequently, the error in velocity divergence for pure ISPH shows rapid spikes. VC-ISPH reduces the error in velocity divergence by maintaining proper particle distribution and preventing erroneous identification of free-surface. It is evident that errors in velocity divergence and kernel calculation in VC-ISPH are not varying significantly with changing free-surface criterion (Fig. 12). 3.4. Scenario-IV Lock-exchange flow for two immiscible fluids is simulated in scenario-IV. To maintain continuity of pressure across interface, pressure gradient needs to be discretized in a different manner compared to dambreak flow. Inter-particle averaged pressure (Eq. 9) is used instead of direct pressure. Free surface criteria is taken as ∇ .r = 1.60 (Bockmann et al., 2012). Effective porosity (η) is taken as unity for multiphase flow. Two cases of lock-exchange flow have been solved using divergence free ISPH and volume conservative ISPH. Lowe et al. (2005) performed an experiment of lock-exchange flow with immiscible fluids in a closed rectangular channel. All sides of the rectangular box are taken as no-slip boundary conditions. Two fluids of different density and equal volume are kept

in the box of length 1.82 m and height 0.2 m (Fig. 13a). The fluid domain is represented by a set of staggered nodes/particles with initial distance 0.0075m. They are separated by a gate of negligible thickness. The ratio (γ ) of density of immiscible fluids is taken as 0.681. Obtained fluid fronts are compared with available experimental data. Front positionand time is non-dimentionalized using following criteria: t ∗ → t (g(1 − γ )/H ) and x∗ → Hx . Both ISPH and VC-ISPH fluid fronts show proper agreement with experimentally obtained lighter fluid front (Fig. 14). Simulated heavier fluid front maintain a constant distance from the experimental data. Instantaneous lifting of gate is considered in the numerical simulations. In experiment, the barrier is lifted with a finite velocity. Thus, heavier fluid gets a headstart compared to the lighter fluid. However, slopes of the obtained fronts have almost identical nature as the experimental front. Similar discrepancy in numerical results with weakly compressible SPH had been observed by Chen et al. (2015). Positions of heavier and lighter fluids for different time instances are shown in Fig. 15. Results obtained from both numerical models are almost identical. Errors in kernel calculation for different time instances are plotted in Fig. 16 for ISPH with and without volume conservation correction. Simultaneous compression and expansion are observed in fluid domain, although total volume of fluids remains unchanged. Errors in kernel calculation and velocity-divergence are compared for both methods (Fig. 17). Using the proposed correction error in particle number density can be minimized to a suitable value. Error in velocity divergence shows a similar trend for both methods. In the second case of lock-exchange flow, top boundary is removed to allow free-surface condition(Fig. 13b). Other parameters remain same as in the first case. Initially, hydrostatic pressure distribution is considered. The gate is removed instantaneously and fluids are allowed to interact with each other. 16288 and 16463 particles are used to simulate lock-exchange flow for closed and free-surface flow cases respectively. Evolution of lock-exchange flow for t ∗ = 1.30 − 6.30 is shown in Fig. 18.

434

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Fig. 19. Comparison of Kernel error in pure and VC-ISPH at t∗ = 1.30(a), 3.30(b), 5.40(c), 6.30(d) of lock exchange flow (Scenario-IVb) . Table 1 Comparison of Computation Time for ISPH and VC-ISPH. Type

Lid-driven cavity Dambreak flow Flow through porous block Lock-exchange flow in closed box Lock-exchange flow with free surface

No of particles

17966 7508 5761 16288 16463

Average Computational Time(s) [for 100 time steps] ISPH

VC-ISPH

180.09 64.62 60.24 144.42 181.06

213.43 75.73 70.59 168.24 215.03

The density profile obtained form purely divergence-free and VC-ISPH is almost same. However, free surface of the flow system in pure ISPH is slightly higher than the results obtained by volume conservative ISPH. Kernel error in both ISPH models are plotted in Fig. 19. Free surface computed by ISPH is dispersed over space, whereas volume conservative algorithm maintains a compact distribution. Errors in kernel and velocity divergence for both ISPH models are in depicted in Fig. 20. VC-ISPH shows an overall improvement of divergence-free velocity criterion over pure ISPH. Free surface criterion plays an important role to maintain particle compactness for pure divergence-free ISPH. It determines the kernel threshold of particles for locating free-surface particles. Due to high pressure gradient at free-surface particles, particles present in the support domain move closer to free-surface particle. Thus, a very high threshold criterion can partly reduce numerical expansion (Fig. 12). In the scenario-III, error in particle number density is reduced to 0.12 from 0.16 for ∇ .r = 1.60 and 1.70 respectively. However, PPE is not directly solved on free surface particles. Instead, a combined form of Dirichlet boundary condition and PPE is solved to implicitly introduce zero pressure on

Percent increase in Computational Time

18.51 17.19 17.18 16.49 18.76

them. This methodology partly reduces error in solenoidal velocity condition at free surface compared to directly putting a reference pressure. More numbers of free-surface particles are identified for a relaxed criterion (∇ .r), which in turn reduces kernel error. However, much variation in velocity-divergence is observed for ISPH. In case of numerical compression, i.e., ∇ .r > 2 (for 2D) the kernel value is unrestricted. Particle density fluctuates widely in pure ISPH based on the nature of the flow and free-surface criteria. Proposed algorithm ensures that kernel value remains close to unity. In the Figs. 4a, 9a, 12a, 17a, 20a maximum error in volume conservative algorithm is within the range of 0.002 − 0.004. For low free-surface criterion, particles near free surface have lesser kernel value because of truncated support domain. The points having lesser kernel value are always included in error calculation. In case of flow passing through porous block, the error has a considerably large value compared to other scenarios. However, particle compaction in VC-ISPH helps to curtail that error. In an immiscible lock-exchange flow, free surface variation over time is relatively low leading to less error in kernel at the beginning. However, the error accumulation is reduced significantly by the proposed corrective measure. Divergence free formulation

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

435

to maintain solenoidal velocity condition and incompressibility of the fluid. The effectiveness of the algorithm is tested for four different 2D numerical examples. With the change in flow domain (ScenarioIII) or lock-exchange (Scenario-IVb) numerical compression is developed in pure ISPH model. Simultaneous compression/expansion in fluid volume is visible in all scenarios. Ekernel and Ediv are used to quantify the numerical errors in density and divergence-free velocity condition. Improvement in Ekernel value (0.12 − 0.16) for VCISPH (compared to ISPH) is observed in case of dam-break and lock-exchange flow. Most visible change (0.3) in volume is found in Scenario-III. Similar improvements in flow inside closed domain (Scenarios-I,IVa) can also be seen with the proposed algorithm. Moreover, error in velocity-divergence (Ediv ) reduces considerably for all the scenarios. The proposed correction can be extended in lateral direction for three-dimensional application. Acknowledgements The authors are grateful to two anonymous reviewers and Associate Editor Prof. Rui M.L. Ferreira for their valuable suggestions. Fig. 20. Comparison of kernel and velocity divergence error in pure and VC-ISPH (Scenario-IVb).

and density invariance are not two mutually exclusive events. Error in divergence free velocity (∇ .u) reduces with decreasing error in particle number density. Solving two simultaneous PPEs in a single time step increases computational cost of the proposed numerical model. All simulations are performed in a computer containing Intel-i5 3.10GHz processor with RAM 4GB. Average CPU runtime for the simulations are given in the Table 1. Computational time mainly depends on number of the particles used in the simulation along with type of particles. Proper balancing of pressure along an interface of multifluid flow may consume high computational time. In a Lagrangian framework, certain time is required for searching particles present in support domain. A linked list algorithm is used for all the simulations used in this study. For a second-order time marching scheme (Cummins and Rudman, 1999; Shadloo et al., 2012), cost of solving PPE may increase upto 19% of total computation time. Computational time may be further reduced by using a first order time stepping and Truly Incompressible SPH (Lee et al., 2008). For divergence-free ISPH simulation, free-surface may get identified at the interior portion leading to velocity fluctuation. Consequently, total computational time may increase compared to VC-ISPH. The proposed scheme does not require any parametric study (Asai et al., 2012; Gui et al., 2015) of weight factor involved in source/sink term or iterative solution of multiple PPEs (Hu and Adams, 2007). Other ISPH methods developed based on density invariance model do not necessarily minimize velocity-divergence error. PPE obtained from density-invariance source/sink term may induce high numerical noise (Xu et al., 2009) in pressure. XSPH may be required to get ordered particle distribution (Keat et al., 2015) for these cases. Proposed correction does not require any stabilization algorithm in form of XSPH. 4. Conclusions A corrective VC-ISPH is proposed for reducing the numerical compressibility induced by traditional divergence-free ISPH algorithm. The framework implements traditional divergence-free algorithm in the first step. In the second step, displacement correction is introduced based on particle density error while maintaining the divergence free velocity field. The algorithm solves two PPEs

Appendix Inside porous media, two phases are present: fluid phase (α ) and solid phase (β ) (Gray, 1975; Gray and Lee, 1977). Thus for a representative volume dV,

dV = dVα + dVβ

(33)

A macroscopic volume averaged quantity < φ > α can be written as,

< φ >α =



1 dV

dV

φγα dv

(34)

Intrinsic volume averaged quantity < φ >α α over α -phase can be written as,

1 < φ >αα = dVα



dV

φγα dv

(35)

γ is a binary function, whose value is dependent on flow phase. dv is microscopic elemental volume. For α -phase, γ = 1. For other phases, γ = 0. Thus, < φ > α can be written as, < φ >α =

dVα < φ >αα dV

(36)

< φ >α = α < φ >αα

(37)

where α denotes the ratio of dVα and dV. Mass-averaged quantity α φ can be written as,



ρφγα dv ργα dv

α

1 dVα dV  1 dVα dV

(38)

α

< ρφ >αα < ρ >αα

(39)

φ = φ =

For incompressible fluid, intrinsic and mass averaging will be α α α α α identical ,i.e. φ =< φ >α α . Thus, < ρφ >α = φ < ρ >α =< φ >α ρα . ρ α is the density of the α -phase. According to Whitaker (2013), Gray (1975) and Gray and Lee (1977), following expressions are used for volume averaging of microscopic equations.

< ∇φ >α = ∇ ( α < φ >αα ) +

1 dVα



sf

ns f φ dS

(40)

where ns f = unit vector normal to sf , s f = interfacial surface. Similarly for a vector A,

< ∇ .A >α = ∇ .( α < A >αα ) +

1 dVα



sf

ns f .AdS

(41)

436

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437

Derivation of Governing equations: Continuity equation for microscopic velocity u can be written as,

∇ .u = 0

(42)

Integrating over the domain, the continuity equation for macroscopic velocity < u > α is written.

< ∇ .u >α = 0

(43)

∇ .( α < u >αα ) = 0

(44)

Momentum equation for microscopic quantity can be written as,

∂ ( ρ ui ) ∂ ( ρ ui u j ) ∂ P ∂τi j + =− + + ρg ∂t ∂xj ∂ xi ∂ x j 

∂u

∂u

where τi j = μ ∂ x i + ∂ x j j i

(45)

. For time derivative of any quantity,

1 ∂A ∂ ( α < A >αα ) < >α = − ∂t ∂t dVα

1 ∂ ( α < P >αα ) ∂ ( α < ui >αα ) ∂ ( α < ui >αα < u j >αα ) + =− ∂t ∂xj ρα ∂ xi R α 1 (∂ α < τi j >α ) R 1 ∂ τi j + + α < g >αα + + (50) ρα ∂xj ρα ρα ∂ x j τR

α α where, ρiαj = ( α (< ui >α α < u j >α − < ui u j >α )).

sf

ns f .(us f A )dS

(46)

∂ ( α < ρ ui >αα ) ∂ ( α < ρ ui u j >αα ) ∂ ( α < P >αα ) + =− ∂t ∂xj ∂ xi  (∂ α < τi j >αα ) 1 + + α < ρ g >αα + ns (τi j − P δi j )dS ∂xj dVα s f f (47)

δi j = Kronecker delta function. The surface integral term  1 n (τi j − P δi j )dS is denoted as the resistive force R indVα s f s f side porous domain. For an incompressible fluid, < ρ >α α = ρα (Eqn. 39). Thus the equation can be written as,

1 ∂ ( α < P >αα ) ∂ ( α < ui >αα ) ∂ ( α < ui u j >αα ) + =− ∂t ∂xj ρα ∂ xi α 1 (∂ α < τi j >α ) R + + α < g >αα + (48) ρα ∂xj ρα For pure fluid flow, Navier-Stokes equation is required to be filtered for Large Eddy Simulation (LES). Any quantity φ can be written as summation of filtered quantity φ and fluctuating term φ  . Following Cheng et al. (2012), a top-hat filter is adopted in the present study, such that,

(49)

Use of this box filter denotes that the spatially averaged quantity will be same as volume-averaged quantity. It is also to be noted that LES is only valid for pores occupied by fluid (not the solid volume). Thus, the equation obtained after volume-averaging is identical to filtered Navier-Stokes Equation with an extra R term. The convective term (< ui u j >α α ) is difficult to model similar to filtered Navier-Stokes equation as it involves knowledge of actual velocity throughout the pores/unfiltered velocity in the pure fluid region.

It

can

be seen that τiRj = Li j + Ci j + Ri j . Lij is the Leonard tensor α α α α ( α (<< ui >α α < u j >α >α − < ui >α < u j >α )). Cij term is for cross-stress tensor, which describes the interaction between large  α  and small scales( α (< ui >α α u j + < u j >α ui )). Rij is sub-grid/subparticle tensor for small-scale interactions ( α (< ui uj >α α )). For the particular filter used here, Leonard tensor becomes zero.

∂ τiRj R 1 ρα + ρα ∂ x j can be written as Resistive force (FR ) term for differ-

ent flow-phases. The surface integral term ( ρRα ) can be physically denoted as the resistive force applied on the fluid particles due to presence of porous structures. Effective porosity (η) is defined as ratio of volume of voids (Vv ) to gross volume (V). Thus η becomes essentially equal to α . For flow through porous media, Darcy velocity (ud ) can be defined as η < u >α α.





us f = velocity of the solid phase. For an immobile solid skeleton, us f = 0. Integrating over domain, the momentum equation for macroscopic scale can be written as,

φ =< φ >αα

Thus the equation can be written as,

α < τij >αα = α μ  α < τij >αα = μ

∂ < ui >αα ∂ < u j >αα + ∂xj ∂ xi

∂ ud i ∂ ud j + ∂xj ∂ xi



(51)

 (52)

The left side can be written as material derivative,

Du d i = Dt

∂ ud i ∂ (ud i < u j >αα ) + ∂t ∂xj

(53)

ρ α is density of fluid at α phase (ρ ).Thus the continuity and momentum equations can be written as,

∇ .ud = 0

(54a)

Du d η = − ∇ P + ηg + ν∇ 2 ud + FR Dt ρ

(54b)

FR term contains a surface-integral term of pressure and viscous drag and a divergence of residual stress tensor. The residual stress tensor is modelled irrespective of the region. References Adami, S., Hu, X., Adams, N., 2012. A generalized wall boundary condition for smoothed particle hydrodynamics. J. Comput. Phys. 231 (21), 7057–7075. http:// dx.doi.org/10.1016/j.jcp.2012.05.005. Akbari, H., 2014. Modified moving particle method for modeling wave interaction with multi layered porous structures. Coastal Eng. 89, 1–19. http://dx.doi.org/ 10.1016/j.coastaleng.2014.03.004. Aly, A.M., Asai, M., 2015. Three-Dimensional Incompressible Smoothed Particle Hydrodynamics for Simulating Fluid Flows Through Porous Structures. Transp. Porous Media 110 (3), 483–502. Asai, M., Aly, A.M., Sonoda, Y., Sakai, Y., 2012. A Stabilized Incompressible SPH Method by Relaxing the Density Invariance Condition. J. Appl. Math. 1–24. http://dx.doi.org/10.1155/2012/139583. Bockmann, A., Shipilova, O., Skeie, G., 2012. Incompressible SPH for free surface flows. Comput. Fluids 67, 138–151. Bonet, J., Lok, T.-S., 1999. Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations. Comput. Meth. Appl. Mech. Eng. 180 (1–2), 97–115. http://dx.doi.org/10.1016/S0 045-7825(99)0 0 051-1. Botella, O., Peyret, R., 1998. Benchmark spectral results on the lid-driven cavity flow. Comput. Fluids 27 (4), 421–433. Chen, Z., Zong, Z., Liu, M., Zou, L., Li, H., Shu, C., 2015. An SPH model for multiphase flows with complex interfaces and large density differences. J. Comput. Phys. 283, 169–188. http://dx.doi.org/10.1016/j.jcp.2014.11.037. Cheng, K., Chun, S., Hweng, H., Ren, T., 2012. Three-dimensional numerical modeling of the interaction of dam-break waves and porous media. Adv. Water Res. 47, 14–30.

G. Pahar, A. Dhar / Advances in Water Resources 96 (2016) 423–437 Cummins, S.J., Rudman, M., 1999. An SPH projection method. J. Comput. Phys. 152 (2), 584–607. http://dx.doi.org/10.1006/jcph.1999.6246. Federico, I., Marrone, S., Colagrossi, A., Aristodemo, F., Antuono, M., 2012. Simulating 2D open-channel flows through an SPH model. Eur. J. Mech. - B/Fluids 34, 35– 46. http://dx.doi.org/10.1016/j.euromechflu.2012.02.002. Gent, M. R. A. V., 1995. Wave Interaction with Permeable Coastal Structures. Ghia, U., Ghia, K., Shin, C., 1982. High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method. J. Comput. Phys. 48 (3), 387–411. Ghimire, B., 2009. Hydraulic Analysis of Free-surfce flows into Highly Permeable Porous Media and its Applications. Kyoto University, The Department of Urban Management. Gingold, R., Monaghan, J., 1977. Smoothed particle hydrodynamics - theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181, 375–389. Gomez-Gesteira, M., Rogers, B.D., Dalrymple, R.a., Crespo, A.J., 2010. State-of-the-art of classical SPH for free-surface flows. J. Hydraul. Res. 48 (sup1), 6–27. http:// dx.doi.org/10.1080/00221686.2010.9641242. Gray, W.G., 1975. A derivation of the equations for multi-phase transport. Chem. Eng. Sci. 30 (2), 229–233. http://dx.doi.org/10.1016/0 0 09-2509(75)80 010-8. Gray, W.G., Lee, P.C.Y., 1977. On the theorems for local volume averaging of multiphase systems. Int. J. Multiphase Flow 3 (4), 333–340. http://dx.doi.org/10.1016/ 0301-9322(77)90013-1. Gui, Q., Dong, P., Shao, S., 2015. Numerical study of PPE source term errors in the incompressible SPH models. Int. J. Numer. Meth. Fluids 77, 358–379. Gui, Q., Shao, S., Dong, P., 2014. Wave Impact Simulations by an Improved ISPH Model. J. Waterw. Port Coastal Ocean Eng. 140, 1–14. Hosseini, S., Manzari, M., Hannani, S., 2007. A fully explicit three-step SPH algorithm for simulation of non-Newtonian fluid flow. Int. J. Numer. Meth. Heat Fluid Flow 17 (7), 715–735. http://dx.doi.org/10.1108/09615530710777976. Hu, X., Adams, N., 2007. An incompressible multi-phase SPH method. J. Comput. Phys. 227 (1), 264–278. http://dx.doi.org/10.1016/j.jcp.2007.07.013. Hu, X., Adams, N., 2009. A constant-density approach for incompressible multi-phase sph. J. Comput. Phys. 228 (6), 2082–2091. Hu, X.Y., Adams, N.A., 2006. A multi-phase SPH method for macroscopic and mesoscopic flows. J. Comput. Phys. 213 (2), 844–861. http://dx.doi.org/10.1016/j.jcp. 20 05.09.0 01. Keat, S., Cheng, N.-s., Xie, Y., Shao, S., 2015. Incompressible SPH simulation of open channel flow over smooth bed. J. Hydro-Environ. Res. 9, 340–353. Khayyer, A., Gotoh, H., 2011. Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230 (8), 3093–3118. Khorasanizade, S., Sousa, J.M.M., 2014. A detailed study of lid-driven cavity flow at moderate Reynolds numbers using Incompressible SPH. Int. J. Numer. Meth. Fluids 76 (10), 653–668. Koh, C., Luo, M., Gao, M., Bai, W., 2013. Modelling of liquid sloshing with constrained floating baffle. Comput. Struct. 122, 270–279. Koshizuka, S., Nobe, A., Oka, Y., 1998. Numerical Analysis of Breaking Waves Using the Moving Particle Semi-Implicit Method. Int. j. Numer. Meth. fluids 769 (April 1997), 751–769.

437

Lee, E., Moulinec, C., Xu, R., Violeau, D., Laurence, D., Stansby, P., 2008. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. J. Comput. Phys. 227 (18), 8417–8436. http://dx.doi.org/10. 1016/j.jcp.20 08.06.0 05. Lowe, R.J., Rottman, J.W., Linden, P.F., 2005. The non-Boussinesq lock-exchange problem. Part 1. Theory and experiments. J. Fluid Mech. 537, 101–124. http://dx.doi. org/10.1017/S0 0221120 050 05069. Lucy, L.B., 1977. A numerical approach to the testing of the fission hypothesis. Astron. J. 82 (12), 1013–1024. Marrone, S., Antuono, M., Colagrossi, A., Colicchio, G., Le Touzé, D., Graziani, G., 2011. δ -SPH model for simulating violent impact flows. Comput. Meth. Appl. Mech. Eng. 200 (13–16), 1526–1542. Monaghan, J., 1994. Simulating Free Surface Flows with SPH. J. Comput. Phys. 110, 399–406. Morris, J.P., Fox, P.J., Zhu, Y., 1997. Modeling Low Reynolds Number Incompressible Flows Using SPH. J. Comput. Phys. 136 (1), 214–226. http://dx.doi.org/10.1006/ jcph.1997.5776. Pahar, G., Dhar, A., 2014. A dry zone-wet zone based modeling of surface water and groundwater interaction for generalized ground profile. J. Hydrol. 519, Part B, 2215–2223. http://dx.doi.org/10.1016/j.jhydrol.2014.09.088. Pahar, G., Dhar, A., 2016. Modeling free-surface flow in porous media with modified incompressible SPH. Eng. Anal. Boundary Elements 68, 75–85, http://dx.doi.org/10.1016/j.enganabound.2016.04.001. Pan, X.-j., Zhang, H.-x., Sun, X.-y., 2012. Numerical simulation of sloshing with large deforming free surface by MPS-LES method. China Ocean Eng. 26 (4), 653–668. http://dx.doi.org/10.1007/s13344- 012- 0049- 6. Liu, P.L.F., Lin, P.Z., Chang, K.A., Sakakiyama, T., 1999. Numerical modeling of wave interaction with porous structures. J. Waterw. Port Coastal Ocean Eng. 125, 322–330. Shadloo, M.S., Zainali, a., Yildiz, M., 2012. Simulation of single mode Rayleigh-Taylor instability by SPH method. Comput. Mech. 51 (5), 699–715. http://dx.doi.org/10. 10 07/s0 0466- 012- 0746- 2. Shao, S., 2010. Incompressible SPH flow model for wave interactions with porous media. Coastal Eng. 57 (3), 304–316. http://dx.doi.org/10.1016/j.coastaleng.2009. 10.012. Shao, S., Gotoh, H., 2005. Turbulence particle models for tracking free surfaces. J. Hydraul. Res. 43 (3), 276–289. Shao, S., Lo, E.Y., 2003. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv. Water Res. 26 (7), 787–800. http://dx.doi.org/10.1016/S0309-1708(03)0 0 030-7. Smagorinsky, J., 1963. General Circulation Experiments with the Primitive Equations. Mon. Weather Rev. 91, 99. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389–396. http://dx.doi.org/10.1007/BF02123482. Whitaker, S., 2013. Diffusion and Dispersion in Porous Media. J. Chem. Inf. Model. 53 (9), 1689–1699. http://dx.doi.org/10.1017/CBO9781107415324.004. Xu, R., Stansby, P., Laurence, D., 2009. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. J. Comput. Phys. 228 (18), 6703–6725. http://dx.doi.org/10.1016/j.jcp.2009.05.032.