Resolving mixing in SPH

Resolving mixing in Smoothed Particle Hydrodynamics

J. I. Read, T. Hayfield and O. Agertz
Institute of Theoretical Physics, University of Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland
Department of Physics and Astronomy, University of Leicester, University Road, LE1 7RH Leicester, UK
Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 16, CH-8093 Zürich, Switzerland
Abstract

Standard formulations of smoothed particle hydrodynamics (SPH) are unable to resolve mixing at fluid boundaries. We use an error and stability analysis of the generalised SPH equations of motion to prove that this is due to two distinct problems. The first is a leading order error in the momentum equation. This should decrease with increasing neighbour number, but does not because numerical instabilities cause the kernel to be irregularly sampled. We identify two important instabilities: the clumping instability and the banding instability, and we show that both are cured by a suitable choice of kernel. The second problem is the local mixing instability (LMI). This occurs as particles attempt to mix on the kernel scale, but are unable to due to entropy conservation. The result is a pressure discontinuity at boundaries that pushes fluids of different entropy apart. We cure the LMI by using a weighted density estimate that ensures that pressures are single valued throughout the flow. This also gives a better volume estimate for the particles, reducing errors in the continuity and momentum equations. We demonstrate mixing in our new Optimised Smoothed Particle Hydrodynamics (OSPH) scheme using a Kelvin Helmholtz instability (KHI) test with density contrast 1:2, and the ‘blob test’ – a 1:10 density ratio gas sphere in a wind tunnel – finding excellent agreement between OSPH and Eulerian codes.

keywords:
Multiphase Smoothed Particle Hydrodynamics, Numerical methods, Monte-Carlo methods

1 Introduction

Smoothed Particle Hydrodynamics (SPH) was first introduced as a tool for studying stellar structure (Gingold & Monaghan 1977; Lucy 1977), but has since found wide application in all areas of theoretical astrophysics (Monaghan, 1992), in engineering (Libersky et al., 1993), and beyond (e.g. Hieber & Koumoutsakos 2008).

Although there are many varieties of SPH, the central idea is to represent a fluid by discrete particles that move with the flow (Monaghan 1992; Price 2005). Typically these particles represent the fluid exactly, though in some variants the fluid is advected on top of the particles (Dilts 1999; Maron & Howes 2003). The key advantages over Eulerian schemes111This does not apply to Lagrangian moving mesh schemes that are Galilean invariant (Springel, 2009). are its Lagrangian nature that makes it Galilean invariant, and its particle nature that makes it easy to couple to the fast multipole method for gravity that scales as (Dehnen 2000; Greengard & Rokhlin 1987). However, SPH has problems correctly integrating fluid instabilities and mixing at boundaries (Morris 1996b; Dilts 1999; Ritchie & Thomas 2001; Marri & White 2003; Agertz et al. 2007). Several different reasons have been suggested for this in the literature so far. Morris (1996b) and Dilts (1999) argue that the problem owes to errors in the SPH gradients that do not show good convergence for irregular particle distributions. Price (2007) argue that the problem owes to the fact that entropies are discontinuous at boundaries, while the densities are smooth. This gives spurious pressure blips at boundaries that drive fluids of different entropy apart. They find that adding thermal conductivity at boundaries to smooth the entropies gives improved mixing in SPH. Wadsley et al. (2008) make a similar argument, phrasing the problem in terms of an inability for SPH particles to mix and generate entropy on the kernel scale. They find that adding a heat diffusion term to model subgrid turbulence gives improved mixing in SPH. Finally, Ritchie & Thomas (2001) suggest that the problem lies in the SPH density estimate. They introduce a new temperature weighted density estimate that is designed to give smoother pressures at flow boundaries, thus combating the spurious boundary pressure blip.

In this paper, we perform an error and stability analysis of SPH in its most general form to understand why mixing fails. In doing this, we show that all of the above authors correctly identified one of two distinct problems with mixing in SPH. The first is an error in the momentum equation identified by Morris (1996b) and Dilts (1999). The second relates to entropy conservation on the kernel scale, as addressed directly by Price (2007) and Wadsley et al. (2008), and indirectly by Ritchie & Thomas (2001). Having identified the problem, we present a new method – Optimised Smoothed Particle Hydrodynamics (OSPH) – that, given sufficient resolution, correctly resolves multiphase fluid flow.

This paper is organised as follows. In §2 and §3, we briefly review standard SPH schemes and introduce our new OSPH scheme. We show that there are two distinct problems with mixing in SPH: the ‘ error’ in the momentum equation, and the ‘local mixing instability’ (LMI), and we show how both can be cured. In §4, we present our implementation of OSPH in the GASOLINE code (Wadsley et al., 2004). In §5, we use a Kelvin Helmholtz instability (KHI) test with density contrast 1:2 and 1:8 to demonstrate mixing in OSPH. We show the effect of turning on each of the OSPH improvements one at a time, arriving at a solution that is in excellent agreement with the Eulerian code RAMSES (Teyssier, 2002). In §6, we use the standard Sod shock tube test to demonstrate that OSPH can successfully model shocks. In §7, we revisit the ‘blob test’ introduced in Agertz et al. (2007), finding excellent agreement between OSPH and the Eulerian code FLASH (Fryxell et al., 2000). Finally, in §8 we present our conclusions.

2 Smoothed Particle Hydrodynamics

In SPH, the fluid is represented by discrete particles that move with the flow. The density of each particle is estimated by a weighted sum over its neighbours:

 ρi=N∑jmjW(|rij|,hi) (1)

where and are the smoothing length and mass of particle and , respectively; we define and similarly for other vectors; and is a symmetric kernel that obeys the normalisation condition:

 ∫VW(|r−r′|,h)d3r′=1 (2)

and the property:

 limh→0W(|r−r′|,h)=δ(|r−r′|) (3)

In the limit (and using ) equation 1 recovers the continuum flow density.

