The Kelvin-Helmholtz instability and smoothed particle hydrodynamics

TL;DR

This paper demonstrates that standard smoothed particle hydrodynamics (SPH) can accurately model the Kelvin-Helmholtz instability and achieve convergence with reference solutions in both linear and non-linear regimes.

Contribution

The study shows that with appropriate physical viscosity and conductivity, SPH solutions converge to reference results, confirming its validity for modeling this instability.

Findings

SPH solutions converge to reference solutions in linear and non-linear regimes.

Physical Navier-Stokes viscosity improves convergence in non-linear regime.

Standard SPH with artificial viscosity can correctly model Kelvin-Helmholtz instability.

Abstract

There has been interest in recent years to assess the ability of astrophysical hydrodynamics codes to correctly model the Kelvin-Helmholtz instability. Smoothed particle hydrodynamics (SPH), in particular, has received significant attention, though there has yet to be a clear demonstration that SPH yields converged solutions that are in agreement with other methods. We have performed SPH simulations of the Kelvin-Helmholtz instability using the test problem put forward by Lecoanet et al (2016). We demonstrate that the SPH solutions converge to the reference solution in both the linear and non-linear regimes. Quantitative convergence in the strongly non-linear regime is achieved by using a physical Navier-Stokes viscosity and thermal conductivity. We conclude that standard SPH with an artificial viscosity can correctly capture the Kelvin-Helmholtz instability.

Tables

Table 1. Table 1: The measured convergence rates, Γ Γ \Gamma , of the ℒ 2 subscript ℒ 2 \mathcal{L}_{2} error as fit to ∝ n x − Γ proportional-to absent superscript subscript 𝑛 x Γ \propto n_{\rm x}^{-\Gamma} .

\br$t$	$\Gamma$
\mr2	0.89
4	0.62
6	0.26
8	0.28
\br

Equations

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla \cdot \mathbf{v},$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{\nabla P}{\rho }-\frac{1}{\rho }\nabla \cdot \mathbf{\Pi },$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}=-\frac{P}{\rho }\nabla \cdot \mathbf{v}+\nabla \cdot (\mathbf{v}\cdot \mathbf{\Pi })+\frac{1}{\rho }\nabla \cdot (\chi \rho \nabla T),$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=\frac{1}{\rho }\nabla \cdot ({\nu }_{\mathrm{c}}\rho \nabla c),$$

$${\Pi }^{ij}=\nu \rho \left(\frac{\partial v^i}{\partial x^j}+\frac{\partial v^j}{\partial x^i}-\frac{2}{3}\frac{\partial v^k}{\partial x^k}{\delta }^{ij}\right),$$

$$\rho =1,$$

$$v_x=v_0\left[\tanh \left(\frac{y-y_1}{a}\right)-\tanh \left(\frac{y-y_2}{a}\right)-1\right],$$

$$v_y=A\sin (2\pi x)\left[\exp \left(-\frac{(y-y_1{)}^2}{{\sigma }^2}\right)+\exp \left(-\frac{(y-y_2{)}^2}{{\sigma }^2}\right)\right],$$

$$P=10,$$

$$c=\frac{1}{2}\left[\tanh \left(\frac{y-y_1}{a}\right)-\tanh \left(\frac{y-y_2}{a}\right)+2\right],$$

$$\text{Re}=\frac{L\Delta v}{\nu },$$

$$s=-c\ln (c)$$

$$S=\int \rho s\mathrm{d}V,$$

$${\mathcal{L}}_2={\left[\sum _a(c_a^{\mathrm{SPH}}-c_a^{D2048}{)}^2\mathrm{d}V\right]}^{1/2},$$

Full text

1 Introduction

The goal of numerical simulation is to recreate real physics. Some numerical methods may yield a more accurate solution, but the true test of any numerical method is whether it produces solutions that are convergent upon a ‘true’ physical answer with increased computational resources.

