Kinetic Schemes for compressible flows

Kinetic schemes are numerical methods for the solution of the equations of fluid dynamics. They are derived from the Boltzmann equation of kinetic theory. Note that we are not talking about solving the Boltzmann equation itself, but only of using it to derive numerical schemes for the macroscopic fluid variables like density, fluid velocity and pressure. Of course the solution of Boltzmann equation is also of practical importance, for example, in highly rarefied situations, where the continuum hypothesis is not valid and it is more appropriate to use a statistical particle model.

Kinetic theory
The Boltzmann equation is an evolution equation for the one-particle velocity distribution function f. The distribution function f is defined such that the numbers of particles in small volume dxdydz around (x, y, z) having velocities in the range (u, v, w), (u + du, v + dv, w + dw) at time t is given by

$$f(t,x,y,z,u,v,w)\textrm{d}x \textrm{d}y \textrm{d}z \textrm{d}u \textrm{d}v \textrm{d}w$$

The total number of particles in the volume dxdydzdudvdw can change because

• some particles move out or move into the physical volume dxdydz because of convection
• some particles move out or move into the velocity volume dudvdw because of collisions

The mathematical expression of the above statement of conservation of particles leads to the Boltzmann equation

$$\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial ... ...partial f}{\partial y} + w\frac{\partial f}{\partial z} = J(f)$$ (1)

where J(f) is the term which gives the rate at which particles move in or out of velocity space due to collisions, and its form depends on the collision model.

The total number of particles in the volume dxdydz is obtained by integrating over all the velocities

$$(\textrm{d}x \textrm{d}y \textrm{d}z)\int\!\!\!\int\!\!\!\int f \textrm{d}u \textrm{d}v \textrm{d}w$$

so that the number density of particles is given by

$$n = \int\!\!\!\int\!\!\!\int f \textrm{d}u \textrm{d}v \textrm{d}w$$

Multiplying this by the particle mass gives the mass density. From now on we will assume that f is already multiplied by the particle mass so that the above equation directly gives the density, i.e.,

$$\rho = \int\!\!\!\int\!\!\!\int f \textrm{d}u \textrm{d}v \textrm{d}w$$

Similarly the momentum density is given by

$$\rho U = \int\!\!\!\int\!\!\!\int uf \textrm{d}u \textrm{d}v \textrm{d}w$$

with similar expressions for the other two components. Finally the energy density is given by

$$\rho E = \int\!\!\!\int\!\!\!\int \frac{u^2 + v^2 + w^2}{2} f \textrm{d}u \textrm{d}v \textrm{d}w$$

where we have assumed that the particles exchange only kinetic energy during collisions. Other forms of energy that may have to considered are rotational and vibrational energy. Writing the conserved vector as

$$\varphi = \left[ \begin{array}{c} \rho \rho U \rho V \rho W \rho E \end{array} \right]$$

and the vector of collisional invariants

$$\Psi = \left[ \begin{array}{c} 1 u v w \frac{u^2 + v^2 + w^2}{2} \end{array} \right]$$

we see that the distribution function satisfies

$$\varphi = \int\!\!\!\int\!\!\!\int \Psi f \textrm{d}u \textrm{d}v \textrm{d}w$$

Since mass, momentum and energy are conserved in a collision, the collision term must satisfy the following identity

$$\int\!\!\!\int\!\!\!\int \Psi J(f) \textrm{d}u \textrm{d}v \textrm{d}w = 0$$

The Maxwell-Boltzmann distribution function which represents the state of thermodynamic equilibrium is given by

$$f = \rho \left( \frac{\beta}{\pi} \right)^{3/2} \exp\left\{ -\beta(u-U)^2 - \beta(v-V)^2 - \beta(w-W)^2 \right\}$$ (2)

and is obtained as the solution of

$$J(f) = 0$$

With the Maxwell-Boltzmann distribution the convective terms give rise to the inviscid fluxes, i.e,

$$\int\!\!\!\int\!\!\!\int u\Psi f \textrm{d}u \textrm{d}v \te... ...ho UV \\ \rho UW \\ (\rho E + p )U \end{array} \right] = F_x$$

and so on.

Hence if we take moments of the Boltzmann equation assuming that the distribution function is Maxwellian, then we obtain the Euler equations of inviscid gas dynamics.

The Maxwellian distribution given above applies to monatomic gases only which have a specific heat ratio of $\gamma=5/3$. In order to handle gases with other values of $\gamma$ we have to introduce an additional internal energy variable $I$ and the Maxwellian in this case reads

$$f = \frac{\rho}{I_o} \left( \frac{\beta}{\pi} \right)^{3/2} ... ...\{ -\beta(u-U)^2 - \beta(v-V)^2 - \beta(w-W)^2 - I/I_o\right\}$$ (3)

Kinetic Schemes
The basic idea of kinetic schemes is the following: since we know the connection between the Boltzmann equation and the Euler equation, and since the Boltzmann equation is a single linear, convection equation (collision term is still non-linear), we first discretize the Boltzmann equation using an upwind scheme; taking moments will then give a discretization of the Euler equations. Let us see this in one-dimension. The Boltzmann equation is

$$\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} = J(f)$$

At time level n the distribution is Maxwellian, so that J(f)=0. The equation is a scalar convection equation and a stable discretization is obtained by using backwards differencing when $u \ge 0$ and forward differencing when $u \le 0$, i.e.,

$$\frac{\bar{f}^{n+1}_j - f^n_j}{\Delta t} + \frac{u + \vert u... ...c{ u - \vert u\vert}{2} \frac{f^n_{j+1} - f^n_j}{\Delta x} = 0$$

This is called the free-flow step. The conserved variable is obtained by taking moments

$$\varphi^{n+1}_j = \int^\infty_0 \int^\infty_{-\infty} \Psi \bar{f}^{n+1}_j \textrm{d}u \textrm{d}I$$

and the distribution function is instantaneously relaxed to the local Maxwellian (collision step)

$$f^{n+1}_j = f(\varphi^{n+1}_j)$$

and this cycle is repeated until the desired time-level is reached. The update equation at the Euler level can be obtained explicitly by performing the integrations,

$$\frac{ \varphi^{n+1}_j - \varphi^n_j }{\Delta t} + \frac{ F^... ...\frac{ F^-(\varphi^n_{j+1}) - F^-(\varphi^n_j) }{\Delta x} = 0$$ (4)

where the split fluxes are given by

$$F^+ = \int^\infty_0\!\!\!\int^\infty_0\!\!\! u \Psi f \textrm{d}u \textrm{d}I$$

$$F^- = \int^\infty_0\!\!\!\int^0_{-\infty}\!\!\! u \Psi f \textrm{d}u \textrm{d}I$$

Kinetic schemes have been found to be very robust; they satisfy the entropy condition and hence do not require any entropy fix. They have been shown to be positivity preserving, i.e., under a CFL restriction the update always yields non-negative values of density and pressure. They do not suffer from the pathologies that plague Riemann solvers like entropy violating shocks, carbuncle phenomenon, kinked Mach stem, etc. They are also trouble-free when solving hypersonic flows and resolve shocks with an accuracy that matches flux differencing schemes. The robustness of kinetic schemes indicates that they have a high dose of numerical viscosity and this is borne out by comparing with Riemann solvers like that of Roe. It can also be seen by looking at the kinetic scheme for a scalar convection equation as will be done later.

It is not necessary to use the Maxwellian distribution only. Any distribution that gives the correct moments can be used and some authors have proposed compact distribution functions. Of course in this case it becomes a purely mathematical approach and the connection with kinetic theory is lost.