The equations of motion for SPH are then derived by discretising the Euler equations – the continuity, momentum and energy equations:

 dρdt=−ρ∇⋅v (4)
 dvdt=−∇Pρ (5)
 dudt=−Pρ∇⋅v (6)

where and are the density, velocity and internal energy per unit mass of the flow, respectively.

The Euler equations can be derived from the Lagrangian for hydrodynamics (e.g. Bennett 2006):

 L=∫(12ρv2−ρu)dV (7)

and in many modern derivations of the equations of motion for SPH, equation 7 is discretised, rather than equations 4-6.

Replacing the volume element with the volume per SPH particle , we obtain (Price 2005):

 L=∑jmj(12v2j−uj) (8)

and the standard SPH equations of motion then follow from the Euler-Lagrange equations:

 dρidt=N∑jmjvij⋅∇iWij (9)
 dvidt=−N∑jmj[Piρ2i+Pjρ2j]∇iWij (10)
 duidt=Piρ2iN∑jmjvij⋅∇iWij (11)

where .

Note that equation 9 is automatically satisfied by time derivative of the SPH density estimate (equation 1). For this reason, equation 1 is often referred to as the integral form of the continuity equation.

The above system of equations are closed by the equation of state:

 Pi=ui(γ−1)ρi (12)

This standard approach to deriving the SPH equations of motion has the advantage that the resulting equations are coherent222Also called consistent (Oger et al., 2007). by construction – that is they are consistent with a Lagrangian. This gives very good conservation properties for the flow. It is also straightforward to calculate the necessary correction terms that arise if the smoothing lengths are a function of space and time (see e.g. Nelson & Papaloizou 1994; Price 2005). We do not include these correction terms in this paper.

However, this standard derivation leads to a scheme that cannot correctly model fluid mixing processes (see §1), which motivates us to move to a more general derivation. Discretising each of the Euler equations separately leads to a free function for each equation: , and , as well as a different smoothing kernel for each. This is the approach we take next in §3. In §3.3 and §3.4, we will then use an error and stability analysis of these more general equations of motion to constrain the new functions and and our new kernels. By choosing these new free functions and kernels such that they minimise the integration error, we will arrive at a new scheme that can, with sufficient resolution, correctly resolve multiphase fluid flow.

3 OPTIMISED Smoothed Particle Hydrodynamics

In the previous section, we presented a standard derivation of the SPH equations of motion. However, this standard derivation leads to a scheme that cannot correctly model fluid mixing processes (see §1). In this section, we move to a more general derivation of the SPH equations of motion. We show that, in general, we have a free function for each of the Euler equations: , and , as well as a different smoothing kernel for each. There is also a freedom in the energy equation in the choice of integration variable (energy or entropy; §3.2). In §3.3 and §3.4, we will then use an error and stability analysis of these more general equations of motion to constrain the new functions and and our new kernels. By choosing these new free functions and kernels such that they minimise the integration error, we will arrive at a new scheme that can, with sufficient resolution, correctly resolve multiphase fluid flow.

3.1 A general derivation of SPH

In general, we have some freedom in how we discretise the Euler equations (equations 4-6) to obtain the equations of motion for SPH (see e.g. Monaghan, 1992; Price, 2005; Rosswog, 2009). The gradients in the Euler equations can be expanded to include a new free function for each equation: , and :

 dρdt=η[v⋅∇(ρη)−∇⋅(ρvη)] (13)
 dvdt=−[Pϕρ2∇(ρϕ)+1ϕ∇(Pϕρ)] (14)
 dudt=Pρ2ζ[v⋅∇(ρζ)−∇⋅(ρvζ)] (15)

In the continuum form, above, and cancel. But in the discrete SPH form, they remain giving a useful additional freedom (Price, 2005):

 dρidt=N∑jmjηiηjvij⋅¯¯¯¯¯Hijrij (16)
 dvidt=−N∑jmj[Piρ2iϕiϕj+Pjρ2jϕjϕi]¯¯¯¯¯Kijrij (17)
 duidt=Piρ2iN∑jmjζiζjvij⋅¯¯¯¯Lijrij (18)

where , and are symmetrised smoothing kernels – one for each Euler equation. Standard SPH (SPH from here on) is a special case of the above with and .

Equation 16 casts the continuity equation in differential form. This is problematic since, in this case, the particles no longer represent the fluid exactly. Instead they represent a moving mesh on which the Euler equations are solved. This leads to the danger that high density regions will contain few particles leading to large errors (Maron & Howes, 2003). For this reason, we use instead a generalised integral form for the continuity equation:

 ρi=N∑jmjηiηj¯¯¯¯¯¯Wij (19)

which, taking the time derivative, gives:

 dρidt=N∑jmjηiηjvij⋅∇i¯¯¯¯¯¯Wij+ϵ (20)

where:

 ϵ=N∑jmj(˙ηiηi−˙ηjηj)ηiηj¯¯¯¯¯¯Wij (21)

and .

This reduces to the continuity equation (equation 16) under the kernel constraint: , and for . The latter can be satisfied by construction if (as is the case for SPH), or if . However, in the continuum limit (, ), and so will vanish with increasing resolution. For this reason, equation 19 gives a valid approximation to the continuity equation for any choice of , with simply contributing an additional error term.

3.2 Energy versus entropy forms of SPH

A final freedom in the equations motion for SPH comes from the energy equation. Equation 18 is the standard energy form of SPH, but there is also an entropy form (Goodman & Hernquist 1991; Springel & Hernquist 2002). Instead of the internal energy, , we evolve a function – the entropy function – that is a monotonic function of the entropy defined by the equation of state:

 Pi=Ai(s)ργi (22)

Away from shocks and in the absence of thermal sources or sinks, is a constant of motion. Thus, taking the time derivative of equation 22 and substituting for equation 12, we recover:

 duidt=Piρ2idρidt (23)

