1 Introduction

  The BGK code distinguishes itself from other hydrocodes in that it has recourse to the physics which generates dissipation, namely the physics of particle collisions. Developed by [7, Prendergast & Xu (1993)], the hydrodynamical scheme invokes a microscopic description of gas flow and it is therefore based on considerations of gas kinetic theory. Recall that kinetic-based hydrocodes rely on the fact that the state of a gas can be described by giving the distribution function of particle velocities $f({\bf x},{\bf u},t)$ at a point in the phase space of a single particle. The codes take advantage of the fact that the quantities of hydrodynamical interest (namely the mass, momentum and energy densities in the gas) are low-order velocity moments of f. f obeys the Boltzmann-equation, which we write as $\frac{Df}{Dt} = {\frac{\delta f}{\delta t}}^{\rm coll}$, where $\frac{D}{Dt}$ is a time rate of change along the trajectory of a single particle moving freely in phase space (under the action of smoothly varying forces, if these are present) and $\frac{\delta f}{\delta t}^{\rm coll}$ is the rate of change of f due to collisions. The classical two-body collision integral is non-linear, and non-local in velocities [2, (Cercignani 1988)]. BGK (from the names of Bhatnagar-Gross-Krook) replaces this integral with the term $\frac{(g-f)}{\tau}$, where g is the Maxwell-Boltzmann distribution function having the same mass, momentum and energy densities as f, and $\tau$ is a relaxation time, which can (ideally) be as small as the mean time between collisions of a particle in the gas [1, (Bhatnagar et al. 1954)]. It is important to notice that f and g have the same mass, momentum and energy densities, but they do not give rise to the same fluxes of these quantities. Fluxes are determined by different moments of the distribution functions. BGK differs from all other kinetic-based codes in several respects: the collisions are active throughout the duration of a time step, and are not imposed as a distinct "equilibration'' process at the end of a timestep. Also, we do not guess the form of the distribution function but solve the BGK equation to find it. The BGK equation is linear in f, which makes it easy to solve if g is known. However g is not known; the "compatability'' conditions that f and g have the same mass, momentum and energy densities would determine g if f were known. Therefore the BGK formulation is also nonlinear and non-local in velocity space, as is the Boltzmann equation. It will be shown that this situation leads to a set of non-linear integral equations for the parameters of g. It might seem that we have made no progress by replacing the Boltzmann equation by the BGK plus compatability conditions; but this is not true, because it will be shown that we need solve for the parameters of g only in the neighborhood of a boundary between computational cells, and for a short time (given by the usual CFL (Courant-Friedrichs-Lewy) condition). Knowledge of the parameters of g is equivalent to knowledge of the mass, energy and momentum densities.

The connection between the BGK equation's microscopic description of gas flow and a macroscopic description has been shown by [2, Cercignani (1988)] for the case of a perfect monatomic gas and by [11, Xu (1993)] for polyatomic gases. Velocity moments of the BGK equation give the Euler equations for negligible particle collision time, $\tau$, the Navier-Stokes equation for small yet non-zero $\tau$and a description of rarefied gas dynamics for large $\tau$.In the Navier-Stokes regime, these derivations also furnish expressions for the shear and bulk viscosity coefficients and the heat conductivity coefficient.

In principle, if the local value of the collision time can be measured as a function of time for a gas, then the BGK scheme may be used to evolve the gas with true physical dissipation parameters. In practice, the resolution of the BGK scheme is limited by the grid resolution. If the grid is not fine enough to resolve a discontinuity, then artificial dissipation must be added to broaden the discontinuity so that it is at least one grid cell thick. Because viscosity and heat conductivity are proportional to $\tau$, the BGK scheme broadens shocks by enlarging $\tau$ at the location of the discontinuities. The expression for the collision time in the BGK scheme therefore contains two terms. One term is the real physical mean collision time and it is chosen according to the desired Reynolds number of the problem. The second term is chosen in such a way that shocks in the flow span at least one grid cell. The latter term tunes the amount of artificial dissipation in the scheme. The notable difference between how the BGK scheme inputs artificial dissipation and how other schemes input artificial dissipation is that the BGK scheme puts it in exactly as if it were real dissipation corresponding to the numerically necessary value for $\tau$.

We emphasize that modelling the real dissipation in an astrophysical object such as that resulting from "turbulent'' viscosity, is still outside the reach of any existing hydrocode. For one thing, values for viscosities in astrophysical objects are highly speculative. Secondly, to solve Navier-Stokes problems with a particular value of the viscosity for an astrophysical object requires grid resolutions which are beyond what is currently achievable. For laboratory scale phenomena and for small enough Reynolds number, the BGK code is successful at modelling the dissipation produced by particle collisions in a fluid. Tests of the Kolmogorov and the laminar boundary layer problem show real dissipative effects [12, (Xu & Prendergast 1994)].

