A well-balanced gas kinetic scheme for Navier-Stokes equations with gravitational potential

TL;DR

This paper introduces a well-balanced gas kinetic scheme for Navier-Stokes equations with gravity, accurately preserving hydrostatic equilibrium and effectively handling small perturbations and long-term simulations.

Contribution

A novel auxiliary variable-based scheme that exactly maintains hydrostatic balance and improves accuracy over traditional methods.

Findings

Preserves hydrostatic equilibrium with high accuracy.

Effectively simulates small perturbations and viscous effects.

Capable of long-time simulations converging to equilibrium.

Abstract

The hydrostatic equilibrium state is the consequence of the exact hydrostatic balance between hydrostatic pressure and external force. Standard finite volume or finite difference schemes cannot keep this balance exactly due to their unbalanced truncation errors. In this study, we introduce an auxiliary variable which becomes constant at isothermal hydrostatic equilibrium state and propose a well-balanced gas kinetic scheme for the Navier-Stokes equations with a global reconstruction. Through reformulating the convection term and the force term via the auxiliary variable, zero numerical flux and zero numerical source term are enforced at the hydrostatic equilibrium state instead of the balance between hydrostatic pressure and external force. Several problems are tested numerically to demonstrate the accuracy and the stability of the new scheme, and the results confirm that, the new scheme can preserve the exact hydrostatic solution. The small perturbation riding on hydrostatic equilibria can be calculated accurately. The viscous effect is also illustrated through the propagation of small perturbation and the Rayleigh-Taylor instability. More importantly, the new scheme is capable of simulating the process of converging towards hydrostatic equilibrium state from a highly non-balanced initial condition. The ultimate state of zero velocity and constant temperature is achieved up to machine accuracy. As demonstrated by the numerical experiments, the current scheme is very suitable for small amplitude perturbation and long time running under gravitational potential.

Equations

