2.2 Hydrodynamics
2.2.2 Particle-based method
As opposed to the grid-based method which solve the hydrodynamical equations in an Eulerian fashion (in the sense that the mesh is stationary), particle-based method is a Lagrangian method where the particles represent the sampling points of the uid and they advect with the uid. The uid equations in the Lagrangian form can be written as
dρ
dt =−ρ∇ ·v, (2.6)
dv
dt =−∇P
ρ , (2.7)
s = constant, (2.8)
where d/dt is the comoving derivative (or convective derivative) and s is the entropy function s ≡ P/ργ (with γ = 5/3 for the adiabatic equation of state). One of the most widely adopted method is the smoothed particle hydrodynamics (SPH). This is the method used in the simulations presented in this thesis.
Figure 2.1: A density slice (left) and the vorticity (right) in an AMR simulation taken from Iapichino et al. (2008). The hierarchical grid structure is superimposed in the right panels.
Smoothed particle hydrodynamics (SPH)
The central concept of SPH is that a uid quantityA at position xcan be estimated by a smoothed kernel weighted integral (or a convolution with the kernel function):
A(x) =˜ Z
A(x0)W(x−x0, h)d3x0, (2.9) where W(x−x0, h) is the kernel function, which is usually chosen to be a Gaussian-like function with a nite truncation radius h, which is called the smoothing length. Fig. 2.2 gives an illustration of the kernel interpolation. So dened, The derivative of a kernel-estimated quantities with respect to x can be written as
∇A(x) =g Z
A(x0)∇W(x−x0, h)d3x0. (2.10) As the derivative of the kernel function is known analytically, spatial discretization like nite-dierence is not required.
The second approximation is to discretize the uid into particles, usually with equal masses. For particle a with an associated densityρa and mass ma, its local volume can be estimated as ma/ρa. The density of particle a can be estimated by summation over the neighboring particles within the kernel at points b :
˜
ρa=X
b
ρbWabmb
ρb =X
b
mbWab, (2.11)
Figure 2.2: An illustration of the kernel interpolation, taken from Springel (2014).
where Wab ≡ W(xa−xb, ha) and ha is the smoothing length at point a. The smoothing length his usually chosen such that the number of particles within the kernel is constant, which implies that h is also a function of x.
At rst sight, it seems straightforward to discretize the equation of motion (Eq. 2.7)
into dva
dt =−1 ρa
X
b
mb
ρb Pb∇aWab. (2.12)
where va is the particle velocity of particle a and Pa is the thermal pressure of particlea. However, this formulation violates the momentum conservation, as the pressure gradient force that particlea exerts on particle b is in general dierent from the force exerted from particleb to particle a, i.e., madva/dt6=mbdvb/dt.
A trick can be applied to remedy this issue by noting the relation
∇P
ρ =∇(P
ρ) + P∇ρ
ρ2 , (2.13)
and insert it into the continuous form of the equation of motion. The discretized equation of motion can be written as1
dva
dt =−X
b
mb(Pa ρ2a +Pb
ρ2b)∇aWab. (2.14)
Eq. 2.14 ensures that the linear momentum in SPH is conserved exactly2 and it is the standard SPH formulation widely used in astrophysical applications. In an adiabatic uid
1Here for simplicity we assumehis constant. In reality wherehvariable, a correction term in Eq. 2.14 has to be introduced to ensure conservation properties.
2Note that there exists an innite number of possible variants in SPH by rearranging the continuous equations as Eq. 2.13, but not all of them preserve the conservation properties. Alternatively, one can derive Eq. 2.14 from a Lagrangian approach which manifestly guarantees conservation properties (Springel
& Hernquist, 2002).
the entropy function of each particle is conserved, and there is no need to evolve s with time. Therefore, both the mass and entropy function are trivially conserved in this SPH formulation.
Problems with SPH
Over the past decades, SPH has become one of the most popular methods in the eld of galaxy formation thanks to its built-in adaptive spatial resolution that can easily achieve a large dynamical range, especially in the situations where gravitational instability is impor-tant. Its Lagrangian nature also means that tracking the gas history is relatively simple compared to grid-based methods. However, it has been noticed that SPH suers from the so-called mixing problem which prevents the uid instabilities (e.g. the Kelvin-Helmholtz instability) (Agertz et al., 2007). Moreover, SPH tends to be very noisy at sub-sonic regimes where the sound speed is comparable or larger than the velocity of the ows (Bauer & Springel, 2012).
The E0 error
A simple error analysis of SPH kernel estimates using the Taylor expansion can shed light on the origin of the numerical diculties of SPH. For a uid quantityA, the kernel estimate at point a is
A˜a=X
b
mb
ρb AbWab =X
b
mb
ρb (Aa+∂iAa(xb−xa)i+O(h2))Wab
=Aa
X
b
mb
ρb Wab+∂iAa
X
b
mb
ρb (xb−xa)iWab+O(h2), (2.15) whereirepresents the Cartesian coordinates and the summation convention is used. There-fore, to be second order accurate, the following conditions have to be satised:
X
b
mb
ρb Wab= 1 and X
b
mb
ρb (xb−xa)iWab = 0. (2.16) Likewise, the gradient of the uid quantity A at point a is estimated by
∂giAa =X
b
mb ρb
Ab∂iWab =X
b
mb ρb
(Aa+∂jAa(xb−xa)j +O(h2))∂iWab
=AaX
b
mb
ρb ∂iWab+∂jAaX
b
mb
ρb (xb −xa)j∂iWab+O(h2), (2.17) where again iand j represent the Cartesian coordinates adopting the summation conven-tion. The partial derivative ∂i is with respect to pointa. The conditions for second order accuracy are
X
b
mb
ρb ∂iWab = 0 and X
b
mb
ρb (xb−xa)j∂iWab=δab, (2.18)
where δab is the Kronecker delta and the partial derivative is with respect to xa.
The conditions Eq. 2.16 and Eq. 2.18 are satised only when the particles are uniformly distributed within the kernel, which can easily be violated in the realistic applications such as strong shocks, contact discontinuities, shearing ows, and/or when a small number of particles within a kernel is adopted. The gradient estimate can be non-zero even in a constant eld A if the particle distribution is non-uniform. These errors are of zeroth order instead of second order when particle distributions are severely non-uniform and they are called the E0-error of SPH (Read, Hayeld & Agertz, 2010), which gives rise to both the mixing problem and the noisy behavior at sub-sonic regimes.
Higher order estimates
One could obtain more accurate gradient estimates by subtracting the rst term in Eq.
2.17
∂^iAa,1st =X
b
mb ρb
(Ab−Aa)∂iWab, (2.19) which is rst order accurate. So dened, the gradient estimate of a constant eld Awould be manifestly zero due to the (Ab −Aa) term, irrespective of the particle distribution.
This is the method commonly used in SPH to estimate the velocity gradients which do not appear in the equation of motion but are only used in the dissipation switches (articial viscosity, see Chapter 3).
One can even go on to pursue a second order accurate gradient by the relationAb−Aa =
∂jAa(xb−xa)j +O(h2) and rewrite Eq. 2.19 as
∂^iAa,1st ≈X
b
mb ρb
∂jAa(xb−xa)j∂iWab. (2.20) This equation can be written as a system of linear equations in a matrix form Ax = b, whereAij ≡P
b mb
ρb(xb−xa)j∂iWab,xj ≡∂jAa, andbi ≡P
b mb
ρb(Ab−Aa)∂iWab. Solving for xj with a matrix inversion one obtains the second-order accurate gradient. In Appendix A we demonstrate how the second-order accurate velocity gradients can be obtained. Higher order accuracy is also possible by further expanding the Taylor series.
The dilemma of SPH
It is therefore tempting to adopt these higher order methods for the pressure gradient in the equation of motion. Indeed, the equation of motion based on a rst-order accurate pressure gradient has been proposed to reduced the E0-error already in Morris (1996).
Unfortunately, the equation of motion that adopts a rst(or higher)-order accurate pressure gradient unavoidably violates the exact momentum conservation and thus leads to even larger errors especially at strong shocks. Therefore, one is forced to adopt the noisy estimate of the pressure gradient in SPH. This error is proportional to Pa/ρa ∝c2s,a wherecs,a is the sound speed at a (as the rst error term is proportional to Aa in Eq. 2.17), which makes SPH particularly noisy in sub-sonic regimes.
Chapter 3
SPHGal: Smoothed Particle
Hydrodynamics with improved accuracy for galaxy simulations
We present the smoothed-particle hydrodynamics implementation SPHGal, which com-bines some recently proposed improvements in Gadget. This includes a pressure-entropy formulation with a Wendland kernel, a higher order estimate of velocity gradients, a mod-ied articial viscosity switch with a modmod-ied strong limiter, and articial conduction of thermal energy. With a series of idealized hydrodynamic tests we show that the pressure-entropy formulation is ideal for resolving uid mixing at contact discontinuities but per-forms conspicuously worse at strong shocks due to the large entropy discontinuities. In-cluding articial conduction at shocks greatly improves the results. In simulations of Milky Way like disk galaxies a feedback-induced instability develops if too much articial viscosity is introduced. Our modied articial viscosity scheme prevents this instability and shows ecient shock capturing capability. We also investigate the star formation rate and the galactic outow. The star formation rates vary slightly for dierent SPH schemes while the mass loading is sensitive to the SPH scheme and signicantly reduced in our favored implementation. We compare the accretion behavior of the hot halo gas. The formation of cold blobs, an artifact of simple SPH implementations, can be eliminated eciently with proper uid mixing, either by conduction and/or by using a pressure-entropy formulation.
This chapter is based on Hu et al. (2014).