This work studies the convergence properties of SPH (smoothed particle hydrodynamics; [1, 2]) on the Kelvin-Helmholtz instability. The Kelvin-Helmholtz instability occurs when two fluids move past each other, causing exponential growth of velocity orthogonal to the flow. This develops into vortices which curl back upon themselves – the hallmark of the instability. As the instability becomes non-linear, secondary instabilities form leading to the two fluids becoming well-mixed along the interface. The Kelvin-Helmholtz instability is inherently a mixing instability.

The ability of SPH, a Lagrangian particle-based method for solving the equations of hydrodynamics [1, 2], to correctly model the Kelvin-Helmholtz instability has been questioned in recent years [3, 4, 5, 6, 7, 8]. Agertz et al [9] concluded that there exists “fundamental differences between SPH and grid methods” in their ability to resolve mixing instabilities on the basis of SPH calculations that showed no growth of the Kelvin-Helmholtz instability.

The abysmal results by Agertz et al [9] were a consequence of their initial conditions which contained discontinuities in density, velocity and internal energy, combined with a lack of numerical treat for contact discontinuities. Discontinuities require special care numerically. For SPH, velocity discontinuities are treated by an artificial viscosity, and density discontinuities are handled naturally by the density summation since it makes no assumption about differentiability [10]. No treatment for the initial contact discontinuity was used by Agertz et al [9]. It has been demonstrated that including an artificial thermal conductivity is effective in capturing the Kelvin-Helmholtz instability and promoting mixing [10, 11].

Discontinuous initial conditions preclude formal convergence studies of the Kelvin-Helmholtz instability, however, even if the discontinuities are treated numerically. All wavenumber modes are unstable to growth when the interface is discontinuous [12], thus numerical noise may seed secondary instabilities that were not present in the initial conditions. Increasing the resolution of the calculation can alter the solution obtained since this permits access to higher wavenumber modes. Assertions have been made that some methods are of higher quality than others since they produce more detailed structure [13, 8], however, it is inappropriate to attribute structure generated from numerical noise as corresponding to resolved substructure [14, 15, 16].

In this work, we study the convergence properties of SPH on the Kelvin-Helmholtz instability. We use the smooth, well-posed initial conditions of Lecoanet et al [16] to obtain a converged solution, avoiding the issues inherent to discontinuous initial conditions. The calculations utilise a Navier-Stokes viscosity and physical thermal conductivity, thereby permitting the calculations to convergence in resolution to a particular solution in the strongly non-linear regime. Lecoanet et al [16] obtained converged solutions for this problem in both the linear and strongly non-linear regime using the finite-volume grid code Athena [17] and the pseudo-spectral code Dedalus.

2 Physical equations and initial conditions

We solve the set of equations given by

[TABLE]

where $\rho$ is the density, ${\bm{v}}$ is the velocity, $P$ is the pressure, ${\bm{\Pi}}$ is the Navier-Stokes stress tensor, $u$ is the internal energy, $T$ is the temperature, and $\chi$ is the thermal diffusivity. A passive scalar, $c$, which we call ‘colour’, is used to quantify the mixing of the two fluid regimes. An ideal equation of state is used, $P=\rho T$, with ratio of specific heats $\gamma=5/3$. The temperature is related to the internal energy according to $T=(\gamma-1)u$. The material derivative is ${\rm d}/{\rm d}t\equiv\partial/\partial t+{\bm{v}}\cdot\nabla$. The Navier-Stokes stress tensor is given by

[TABLE]

with shear viscosity $\nu$. The colour is passively advected with the flow, but includes a ‘physical’ diffusion term analogous to thermal conductivity. The evolution of the Kelvin-Helmholtz instability in the non-linear regime is strongly sensitive to the dissipation, both numerical and physical. In order to enforce one particular solution, the dissipation is made resolution independent through the physical dissipation terms.

The initial conditions are two-dimensional and given by

[TABLE]

where $a=0.05$, $\sigma=0.2$, $A=0.01$, $v_{0}=1$, $y_{1}=0.5$ and $y_{2}=1.5$. The calculations are performed in a periodic box of size $x\in[0,L]$ and $y\in[0,2L]$ with $L=1$. The Reynolds number used is Re=$10^{5}$, defined according to

[TABLE]

where $\Delta v=2v_{0}$. This yields $\nu=2\times 10^{-5}$, and for these calculations $\nu=\chi=\nu_{\rm c}$.