by construction. Schemes that obey equation 23 are called thermodynamically consistent.

In practice, we find – for the tests presented in this paper – that the energy and entropy forms of SPH give near-identical results, provided that equation 23 is satisfied (for adiabatic flow). We use the thermodynamically consistent energy form throughout this paper. This gives us the constraints: and , which we apply from here on. We also use , as in standard SPH. This is not a formal requirement, but ensures that coherence is recovered in the limit of constant density.

3.3 Errors: choosing the free functions

In this section, we perform an error analysis of the generalised equations for SPH (equations 17, 18 and 19) derived in §3.1. We will then choose our free functions , and so that these errors are minimised.

3.3.1 Error analysis

We assume that the pressure and velocity of the flow are smooth. In this case, we can Taylor expand to give:

 Pj≃Pi+hxij⋅∇iPi+O(h2) (24)

and

 vj≃vi+h(xij⋅∇i)vi+O(h2) (25)

where , and we have assumed a constant smoothing length .

Substituting equations 24 and 25 into the continuity and momentum equations gives:

 dρidt≃−ρi(Ri∇i)⋅vi+ϵ+O(h) (26)

and

 (27)

where is a dimensionless error vector given by:

 (28)

and and are dimensionless error matrices given by:

 Ri=N∑jmjρjfijSij;Vi=N∑jmjρjg−1ijSij (29)

with:

 Sij=1x∂¯¯¯¯¯¯Wij∂x⎛⎜ ⎜ ⎜⎝x2ijxijyijxijzijyijxijy2ijyijzijzijxijzijyijz2ij⎞⎟ ⎟ ⎟⎠ (30)

where ; ; ; ; and .

The accuracy of the continuity equation (26) is given by the extent to which (see equation 21) and , the identity matrix. The accuracy of the momentum equation (27) is given by the extent to which and . (The energy equation behaves similarly to the continuity equation with .)

3.3.2 Minimising errors: the continuity equation

Let us consider how accurately equation 26 approximates its Euler equation equivalent (equation 4). First, consider standard SPH where and by construction. Typically in the literature, the error is calculated only in the continuum limit (; ; see e.g. Price (2005)). In this case, the sums become integrals, and (using ) we obtain terms like:

 limN→∞R33(x)=∫Vd3x′f(x,x′)(z−z′)2|x−x′|∂W∂x (31)

and

 limN→∞R12(x)=∫Vd3x′f(x,x′)(x−x′)(y−y′)|x−x′|∂W∂x (32)

where the notation refers to element in the matrix .

If we assume smooth densities, then we can Taylor expand also to obtain:

 f=ρ(r′)ρ(r)≃1+h(x−x′)ρ⋅∇ρ+O(h2) (33)

and we see that, by symmetry of , to . In fact, Taylor expanding to an order higher than above, it is straightforward to show that the whole continuity equation is accurate to in the limit (see e.g. Price 2005). A similar argument applies to the other SPH equations of motion and leads to the conclusion that SPH is accurate to . However – and this is a key point – this formal calculation is only valid for smoothly distributed particles in the limit . In practical situations, where we have a finite number of particles within the kernel and these are not perfectly smoothly distributed, the leading order errors in the continuity equation appear at and are contained within the matrix . We will quote orders of error from here on in this finite particle limit.

We can think of each term of as a finite sum approximation to a dimensionless integral that should be either 0 (for the off diagonal terms), or 1 (for the diagonal terms). For smooth particle distributions, this approximation is a good one since , while gives a good estimate of the volume of each particle within the kernel. However, if the particles are distributed irregularly on the kernel scale – for example at a sharp density step – then can grow arbitrarily large, while becomes a poor volume estimate. We will demonstrate this in §5.

We can improve matters by choosing , which fixes always. However, the integral form of the continuity equation then becomes:

 ρi=N∑jρiρjmj¯¯¯¯¯¯Wij (34)

which must be solved iteratively and is not guaranteed to converge. Worse still, is now no longer zero and contributes an additional error.

Ritchie & Thomas (2001) present an interesting solution to this dilemma. If the pressures are approximately constant across the kernel () then, for the energy form of SPH (see equation 12 and §3.2), and equation 34 is well approximated by the integral continuity equation:

 ρi=N∑jujuimj¯¯¯¯¯¯Wij (35)

This can be solved without the need for iteration.

The above suggests that we use . Thermodynamic consistency then requires that we set (see §3.2).

There may be some advantage, however, to using the entropy form of SPH. In this case, the equation of state is given by equation 22. For approximately constant pressure across the kernel, we now have that , and the integral continuity equation becomes:

 ρi=N∑j(AjAi)1γmj¯¯¯¯¯¯Wij (36)

This has the advantage that, in the absence of shocks or thermal sources/sinks, and so the error term by construction (see equation 21). In practice, however, we find no appreciable difference between the energy and entropy forms of SPH for the tests presented in this paper. This suggests that is not a significant source of error.

Equations 35 and 36 retain the desirable integral form for the density, while giving significantly improved error properties. They also have a second important advantage that we discuss in §3.5. We refer to equation 35 as the ‘RT’ density estimator for the energy form of SPH; and equation 36 as the RT density estimator in the entropy form.

3.3.3 Minimising errors: the momentum equation

The momentum equation (equation 27) is more problematic than the continuity equation. Its accuracy is governed not only by the extent to which , but primarily by the leading term that should vanish.

First, consider the situation in standard SPH where and . As for the continuity equation, in the continuum limit (), since is antisymmetric. However, this analysis is only relevant if the particles are smoothly distributed on the kernel scale. For irregularly distributed particles, can grow arbitrarily large, while is not guaranteed to be a good volume estimate. In such situations, contributes a significant error. Worse still, moving to higher resolution is not guaranteed to help. In order for the SPH integration to converge as , we require that shrinks faster than . This requires some care in making sure that does not shrink too fast as the number of particles is increased.