Thus far the BGK code has been extensively tested without the inclusion of gravity. Several papers document the tests which verify its accuracy and robustness as both an Euler and Navier-Stokes solver in 1D and 2D Cartesian geometries and for two-dimensional adaptive unstructured grids. The list of one-dimensional Euler tests which have been performed with the BGK scheme includes: the Roe, Sod, Lax-Harten, Woodward-Colella, and Sjögreen tests for subsonic and supersonic expansion [7, (Prendergast & Xu 1993]; [11, Xu 1993]; [13, Xu et al. 1995)]. The results from these one-dimensional test cases are that the BGK scheme produces shock fronts which are typically one to two cells wide, and contact discontinuities which are slightly broader. This is competitive with high resolution codes which do not employ regridding in the neighborhood of a shock front. The BGK scheme also exhibits negligible under and over shooting even at strong shock fronts. The other notable features in the tests with rarefaction waves are the sharp corners at the junctions between the rarefaction waves and the undisturbed, uniform regions.

Many codes have trouble treating a low density region with a flow gradient. The BGK scheme successfully handles low density regions because it satisfies both an entropy condition and a positivity condition [14, (Xu et al. 1996)]. In the BGK scheme collisions relax the gas toward local thermodynamic equilibrium states and this relaxation process is accompanied by an increase in entropy. Godunov-type schemes on the other hand typically demand an entropy fix when they encounter strong rarefaction waves (cf. the Sjögreen test) otherwise they produce unphysical rarefaction shocks. Usually the fix is an addition of artificial viscosity. Because the BGK scheme naturally satisfies the entropy condition, it simply cannot generate these unphysical phenomena. In satisfying the positivity condition, the BGK scheme avoids producing states with negative density or internal energy. The Roe test is conducive to the creation of these non-physical states but the BGK scheme does not encounter them [7, (Prendergast & Xu 1993)]. The same cannot be said of the performance of conservative finite difference schemes and codes employing Riemann solvers (e.g. Roe's approximate Riemann solver).

The list of two-dimensional Euler test cases which have been performed with the BGK scheme includes: (a) uniform Mach 3 flow in a tunnel with a forward facing step (the Emery test). The BGK code is not modified in any way near the step to treat the flow past it and its corner, which is a singular point. With the BGK scheme expansion shocks never emerge from the corner [7, (Prendergast & Xu 1993]; [11, Xu 1993]; [13, Xu et al. 1995)], (b) double Mach reflection in supersonic flow over a wedge, (c) the diffraction of a strong shock (${\rm Mach} = 5.09$) around a corner. Test cases (b) and (c) are further examples of the BGK scheme succeeding at simulating flow without having to summon any detection algorithms or entropy fixes. According to [13, Xu et al. (1995)], the original Godunov scheme, the Roe scheme without the entropy fix, and the Osher scheme could produce a rarefaction shock at the corner in test case (c); (d) flow around an impulsively started cylinder. When applied to this problem, many schemes either fail or have severe problems in maintaining positive pressure and density in the near vacuum low pressure and low density region created behind the cylinder. The BGK scheme seems to be able to preserve positivity without ad hoc fixes, and to reach a steady state solution for the problem [14, (Xu et al. 1996)].

To test its performance as a Navier-Stokes solver and to see if it has real viscosity effects, the BGK scheme has been applied to the laminar boundary layer problem and the Kolmogorov problem [11, (Xu 1993)]. The laminar boundary layer problem models the flow of gas above a flat plate. Even on coarse grids (e.g. (32 X 16), (16 X 8)), the BGK scheme impressively recovers the Blasius profile. In the Kolmogorov problem, a one-dimensional sinusoidal velocity field is imposed in a uniform density and isothermal fluid. The BGK code fulfills the Navier-Stokes prediction which is that the shape of the fluid's velocity profile remains unchanged while the amplitude of the velocity decreases in such a way that the fluid's kinetic energy decays exponentially. The agreement between the theoretical viscosity coefficient and the numerical viscosity coefficient (deduced from the measured decay rate) is excellent.

To improve the resolution of physical discontinuities occurring in complicated flows, a version of the BGK scheme on a two-dimensional adaptive unstructured grid has now been developed [6, (Kim & Jameson 1998)].

Most astronomical applications of a hydrocode require a consideration of gravity. To this end we have been developing the BGK scheme. Recently the BGK scheme has been used for a cosmological simulation [10, (Xu 1997)] but this is prior to results showing the long-term stability of the BGK scheme with gravity and the BGK scheme's convergence to the equilibrium state with gravity. In this paper we give such results. In addition to the incorporation of gravity, we have modified the geometry of the Eulerian grid to axisymmetric cylindrical coordinates. Changing the geometry of the grid has effects similar to the effects of adding gravity to the BGK scheme. It gives rise to source terms in the hydrodynamic equations. For the purpose of simplicity, in this paper we describe modifications to the BGK method when a time-independent gravitational potential is incorporated into a BGK scheme in Cartesian coordinates (Sect. 2). We present results however for both gas flow on a one-dimensional Cartesian grid and on an axisymmetric cylindrical grid (Sect. 3).