3 Numerical method

We use SPH to solve equations 1–4. We use a standard formulation of SPH where the density is computed through summation over neighbouring particles, and the smoothing length and density are self-consistently obtained through iteration [18, 19]. Artificial viscosity is used with the Morris & Monaghan limiter [20]. The calculations use the septic (M8) spline, a high-order kernel from the same family as the cubic and quintic B-splines [21].

The Navier-Stokes viscosity is solved using a two-first derivatives implementation, similar to [22, 23, 24]. The Navier-Stokes stress tensor (equation 5) is computed first, then the corresponding accelerations and heating are computed using the result. This implementation exactly conserves energy and momentum. The thermal conductivity and colour diffusion are computed directly using a second derivative in the manner of Brookshaw [25] and Cleary & Monaghan [26]. Total energy and colour are preserved with these schemes.

Calculations are performed using $256\times 592$, $512\times 1184$, $1024\times 2364$ and $2048\times 4728$ particles arranged on triangular lattices, labelled as $n_{\rm x}=256$, 512, 1024 and 2048, respectively. The $n_{\rm x}=2048$ calculation required $\sim 70~{}000$ cpu-hours of computational time.

4 Results

The goal is to obtain convergence of SPH calculations on the Kelvin-Helmholtz instability test of Lecoanet et al [16]. In absence of an analytic solution, convergence is measured by comparing SPH results to the D2048 solution obtained using the pseudo-spectral code Dedalus at a resolution of $n_{\rm x}=2048$ [16].

Figure 1 shows the colour field at $t=2$, 4, 6 and 8 for the $n_{\rm x}=2048$ SPH calculations alongside the D2048 reference solution. The structures in both solutions are remarkably similar at all times. In each case, the seeded mode forms a large singular curl ($t=2$), which continues to wind ($t=4$), leading to the two fluids becoming well mixed along the interface at late times ($t=6$ and 8). One minor difference is the degree of winding of the inner tip of the curl, which is not as tight for the SPH calculations compared to the reference solution. Several of the spurs at in the $t\geq 6$ snapshots are also of slightly different lengths. Despite these small differences, the SPH solution is in close proximity to the reference solution.

Figure 2 shows the growth of the amplitude of the seeded mode in the linear regime, obtained in a manner according to McNally, Lyra & Passy [15]. The measured growth rate of the mode is converged to $\propto\exp(\pi t)$. It is difficult, however, to obtain an analytic estimate of the growth rate for this problem. An incompressible fluid with discontinuous interfaces should have a growth rate that is $\propto\exp(2\pi t)$ [12]. A smoothed velocity interface will reduce the growth rate by approximately 20% [27]. It is reasonable to expect that the growth rate for this problem should be even slower than that. These calculations are for a compressible fluid, not incompressible, and additionally smooth the initial velocity perturbation used to seed the instability. Thus, while no analytic estimate of the growth rate of this mode is available, the results are at least consistent with the estimates available acting as upper bounds. The growth rate is converged even for our lowest resolution of $n_{\rm x}=256$ particles.

The degree is mixing is quantified by an entropic quantity. Defining

[TABLE]

as the colour entropy, the total colour entropy is calculated by the volume integral

[TABLE]

which may computed in SPH by the summation $\sum_{a}m_{a}s_{a}$. The total colour entropy only increases when $\nu_{\rm c}>0$. Figure 3 shows the total colour entropy as a function of time, with black circles the values of the reference solution at $t=2$, 4, 6 and 8. The shape of the total colour entropy curve matches the reference solution (see also [16]), and is converging towards the reference data points as the resolution improves.

Formal convergence of our calculations is assessed through the $\mathcal{L}_{2}$ error. The SPH particles for each calculation are interpolated to a $2048\times 4096$ grid and compared grid-cell to grid-cell with the reference solution. The error is computed according to

[TABLE]

where ${\rm d}V=2048^{-2}$ is the volume of each grid-cell, and $c^{\rm SPH}$ and $c^{\rm D2048}$ are the colour fields of the gridded SPH data and reference solution, respectively.

