Incompressible–compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH)

Accepted Manuscript Incompressible–compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH)

S.J. Lind, P.K. Stansby, B.D. Rogers

PII: DOI: Reference:

S0021-9991(15)00820-7 http://dx.doi.org/10.1016/j.jcp.2015.12.005 YJCPH 6283

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

16 March 2015 26 October 2015 7 December 2015

Please cite this article in press as: S.J. Lind et al., Incompressible–compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH), J. Comput. Phys. (2015), http://dx.doi.org/10.1016/j.jcp.2015.12.005

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.

2

Incompressible-compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH)

3

S. J. Lind, P. K. Stansby, B. D. Rogers

1

4

1 School

5

Abstract

of Mechanical, Aerospace and Civil Engineering, The University of Manchester, Manchester, UK, M13 9PL

A new two-phase incompressible-compressible Smoothed Particle Hydrodynamics (SPH) method has been developed where the interface is discontinuous in density. This is applied to water-air problems with a large density difference. The incompressible phase requires surface pressure from the compressible phase and the compressible phase requires surface velocity from the incompressible phase. Compressible SPH is used for the air phase (with the isothermal stiffened ideal gas equation of state for low Mach numbers) and divergence-free (projection based) incompressible SPH is used for the water phase, with the addition of Fickian shifting to produce sufficiently homogeneous particle distributions to enable stable, accurate, converged solutions without noise in the pressure field. Shifting is a purely numerical particle regularisation device. The interface remains a true material discontinuity at a high density ratio with continuous pressure and velocity at the interface. This approach with the physics of compressibility and incompressibility represented is novel within SPH and is validated against semi-analytical results for a two-phase elongating and oscillating water drop, analytical results for low amplitude inviscid standing waves, the Kelvin-Helmholtz instability, and a dam break problem with high interface distortion and impact on a vertical wall where experimental and other numerical results are available.

6

1. Introduction

7

Multi-phase flows are common, in fact quite general, in environmental and industrial processes. Broadly these may be modelled as continuous problems where phases are mixed (e.g. oil-water homogenisation [36], sediment transport [18]) or interface problems where phases are distinct and interact at the interface (e.g. gas-assisted injection moulding [21], liquid jet breakup [40]). In some cases flows start as interface problems but as mixing occurs at the interface they become effectively continuous, at least locally. Air entrainment, perhaps due to wave breaking, is an obvious example. We consider here two-phase interface problems where the interface remains distinct and the density difference is high, e.g. air and water, and where one phase may be considered incompressible. The interface is transient and may become highly distorted and interconnected. Such problems have been tackled with mesh-based methods using periodic (or adaptive) re-meshing or additional phase tracking functions [40]. However, these approaches can be time-consuming to implement and prone to errors in surface representation [50] or mass conservation [34].

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

Lagrangian particle methods are ideally suited to multi-phase flows as they are naturally surface or interface following and conservative: here we employ the smoothed particle hydrodynamics (SPH) method. SPH has been previously applied to such problems either in (weakly) compressible form or incompressible form, but not a combination of the two. Here we impose the incompressible divergence-free condition for one phase while the other phase is compressible. This requires a new interface treatment which is the subject of this paper. The air phase of interest here is weakly compressible, but, in SPH, this term generally refers to the form with a stiff equation of state to represent incompressible flow (originally due to Cole [8], but often attributed to Tait [46]). Here we refer to compressible SPH or CSPH to distinguish the method as the physics of compressibility is ∗ Corresponding

author: S. J. Lind, [email protected]

Preprint submitted to Elsevier

December 21, 2015

28

represented through the isothermal gas equation appropriate for low Mach numbers.

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

The majority of multi-phase SPH simulations are based on the weakly compressible SPH (WCSPH) formulation, which originates from Lucy [25] and Gingold & Monaghan [13]. Monaghan & Kocharyan [28] presented one of the first extensions of SPH to multi-phase flows, but applied to astrophysical dusty gas simulations. The work of Colagrossi & Landrini [7] was one of the first SPH studies to consider water-air flows at realistic density ratios of 1 to 1000. They used WCSPH with gradients recast in terms of volume (not density) to avoid issues around the density discontinuity at the water-air interface. The method also required an empirical adhesion term that physically resembled surface tension and maintained interface integrity. Results for bubble and dam-break problems showed good agreement with other numerical methods and experimental results. This work was extended by Grenier et al. [15] who developed a multifluid SPH method derived from Lagrangian variational principles. Surface tension is included explicitly through the Continuum Surface Force approach, in addition to an empirical repulsive term in the pressure gradient that minimises the fragmentation of the interface (this term bears semblance to the cohesion force introduced by [7]). Both [7] and [15] demonstrated good results for their chosen test cases, but both methods are based on WCSPH and are restricted to using non-physical speeds of sound. In both cases, the air phase requires a speed of sound typically 10 times larger than the water phase, when, in reality, the opposite should occur. Recent work by Monaghan & Rafiee [31] presented a similar approach to modelling multi-phase flows with large density ratios (up to 1000:1). The method is attractive in its simplicity and, in a similar manner to [15], includes a repulsion term that acts between fluids in different phases. However, a larger sound speed is still required in the less-dense fluid.

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

The inconsistency of non-physical sound speeds is avoided if the flow is such that incompressibility can be safely assumed in both phases. In truly incompressible SPH (ISPH) formulations incompressibility is enforced through a projection method, not approximated through large values of sound speed as in WCSPH. First proposed by Cummins & Rudman [11], ISPH solves a pressure Poisson equation (PPE) for the fluid pressure, and, provided particle distributions are near uniform, pressure predictions are smooth, accurate and noise-free [49, 24, 45]. In contrast WCSPH pressure predictions are known to be inaccurate and oscillatory [20], unless appropriate smoothing or diffusive terms are applied [27]. As with WCSPH, ISPH has been readily extended to multi-phase flows by a number of authors, notably Hu & Adams [16, 17] who utilise a form of spatial derivative that enables discontinuities across phases to be handled naturally. Incompressibility is enforced in both phases through a projection method and the results presented for droplet deformation and the Rayleigh-Taylor instability showed good agreement with reference solutions. However, the density iterations used in [16] and the multiple solutions of a PPE per time step in [17] can make both these approaches quite time consuming [48]. Multi-phase incompressible flows have also been considered by Shao [43] where coupled and decoupled modelling approaches are compared. The coupled approach utilises a single Poisson equation for both fluids, while the decoupled approach employs a separate PPE in each phase with appropriate normal and shear stress conditions at the interface. The decoupled approach was found to perform better for larger density ratios, but the ratios considered were still relatively small (the largest being 1:1.3). Nevertheless, the potential of multi-phase ISPH has seen its application to a number of complex flow problems including multi-phase Newtonian and viscoelastic flows [51] and thermocapillary flow with appropriate heat transfer and Marangoni force sub-models [47].

70 71 72 73 74 75 76 77 78 79

Despite the improvements ISPH multi-phase flows offer with regard to accurate pressure prediction, its application is still excluded from a class of problems where air retains a significant amount of compressibility alongside an incompressible water phase. Water-air flows of practical interest in coastal and offshore engineering (such as wave breaking, wave impact, and sloshing) typically occupy flow regimes where air should be modelled as compressible, but water remains incompressible. The need for an accurate description of such flows with a careful interface treatment motivates this study. Two-phase incompressible-compressible schemes have been developed for particle methods recently (see Khayyer et al. [19], for example), however these approaches utilise projection methods in both air and water, requiring a smoothing of density across the interface in order to solve the governing PPE. In contrast, 2

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