A density step is an extreme example of an irregular particle distribution, and this suggests that the error is at least in part responsible for SPH’s failure to correctly model mixing processes between different fluid phases. We demonstrate this in §5.

There are three key problems with ensuring that will shrink with increasing resolution. The first is the function . In SPH, this is the ratio which is large when there are large density gradients. We can significantly improve on this we choose our free function . In this case, we have by construction, and no longer contributes to the error even for large density gradients across the kernel333It is interesting to note that other work in the literature has also found that is the preferred choice if density gradients are large (Oger et al. 2007; Price 2005; Ritchie & Thomas 2001; Dilts 1999). But we could not find a detailed proof similar to that presented here. Interestingly, Marri & White (2003) find empirically that gives the best performance for multiphase flow in their tests. Our analytic results here suggest that this is not the optimal choice, though perhaps the inclusion of cooling and/or other physics makes a difference.. The second problem relates to kernel scale smoothness. If particles clump or band on the kernel scale, then we will have poor kernel sampling and will not approach its integral limit even at very high resolution. Ensuring that this does not happen means ensuring that our OSPH scheme is stable to perturbations. We discuss this next in §3.4. The third and final problem is the volume estimate of each particle . This will be poor if the particles are irregularly distributed on the kernel scale (for example at a density step) leading to a large error. We discuss this further in §3.5.

The choices ; and the kernel constraints are the first important ingredients in our OSPH scheme. These choices mean that we are no longer coherent444Note that it is possible to construct pseudo-coherent versions of OSPH using for the energy form, or for the entropy form. Introducing ‘grad ’ terms as in Nelson & Papaloizou (1994), such schemes can then be made to conserve energy exactly in the limit of constant timesteps. However, they are only truly coherent up to the approximation that the error in the continuity equation is small (see equation 21). Nonetheless, it would be interesting to explore such schemes in future work., but this only introduces tolerable errors in the energy conservation (Hernquist & Katz, 1989).

3.4 Stability: the choice of kernel function

In §3.3, we used an error analysis of the generalised SPH equations of motion to show that the dominant source of error in SPH is in the momentum equation – the error. We showed that choosing the free functions ; and the kernel constraints should minimise both this error and errors in the continuity equation, and we called these choices OSPH.

In OSPH, provided the particles are regularly distributed on the kernel scale, we can make arbitrarily small simply by increasing the neighbour number. However, if the particles are irregularly distributed, can shrink very slowly with increasing resolution. In this section, we show that for large neighbour number the cubic spline kernel typically used in SPH calculations is unstable to both particle clumping (§3.4.1) and particle banding (3.4.2), and we derive a new class of kernels that are stable to both even for large neighbour number. In §5, we show that these new kernels give significantly improved performance at no additional computational cost.

3.4.1 The clumping instability

The clumping instability555Also called the tensile instability. can be derived from a linear 3D stability analysis of the OSPH equations of motion. Following Morris 1996a and Morris 1996b, we imagine a lattice of equal masses of equal separation with initial density and pressure . We perturb these with a linear wave of the form:

 xi=x0,i+aexp[i(k⋅x0,i−ωt)] (37)
 ρi=ρ0+Dexp[i(k⋅x0,i−ωt)] (38)
 Pi=P0+c2sDexp[i(k⋅x0,i−ωt)] (39)

and similar for particle , where is the sound speed assuming an adiabatic equation of state ( gives an isothermal equation of state), is the amplitude of the perturbation and is the wave vector.

To simplify the analysis, we assume that we have a lattice symmetry such that for every displacement vector to a neighbour, there is also one at . Then, plugging equations 37, 38 and 39 in to 17, discarding terms higher than first order, and connecting to through the continuity equation666We use here the full OSPH continuity equation in the entropy form (equation 16 with ). However, for plane waves on a constant density lattice and so this is identical to the SPH continuity equation with . (), we obtain the 3D OSPH dispersion relation777This is actually identical to the SPH dispersion relation derived under the same assumptions in 3D (Morris, 1996b).:

 ω2a = [2mP0ρ20∑jH(¯¯¯¯¯¯W0,ij)(1−cosk⋅x0,ij)+ (40) (γ−2)m2P0ρ30(qi∧qi)]⋅a

where is the outer product of , is the Hessian888Recall that the outer product of two vectors is a matrix, while the Hessian is a square matrix of second order partial derivatives. of :

 Haa=∂2¯¯¯¯¯¯W∂xa2=d2¯¯¯¯¯¯W(r)dr2xa2r2+d¯¯¯¯¯¯W(r)dr1r(1−xa2r2) (41)
 (42)

and is given by:

 qi=∑jsink⋅x0,ij∇i¯¯¯¯¯¯W0,ij (43)

Our scheme is stable if . It is also desirable for the numerical sound speed to equal the true sound speed: .

In SPH it is typical to use the cubic spline (CS) kernel given by:

 W=8πh3⎧⎪ ⎪⎨⎪ ⎪⎩1−6x2+6x30

where is the distance from the centre of the kernel in units of the smoothing length.

In Figure 1, we show contours of as a function of wavenumber and the smoothing length in units of the inter-particle spacing , for the CS kernel. We assume an adiabatic equation of state with . From left to right the plots show and . The three rows show the longitudinal wave and the two transverse waves for these orientations.

From Figure 1, it is clear that the CS kernel in 3D is unstable to longitudinal waves for , and very unstable to transverse waves. The unstable longitudinal waves drive the clumping, or tensile instability that causes particles to clump on the kernel scale (Schuessler & Schmitt 1981; Thomas & Couchman 1992; Herant 1994; Morris 1996a; Monaghan 2000).