Figure 4 shows the $\mathcal{L}_{2}$ error of the SPH calculations with respect to the reference solution at $t=2$, 4, 6 and 8. Table 1 lists the measured convergence rates fit to $\propto n_{\rm x}^{-\Gamma}$. The errors in the SPH calculations may be contrasted to Athena calculations in Figure 2 of [16]. The maximum errors of the SPH calculations are of order $10^{-1}$, whereas the A1024 and A2048 calculations have errors on the order of $\sim 10^{-3}$ and $\sim 10^{-4}$, respectively, for $t>4$.

The error reduces linearly with respect to resolution at $t=2$, in line with the expected rate of convergence from the artificial dissipation. Even for the highest resolution calculation, the dissipation of kinetic energy from the artificial dissipation remains comparable to the dissipation by the Navier-Stokes viscosity. The convergence rate becomes sub-linear in the non-linear regime ($t\geq 4$), caused primarily by errors in the pressure gradient, which scale as $\mathcal{O}(1)$ with respect to resolution. This may be contrasted to Athena calculations, which have a measured rate of convergence that is between second to third order [16]. Importantly, the error of the SPH calculations is reducing at all times once the resolution is greater than $n_{\rm x}=512$ particles.

Reducing the error in the pressure gradient can be done by using higher-order smoothing kernels. For this problem, it was found that to obtain high-quality results necessitated the use of the septic spline, a high-order smoothing kernel from the B-spline family of kernels [21]. The septic spline is the high-order relative of the commonly used cubic and quintic splines. Using a high-order kernel is not an intrinsic requirement to activate the Kelvin-Helmholtz instability, but rather is needed to capture the amplitude of the initial velocity perturbation in the initial conditions and reduce spurious growth of other modes.

Figure 5 show the colour field at $t=4$ for calculations using the cubic spline through to nonic spline. For the cubic and quartic spline, the evolution of the vortex is significantly altered with respect to the reference solution (c.f. Figure 1). The resemblance of the SPH colour fields to the reference solution improves as the quality of smoothing kernel improves. The difference between the septic and nonic spline results is negligible, thus we conclude that the kernel bias is a sub-dominant source of error when these high-order kernels are used. The calculations examined in this work have used the septic spline.

5 Summary

We have performed SPH calculations of the Kelvin-Helmholtz instability. We used the smooth, well-defined initial conditions of Lecoanet et al [16], thereby avoiding convergence issues with discontinuous initial conditions. A Navier-Stokes viscosity and thermal conductivity were included to force the evolution of the instability in the non-linear regime. A passive scalar, referred to as ‘colour’, was added to the two fluids to measure the degree of mixing. A high-resolution calculation using the pseudo-spectral code Dedalus provided a reference solution to which results were compared.

The SPH results qualitatively matched the reference solution. The two fluids formed a curled vortex along the interface, which continued to spiral producing substructure and mixing of the two fluids. The substructure generated within the vortex matched between the SPH calculations and the reference solution well into the late non-linear regime.

Quantitative numerical convergence was measured for SPH on the Kelvin-Helmholtz instability, albeit with only linear convergence in the $\mathcal{L}_{2}$ error in the linear regime of the instability ($t=2$), and sub-linear convergence in the strongly non-linear regime. The degree of mixing was measured by defining an entropic quantity for the ‘colour’ field, and the total entropy of the colour field increased monotonically in agreement with the curve from the reference solution.

Importantly, convergence has been demonstrated using the standard form of SPH with an artificial viscosity, of the kind that has been used for decades [28, 18]. No alternative SPH formulations, modifications, or hacks were employed. The only requirement to achieve the results presented here was a high-order smoothing kernel, which was needed in order to capture the initial amplitude of the velocity perturbation in the initial conditions and to curb growth of other modes from numerical noise. We conclude that the rumours of the ‘fundamental flaws’ of SPH have been grossly exaggerated.

\ack

This work was motivated by attendance at AstroNum 2017 in Saint-Malo, France. Thank you to Daniel Lecoanet for providing the reference Dedalus solution, and to Daniel Price for useful discussions on physical dissipation terms in SPH.