this paper presents a novel two-phase incompressible-compressible SPH method (ICSPH) for water-air flows, that utilises weakly compressible SPH in the air phase and projection-based incompressible SPH in the water phase. The water phase imposes kinematics on the air phase at the air-water interface while the air phase imposes pressures on the water at the interface. This interface treatment allows a true material discontinuity to be maintained in the flow, with the correct step change in density. To the authors’ knowledge, this coupling of CSPH and ISPH is the first of its kind. This paper presents analytical and fundamental validation cases of the ICSPH method, and it supports a parallel paper that focuses on a specific and challenging application in offshore engineering on wave slam on structures [22]. Given our interest in air-water flows, throughout this paper density ratios are maintained at 1000:1. Four test cases are considered here. The first test compares ICSPH with semi-analytical solutions for the deformation of a two-fluid elliptical drop. The second test considers low amplitude standing wave oscillation for which an analytical expression for the frequency exists [10]. Both tests one and two are standard (semi) analytical comparisons validating the method and air-water coupling. The third test case considers the Kelvin-Helmholtz instability while the final case is a two-phase dam break, where comparisons are made against existing numerical (level-set, boundary element and WCSPH [14, 9, 7]) and experimental results for the impact pressure [52].

96

100

The manuscript is structured as follows: Section 2 details the two-phase ICSPH numerical method and the air-water coupling. Section 3 then presents and discusses the numerical results for each of the aforementioned tests, with relevant experimental and analytical comparisons. Conclusions are made in Section 4.

101

2. The Numerical Method

102

2.1. The Governing Equations Air-water flows are considered in this study, and both phases obey the governing equations of a viscous Newtonian fluid. Namely, the conservation of momentum,

97 98 99

103 104

du 1 = − ∇p + ν∇2 u + f , dt ρ 105 106

(1)

(which retains the same form in both the compressible and incompressible phases), and the conservation of mass,  0, in the incompressible water phase; (2) ∇·u= 1 dρ − ρ dt , in the compressible air phase.

112

The symbols u, p, ρ, ν, and f denote the fluid velocity, pressure, density, constant kinematic viscosity, and constant gravity body force, respectively. The density and viscosity are assigned different values appropriate to each phase (initially: ρa = 1 kgm−3 (air), ρw = 1000 kgm−3 (water), νa = 1 × 10−5 m2 s−1 , and νw = 1 × 10−6 m2 s−1 ). Note that the dilatational/bulk viscosity terms present in the general expression of the compressible momentum equation are neglected in this study. At the interface between the two phases, surface tension effects are also neglected.

113

2.2. The SPH Method

114

In SPH, a variable A at a point r is approximated by a convolution product of A with a weight function, ωh (| r − r |), called a smoothing kernel, with a characteristic or smoothing length h, and is written as  A(r) ≈ A(r )ωh (| r − r |) dr , (3)

107 108 109 110 111

115

Ω

116 117

where Ω is the supporting domain. When discretised over surrounding Lagrangian fluid point masses, or particles, the interpolation is approximated as

3

A(ri ) ≈



Vj Aj ωh (rij ),

(4)

j

123

where Vj is the volume of particle j, and rij is a distance vector, of magnitude rij , between particle i and j. Hereafter ωh (rij ) = ωh (| ri − rj |) will be simply written as ωij . In this paper a quintic spline kernel, continuous to the fifth derivative [32] is used for all cases. A smoothing length of h = 1.3dx is typically used, where dx is the initial particle spacing. To aid the description of the method, the following notation is introduced: let W be the set of all water particles and A the set of all air particles. Accordingly, S = W ∪ A is the set of all SPH particles.

124

2.3. The Incompressible Water Phase

125

2.3.1. SPH Operators In this study the incompressible phase is discretised in exactly the same way as in Xu et al. [49] and Lind et al. [23] using ISPH with particle shifting. The discretised gradient and divergence operators take the form  ∇φi ≈ − Vj (φi − φj )∇Wij , (5)

118 119 120 121 122

126 127 128

j∈W 129 130 131 132 133

for a general variable φ and normalised kernel gradient ∇Wij . This choice of gradient operator is preferred over others due to its demonstrated accuracy in ISPH [23] when used in combination with kernel gradient normalisation [33, 3]. A similar operator is used for calculating the divergence. Kernel gradient normalisation is used only in Eqn. (5) and only in the incompressible phase, with ∇Wij = L(r)∇ωij where,   ∂ω Vj (xj − x) ∂xij  L(r) = ∂ω Vj (yj − y) ∂xij

134

−1  ∂ω Vj (xj − x) ∂yij  . ∂ω Vj (yj − y) ∂yij

(6)