The clumping instability is a problem because it means that increasing the neighbour number will not give improved sampling of the kernel, and the error will remain large. However, the situation is dramatically improved if we add a constant central core to the kernel gradient . This gives a constant force term at the centre of the Kernel that physically prevents clumping. We choose a kernel that is maximally similar to the CS kernel, while obeying , where is the core size. This leads us to the Core Triangle (CT) kernel:

 W=Nh3⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(−12α+18α2)x+β0

where , , and the core size is fixed at by the requirement that be continuous.

Figure 2 shows stability plots for the CT kernel. The CT kernel has greatly improved stability for the longitudinal waves (top row) compared to the CS kernel and should give significantly improved performance for large neighbour number. We demonstrate this in §5.

Note that for all of the kernels we use in this paper, we consistently apply the kernel for the density estimate and its gradient for the energy and momentum equations. For small neighbour number, the central triangle in the CT kernel will degrade the quality of the density estimate. However, in this paper we typically use large neighbour numbers (). In this case, very few particles sample the inner regions of the kernel and the bias introduced in the density is negligible. (The quality of the density estimate in OSPH can be seen in the Sod shock tube test in §6.) We found in tests that retaining the CS kernel just for the density estimate gives near-identical results.

3.4.2 The banding instability

The clumping instability is a result of unstable longitudinal waves. A related instability – the banding instability – is a result of unstable transverse waves. For both the CT and CS kernels, there are broad bands of instability to transverse waves (see Figures 1 and 2). If the neighbour number is carefully chosen to lie in a stable region, banding will not occur. However, banding can still be excited at boundaries if changes there, moving into an unstable region.

For both the CS and CT kernels, it is difficult to find a suitable neighbour number for which the kernel is stable to all transverse and longitudinal modes. This suggests hunting for an even more stable kernel. A full search is beyond the scope of this paper. Here, we present a simple class of kernels that improve stability by moving to higher order (Morris, 1996b). Following Price (2005), we generalise our CT kernel to order to obtain the following class of kernels that we call the High Order Core-Triangle (HOCT) kernels:

 W=Nh3⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Px+Q0

where:

 A=1−β2γnk−3(γ2−β2) (47)
 B=−1+Aγnk−1βnk−1 (48)
 P=−nk(1−α)nk−1−nkA(γ−α)nk−1−nkB(β−α)nk−1 (49)
 Q=(1−α)nk+A(γ−α)nk+B(β−α)nk−Pα (50)

and and are calculated numerically for a given choice of . Continuity requires that solves the equation:

 0=(1−α)nk−2+A(γ−α)nk−2+B(β−α)nk−2 (51)

where and are free parameters. In this paper we choose . Other choices, and indeed other high-order kernels, may give better results than those presented here. We tabulate values for and as a function in Table LABEL:tab:hoctkern. Notice that the core size decreases with .

Stability plots for the HOCT4 kernel (with ) are given in Figure 3. Notice the improvement over the CT kernel, particularly for the transverse waves. There are two bands where the kernel is fully stable to both longitudinal and transverse waves on a lattice: 96 neighbours and 442 neighbours, corresponding to and , respectively. We use the latter choice since this also gives very low . The CT kernel also has a stability band for , but this is narrower than for the HOCT4 kernel, while the HOCT4 kernel with this many neighbours gives better spatial resolution. (It is important to realise that the smoothing length for different kernels takes on a different meaning in terms of spatial resolution. We suggest a resolution criteria based on the numerical sound speed versus the true sound speed for longitudinal waves: . Spatial scales are well resolved if . By this definition, our choice of 442 neighbours () for the HOCT4 kernel gives a very similar spatial resolution to 128 neighbours () for the CT kernel.)

Note that our stability analysis only applies for particles arranged on a lattice. Hexagonal close-packed particles, randomly arranged particles, and indeed boundaries may have different preferred stability regions. A full analysis is beyond the scope of this present work.

The banding instability is not as problematic as the clumping instability for the tests we present in §5 and §7. Unlike the clumping instability, it does not seem to (directly) play a major role in preventing mixing from occurring in SPH (see §5.2.2).

3.5 The local mixing instability & RT densities

Our error analysis in §3.3 missed one very important error term. This is because the Taylor expansion assumed that both the pressures and velocities in the flow are smooth. Unfortunately, in SPH at sharp boundary this is not the case. The reason for this is easiest to understand using the entropy form of SPH, as follows (similar arguments apply also for the energy form).

Imagine a density step of ratio initially in pressure equilibrium, such that the entropy function (equation 22) is given by . Now imagine that we perturb the boundary very slightly by pushing a low density particle towards it. The particle’s entropy is conserved, but its density increases very rapidly proportional to . This leads to an increase in pressure: , where is some constant that depends on the perturbation size and the kernel. On the other side of the boundary, if we push a high density particle towards the low density region, however, its density will rapidly decrease giving a decrease in pressure: . This drives us towards a pressure discontinuity at the boundary which drives an associated error in the momentum equation. It can be thought of as a fundamental result of particles trying to mix on the kernel scale, but being unable to as a result of entropy conservation. We call this the local mixing instability (LMI).

Although not phrased in terms of the LMI, the LMI is a recognised problem in the literature and there are essentially two classes of solution. We can generate entropy at the boundary to give smooth entropies and therefore smooth pressures, as in (Wadsley et al., 2008) and Price (2007); or we can try to obtain sharper densities that are consistent with the discrete entropies. This is the approach adopted by Ritchie & Thomas (2001), and the approach we take in this paper. The key advantage of sharpening the densities is that we do not need to specify a subgrid mixing model.

The sharper densities we require are exactly what we get from the density estimate given in equation 36, and originally proposed by Ritchie & Thomas (2001) – the ‘RT’ density estimate. Consider the perturbation discussed above, but now using the RT density estimate999We use here the entropy form given in equation 36 since we use the entropy form of SPH in this analysis. If instead, we use the energy form of SPH then we should use instead equation 35.. A low density particle which has half of its kernel in the high density phase (an extreme example), will have a density:

 ρlow=Nl∑jmj¯¯¯¯¯¯Wij+Nr∑j(AjAi)1γmj¯¯¯¯¯¯Wij (52)