$$\displaystyle\left\{\begin{array}[]{rcc}\rho_{t}+(\rho U)_{x}&=&0,\\ (\rho U)_{t}+(\rho U^{2}+p)_{x}&=&0,\\ \displaystyle(\frac{1}{2}(\rho U^{2}+\frac{2}{\gamma-1}p))_{t}+(U\frac{1}{2}(\rho U^{2}+\frac{2\gamma}{\gamma-1}p))_{x}&=&0,\end{array}\right.$$

$$\rho ={\rho }_{ref},U=U_{ref},p=p_{ref}$$

$$\displaystyle\left\{\begin{array}[]{rcc}h_{t}+(hU)_{x}&=&0,\\ (hU)_{t}+(hU^{2}+\frac{1}{2}gh^{2})_{x}&=&-ghB_{x},\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{rcc}\rho_{t}+(\rho U)_{x}&=&0,\\ (\rho U)_{t}+(\rho U^{2}+p)_{x}&=&-\phi_{x}\rho,\\ \displaystyle(\frac{1}{2}(\rho U^{2}+\frac{2}{\gamma-1}p))_{t}+(\frac{1}{2}U(\rho U^{2}+\frac{2\gamma}{\gamma-1}p))_{x}&=&-\phi_{x}(\rho U),\end{array}\right.$$

$$\rho ={\rho }_{ref}\exp (-\frac{\varphi }{R{T}_{ref}}),U=0,T=T_{ref}.$$

$$\rho ={\rho }_{ref}\exp (-\frac{\varphi }{\alpha }).$$

$$\displaystyle\left\{\begin{array}[]{rcl}\displaystyle(\rho U)_{x}&=&0,\\ \displaystyle(\rho U^{2}+p)_{x}&=&-\phi_{x}\rho,\\ \displaystyle(\frac{1}{2}U(\rho U^{2}+\frac{2\gamma}{\gamma-1}p))_{x}&=&-\phi_{x}(\rho U).\end{array}\right.$$

$$\rho {\varphi }_x$$

$$(\rho U^2+p-\rho \alpha {)}_x+\rho {\alpha }_x(1+\frac{\varphi }{\alpha })=0.$$

$$(p-\rho \alpha {)}_x+\rho {\alpha }_x(1+\frac{\varphi }{\alpha })=0.$$

$$\frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}-{\varphi }_x\frac{\partial f}{\partial u}=\frac{g-f}{\tau },$$

$$g$$

$$\mathbf{W}=⟨\mathbf{\psi }f⟩$$

$$⟨f⟩={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }fdud\xi .$$

$$⟨u^n{\xi }^m{f}_u⟩$$

$$⟨\mathbf{\psi }(g-f)⟩=0.$$

$$\displaystyle\langle\psi f\rangle_{t}+\langle u\psi f\rangle_{x}-\phi_{x}\langle\psi f_{u}\rangle=0\quad\Rightarrow\left\{\begin{array}[]{rcl}\rho_{t}+(\mathcal{F}^{\rho})_{x}&=&0\\ (\rho U)_{t}+(\mathcal{F}^{\rho U})_{x}+\phi_{x}\rho&=&0\\ (\rho E)_{t}+(\mathcal{F}^{\rho E})_{x}+\phi_{x}(\rho U)&=&0\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{rcl}\displaystyle\rho_{t}+(\mathcal{F}^{\rho})_{x}-\mathcal{S}^{\rho}&=&0,\\ \displaystyle(\rho U)_{t}+(\mathcal{F}^{\rho U}-\rho\alpha)_{x}-\mathcal{S}^{\rho U}&=&0,\\ \displaystyle(\rho E)_{t}+(\mathcal{F}^{\rho E})_{x}-\mathcal{S}^{\rho E}&=&0.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{lll}\displaystyle\mathcal{F}^{\rho}=\langle uf\rangle,&\mathcal{F}^{\rho U}=\langle uuf\rangle,&\mathcal{F}^{\rho E}=\langle u\frac{1}{2}(u^{2}+\xi^{2})f\rangle,\\ \displaystyle\mathcal{S}^{\rho}=0,&\mathcal{S}^{\rho U}=-\rho\alpha_{x}(1+\frac{\phi}{\alpha}),&\mathcal{S}^{\rho E}=-\phi_{x}(\rho U)\end{array}\right.$$

$$f=g-\tau (g_t+u{g}_x-{\varphi }_x{g}_u)+O({\tau }^2).$$

$${\mathbf{W}}_i^{n+1}={\mathbf{W}}_i^n-\frac{1}{\Delta x}(F_{i+1/2}-F_{i-1/2})+S_i,i=1,2,...,n$$

$$\displaystyle F=\int_{0}^{\Delta t}\left(\begin{array}[]{c}\mathcal{F}^{\rho}\\ \mathcal{F}^{\rho U}-\rho\alpha\\ \mathcal{F}^{\rho E}\end{array}\right)dt=\int_{0}^{\Delta t}\left(\langle u\bm{\psi}f\rangle+\left(\begin{array}[]{c}0\\ \displaystyle-\rho_{0}\alpha_{0}\\ \displaystyle 0\end{array}\right)\right)dt,$$

$$\displaystyle S=\frac{1}{\Delta x}\int_{0}^{\Delta t}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left(\begin{array}[]{c}\mathcal{S}^{\rho}\\ \mathcal{S}^{\rho U}\\ \mathcal{S}^{\rho E}\end{array}\right)dxdt=\frac{-1}{\Delta x}\int_{0}^{\Delta t}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left(\begin{array}[]{c}0\\ \displaystyle\rho_{0}\alpha_{x}(1+\frac{\phi}{\alpha_{0}})\\ \displaystyle-\phi_{x}(\rho U)\end{array}\right)dxdt.$$

$$f(t)=g_0-\tau (g_t+u{g}_x-{\varphi }_x{g}_u)+g_{t}t,$$

$$⟨\psi (f-g)⟩=\tau ⟨\psi (g_t+u{g}_x-{\varphi }_x{g}_u)⟩=0.$$

$$⟨\psi g_t⟩=-⟨\psi (u{g}_x-{\varphi }_x{g}_u)⟩,$$

$$\displaystyle\Delta\bm{W}^{*}=\left(\begin{array}[]{c}\Delta\rho\\ \Delta(\rho U)\\ \Delta(\rho E)\end{array}\right)=-\Delta t\left(\begin{array}[]{c}\langle ug_{x}\rangle\\ \langle uug_{x}\rangle+\phi_{x}\rho_{0}\\ \langle\frac{1}{2}u(u^{2}+\xi^{2})g_{x}\rangle+\phi_{x}\rho_{0}U_{0}\end{array}\right).$$

$$g^{\ast }=\mathcal{M}({\mathbf{W}}^{\ast })=\mathcal{M}({\mathbf{W}}_0+\Delta {\mathbf{W}}^{\ast }).$$

$$g_t=\frac{g^{\ast }-g_0}{\Delta t}.$$

$$\alpha =\frac{\varphi }{\ln {\rho }_{ref}-\ln \rho }.$$

Full text

well-balanced, source term, gravity, gas kinetic scheme

I Introduction

Gravity is involved in many physical problems, including astrophysical problems like core-collapse supernova, atmospheric motions on planet, smoke stratification in compartment fires etc. In order to understand these phenomena and make reliable prediction, conservation laws with gravitational force are invoked in the form of partial differential equations. However, numerical simulations of these systems are not easy from the following aspects: (1) for a long time evolution, the truncation error will accumulate and dramatically affect the final solution of an isolated gravitational system Tian et al. (2007); Xing and Shu (2013); (2) to predict small perturbation, say, numerical weather prediction and climate modeling, the truncation error will mask the small perturbations on the top of stationary solution Xing and Shu (2013); Chandrashekar and Klingenberg (2015).

These two problems can be attributed to unbalanced discretization of the convection term and gravitational force term. Consider a fluid system under gravity governed by the Euler equations. The fluid system possesses a stationary state known as the hydrostatic equilibrium state in which the gravitational force is exactly balanced by the pressure gradient. However, in conventional numerical schemes, the gravitational force term and the convection term are discretized separately, thereby, the truncation errors cannot cancel each other. As a result, the conventional numerical schemes are not able to preserve the hydrostatic equilibrium state, and can induce unacceptable spurious motions Botta et al. (2004). In this context, a numerical scheme which ensures the hydrostatic balance exactly on discrete level is termed well-balanced.

Many well-balanced schemes has been developed, especially for shallow water equations Audusse, Bristeau, and Perthame (2000); LeVeque (1998); Xing and Shu (2005). But the techniques developed for shallow water equations seem not easy to be implemented in the Euler equations with gravitational force term. This problem bothers the CFD community for a long time.

Botta et al.Botta et al. (2004) developed a well-balanced finite volume method for the Euler equations, using a discrete Archimedes principle to express the gravity source term as the cell surface integral of the reconstructed hydrostatic pressure.With the help of kinetic theory, Xu et al.Xu, Luo, and Chen (2010) proposed a kind of well-balanced scheme for the Euler equations in 2010. The gravitational potential is approximated as a step function inside each cell, and the amount of particle penetration and reflection from the cell interface is evaluated according to the incident particle velocity and the strength of the potential barrier at the cell interface. This scheme can maintain hydrostatic equilibrium state exactly if the numerical integration used in kinetic flux is evaluated accurately.

In recent years, the development of well-balanced schemes for the Euler equations has gained more attention. In 2012, Xing and Shu Xing and Shu (2013) developed a special source term discretization so that the resulting WENO scheme balances the zero-velocity and constant temperature steady state solutions to machine accuracy, and at the same time maintains the high order accuracy and essentially non-oscillatory property for general solutions. Unlike the kinetic scheme with step potential assumption, their scheme is free of analytical or numerical integration. Käppeli and Mishra Käppeli and Mishra (2014) developed a novel reconstruction of the enthalpy which is based on local constant entropy assumption. After that they proposed a more general pressure reconstruction using a local analytical integration of hydrostatic equation, and demonstrated the efficiency of their well-balanced schemes for a broad set of astrophysical scenarios with several types of equation of state Käppeli and Mishra (2016). Ghosh and Constantinescu Ghosh and Constantinescu (2015, 2016) extended Xing and Shu’s work to more general flows encountered in atmospheric simulations. Besides an isothermal equilibrium, the proposed well-balanced scheme can hold for many other hydrostatic equilibrium states. Chandrashekar and Klingenberg Chandrashekar and Klingenberg (2015) also proposed a well-balanced scheme which involves a specific combination of source term discretization similar to Xing and Shu’s, and requests that numerical flux exactly resolves stationary contacts. The scheme is able to preserve isothermal and polytropic stationary solutions up to machine precision. Chertock et al. Chertock et al. (2018) proposed a global flux which reflects the accumulating effect of the gravitational potential along a direction and developed a well-balanced scheme for the Euler equations.

For the Navier-Stokes equations, the well-balanced scheme is rarely reported except the gas kinetic scheme(GKS) Luo, Xu, and Liu (2011). As early as last century, Slyz and Prendergast Slyz and Prendergast (1999) incorporated a time-independent gravitational potential into the gas kinetic scheme and guaranteed conservation of total (kinetic+internal+gravitational) energy, but did not consider the gravity in numerical flux. Tian et al Tian et al. (2007) introduced the gravitational source term into numerical flux through the Chapman-Enskog expansion of the BGK equation. Although many attempts have been made in the GKS framework, the GKS at that time were only of second-order accuracy with respect to holding the hydrostatic equilibrium state, and could not maintain the zero velocity and constant temperature to the machine zero. In 2011, Luo et alLuo, Xu, and Liu (2011) followed Xu’s idea, used a piecewise constant function inside each cell to represent gravitational potential and the physical mechanism of particle transport across the potential barrier is explicitly used in the flux evaluation. The proposed symplecticity-preserving gas kinetic scheme was proofed to be well-balanced for the Navier-Stokes equations (WB-NS, for short). However, the disadvantages of the symplecticity-preserving gas kinetic scheme are also obvious: (1) the step piecewise representation of potential restricts the accuracy of numerical solution at a very low level; (2) potential barrier at cell interface complicates the computation of numerical flux.

The present study mainly focuses on the development of well-balanced gas kinetic scheme for the Navier-Stokes equations. We are going to introduce an auxiliary variable which allows us to develop a simple and computational efficient well-balanced gas kinetic scheme under arbitrary potential function. Through the evolution towards the equilibrium state, we will show the difference between the equilibrium states of the Euler equations and the Navier-Stokes equations.

The remaining part is organized as follows: section 2 briefly introduces the auxiliary variables and modified the source terms in hydrodynamic equations; section 3 introduces the numerical flux calculation and the interpolation of the auxiliary variables; section 4 discusses the well-balanced property and the convergence towards isothermal hydrostatic equilibrium state through several numerical tests; section 5 concludes this study.

II Auxiliary variable and modified source terms

II.1 Two strategies to eliminate the truncation errors

Truncation error is inevitable in many numerical schemes. In order to eliminate the truncation error, Xing and Shu Xing and Shu (2013) proposed a strategy in which the convection term and the source term are discretized by the same difference scheme. Although every single discretization generates truncation error, identical difference scheme guarantees the truncation error cancels each other. This strategy can be labeled as generating and eliminating.

On the contrary, another strategy is to prevent the occurrence of truncation error. As we know, zero truncation error only occurs under special circumstance, say, the discretization of a constant function. For example, consider the one dimensional Euler equations without external force,

[TABLE]

where $\rho$ represents the density, $U$ represents the velocity, $p$ denotes the pressure, and $\gamma$ denotes the specific heat ratio. The corresponding hydrostatic equilibrium state is trivial as the density, velocity and pressure are uniform everywhere.

[TABLE]

In fact, every available numerical schemes are well-balanced for this hydrostatic equilibrium state.

Another example is the hydrostatic solution for the shallow water equations,

[TABLE]

where $h$ represents the water height above this bottom, $U$ is the velocity, $g$ is the gravitational constant, and $B$ is the bottom elevation. The top surface is at $B+h$ which is a constant at hydrostatic equilibrium. By taking advantage of constant surface elevation ($B+h$), a kind of well-balanced scheme is developed Xu (2002); Zhou et al. (2001).

However, constant function will not happen naturally in most cases. Consider the Euler equations under an external potential,

[TABLE]

where $\phi$ is the external potential. If ideal gas state equation is adopted, $p=\rho RT$, where $R$ is the gas constant and $T$ is the temperature, the isothermal hydrostatic equilibrium state is

[TABLE]

Since the exponential function cannot be accurately approximated by numerical discretizations based on Taylor expansion (polynomial), many numerical schemes cannot exactly hold this hydrostatic equilibrium state.

II.2 Auxiliary variable

Inspired by Zhou et al.’s workZhou et al. (2001), we adopt the second strategy to develop a well-balanced scheme for the Navier-Stokes equations with external force. To do so, we first introduce an auxiliary variable, $\alpha$, which varies with location and satisfies the following equation,

[TABLE]

where $\rho_{ref}$ is a constant reference density. $\alpha$ can be regarded as an analogue of $RT$ based on dimensional analysis. More importantly, when the gas system rests at hydrostatic equilibrium state, the temperature is uniform throughout the space for isolated systems. Therefore, $\alpha$ is also a constant if $\rho_{ref}$ is chosen properly.

In numerical scheme, the interpolation of density can be replaced by the interpolations of $\alpha$ and potential function $\phi$ through Eq.(14). Therefore, we can circumvent the interpolation of the exponential function, and only deal with a constant function at the hydrostatic equilibrium state.

II.3 Modified source terms

The stationary Euler equation with external potential is written as follows,

[TABLE]

It is obvious that the mass and energy equations will be satisfied at hydrostatic equilibrium state because of zero velocity. But the flux of momentum is nonzero and balanced by the external force. Actually, the momentum equation is the only obstacle for numerically holding the hydrostatic equilibrium state.

If the external force term are discretized directly in the momentum equation, the nonlinear distributed density and pressure will introduce nonzero truncation error and then induce nonzero velocity or oscillate the solution. In order to remove the truncation error, we propose a new strategy to calculate the source term. Substitute Eq.(14) into the source term in the momentum equation,

[TABLE]

Then reformulate the steady momentum equation as follows,

[TABLE]

At the hydrostatic equilibrium state, the velocity is zero, thereby the Eq.(20) becomes,

[TABLE]

The first term is solely composed of the static pressure and $-\rho\alpha$ which can be taken as a negative pressure ($-p$) and will completely cancel the static pressure. As a result, the first term on left hand side of the Eq.(21) is zero. The consequence is that the discretization of the first term will not introduce any truncation error at the hydrostatic equilibrium. The last term on the left hand side of the Eq.(21) also vanishes since $\alpha$ is a constant. Therefore the hydrostatic equilibrium state is exactly held by the above equation, as long as the discretization can exactly approximate the zero derivative of the constant value of $\alpha$, which can be easily fulfilled by most of numerical schemes. One thing we should emphasize is that, the modified momentum equation will lead to non-conservative scheme, which is discussed in Appendix B. In the following section, we will use the auxiliary variable and modified momentum equation (Eqs.(14,19)) to construct a well-balanced scheme for the Navier-Stokes equations.

III Well-balanced gas kinetic scheme with external force

III.1 BGK equation with external force

Consider the dimensionless BGK equation under external potential force in one dimensional space,

[TABLE]

where $\tau$ denotes the relaxation time, $u$ denotes the particle velocity, $f=f(x,u,\xi,t)$ represents the velocity distribution function, $g=g(x,u,\xi,t)$ represents the corresponding equilibrium state, the Maxwellian distribution function, which can be expressed as follows,

[TABLE]

where $\xi=(\xi_{1},\xi_{2},...,\xi_{k})$ denotes the effective internal freedom and $k$ is the degree of effective internal freedom ($\gamma=\frac{k+3}{k+1}$ in one dimensional problem). The macroscopic variables can be derived by taking the moments of the microscopic distribution function,

[TABLE]

where $\rho E=\frac{1}{2}\rho(U^{2}+(k+1)RT)$ denotes the total energy. The symbol $\langle f\rangle$ is defined as,

[TABLE]

The moments of $f_{u}$ can be expressed by lower order moments of $f$,

[TABLE]

The conservation of collision term requires,

[TABLE]

Taking moments of the BGK equation, we have,

[TABLE]

Using the auxiliary variable to reformulate the above equations, we have,

[TABLE]

The flux and the source term can be explicitly expressed as follows,

[TABLE]

In this study the Chapman-Enskog expansion is adopted in order to solve the Navier-Stokes equations,

[TABLE]

The second term on the right hand side corresponds to the Navier-Stokes constitutive relationship Chapman and Cowling (1970).

III.2 GKS flux with external force

Consider the conservation law of $\bf{W}$ in a one dimensional control volume $\Delta x$ during time interval $\Delta t$,

[TABLE]

where $F$ denotes the numerical fluxes and $S$ represents the numerical source term during a time step. As the modified equations present, the flux term and source term in this study become,

[TABLE]

[TABLE]

The local approximate solution at cell interface is,

[TABLE]

where $g_{0}$ defined at a cell interface is the equilibrium state at the beginning of the time step. The spatial derivative of the equilibrium state is calculated from the derivative of the macroscopic variables (Appendix A). The moments involved with $g_{u}$ can be explicitly calculated by Eq.(35). Because of the conservation of the collision term, we have the following equation,

[TABLE]

The time derivative of conservative variables $\langle\psi g_{t}\rangle$ can be derived,

[TABLE]

Assume that $g_{x}$ and $\phi_{x}$ are constant during the time step and conservative variables can be expressed as the Taylor expansion in terms of time. Integrate Eq.(65) over the time step, and only retain leading order terms up to $O(\Delta t)$,

[TABLE]

where subscript ”0” denotes initial time at the beginning of the time step. Then an intermediate equilibrium state at the end of the time step can be constructed,

[TABLE]

This intermediate state is not the new state for the next time step, but only used to estimate the numerical time derivative,

[TABLE]

In fact, this procedure can be seen as the use of the Euler equations to predict the time derivative. As all the terms in Eq.(63) are derived, the numerical fluxes (Eq.(55)) and the numerical source term (Eq.(62)) can be calculated.

III.3 Discretization

Suppose the computational domain is uniformly discretized by $n$ cells, and the cell size is $\Delta x$. The variables defined at cell $i$ are denoted by their subscript $i$ and the variables defined at the cell interface between $i$ and $i+1$ cells are denoted by subscript $i+1/2$.

III.3.1 Reference density

Reformulate Eq.(14) as follows,

[TABLE]

Considering the above equation and Eq. (14), it requires that $\rho_{ref}\neq\rho$ and $\phi\neq 0$. Therefore, we define the reference density as follows,

[TABLE]

where $\max(\rho)$ denotes the cell index whose density is maximum in the computational domain. Furthermore, a positive value is added to the potential function to ensure its positivity throughout the entire computational domain. $\rho_{ref}$ is calculated at the beginning of every time step globally.

III.3.2 Auxiliary variables

Instead of using the original variables in Eq. (13),

[TABLE]

we perform interpolation with the following set of variables in which $\alpha$ instead of $\rho$ will be constant at the hydrostatic equilibrium state,

[TABLE]

Linear interpolation and central difference are adopted to approximate the interfacial values and corresponding derivatives respectively, that is,

[TABLE]

Then conservative variables and their the derivatives are derived from $\bar{\bm{v}}$ and $\frac{\partial\bar{\bm{v}}}{\partial x}$ by the transformation and the chain role in calculus. The quantities at cell interface (denoted by subscript ”$i+1/2$”) can be calculated via the interpolated values and derivatives.

The spatial integral of the source term is approximated by the trapezoidal rule of two interfacial values,

[TABLE]

This procedure guarantees identical discretization is adopted for the flux and the source term, which is crucial for WB scheme as mentioned by Xing and Shu Xing and Shu (2013).

III.4 Well-Balanced property

Under isothermal hydrostatic equilibrium state ($\alpha=\alpha_{iso}=RT_{iso},U=U_{iso}=0,T=T_{iso},\phi$), the discretization introduced in last subsection will exactly reproduce the isothermal quantities, $\alpha_{iso},U_{iso},T_{iso}$ and their zero derivatives. It is easy to verify the following identities (Appendix A),

[TABLE]

As a result, the time derivative, $g_{t}$ (Eq.(65)) will completely vanish at the hydrostatic equilibrium state. Hence, the numerical fluxes (Eq.(55)) become,

[TABLE]

and the interfacial value of source term is zero too,

[TABLE]

so is the volume integral (Eq.(81), $S_{i}=0$). The numerical source term and numerical flux term are both zero exactly. Therefore, the isothermal hydrostatic equilibrium state can be exactly held by the proposed numerical scheme regardless of the shape of $\phi$.

IV Numerical results

IV.1 Adiabatic Boundary condition

The adiabatic boundary condition is adopted in order to exclude the effects of external heating when simulating the converging process towards the hydrostatic equilibrium state. Under this circumstance, the entire computational domain becomes an isolated system. We use ghost cell technique to realize adiabatic boundary condition. The density, temperature and potential functions in the ghost cell at the boundary are assigned the same as those of the direct neighboring cell in the flow field, and the fluid velocity is set to be the opposite of that in the neighboring cell in the flow field,

[TABLE]

where subscript ”$g$” represents the ghost cell, and ”$f$” represents the flow field. Thus no mass and no heat penetrate the solid boundaries. In the following numerical simulations, all the boundary conditions are adiabatic boundary condition, specific heat ratio of gas is 1.4 and gas constant is 1.0 if not specified. In most cases, the sound speed of the initial condition is $c_{s}=\sqrt{1.4}$; the length of computational domain is $L=1$; so the sound crossing time defined as $\tau=2L/c_{s}$ is about $1.69$.

IV.2 Maintaining the isothermal hydrostatic equilibrium state

The well-balanced property of the present scheme is demonstrated firstly. We choose three different potentials to verify our code, and stop the simulations at time $t=2.0$ to check whether the hydrostatic equilibrium state is kept. The potential functions are given as follows,

[TABLE]

The initial condition is given by Eq.(13) with $\rho_{ref}=1$ and $T_{ref}=1$. The computational domain is uniformly divided into 100 cells. Adiabatic boundary condition is adopted at the two ends of computational domain.

As a comparison, a primary scheme is employed for the simulation of $\phi_{1}(x)=x+1$ case, in which the conservative variables ($\rho,\rho U,\rho E$) and original source term are adopted, and the spatial integral of the Eq.(62) is also calculated by trapezoidal rule. This scheme is referred as ”nWB” hereafter. As shown in Fig. 2, the velocity cannot stay at zero and oscillates at the boundaries, and the temperature also deviates from the equilibrium condition, which implies the scheme is not well-balanced.

Then, auxiliary variables $\bar{\bm{v}}_{i}$ are used for the interpolation and the Eq.(81) is adopted for the spatial integral of the source term. As we expected, the modified equations and the new interpolation technique work perfectly. In Fig. 3 the hydrostatic equilibrium states are maintained up to the machine accuracy for all the potential functions, which verifies the well-balanced property.

IV.3 Propagation of perturbation on an isothermal hydrostatic equilibrium state

Well-balanced property is of great importance for simulating the propagation of small perturbation riding on an hydrostatic equilibrium state. We simulate a test case from the referenceLeVeque and Bale (1999) to demonstrate the efficiency of the present scheme. Consider an ideal gas with $\gamma=1.4$ staying initially at an isothermal hydrostatic equilibrium state,

[TABLE]

for $x\in[0,1]$. Then the initial pressure is perturbed by

[TABLE]

where $\eta$ is the amplitude of the perturbation. The potential is given as follows,

[TABLE]

The computation is conducted with $100$ uniform cells in the whole domain and stops at time $t=0.25$. With $0<\eta\ll 1$, the initial perturbation splits into two waves spreading on both sides.

IV.3.1 Inviscid flow

The Euler equations are considered firstly. The initial conditions are shown in Fig. 4(a). The pressure is disturbed by a moderate small perturbation, $\eta=10^{-5}$. The benchmark solution is derived by the present well-balanced scheme with 24300 uniform cells. As shown in the figure 4(b), the green dash dot line derived by WB scheme is very close to the benchmark solution, while the result derived by nWB scheme deviates from the benchmark solution. Then the convergence study is conducted by refining the mesh. And the error is defined as the $L_{1}$ norm of the deviation from the benchmark solution. Figure 4(c) shows that, the convergence rate is 1.9974 for the WB scheme, and is 2.0050 for the nWB scheme. Although both of them are of second order spatial accuracy, the WB scheme is much more accurate than the nWB scheme on coarse mesh.

Then we tested two smaller amplitudes of the perturbation, $\eta=10^{-13}$ and $10^{-14}$. Since no other numerical results has been reported for such small perturbation, we propose a self-evaluation procedure to assess the present well-balanced scheme. As the perturbation is very small, the flow system will respond linearly, namely, if normalize the numerical results by their amplitude of the initial perturbation, the rescaled numerical results will collapse to a single curve. Figure 5 shows that the normalized solutions for $\eta=10^{-5},10^{-13}$ coincide with each other, but the normalized solution for $\eta=10^{-14}$ deviates from others on the level of $10^{-16}$. In fact, the rounding error is about $10^{-16}$ if double-precision is adopted in the computation.

IV.3.2 Viscous flow

As we claimed the well-balanced scheme for the Navier-Stokes equations, the small perturbation propagating in viscous flow ($\nu=0.01$) is also simulated. Two amplitudes of the perturbation, $\eta=10^{-5},10^{-13}$, are considered and the results are also compared with the benchmark solution derived by the present WB scheme with 7290 uniform cells. As shown in Fig. 6, the final amplitudes are smaller than that in the inviscid fluid. More importantly, the normalized solutions are also identical to each other, which means the present scheme for viscous term discretization is well-balanced. Otherwise, the truncation error will pollute the numerical results for smaller perturbation more severely, and separate two profiles in Fig.6. Similar to the inviscid case, the convergence rate is 2.1519 for the WB scheme, and is 1.9995 for the nWB scheme. This challenging test demonstrated that our scheme can predict very accurate numerical results for small perturbations up to the machine zero for the Navier-Stokes equations.

IV.4 Evolution towards the isothermal hydrostatic equilibrium state

Since the Euler equations possess many equilibrium states, a sort of well-balanced scheme has been proposed for a variety of equilibrium states Chandrashekar and Klingenberg (2015). On the other hand, without heat conduction, the fluid system governed by the Euler equations allows temperature stratification, and cannot determine which equilibrium state the system will eventually stay at. As a result, the previous Euler equations’ well-balanced schemes generally did not test the convergence towards isothermal hydrostatic equilibrium state. As a contrast, the Navier-Stokes equations only allow the fluid system eventually converge to the isothermal hydrostatic equilibrium state for an isolated system. Therefore, the evolution towards the isothermal hydrostatic equilibrium state distinguishes the well-balanced scheme for the Navier-Stokes equations from the well-balanced scheme for the Euler equations.

To test the convergence process, the three potential functions (Eq.(93)) are employed again, and the gas properties are $\nu=0.005$, $\gamma=1.4$ and $\text{Pr}=1$. The initial condition is given as follows,

[TABLE]

Notice that this initial condition deviates far from the equilibrium state regarding anyone of the three potential functions. The computational domain is uniformly divided into 100 cells. The simulations stop at $t=1000$. Fig. 7 shows the velocity and temperature solution at $t=40,1000$ under the three potential functions. It can be seen that the system converges to the stationary isothermal hydrostatic equilibrium state for all three potential functions. Even for the most complex sine function, the final velocity is less than $10^{-14}$.

By contrast, Fig. 8 shows gas evolution at three different times with $\phi_{1}$, zero viscosity and zero thermal conductivity. It can be seen that the vital movement of gas lasts for long time. Since the light gas with high temperature will move towards the end of low potential, the high temperature gas will accumulate at the right end if the thermal conduction is absent. As a result, the temperature profile becomes very steep, and numerical solution oscillates due to the central difference in the scheme. Fig. 8 indicates that the Euler equations are inappropriate for long time simulation because of the lack of the physical dissipation mechanism.

IV.5 2D Rayleigh-Taylor instability

In this subsection, two dimensional two-layer flow under gravity is simulated to validate the proposed scheme. A radially symmetric potential is given in polar coordinate as follows,

[TABLE]

where $r$ is the radius, and $\theta$ is the azimuth in polar coordinate. Two isothermal equilibrium states are assigned to the inner layer ($r<r_{i}$) and outer layer ($r\geq r_{i}$) respectively,

[TABLE]

[TABLE]

where $\rho_{0}=0.1$, $\Delta\rho=0.1$, $a=1.0$, $b=1.0$, $T_{0}=0.3$ and $R=1.0$. An interface located at $r=r_{i}$ separates the two layers, and is twisted so as to make the cold fluid (denser at the interface) penetrate into the hot fluid (lighter at the interface),

[TABLE]

where $r_{0}=0.6$, $\eta=0.02$, $\kappa=20$. Note that, there are varied numerical configurations in previous literature LeVeque and Bale (1999); Tian et al. (2007); Luo, Xu, and Liu (2011); Chandrashekar and Klingenberg (2015). In the following numerical simulations, the simpler configuration (Eqs.(100,104,105)) is adopted and a complete set of parameter including kinematic viscosity is presented. The simulations are conducted on a uniform grid covering $[-1,1]\times[-1,1]$ with $200\times 200$ cells.

In Fig. 9, the density contour is presented and a temperature contour line ($T=0.23$) which separates the hot and cold fluids is also provided to illustrate the spike structure more clearly. Theoretically, if polar grid system is employed, the numerical solution will be symmetric. However, the projection onto the Cartesian mesh can be regarded as a force disturbing the symmetry. This is the origin of asymmetrical pattern in this test problem. On the other hand, the viscous term can be regarded as a resistance to the disturbing force. As shown in Fig. 9, the flow pattern of the RT instability is highly dependent on the viscosity. When $\nu$ is small (Fig.(9(a)), the disturbance from the projection onto the Cartesian mesh breaks the symmetry of the solution. As $\nu$ increases (Fig.(9(b))), the viscous effect suppresses the disturbance and stabilizes the mushroom-shaped spikes and maintains the symmetry much better. Further increasing $\nu$ (Fig.(9(c,d))), the mushroom-shaped structures become smooth and further fade away with large viscosity. In fact, all these three types of flow pattern were reported in the literatureLeVeque and Bale (1999); Tian et al. (2007); Luo, Xu, and Liu (2011); Chandrashekar and Klingenberg (2015), because the numerical dissipation was implicitly implemented in their solvers of the Euler equations.

We further refined the meshes, and present a high resolution result (Fig.(10)) on a $800\times 800$ uniform mesh. It can be found that, since the disturbing force becomes smaller on the refined mesh, the flow pattern keeps its symmetry more easily even with small viscosity ($\nu=0.0001$). According to this result, the numerical solutions on $200\times 200$ mesh are far from convergency.

V Conclusion

In this study, we introduced an auxiliary variable which becomes constant at isothermal hydrostatic equilibrium state and reformulated the source term in the Navier-Stokes equations into a convenient form. A part of the original source term is merged into the flux term, and the remaining source term becomes zero at isothermal hydrostatic equilibrium state.

Based on the reformulated source term, we proposed a second-order well-balanced gas kinetic scheme for the Navier-Stokes equations. Through the global reconstruction of the auxiliary variables, the numerical fluxes and numerical source term vanish simultaneously when approaching the hydrostatic equilibrium state, which guarantees the well-balanced property. The new scheme has no assumption of the potential function, and hence, is simple and computationally efficient compared to symplecticity-preserving GKS.

Several test cases were presented to demonstrate the accuracy and the stability of the new scheme. The one-dimensional hydrostatic equilibrium can be exactly held up to machine accuracy by the proposed scheme. The results for wave propagation riding on a hydrostatic equilibria has shown a significant gain in accuracy with current scheme. The small perturbation only several orders greater than the machine accuracy still survived as the simulation proceeded. Moreover, the linear response with viscous effect was also predicted accurately.

Another featured property of well-balanced scheme for the NS equations is the capability of simulating the evolution towards the hydrostatic equilibrium state. It requires that the physical dissipation and heat transfer, which are missing in the Euler equations, must be properly represented in the scheme for long period simulation. And the proposed scheme correctly predicted an isothermal hydrostatic equilibrium state after a long running time.

In summary, a well-balanced gas kinetic scheme for the Navier-Stokes equations is proposed and validated through several challenging numerical problems. The evolution towards the isothermal hydrostatic equilibrium from highly non-equilibrium state is simulated and the final equilibria is achieved. The current scheme is capable of making accurate prediction for small amplitude perturbation and long time running.

Acknowledgments This work was supported by the National Science Foundation of China (11602091, 91530319) and the National Key Research and Development Plan (No. 2016YFB0600805).

Appendix  A GKS formula

Let $g^{\prime}$ be the normalized equilibrium distribution function,

[TABLE]

A.1 Moments of the Maxwellian distribution function

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.2 Derivative of the Maxwellian distribution function

[TABLE]

where $s$ represents the space coordinate or time coordinate. Since $u$ and $\xi$ are independent of coordinate $s$, the derivatives of the Maxwellian distribution function can be expressed by the derivative of macroscopic variables, say, $\frac{\partial\rho}{\partial s}$, $\frac{\partial U}{\partial s}$ and $\frac{\partial T}{\partial s}$, or other set macroscopic variables,

[TABLE]

A.3 The identities at isothermal hydrostatic equilibrium state

The isothermal hydrostatic equilibrium state is represented by Eq.(13), at which the velocity is zero, and $p_{x}+\rho\phi_{x}=0$. Then, we have,

[TABLE]

Appendix  B Non-conservative discretization

As aforementioned, a portion of force term is treated as a flux term in the modified momentum equation. Actually, the proposed scheme is presented in a non-conservative form which may lead to wrong shock speed. That’s why we only consider the smooth flow in this study. However, the reader might curious about the performance of the proposed scheme on simulating the supersonic flow with shocks.

In order to estimate the effect of our non-conservative scheme, we compare the numerical results of the Sod shock tube without external force. Two numerical schemes, the proposed non-conservative scheme (labeled as WB-GKS) and a conserved scheme (labeled as GKS), are employed. The vanLeer slope limiter is adopted for the interpolation of $\alpha,U,T$. The other discretizations keep unchanged. For more details of high speed GKS, please refer to Xu (2001); Chen, Xu, and Li (2016). The initial conditions are given as follows,

[TABLE]

The kinematic viscosity is $0.0$. The numerical results are shown in figure 11. The overall results are good. It can be seen that the shock speed predicted by our non-conservative scheme is almost identical with result predicted by conserved GKS. However, there is a little defect near the rarefaction wave as shown in figure 11(c).

Then a linear gravitational potential, $\phi=x+1$, is applied. The kinematic viscosity is $0.0$ and $0.0001$ for the Euler equations and the NS Equations respectively. The results at $t=0.2$ and $t=35000$ are presented in figure 12. The system eventually return to a quiescent isothermal equilibrium state as expected.

Although, Sod shock tube problem is simulated here, we strongly recommend limiting the application of the proposed scheme to low-speed continuous flow.