The viscous and Laplacian operator employed is that suggested by Morris et al. [32]: (∇ · μ∇u)i ≈

 mj (μi + μj )rij · ∇ωij uij . 2 + η2 ) ρj (rij

(7)

j∈W 135 136 137 138

139 140 141 142 143 144

Here η is a very small parameter (O(10−5 dx)) included to mitigate computational issues around the singularity at rij = 0. Other forms of the Laplacian exist and can offer improvements in accuracy, especially near free boundaries (see [41, 26], for example). In this study however, the gains in accuracy are small and the simpler (and more computationally efficient) Morris operator is retained. 2.3.2. The Projection Method The projection method [6] was first applied by Cummins & Rudman [11] to the SPH method in order to maintain a divergence-free, incompressible velocity field. A standard first-order time marching scheme is applied, with the density and mass of the incompressible particles taken as constant. Higher order temporal schemes are available for ISPH (e.g. [2]), but are not considered here. Firstly, particle positions, rni , are advected with velocity uni to positions r∗i , r∗i = rni + tuni .

145 146

An intermediate velocity u∗i is then calculated at the position, r∗i , based on the momentum equation without the pressure gradient term,  u∗i = uni + ν∇2 uni + fin t.

147

(8)

(9)

The pressure at time n + 1 can then be obtained from the pressure Poisson equation (PPE), written as

4

1 2 n+1 1 ∇ pi ∇ · u∗i . = ρ t 148 149

(10)

An application of the Laplacian operator (7) and gradient operator (5) results in the following discretised form of Eqn. (10),  Vj (pn+1 − pn+1 )rij · ∇ωij 1  i j 2 = Vj (u∗j − u∗i ) · ∇ωij . 2 2 ρ rij + η Δt

j∈W 150 151 152

Equation (11) forms a linear system for pi that can be solved using an iterative solver (the stabilised bi-conjugate gradient method is used). The desired divergence-free velocity at time n + 1, un+1 , then i results from the projection of u∗i : un+1 = u∗i − i

153

t n+1 ∇pi . ρ



155

156 157 158 159 160 161 162 163

(12)

Finally, the particle positions are advanced in time, rn+1 = rni + t i

154

(11)

j∈W

un+1 + uni i 2

 .

(13)

For further details on the projection method employed herein, including comparisons between different SPH projection approaches and iterative solvers, readers are referred to [48]. 2.3.3. Fickian Shifting Algorithm In the ISPH algorithm presented above, particles will move along streamlines when the Lagrangian Navier-Stokes equations are solved accurately. Therefore, near a stagnation point, for example, particle clustering cannot be avoided. Since the first derivative of the kernel function goes to zero at the particle origin, this particle clustering can persist, increase numerical error, and jeopardise the stability of the simulation. Therefore, to stabilise the simulation, following the approach of Xu et al. [49], after the pressure and velocity fields are calculated and the particles are advanced, the particles are then shifted slightly and the hydrodynamic variables are corrected by the Taylor series approximation, 2 φi = φi + (∇φ)i · δrii + O(δrii  ).

(14)



164 165 166 167 168

Here φ is a general variable, i and i are the particle’s old position and new position respectively, and δrii is the vector between the particle’s new position and its old position. Following Lind et al. [23], the magnitude and direction of the position shift is governed by Fick’s law of diffusion. This acts to shift the particles from regions of high concentration to regions of low concentration, thereby ensuring a relatively equispaced particle distribution. The original mathematical statement of Fick’s law is J = −D ∇C,

169 170 171 172

where J is the flux, C is a concentration, and D is a diffusion coefficient. Assuming that the flux, i.e. the number of particles passing through a unit surface in unit time, is proportional to the velocity of the particles (and of the same order of magnitude), a particle shifting velocity, vs , and subsequently a particle shifting distance, δrs = vs Δt, can be found. Consequently, δrs = −D∇C.

173 174 175

(15)

(16)

Note that D is a new diffusion (or shifting) coefficient which absorbs the constants of proportionality arising between equations (15) and (16). A measure of the particle concentration comes simply from the sum of the kernel function,  Ci = Vj ωij , (17) j∈W

5

176

while a measure of the concentration gradient is given by  ∇Ci ≈ Vj (1 + fij )∇ωij .

(18)

j∈W 177

The term fij is the tensile instability term defined by Monaghan [30], where n

ωij fij = R , ω(dx)

(19)

187

with constants R = 0.2 and n = 4. This term is required to alleviate problems associated with a vanishing kernel gradient when approaching particles are virtually on top of each other. The magnitude of the diffusion coefficient D can be initially estimated through a von Neumann stability analysis of the advection-diffusion equation to give an approximate upper limit of D = 0.5h2 . By combining this estimate with a local CFL condition, Skillen et al. [45] defined a particle-wise diffusion coefficient, D = hΔt|ui |, which is more effective for very violent wave impacts, but it is not used in this study. As in its original application [23], shifting is applied in the incompressible water phase as if the computation were entirely single phase and the air not present (a free-surface shifting correction is therefore applied). This ensures a well-defined water-air boundary as shifting across the interface is restricted and any artificial mixing of phases is minimised.

188

2.4. The Compressible Air Phase

189

2.4.1. SPH Operators Within the compressible phase, the governing equations are solved as in the weakly compressible SPH (WCSPH) approach, first developed by Lucy [25] and Gingold & Monaghan [13]. Investigations have shown that the locally conservative and coherent combination of SPH gradient operators recommended by Colagrossi & Landrini [7] and [33] is the most effective discretisation for the air phase. Unlike the incompressible phase, the gradients in the compressible phase are not normalised (neither is zeroth-order consistency enforced) in order to maintain pairwise particle antisymmetry and associated momentum conservation properties. In the context of WCSPH, exact conservation is often recommended over kernel corrections for long-term stability and accuracy [38]. However, the stabilising benefits of conservative gradients are less apparent in ISPH due to the constant particle density and absence of local compressible pressure increases that act to self-regulate particle distributions. Indeed, the necessity for additional particle regularisation (shifting) in ISPH means that kernel corrections may be applied more readily in the incompressible phase with no detriment to long-term stability (as is the case here). A numerical example is presented in Section 3.1.2 that demonstrates the above and supports the choice of pressure gradient approximation used in either phase.

178 179 180 181 182 183 184 185 186

190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207

So, for a particle i ∈ A, the conservative and uncorrected pressure gradient for the compressible momentum equations is given by  ∇pi = (pi + pj )Vj ∇ωij , (20) j∈S

208

while the discretised conservation of mass becomes,  dρi = ρi (ui − uj ) · ∇ωij Vj . dt

(21)

j∈S

209 210 211

In this compressible phase the volume Vj is permitted to vary. The governing equations in the compressible phase are updated in time using a standard second-order two-step method, which, for equation (21), reads as

6

ρn+1 i 212 213

214 215 216 217

=

ρni



n

n−1  3 dρi 1 dρi , + Δt − 2 dt 2 dt

with an equivalent expression used for the velocity update. Particle positions are then updated as in Eqn. (13). 2.4.2. The Equation of State The air flows considered in this paper are of low Mach number, M a ≤ 0.1, satisfied for all test cases. To a good approximation, therefore, the air phase may be considered isothermal (assuming no external heat transfer) [35] with the stiffened ideal gas equation, p = c2 (ρ − ρ0 ),

218 219 220 221

222 223

224 225 226 227 228 229 230 231

(22)

(23)

providing a physical thermodynamic relation between pressure and density [39]. The constant c is the speed of sound while ρ0 = 1 kgm−3 is the initial air density. As mentioned in the Introduction, the air flows in this paper may rightly be termed weakly compressible as Mach numbers are small and thermal effects are negligible. 2.4.3. Diffusion and Artificial Viscosity We use the conventional artificial viscosity term [29] to provide numerical stability in the air phase:

 uij · rij ci + cj . (24) Πij = −αh 2 rij + η 2 ρi + ρj An artificial viscosity coefficient between 0.01 ≤ α ≤ 0.1 is used in this paper, depending on the test case (a table providing all numerical and physical parameters used is presented for each test case). For the α values considered here, the small physical laminar viscosity tends to have negligible influence in the presence of artificial viscosity. In some cases it was also found to be beneficial to employ the Fickian shifting algorithm in the compressible phase, as well as the incompressible phase (for details please refer to the relevant parameter summary table for each test case). As is consistent with WCSPH schemes that require accurate pressure predictions, every twenty to thirty time steps a Shepard filter [44] is also applied to the density field. The filtered density, ρ˜i , is given by   j∈A ρj Vj ωij ρ˜i = ρj V j ω ˆ ij =  . (25) j∈A Vj ωij j∈A

232 233 234 235 236

237 238 239 240 241 242 243 244

2.5. Coupling the Two SPH Methods Although both phases employ SPH methods to solve their respective governing equations, the two SPH approaches are quite distinct numerically and must be coupled in a mathematically viable and physically correct manner. The manner in which the two phases are coupled is described below and illustrated in Figure 1. 2.5.1. Continuity in Pressure Formally, the normal components of the Cauchy stresses in either fluid are continuous across the fluidfluid interface [1]. However, the problems of interest are water-air flows where the role of fluid viscosity and surface tension is small. Therefore, the pressure (the isotropic part of the Cauchy stress) is taken to be continuous across the interface. Hence, in the two-phase ICSPH method, the compressible phase provides a pressure boundary condition at the interface for the incompressible phase (Fig. 1(a)). Not only does this ensure continuity in pressure at the interface, but it also provides the necessary Dirichlet boundary condition for the solution of the pressure Poisson equation (Eqn. (10)). The pressure values

7

(a)

(b)

Figure 1: An illustration of the manner in which the two phases are coupled: (a) Compressible particles (white) assign a pressure to incompressible interface particles, and thereby provide the necessary pressure boundary condition. (b) Incompressible particles (blue) are used in velocity calculations in the momentum and mass equations in the compressible phase, and thereby provide a velocity boundary condition.

245 246

of neighbouring compressible particles are assigned to the interface incompressible particles, i, using the normalised SPH interpolation,  Pi = Pj ω ˆ ij Vj , (26) j∈A

247 248 249 250 251 252

253 254 255 256 257 258 259 260 261 262 263 264

where the summation is over surrounding compressible particles only and ω ˆ is defined by Equation (25). To determine which incompressible particles are interface particles, the criterion proposed by Lee et al. [20] is used, whereby interface particles are those whose position vector ri satisfies ∇ · ri < 1.6 in two dimensions. The pressure of interface particles can be enforced in the PPE by simply inserting the desired values into the right-hand side of the linear system (Eqn. (11)) and adjusting the relevant rows in the coefficient matrix before solution. 2.5.2. Continuity in Velocity Across the interface, the normal and tangential velocity components should also be continuous (due to viscosity). Consequently, the coupling is completed by using the incompressible phase to provide a velocity boundary condition for the compressible phase. Unlike the incompressible phase, a pressure boundary condition is not formally required for the compressible air as pressure information is already known explicitly from the equation of state. Instead, it is the velocity boundary condition at the interface that provides new flow information for the air phase. At the discrete level the boundary condition is implemented as in Fig. 1(b): incompressible particles are used to complete the kernel support of nearby air particles and act as a (moving) boundary for the governing equations in the air phase. Note that, in completing kernel support, incompressible particles provide a pressure (as well as velocity) to improve accuracy in compressible pressure gradient calculations.

265 266 267 268 269

270 271 272

The exchange of kinematic and dynamic information across the interface in this manner results in a fully coupled system with a well-defined material discontinuity. There is no smoothing of density across the interface (as employed in volume fraction SPH methods [28]) and realistic density ratios can also be maintained. 2.6. Time-Stepping and Solid Boundary Conditions Generally, in SPH the time step sizes are determined through the Courant-Friedrichs-Lewy (CFL) condition, 8

Δt = Crh/U, 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293

(27)

where the Courant number, Cr ≤ 1. Incompressible SPH and compressible SPH utilise different characteristic velocities, U , in their formulations: the speed of sound, c, is used in the compressible case, while a characteristic kinematic flow velocity, Uk , is used in the incompressible case. As a consequence, the time steps for incompressible SPH may be an order of magnitude larger than those for compressible SPH, depending on the sound speed (typically c ≥ 10Uk for WCSPH [29]). Therefore, substantial gains in efficiency can be achieved by employing separate time step sizes in either phase. Here, the Courant number in the compressible phase is taken as Cr = 0.1 for all cases. In line with observations made in Lee et al. [20], the incompressible time step, Δtinc , is then chosen to be ten times that of the calculated compressible time step (i.e. Δtinc =10Δtcom = 10Crh/c). The equations governing the incompressible phase are, therefore, updated once every 10 compressible time steps. The use of separate time-stepping schemes is widely-used in multi-scale simulations where different time scales exist (e.g. [12]), and one process may look stationary (for a short time) compared to a rapidly evolving second process. The problem here is analogous as there is a need to resolve the high frequency acoustic transients in the compressible phase which evolve on a faster temporal scale than any flow features in the incompressible phase; in other words the incompressible phase looks to be stationary for short intervals, Δtinc . From a mathematical viewpoint, the separate time-step approach employed here is also entirely consistent with the global first order nature of the method, as at any point interfacial boundary information is also accurate to O(Δtinc ). In Section 3.1.2 a numerical example is provided supporting this approach and the choice of incompressible time step.

294

296

Note that, when required, solid wall boundary conditions are imposed using the mirror particle approach [11].

297

3. Numerical Results

298

3.1. Two-Layer Elliptical Drop Deformation To validate the pressure-velocity coupling across the air-water interface, we consider the semi-analytical solution for deformation in a two fluid layer elliptical drop [31]. An inner ellipse of dense fluid (water), with a major and minor axis values of a and b, respectively, is surrounded by an air layer occupying the space between the inner ellipse and an outer elliptical surface with major and minor axis values A and B, respectively. Note that a, b, A, B are all time dependent, and initially A(0) = B(0) = 1 and a(0) = b(0) = 0.5 (so the two fluids initially occupy two concentric circles, as seen in Fig. 2). The analytical solution assumes both phases are inviscid and incompressible. The latter condition can be accurately represented in the air phase with a suitably large speed of sound. A value of c = 25m/s is used here. Accordingly, note that this test is a measure of the coupling of the two phases, not the role of air compressibility. For this test case, Table 1 summarises the key physical and numerical parameter values used.

295

299 300 301 302 303 304 305 306 307 308 309

310 311

312

3.1.1. Drop Elongation Given an initial incompressible velocity field of the form u =

σ(t)x,

(28)

v

−σ(t)y,

(29)

=

the inviscid governing equations for both phases can be reduced to  2 2 dσ 2 B −A =σ , dt (B 2 + A2 ) 9

(30)

Figure 2: Initial set-up for the two-layer drop oscillation test case, with the dense inner fluid (water) surrounded by a layer of air (shaded area).

Numerical Parameter dx (m) Δt/Δtcom Artificial viscosity, α Shepard filter period Shifting coefficient, D Physical Parameter ρ (kg/m−3 ) ν (m2 s−1 ) c (m/s)

Water phase value 0.01 10 0.25h2

Air phase value 0.01 1 0.01 20 time steps 0

1000 1 × 10−6 -

1 1 × 10−5 25

Table 1: Summary of numerical and physical parameter values used for two-layer elliptical drop deformation.

313 314 315

with an identical equation for the inner fluid, where A and B are replaced by a and b, respectively [31]. Initially, σ is chosen to be σ(0) = 0.5. Values of A, B, a and b can be found from velocity considerations at the interface and free surface. Namely, if u = σx then dA = σA (31) dt on the outer ellipse surface. Similar relations hold for the remaining major/minor axis values. Equations (30) and (31) can be integrated numerically (rapidly and to very high accuracy) using standard numerical integration techniques (here a second order predictor-corrector scheme is used with a time step size of order 10−8 ). The semi-analytical solution, A, quickly asymptotes to a straight line. Figure 3 shows the semi-analytical solution (from Eqns. (30) and (31)) compared with results from ICSPH for the size of the outer major axis of the drop, A. The agreement is excellent, indicating that the pressure-velocity air-water coupling is working well. Table 2 shows the time averaged major axis value, Aav , (averaged over the first 3s) and the associated percentage relative error for different particle spacings. In all cases, the relative error is small and convergence is indicated. The convergence rate is around linear, which is u=

316 317 318 319 320 321 322 323 324

10

Figure 3: Variation of the outer major axis, A, with time for a stretching elliptical drop. The line denotes the semi-analytical results while the squares are predictions from ICSPH.

325

typical for the cases studied in this paper. dx (m) 0.04 0.02 0.01

Aav (m) 4.932 4.936 4.953

Rel. Error (%) 0.29 0.23 0.14

Table 2: Numerical results for the time-averaged position of the major axis, Aav , of a stretching elliptical drop. The semi-analytical value is 4.947m (3dp).

326 327 328 329 330 331 332

333 334 335 336 337 338 339 340

3.1.2. Drop Oscillation The introduction of an appropriate conservative body force (or potential) to the inviscid governing equations can result in the drop oscillating about its equilibrium position, with extension along both minor and major axes. The oscillations are reminiscent of those due to surface tension, making this a convenient test of drop deformation for flows where surface tension is negligible. With the introduction of a body force, f = Ω2 (r − r0 ), into the governing equations (1), the simplified equation governing motion of the ellipse is  B 2 − A2 dσ  2 2 = σ +Ω . (32) dt (B 2 + A2 ) A value of Ω2 = 1 produces the drop oscillations observed in Figure 4. Despite the large (unsmoothed) density difference, particle distributions remain well-behaved around a welldefined interface - see Figure 5, a close up of the region highlighted in Fig. 4(d). Clearly, the particle shifting has worked well to inhibit any artificial diffusion and dispersion of particles into either phase (which can readily and erroneously occur in particle approaches, and has been treated in the past by adding adhesion-like terms to the momentum equations [7, 31]). In preventing numerical dispersion of particles into a different phase, some clustering is evident at the interface. It is important to highlight that this is a superficial 11

341 342 343

effect that has little impact on the accuracy and stability of the simulation (as evidenced by the results in this paper), and, in fact, is quite beneficial in providing a clearly defined and well-resolved interface.

344 345 346 347 348 349 350 351 352 353 354

Figure 6 shows the pressure and vertical velocity contours over the whole drop (through both phases) at t = 2.2s. Note that, for brevity, horizontal velocity contours are not shown, but behave identically to the vertical component given the symmetry of the problem. As required, both pressure and velocity variables are continuous across the interface, supporting the validity of the interface treatment and coupling approach. There is a small amount of noise in the velocity value around (3.1, 3.5), but this is clearly a very localised disturbance within the compressible phase that has little impact on global accuracy. Figure 7 compares the value of the outer major axis, A, with the semi-analytical solution from (32) over time. Once again, the agreement is excellent, with drop deformation well-captured by ICSPH for relatively long times (at least up to 6s).

(a) t = 0

(b) t = 1.1s

(c) t = 2.2s

(d) t = 3.5s

Figure 4: Oscillating drop particle positions at select times. The inner fluid is water and the outer fluid is air (in the physical density ratio 1000:1). 355

To support the choice of pressure gradient approximation used in either phase, the effect 12

Figure 5: Particles near the water-air interface inside the boxed region in Figure 4(d). The density ratio is 1000:1. The interface remains distinct with variables unsmoothed.

356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384

of different pressure gradient formulations on the solution will now be explored. As Fig. 7, Fig. 8 plots the variation of the outer major axis A with time but for three different combinations of pressure gradient approximation. Option A uses Eqn. (5) in both phases with kernel gradient correction [3] (with all other parameters as Table 1); option B utilises the conservative operator, Eqn. (20), in both phases and without any conservation disrupting corrections. Results are plotted alongside the current approach (a combination of option A and B in the water and air phase, respectively) and the semi-analytical solution. Clearly, the current approach yields the best results with stable and accurate solutions predicted. Option A corrects gradients in the air phase and thereby disrupts exact momentum conservation, causing considerable error in the air phase, particularly at the outer free surface as shown in Fig. 9(a). Application of shifting in the air phase will improve the solution, but given that accurate (and conservative) solutions can be obtained through careful choice of the pressure gradient alone (Eqn. (20)), the shifting process may remain an optional numerical regularisation device for the air phase. As mentioned in Section 2.4.1, while beneficial to the compressible air phase, the conservative form (Eqn. (20)) is not as effective for ISPH. Option B uses the conservative form for both phases, with no shifting in either phase due to the operator’s self-particle-regulating properties [38]. Good agreement with the analytical solution is observed until around t = 1.9s, when the solution again crashes. This time, however, instability arises at the surface of the incompressible phase (Fig. 9(b)). As mentioned, the conservative form alone is insufficient for particle self-regulation in ISPH (essentially due to the accurate pressure field forcing the particles to follow their Lagrangian paths more precisely) and additional particle regularisation (shifting) is required. Indeed, if additional regularisation is necessary then ISPH may readily benefit from non-conservative (non-regulating) gradients and corrections that guarantee first order consistency, as is the case here. As mentioned in Section 2.6, the incompressible phase uses a time step that is a factor of 10 times larger than that of the compressible phase (in the interests of computational efficiency). To demonstrate that this time step factor is of an appropriate magnitude, Figure 10 shows the (semi) analytical solution for the oscillation of the outer major axis 13

(a) Pressure, p

(b) Vertical velocity, v

Figure 6: Pressure and vertical velocity contours across the two-fluid drop (in the physical density ratio 1000:1) indicating continuity across the interface (t = 2.2s).

394

compared with ICSPH predictions for different incompressible time step sizes. Even a two order of magnitude increase in the incompressible time step size over the compressible (with Δtinc = 100Δtcom ) results in good agreement with the analytical solution (the smaller Δtinc is, of course, more accurate). Indeed, the CFL condition remains satisfied in both phases and, as interface information remains accurate to O(Δtinc ) and time steps remain small (with Δtcom ∼ 10−5 s here), the simulation remains stable and quite accurate. Nevertheless, a choice of Δtinc = 10Δtcom will continue to be used throughout this paper, as a compromise between temporal and interfacial accuracy against computational efficiency. For all the cases considered here, a choice of Δtinc = 10Δtcom results in a 15-20% reduction in total CPU time with little or no impact on accuracy.

395

3.2. Standing Gravity Waves

396

As depicted in Figure 11, the second test case considers an oscillating standing wave comprised of some dense fluid (water), positioned beneath a body of lighter fluid (air), in the density ratio 1000:1. No penetration conditions are imposed at the top and bottom of the domain, while periodicity is enforced at the side walls. For small amplitudes and viscosities, the frequency of the oscillation of the fluid-fluid interface (initially of the form y = a cos(kx)) can be determined analytically [10] from

385 386 387 388 389 390 391 392 393

397 398 399 400

f2 = 401 402 403 404 405 406 407 408 409 410 411

1 1 ((ρw − ρa )gk) tanh (kD) . 4π 2 (ρw + ρa )

(33)

Here k = 2πm−1 , D = 0.5m, and the length of the periodic domain is L = 1m. The time period, T , can be determined from Eqn. (33) straightforwardly. For this particular test case, Table 3 presents the numerical and physical parameter values used. Table 4 shows numerical results obtained at the third half period, including the predicted values of 3T /2 and the maximum y-position (ymax ) of the interface. The initial amplitude of the wave is a = 0.025m. Results are compared with the analytical values of 3T /2 (from Eqn. (33)) and ymax . Rapid convergence in 3T /2 is observed initially (second order, approximately) before the error remains fixed at around 1% for dx = 0.01m and lower. This limiting error arises due to physical non-linear wave effects in the numerical simulation not considered in the linear analytical solution. The calculated value of ymax converges uniformly (approximately linearly), with relative error reducing to 3.0% at dx = 0.005m. Table 5 presents results for the same test case and parameter values, but with a larger initial wave amplitude

14

Figure 7: Variation of the outer major axis, A, with time for an oscillating elliptical drop. The line denotes the semianalytical results while the squares are predictions from ICSPH.

Figure 8: Variation of the drop outer major axis, A, with time, calculated using different pressure gradient formulations. Option A applies Eqn. (5) everywhere (with correction); Option B uses Eqn. (20) everywhere (no corrections).

15

(a) t = 0.03s

(b) t = 1.85s

Figure 9: Two-layer drop particle distributions for different pressure gradient formulations (in combinations not used in this work): (a) Option A (Eqn. (5) based everywhere), (b) Option B (Eqn. (20) based everywhere). Times given are those at which the solution crashes.

Numerical Parameter dx (m) Δt/Δtcom Artificial viscosity, α Shepard filter period Shifting coefficient, D Physical Parameter ρ (kg/m−3 ) ν (m2 s−1 ) c (m/s)

Water phase value See Table 4 10 0.25h2

Air phase value See Table 4 1 0.05 20 time steps 0.25h2

1000 1 × 10−6 -

1 1 × 10−5 10

Table 3: Summary of numerical and physical parameter values used for standing gravity wave and Kelvin Helmholtz test cases.

412 413 414 415 416 417 418

419 420 421 422 423 424 425 426 427

of a = 0.05m. Non-linear effects now have a greater role in the simulation, and measurements of 3T /2 against linear theory differ by around 3.0%, regardless of discretisation. The numerical value of 3T /2 does converge, however, to a fully non-linear approximation of 1.24s (to two decimal places). The value of ymax still converges uniformly (but more rapidly than the case a = 0.025m with a rate closer to one) to an error as small as 0.9% for dx = 0.005m. Indeed, results for ymax are improved for a = 0.05m due to an increase in relative resolution around the interface (there are a greater number of particles per wave height). 3.3. Kelvin-Helmholtz Instability The Kelvin-Helmholtz instability provides a challenging test case for two-phase flows. Although usually studied as a fully incompressible problem, it remains a valid test of the ICSPH method and the interface treatment if the speed of sound in the air phase is suitably large and flow velocities are small. As before, we chose c = 10 m/s to yield a Mach number of 0.05 for this problem. In a similar manner to Price [37], the initial set-up for this case is given in Figure 12. The water phase occupies a central rectangular region, [0, 1] × [0.25, 0.75]m2 , surrounded by air. Periodicity is applied in the x-direction, with a no penetration condition imposed at the top and bottom boundaries using mirror particles. The water and air phase are initialised to move horizontally, but in opposite directions, with Uw = −Ua = 0.5 m/s. The physical 16

Figure 10: Variation of the drop outer major axis, A, with time, for different incompressible time step sizes.

Figure 11: Standing gravity wave initial set up.

17

dx (m) 0.04 0.02 0.01 0.005

3T /2 (s) 1.107 1.187 1.193 1.217

Rel. Error (%) 8.0 1.4 0.9 1.0

ymax (m) 0.473 0.492 0.501 0.509

Rel. Error (%) 9.8 6.2 4.6 3.0

Table 4: Numerical measurements for the third half period, 3T /2, and maximum y-position of the interface (at t = 3T /2) for a standing gravity wave. The initial amplitude is a = 0.025m.

dx (m) 0.04 0.02 0.01 0.005

3T /2 (s) 1.167 1.221 1.240 1.244

Rel. Error (%) 2.9 1.4 3.0 3.2

ymax (m) 0.494 0.532 0.541 0.545

Rel. Error (%) 10.0 3.3 1.7 0.9

Table 5: Numerical measurements for the third half period, 3T /2, and maximum y-position of the interface (at t = 3T /2) for a standing gravity wave. The initial amplitude is a = 0.05m.

428 429 430 431 432

and numerical parameter values used for this case are identical to those of the standing gravity wave, as summarised in Table 3; importantly, however, gravity is now omitted. In the absence of gravity and surface tension (as is the case here), hydrodynamic instability will always develop at the interface for a non-zero velocity difference in the two fluids [4]. Linear stability analysis predicts exponential growth in the magnitude of any applied perturbation to the primitive variables, δv. Namely, δv ∼ exp(ωt),

433

where ω=

434 435

436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452

2π(ρw ρa )1/2 |Uw − Ua | . λ(ρw + ρa )

(34)

(35)

Here λ is the wavelength of the initial perturbation. As in Price [37], to start the simulation we introduce a vertical velocity perturbation around the interface according to  a sin(2πx/λ) if |y − 0.25| < 0.025, δv = a sin(−2πx/λ) if |y − 0.75| < 0.025, where a = 0.025 m and λ = 1/6 m. Figure 13 shows the subsequent development of the interface disturbance at different times for the density ratio 1:1000 (dx = 0.01m). Note that, to the authors’ knowledge, no other SPH method has considered this flow at this density ratio. Characteristic KelvinHelmholtz flow structures are observed, with periodically arranged branches of the denser fluid growing rapidly over time. As is observed for larger density ratios (ρw /ρa  10), vortex rolls do not readily form (see Shadloo & Yildiz [42], for example, who observed similar structures to Fig. 13(c) at ratio 1:10). Figure 14 compares the growth of the maximum position of the interface disturbance, ηmax , in time with the theoretical predictions from linear stability theory. Results are presented for the density ratio 1:1000 and an additional test case where ρw /ρa = 10. For both density ratios, the agreement with linear theory in the early stages of the flow (t < 0.1s) is good. Figure 15 provides a close-up of Figure 14 for t ≤ 0.1s, where the difference in predictions of the interface position is always less than one quarter of a particle spacing. The discrepancy between theoretical and numerical predictions increases in the later stages, although this is expected as the linear perturbation analysis becomes increasingly invalid. However, quantitatively similar behaviour is still observed, with ICSPH and linear theory predicting similar rates of disturbance growth, appropriate to the density ratio. In particular, more rapid perturbation growth is seen in the 1:10 case (which reaches ηmax = 0.8m at approximately t = 0.3s), compared to the 1:1000 case (which reaches ηmax = 0.8m at approximately t = 1.0s).

18

Figure 12: The initial domain set up for the Kelvin-Helmholtz instability test case. Periodicity is imposed in the x-direction.

453

3.4. Two-Phase Dam Break

454

In this section the ICSPH predictions for a standard dam break flow are presented and compared with other numerical methods (including boundary element [14], a second-order finite difference ENO level-set scheme [9] and two-phase WCSPH [7]) and experimental data [52]. The flow is well approximated by a single phase SPH model up to the point of overturning, where trapped air pockets (with associated compressibility) can then play a significant role in dynamics. Therefore, this case provides a test of the interface conditions for a range of air-water interactions.

455 456 457 458 459

Numerical Parameter dx (m) Δt/Δtcom Artificial viscosity, α Shepard filter period Shifting coefficient, D Physical Parameter ρ (kg/m−3 ) ν (m2 s−1 ) c (m/s)

Water phase value 0.015 10 0.25h2

Air phase value 0.015 1 0.1 30 time steps 0.25h2

1000 1 × 10−6 -

1 1 × 10−5 50

Table 6: Summary of numerical and physical parameter values used for the two-phase dam break case.

460 461 462 463 464 465 466 467

A column of water is surrounded by air in a domain sized [0, 5.366]m×[0, 3]m. Table 6 presents the numerical and physical parameters used here. For a water column of dimension [0, 1]m×[0, 1]m, Figure 16 shows the evolution of the water front (the toe position) with time and comparisons with a single phase boundary element method (as in [7]). The agreement is excellent and demonstrates the single phase nature of the flow in this early stage. For a slightly different water column geometry ([0, 2]m×[0, 1]m), Figure 17(a) compares the flow profiles of ICSPH and the single-phase boundary element method (BEM) at time t ≈ 1.9s [14, 7]. The air phase is not pictured for the purpose of clarity. As expected the role of the air phase is still minimal at this point in the flow, and the single and two-phase results for the 19

(a) t = 0.52s

(b) t = 0.91s

(c) t = 1.04s Figure 13: Development of the interface disturbance in the water phase at select times. The surrounding air phase is not shown for clarity.

468 469 470 471 472 473 474 475 476 477 478 479 480 481 482

water surface are in reasonable agreement. Furthermore, the pressure predictions from incompressible SPH remain characteristically smooth and noise-free, highlighting the efficacy of the compressible air phase in prescribing a suitable dynamic (pressure) boundary condition at the water-air interface. Figure 17(b) compares ICSPH results with a two-phase level-set method [9, 7] for the same dam break at the later time of t ≈ 2.16s. Here, there is a notable difference in the flow profiles: the rebounding liquid jet forms at approximately the same point in both methods (x ≈ 3.8m), but is more upwardly directed and fragmented for ICSPH. The predictions for the size and the shape of the trapped air pocket also differ. However, it is important to note that comparisons are being made against a two-phase level-set model which assumes incompressibility in both phases. Indeed, one would expect predicted air pocket shapes to be smaller in ICSPH as air compression and changes in volume are permitted. Indeed, the modelling assumption that the air pocket in the level-set method will not change volume is a potentially unphysical feature. The results from the WCSPH method [7] agree well with those of the incompressible level-set method, but the inability to choose physically meaningful speeds of sound means that the air phase in the WCSPH method is more incompressible than the water, raising similar issues over the form of the surface profile and air pocket. In an attempt to discern which model provides a better physical

20

Figure 14: Comparison between ICSPH (solid line) and linear stability (dashed line) predictions for the maximum interface disturbance position with time. The red lines are results for the density ratio 1:10, while the black lines denote the ratio 1:1000. Colour figures are available online.

Figure 15: A close-up of Fig. 14 for t ≤ 0.1s comparing ICSPH (solid line) and linear stability (dashed line) results for the maximum interface disturbance position with time. The red lines are results for the density ratio 1:10, while the black lines denote the ratio 1:1000. Colour figures are available online.

21

Figure 16: Comparison between ICSPH and BEM [14, 7] predictions of the evolution of the water front (or toe) in a dam break. The line denotes ICSPH results while the squares are from the BEM study [14, 7].

498

representation of this flow, Figure 18 presents pressure measurements for ICSPH, the two-phase WCSPH method [7], and experimental results from Zhou et al. [52]. The experimental pressure transducer occupies a circle of diameter 0.09m, centred at 0.16m above the bed, on the vertical wall opposite the initial water column. The dimensions used in the simulation are now scaled to match experiment, where the water column occupies [0, 1.2]m×[0, 0.6]m. Incompressible-compressible SPH compares well with experiment up to t ≈ 1.5s, and, unlike WCSPH, ICSPH captures the slight overshoot after the initial impact at t = 0.7s. There is a slight delay and notable over-prediction of the secondary pressure peak at t = 1.5s, which is also observed for the two-phase WCSPH method (which predicts a pressure peak of a similar magnitude). This disagreement is likely a consequence of the 2D assumption in both methods; the secondary pressure peak arises due to the water impact of the plunging front in Fig. 17(a), and such breaking wave-like events (with associated air entrapment) are well-known to be three dimensional [5]. Importantly, in these later stages, the WCSPH results demonstrate large amplitude high frequency pressure oscillations which are not observed in ICSPH or experiment. These oscillations likely result from an overly-incompressible air pocket undergoing rapid expansion and contraction when entrapped, suggesting that ICSPH (with its physical descriptions of compressibility in either phase) may provide the more representative predictions for the flow at this stage.

499

4. Conclusions

500

This paper has presented a novel two-phase incompressible-compressible (water-air) SPH numerical method (ICSPH) that couples compressible SPH with a truly incompressible (projection-based) SPH approach. Crucially, the interfacial kinematic and dynamic boundary conditions enable a true material discontinuity to be maintained, with the physically correct step change in density for air-water flows. The method has been validated using a semi-analytical solution for droplet oscillation, a low-amplitude standing wave analytical solution, the Kelvin-Helmholtz instability, and numerical and experimental results from the literature for a benchmark dam break case. The results from the method are promising: high levels of accuracy are obtained in the (semi) analytical comparisons and convergence is demonstrated. Pressure predictions for the dam break flow agree well with experimental data up to the point of impact

483 484 485 486 487 488 489 490 491 492 493 494 495 496 497

501 502 503 504 505 506 507 508

22

(a)

(b)

Figure 17: Comparisons of the fluid profile of the incompressible phase with other numerical studies. (a) Comparison with a BEM study [14, 7] (black squares) at t ≈ 1.9s (b) Comparison with a two-phase level-set study [9, 7] (black squares) at t ≈ 2.16s. For clarity, the compressible air phase is not shown.

511

of the overturning front, when the flow becomes three dimensional. However, subsequent pressure predictions do show a reduction in the spurious high frequency pressure oscillations observed with two-phase WCSPH, likely due to a more physical representation of the air pocket compressibility.

512

5. Acknowledgements

513

Support from EPSRC grants EP/H018638/1 (2010-2012) and EP/L014661/1 (2014-2017) is gratefully acknowledged.

509 510

514

515

516 517

[1] Batchelor, G. K. 1967 An introduction to Fluid Dynamics. Cambridge University Press. [2] Bøckmann, Arne, Shipilova, Olga & Skeie, Geir 2012 Incompressible {SPH} for free surface flows. Computers & Fluids 67, 138 – 151.

519

[3] Bonet, J. & Lok, T. S. L. 1999 Variational and momentum preservation aspects of smooth particle hydrodynamics formulations. Computer Methods in Applied Mechanics Engineering 180, 97–115.

520

[4] Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. New York: Dover.

521 522

[5] Chang, K.-A. & Liu, P. L.-F. 1998 Velocity, acceleration and vorticity under a breaking wave. Phys. Fluids 10, 327.

523

[6] Chorin, A. J. 1968 Numerical solution of the Navier Stokes equations. J. Math. Comp. 22, 745–762.

524 525

[7] Colagrossi, A. & Landrini, M. 2003 Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J. Comput. Phys. 191, 448–475.

526

[8] Cole, R. H. 1948 Underwater Explosions. Princeton University Press.

527

[9] Colicchio, G., Landrini, M. & Chaplin, J.C. 2003 Level-set modelling of the air-water flow generated by a surface piercing body. in: Proc. 8th Int. Conf. on Numerical Ship Hydrodynamics, Korea .

518

528 529

530

[10] Crapper, G. D. 1984 Introduction to Water Waves. Ellis Horwood. 23

Figure 18: Pressure measurements on wall for the dam break experiment described in Zhou et al. [52] compared with predictions from ICSPH and a two-phase WCSPH method [7]. The data used here has been digitised from the literature.

531

[11] Cummins, S. J. & Rudman, M. 1999 An SPH projection method. J. Comput. Phys. 152, 584–607.

532

[12] Ganine, V., Hills, N. J. & Lapworth, B. L. 2013 Nonlinear acceleration of coupled fluidstructure transient thermal problems by anderson mixing. International Journal for Numerical Methods in Fluids 71 (8), 939–959.

533 534

535 536

537 538

539 540 541

542 543

544 545

546 547

548 549 550

551 552 553

[13] Gingold, R.A. & Monaghan, J.J. 1977 Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181, 37589. [14] Greco, M., Landrini, M. & Faltinsen, O.M. 2004 Impact flows and loads on ship-deck structures. J. Fluids and Structures 19, 251–275. [15] Grenier, N., Antuono, M., Colagrossi, A., Touze, D. Le & Alessandrini, B. 2009 An hamiltonian interface SPH formulation for multi-fluid and free surface flows. J. Comput. Phys. 228, 8380-8393. [16] Hu, X. & Adams, N. 2007 An incompressible multi-phase SPH method. J. Comput. Phys. 227, 264–278. [17] Hu, X. Y. & Adams, N. A. 2009 A constant-density approach for incompressible multi-phase SPH. Journal of Computational Physics 228, 2082–2091. [18] James, S. C., Jones, C. A., Grace, M. D. & Roberts, J. D. 2010 Advances in sediment transport modelling. Journal of Hydraulic Research 48 (6), 754–763. [19] Khayyer, A., Gotoh, H., Ikari, H. & Tsuruta, N. 2013 A novel error-minimizing scheme to enhance the performance of compressible-incompressible multiphase projection-based particle methods. 8th international SPHERIC workshop, Trondheim, Norway pp. 68–73. [20] Lee, E.-S., 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, 8417–8436. 24

554 555

556 557 558

559 560 561

[21] Li, Q., Ouyang, J., Yang, B. & Jiang, T. 2011 Modelling and simulation of moving interfaces in gas-assisted injection moulding process. Applied Mathematical Modelling 35 (1), 257 – 275. [22] Lind, SJ, Stansby, PK, Rogers, BD & Lloyd, PM 2015 Numerical predictions of water– air wave slam using incompressible–compressible smoothed particle hydrodynamics. Applied Ocean Research 49, 57–71. [23] Lind, S.J., Xu, R., Stansby, P.K. & 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. J. Comput. Phys. 231 (4), 1499 – 1523.

563

[24] Lind, S. J., Stansby, P. K. & Rogers, B. D. 2012 A multiphase incompressible-compressible smoothed particle hydrodynamics method. Proc. 7th SPHERIC workshop pp. 324–332.

564

[25] Lucy, L. 1977 Numerical approach to the testing of the fission hypothesis. Astron. J. 82, 1013–1024.

565

´ , F., Antuono, M., Gonza ´ lez, L. M. & Colagrossi, A. 2011 Theoretical analysis of the [26] Macia no-slip boundary condition enforcement in SPH methods. Progress of Theoretical Physics 125 (6), 1091–1121.

562

566 567

568 569 570

[27] Marrone, S., Antuono, M., Colagrossi, A., Colicchio, G., Touze’, D. Le & Graziani, G. 2011 Delta-SPH model for simulating violent impact flows. Computer Methods in Applied Mechanics and Engineering 200, 1526–1542.

572

[28] Monaghan, J.J. & Kocharyan, A. 1995 SPH simulation of multi-phase flow. Computer Physics Communications 87, 225–235.

573

[29] Monaghan, J. J. 1994 Simulating free surface flows with SPH. J. Comput. Phys. 110, 399–406.

574

[30] Monaghan, J. J. 2000 SPH without a tensile instability. J. Comput. Phys. 159, 290–311.

575

[31] Monaghan, J. J. & Rafiee, Ashkan 2013 A simple SPH algorithm for multi-fluid flow with high density ratios. International Journal for Numerical Methods in Fluids 71 (5), 537–561.

571

576

577 578

579 580

[32] Morris, J. P., Fox, P. J. & Zhu, Y. 1997 Modelling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136, 214–226. [33] Oger, G., Doring, M., Alessandrini, B. & Ferrant, P. 2007 An improved SPH method: towards higher order convergence. J. Comput. Phys. 225, 1472–1492.

582

[34] Olsson, Elin & Kreiss, Gunilla 2005 A conservative level set method for two phase flow. Journal of Computational Physics 210 (1), 225 – 246.

583

[35] Panton, R.L. 2005 Incompressible Flow . Wiley.

584

[36] Paul, E. L., Atiemo-Obeng, V. & Kresta, S. M. 2004 Handbook of Industrial Mixing: Science and Practice. New Jersey: John Wiley & Sons.

581

585

586 587

588 589

590 591

592 593

[37] Price, D. J. 2008 Modelling discontinuities and Kelvin-Helmholtz instabilities in SPH. J. Comput. Phys. 227, 10040–10057. [38] Price, Daniel J. 2012 Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics 231 (3), 759 – 794. [39] Saurel, R. & Abgrall, R. 1999 A simple method for compressible multifluid flows. SIAM Journal on Scientific Computing 21 (3), 1115–1145. [40] Scardovelli, R. & Zaleski, S. 1999 Direct numerical simulation of free-surface and interfacial flow. Ann. Rev. Fluid Mech. 31, 567–603. 25

594 595

596 597 598

599 600

601 602

[41] Schwaiger, H. F. 2008 An implicit corrected SPH formulation for thermal diffusion with linear free surface boundary conditions. Int. J. Numer. Meth. Engng. 75, 647–671. [42] Shadloo, M. S. & Yildiz, M. 2011 Numerical modeling of kelvin–helmholtz instability using smoothed particle hydrodynamics. International Journal for Numerical Methods in Engineering 87 (10), 988–1006. [43] Shao, S. 2012 Incompressible smoothed particle hydrodynamics simulation of multifluid flows. Int. J. Numer. Meth. Fluids 69, 1715–1735. [44] Shepard, D. 1968 A two-dimensional interpolation function for irregularly spaced data. Proceedings of the 23rd ACM national conference, New York, USA pp. 517–524.

606

[45] Skillen, A., Lind, S., Stansby, P.K. & Rogers, B.D. 2013 Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised fickian smoothing applied to body-water slam and efficient wave-body interaction. Computer Methods in Applied Mechanics and Engineering 265, 163–173.

607

[46] Tait, P. G. 1888 Physics and chemistry of the voyage of h.m.s. challenger. HMSO, London 2 (4).

608

[47] Tong, M. & Browne, D. J. 2014 An incompressible multi-phase smoothed particle hydrodynamics (SPH) method for modelling thermocapillary flow. International Journal of Heat and Mass Transfer 73 (0), 284 – 292.

603 604 605

609 610

611 612

613 614

615 616

617 618 619

620 621 622

[48] Xu, R. 2010 An improved incompressible smoothed particle hydrodynamics method and its application in free-surface simulations. PhD thesis, University of Manchester. [49] Xu, R., Stansby, P. K. & Laurence., D. 2009 Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. J. Comput. Phys. 228, 6703–6725. [50] Yokoi, K. 2007 Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm. J. Comput. Phys. 226, 1985–2002. [51] Zainali, A., Tofighi, N., Shadloo, M.S. & Yildiz, M. 2013 Numerical investigation of newtonian and non-newtonian multiphase flows using ISPH method. Computer Methods in Applied Mechanics and Engineering 254, 99–113. [52] Zhou, Z.Q., Kat, J.O. De & Buchner, B. 1999 A nonlinear 3-D approach to simulate green water dynamics on deck. In Proc. 7th Int. Conf. Num. Ship Hydrod. (ed. J. Piquet), , vol. 15, pp. 5.1–1. Nantes.

26