where is the number of particles in the low density region, is the number in the high density region, and we have used the fact that the ratio for the low density region.

If the simulation is adiabatic and started in pressure equilibrium, then for the high density region , and since the high density particles sample the kernel times more often than the low density particles, we recover:

 ρlow=N∑jmj¯¯¯¯¯¯Wij (53)

which is identical to a particle in the low density region. A similar derivation applies to a high density particle at the boundary. Thus, the RT density estimate ensures that densities remain sharp. It is straightforward to show that it also ensures the pressures are single valued throughout the flow. Substituting the RT density estimator (equation 36) into equation 22, we obtain:

 Pi=[N∑jmjA1γj¯¯¯¯¯¯Wij]γ (54)

Notice that the entropy function now appears inside the sum, whereas in standard SPH it would appear as outside of the sum. This difference ensures that the will be everywhere single valued throughout the flow – even at boundaries. A similar derivation can be made for the energy form of SPH, in which case we should use the density estimate given in equation 35.

The RT density estimate is robust to particle mixing on the kernel scale and should lead to a dramatically reduced LMI. We demonstrate this in §5. Furthermore, the RT density estimate ensures that our error analysis in §3.3 is valid by construction since it ensures smooth pressures (recall that we assumed that both the pressures and the velocities were smooth, but not the densities). And, since the RT density estimate leads to sharper densities, it gives improved volume estimates for the particles. This suggests that we can expect the RT density estimate to reduce at boundaries. We demonstrate this also in §5.

Note that the RT density estimate is chosen to ensure single valued pressures throughout the flow. However, when extracting results from a simulation, it is the positions of the particles themselves that describe the state of the fluid. This suggests using the density estimate in equation 1 for calculating the observable flow density, rather than the RT density estimate. This is the approach we adopt in this paper, though the difference is negligible.

4 Implementation

We implemented OSPH in the GASOLINE code (Wadsley et al., 2004), a parallel implementation of TreeSPH that uses a fixed number of smoothing neighbours101010Allowing for varying neighbour number is needlessly dissipative (Nelson & Papaloizou 1994; Attwood et al. 2007)., and a standard prescription for the artificial viscosity as in Gingold & Monaghan (1983) with , , controlled with a Balsara switch (Balsara, 1989). We used variable timesteps controlled by the Courant time with a Courant factor of .

The improved stability and error properties of OSPH motivate a full re-examination of the standard SPH artificial viscosity. This is beyond the scope of this present work. However, we note that the improved stability in OSPH means that particles better follow characteristics of the flow, while the gradients in the Balsara switch will be less noisy. Both of these effects should act to decrease the viscosity in regions of steady flow. (Note that all numerical schemes carry numerical viscosity, whether it is manifested through limited resolution or artificial shock-capturing viscosity. Indeed, these viscous terms are vital for successfully modelling shocks.) In §5, §6 and §7, we show that our OSPH results agree very well with analytic expectations, and with the results from Eulerian codes. This suggests that the viscosity prescription in OSPH is not a significant source of error. Certainly it is not responsible for SPH’s inability to model mixing processes.

5 The Kelvin-Helmholtz instability

In this section, we use a 1:2 and 1:8 density ratio shearing fluid simulation to test mixing in OSPH. We use the naming convention XSPH-K-N, where X denotes the variety of SPH, K the choice of kernel, and N the neighbour number (see Table 2).

5.1 Numerical set-up

A Kelvin-Helmholtz instability (KHI) occurs when two shearing fluids are subjected to an infinitesimal perturbation at the boundary layer. The result of the perturbation is a linearly growing phase in which the layers start to interpenetrate each other, progressively developing into a vortex in the non-linear phase that mixes the two fluid layers. The growth-rate of the instability is in general a complicated function of the shear velocity, fluid densities, compressibility, interface thickness, gravity, viscosity, surface tension, magnetic field strength etc. In this test, we are only interested in the behaviour of inviscid, incompressible (i.e. with bulk motions very much less than the sound speed) perfect fluids neglecting gravity. In this case, the linear growth rate of the KHI is (Chandrasekhar 1961):

 w=k(ρ1ρ2)1/2v(ρ1+ρ2), (55)

where is the wavenumber of the instability, and are the densities of the respective layers and is the relative shear velocity. The characteristic growth time for the KHI is then:

 τKH≡2πw=(ρ1+ρ2)λ(ρ1ρ2)1/2v. (56)

This is a particularly challenging test for SPH/OSPH because the velocity due to particle noise can approach the sound speed which can wash out the physical velocity perturbation relevant for this test.

For the simulations, we set up the problem in three dimensions using a periodic thin slab defined by , and . The domain satisfied:

 ρ,T,vx={ρ1,T1,v1|y|<0.25ρ2,T2,v2|y|>0.25 (57)

The density and temperature ratio were , ensuring that the whole system was pressure equilibrium. The two layers were given constant and opposing shearing velocities, with the low density layer moving at a Mach number and the dense layer moving at . The density ratios considered in this work are small which assures a subsonic regime where the growth of instabilities can be treated using equation 56 (Vietri et al., 1997).

To trigger instabilities, velocity perturbations were imposed on the two boundaries of the form:

 vy = δvy[sin(2π(x+λ/2)/λ)exp(−(10(y−0.25))2) (58) −sin(2πx/λ)exp(−(10(y+0.25))2)]

where the perturbation velocity and is the wavelength of the mode.

Equal mass particles were placed in lattice configurations to satisfy the setup described above. To satisfy pressure equilibrium everywhere, in TSPH the temperatures were adjusted at boundaries to be coherent with the smoothed density step measured by equation 1. This was not done for the OSPH simulations since these sharpen the densities using the discrete initial temperatures.

The low density region was set up using 256 particles in the -direction and the appropriate number of particles in the other dimensions to satisfy a fixed inter-particle distance. The high density region was created in the same way with 320 particles in the -direction. We adopted a periodic simulation domain.

The RAMSES simulation used the same numerical set-up as described above, but in 2D rather than in a thin slab. We performed the simulation using the LLF Riemann solver (Toro, 1999) on a fixed Cartesian grid. The LLF solver is rather diffusive and is used in order to suppress the growth of undesirable small scale KHIs arising from grid irregularities.

We note that all numerical schemes carry numerical viscosity, whether it is manifested through limited resolution or artificial shock-capturing viscosity. A detailed study of this effect on the KHI and the relation to physical viscosity is beyond the scope of this paper.

5.2 Results

Figure 4 shows our results for the KHI test (density ratio ) at modelled with SPH, TSPH and OSPH, using three different kernels: CS, CT and HOCT4, and different neighbour numbers as marked on each plot (see also Table 2). From left to right, the panels show, in a slice of width about the z-axis: density contours of the simulation box, a zoom in on the particle distribution around on of the rolls, the magnitude of the error (equation 28) as a function of , and the pressure as a function of in a slice of width about the x-axis.

5.2.1 The clumping instability

Using the standard CS kernel, SPH-CS-128 (top row, Figure 4) and TSPH-CS-128 (second row) give poor results that improve very slowly with increasing neighbour number. This can be seen both in the lack of strong evolution on the boundary, and in the large error, even for 128 neighbours. TSPH-CS-128 gives slightly better results than SPH-CS-128, showing the first beginnings of a KHI roll, but both are in poor agreement with the RAMSES results (bottom row).

The reason for the poor performance in both SPH-CS-128 and TSPH-CS-128 is the clumping instability (§3.4.1). Particles gather together on the kernel scale, giving poor kernel sampling, and poor associated error. This can be seen in the particle distribution for SPH-CS-128 and TSPH-CS-128 (second row, Figure 4) which show visible holes and over-densities in the particle distribution. Using instead the CT kernel introduced in §3.4.1, the results improve dramatically (third row, Figure 4). Now the errors reduce for increasing neighbour number (see Appendix A). With 128 neighbours, we successfully resolve a KH roll up to with the correct growth time.

It has been noted previously in the literature that putting a small core inside a cubic spline kernel suppresses the clumping instability (Thomas & Couchman 1992; Herant 1994), though its importance for modelling multiphase flow was not realised. Alternative fixes include adding an negative pressure term (Monaghan, 2000), which in tests we find works also. However, we prefer changing the kernel to introducing new forces since we may then still estimate our errors through .

5.2.2 The banding instability

In addition to the clumping instability, there is also an instability to transverse waves – the banding instability (§3.4.2). For the KHI tests we present here, the banding instability occurs only on the boundary and appears to be relatively benign. This is shown in Figure 5, that shows a zoom in on the boundary at for TSPH-CT-128, TSPH-HOCT4-442 and OSPH-HOCT4-442. The TSPH-CT-128 simulation has a kernel and neighbour number combination that are unstable to transverse waves (see Figure 2), and banding is clearly visible on the boundary. However, TSPH-HOCT4-442 should be stable to transverse waves, yet the banding persists. Only in our full scheme, OSPH-HOCT4-442, is the banding is gone.

To understand the above results, we ran an additional test that we omit for brevity – TSPH-HOCT4-96. This simulation showed little boundary evolution because the low neighbour number and associated large significantly damped the KHI. However, interestingly, there was no banding observed on the boundary (recall that 96 neighbours for the HOCT4 kernel should be stable to both transverse and longitudinal wave perturbations).

Taken together, our results suggest that the observed banding at the boundary is a result of a transverse wave instability driven by the local mixing instability (LMI; §3.5). Where there is little evolution at the boundary and the kernel is chosen to be stable to transverse waves, the banding disappears, as was the case for our extra TSPH-HOCT4-96 simulation. Where there is strong evolution at the boundary, as was the case for TSPH-HOCT4-442, the LMI drives banding irrespective of the choice of kernel. Only in our full scheme, OSPH-HOCT4-442, where the LMI is cured and the kernel is stable to transverse waves, is the banding cured.

5.2.3 The E0 error

Away from boundaries, the error in TSPH decreases with the neighbour number, as expected for smooth flow (see Appendix A). However, on the boundary the error grows by 2-3 orders of magnitude. Increasing the neighbour number does result in better long-term evolution, but the results improve very slowly. This is shown in Figure 6. Notice that TSPH-HOCT4-442 resolves two wraps of the KH roll at , whereas TSPH-CT-128 only manages one. However, even in TSPH-HOCT4-442, the long-term evolution eventually degrades. By , the results are ‘gloopy’, rather similar to simulations that explicitly model fluid surface tension (see e.g. Herrmann 2005).

The poor on the boundary is the result of a poor volume estimate for each particle (see §3.3). However is not solely responsible for the gloopy behaviour. There is a second problem – similar to a numerical surface tension term – that needs to be solved extra to minimising . This is the local mixing instability (LMI) error (§3.5).

5.2.4 The local mixing instability error

The right panels of Figure 4 show the pressure as a function of in a slice of width about the z-axis and width about the x-axis. In SPH and TSPH  there is a clear pressure discontinuity on the boundary. This is caused by the local mixing instability (LMI) discussed in §3.5.

Notice that the pressure blip is larger in TSPH than in SPH, yet the KHI roll progresses further in TSPH than in SPH. This apparent paradox is the result of the improved performance in TSPH. As the KHI roll progresses in TSPH, particles are pushed closer to the boundary making the LMI worse and increasing the pressure blip. In SPH, there is a larger gap at the boundary due to the larger surface tension error. This leads to less evolution and a smaller associated pressure blip. We will see a similar effect occurring in the blob test in §7.

As discussed in §3.5, the LMI should be cured by the RT density estimate (equation 35). This is shown in Figure 4, third row which shows the results for our full OSPH scheme. With the RT densities, the pressure at the boundary has a much smaller blip, while is reduced by over an order of magnitude. This latter effect occurs since the RT densities also give an improved volume estimate for each particle (see §3.3 and §3.5). The long term evolution is now in excellent agreement with the RAMSES results (compare Figure 4 third and bottom rows).

Although OSPH gives significantly improved results as compared with SPH, the scheme is numerically expensive. Simulations with larger density gradients require very high resolution. This is shown in Figure 6, second from bottom row. This shows the long term evolution of a KHI test with density ratio in OSPH. The solution should be similar to the simulation, but it is not. The ‘gloopy’ behaviour indicative of large surface tension errors has returned. Further increasing the neighbour numbers would reduce this problem, but at increased numerical cost. We will address this issue in future work (Hayfield & Read in prep.).

6 The Sod shock tube

Before we embark on the blob test in §7, it is worth checking that our new OSPH scheme can still correctly resolve shocks. To test this, we use a standard Sod shock tube test (Sod, 1978).

The Sod shock tube consists of a 1D tube on the interval with a discontinuous change in properties at designed to generate a shock. The left state is described by , , , and the right state by , , , where and are the density, pressure and velocity along the axis. We use an adiabatic equation of state with . The subsequent evolution of the problem has a self-similar analytic solution that has a number distinct features which quite generally test a code’s conservation properties, artificial viscosity, ability to handle nonlinear waves, and shock resolution.

Figure 7 shows the results for the Sod shock tube test at time in SPH (top) OSPH (bottom). Since we are primarily concerned with the 3D performance of the code, the test was performed in 3D on the union of a lattice on the left, with a lattice on the right, giving a 1D resolution of 450 points. We use 442 neighbours for this test in both SPH and OSPH to ensure that any difference is not simply due to improved kernel sampling in OSPH.

For SPH, the only strong disagreement with the analytic solution is in the pressures that have a blip at , and the temperatures that overshoot at . The former feature is due to the LMI (see §3.5 and §5.2.4). The latter feature is seen in all SPH Sod shock tube tests and results from the well-known ‘wall heating’ effect (Noh, 1987). This is an error due to the artificial viscosity prescription and is beyond the scope of this present work.

For OSPH, the results are even better than for SPH. The pressure blip is now gone, while the temperature overshoot at is reduced. Only the velocities appear to be worse, with some remaining dispersion at . This owes to the jump in density at this point, and the associated jump in . This gives a force error at the discontinuity which introduces some dispersion into the velocities. In SPH this cannot occur since the LMI causes a pressure blip at the boundary that prevents mixing. We will discuss this issue further in a forthcoming paper (Hayfield & Read in prep.).

7 The blob test

The KHI test presented in §5 is a worst-case scenario for OSPH, since it has a pure adiabatic sharp boundary. For many practical situations, boundaries will be less sharp, while physical entropy generation due to shocks and/or cooling will suppress the LMI. We give a practical example of this in this section using the blob test described in Agertz et al. (2007). A spherical cloud of gas of radius is placed in a wind tunnel with periodic boundary conditions. The ambient medium is ten times hotter and ten times less dense than the cloud so that it is in pressure equilibrium with the latter. We refer to the initial density contrast between the cloud and the medium as . The wind velocity () has an associated Mach number . This leads to the formation of a bow shock after which the post-shock subsonic flow interacts with the cloud and turns supersonic as it flows past it.

The blob test is useful for investigating how different hydrodynamics codes model astrophysical processes important for multiphase systems, such as shocks, ram-pressure stripping and fragmentation through KH and Rayleigh-Taylor (RT) instabilities. As , an approximate timescale of the cloud destruction is that of the full growth of the largest KH mode, i.e. the wavelength of the cloud’s radius. This can be obtained by considering the post-shock flow on the cloud and its time-dependence as the shock weakens and the cloud is accelerated. A full analysis of this test is presented in Agertz et al. (2007) and gives where is the crushing time and the velocity refers to the streaming velocity in the reference frame of the cloud. After this time, the cloud is expected to show a more complicated non-linear behaviour leading to disruption. The original blob test was initialised in a glass-like configuration obtained using simulated annealing using a standard SPH code. Since we now use OSPH rather than SPH, we must set up new ICs for the blob. Our new IC set up is described in detail in Appendix B. Unlike the previous blob test, where perturbations were seeded by random noise in the particle distribution, here we deliberately seed an inward growing mode on the front surface of the blob. This makes comparison between OSPH and FLASH simpler, since then the morphology of the blob is less dependent on small scale numerical noise.

The results are presented in Figure 8, where we compare SPH, TSPH and OSPH with increasing resolution with similar results from the Eulerian code FLASH (Fryxell et al., 2000). The SPH results (top panels) are similar to those presented in Agertz et al. (2007). The blob is squashed by the shock, but does not break up. There are no visible surface KHI or Rayleigh-Taylor instabilities. TSPH (second row) gives significantly improved results. The central depression is now resolved and the blob is mostly destroyed by . However, the density remains clumpy as compared to the FLASH simulation (bottom row). Our full OSPH scheme (third row) gives excellent agreement with the FLASH results. There are clear surface KHI and RT instabilities and the blob breaks up fully by . The precise details of the break up in FLASH and OSPH are different. However, these differences are smaller than those we observed between FLASH simulations of varying resolution. They are caused by the non-linear break up of the blob that is affected by resolution-scale perturbations.

Figure 9 shows and the pressure blips for the blob test in SPH (left), TSPH (middle) and OSPH (right) at . In SPH and TSPH, the two fluid phases (marked by the black and grey solid circles) remain well separated at all times. In both cases the pressure distribution shows discontinuities. By contrast, in OSPH the fluids are already mixed at