key: cord-0045912-3rg0qcb5
authors: Koellermeier, Julian; Samaey, Giovanni
title: Projective Integration for Moment Models of the BGK Equation
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50433-5_25
sha: 2118bc7f6bf90299b3a683d64a1668ac2bc9f8c7
doc_id: 45912
cord_uid: 3rg0qcb5

In this paper, we apply the projective integration method to moment models of the Boltzmann-BGK equation and investigate the numerical properties of the resulting scheme. Projective integration is an explicit, asymptotic-preserving scheme that is tailored to problems with a large spectral gap between fast and slow eigenvalues of the model. A spectral analysis of the moment model shows a clear spectral gap and reveals the multi-scale nature of the model. The new scheme overcomes the severe time step constraint of standard explicit schemes like the forward Euler scheme by performing a number of inner iterations and then extrapolating the solution forward in time. The projective integration scheme is non-intrusive and yields fast and accurate solutions, as demonstrated using a 1D shock tube test case. These observations open up many possibilities for further use of the scheme for high-resolution discretizations and different collision models.

The Boltzmann BGK equation is widely used to model applications in kinetic theory, such as rarefied gases [15, 19] . It is a high-dimensional equation for the evolution of the particle density function that describes the distribution of particles in position-velocity phase space. An efficient discretization of the BGK equation is necessary due to its high-dimensionality. A simple mesh-based discretization in the velocity dimensions, the so-called Discrete Velocity Method (DVM) typically requires a large number of variables and thus leads to a large system of equations [1, 15] . Another approach is the discretization by means of moments [19, 23] . Moments are integrals of some quantities of interest with respect to the distribution functions over the velocity dimensions. These moment models lead to smaller PDE systems and allow for more physical insight.

Both for the moment models and the DVM, the resulting PDE models are in balance law form and typically contain a stiff right-hand side governed by a smallness parameter . For small values of , which correspond to flow conditions close to the continuum regime, suitable time-stepping schemes need to be applied to overcome the severe time step constraint posed by the stiff right-hand side.

Implicit time stepping methods can in principle be used with a very large time step, but they are costly because they require the solution of a large non-linear system of equations in each time step. As they use the information of all spatial grid points for the update of each single unknown, this method may not reflect the properties of the underlying BGK equation, which is a hyperbolic transport problem describing a directed transport of information in our application.

Explicit time stepping methods on the other hand are fast and straightforward to apply but typically suffer from a severe time step constraint when the stiff right-hand side function features multiple scales. The stability condition then poses a more severe restriction to the time step size than the typical CFL constraint, which might otherwise be sufficient to acquire a desired level of accuracy.

Several methods have been developed to overcome these problems. For example, IMEX Runge-Kutta methods split the right hand side in a stiff and a non-stiff part and apply the respective explicit and implicit solvers to each of the parts to allow for an efficient solution [5, 16] . However, they need to be tailored to the specific form of the right-hand side. In the context of discrete velocity methods, we refer to [22] for an efficient method to evaluate the full Boltzmann operator and [20, 21] which both make use of a splitting into a relaxation and a transport stages for near-equilibrium flows.

In this paper we will use projective integration, a fully explicit, potentially high-order scheme that is asymptotic preserving [14] . The asymptotic preserving property means that the time step constraint of the projective integration scheme does not depend on the stiffness parameter in the limit where tends to zero [8] . Instead, a standard CFL-type time step constraint will be sufficient to allow for accurate and fast solutions. The scheme is furthermore non-intrusive and can be used with different collision models and spatial discretizations.

In this paper we present the first application of a first order scheme using projective integration for moment models. We thus combine the efficient discretization in velocity space using moment models with the asymptotic preserving discretization in time using projective integration. This is a first step towards high-order schemes in future work.

The rest of this paper is organized as follows: in Sect. 2 we present a fully non-linear, hyperbolic moment model and a linearized version thereof. In Sect. 3 we analyse the spectrum of the linearized, semi-discrete system to determine the parameters of the projective integration scheme in the following Sect. 4. Simulation results for a shock tube test case are presented in Sect. 5. The paper ends with a short conclusion.

The Boltzmann BGK equation is used to model rarefied gases. It describes the evolution of the particles' mass density distribution function f (t, x, c), where

x, c ∈ R d denote the physical space and the velocity space, respectively. For the rest of this paper we consider the 1D case (d = 1) so that the kinetic equation

The right-hand side operator models collisions. Several collision models are possible, but we use the simplified BGK model [2] , which describes a relaxation with relaxation time ∈ R + towards the local Maxwellian f M (t, x, c) given by

The macroscopic quantities density ρ(t, x), bulk velocity u(t, x) and temperature θ(t, x) are so-called moments of the distribution function f (t, x, c) in velocity space and they are computed via integration over velocity space

For vanishing , Eq. (1) only allows solutions that are in the kernel of the righthand side collision operator. These solutions are given by the standard compressible Euler equations for ρ, u and θ. The dynamics of the kinetic equation is governed by the relaxation speed with propagation speeds proportional to 1 , as we will see in Sect. 5, whereas the macroscopic quantities might evolve on a much slower time scale. We are interested in numerical schemes that allow for an efficient solution in the limit → 0.

Equation (1) is difficult to solve because of the additional microscopic velocity dimension. A direct discretization of the microscopic velocity variable, the socalled Discrete Velocity Method (DVM) is possible, but leads to a large system of equations [1, 15] . A more efficient choice of variables is the discretization via moments leading to a smaller set of moment equations [23] .

The distribution function f is therefore expanded in a series of basis functions multiplied with basis coefficients around some equilibrium state. In this paper the expansion will be performed around local equilibrium (2) first, before we consider a linearized version in the next section. In [7] the expansion uses a series of Hermite basis functions

for coefficients (also called moments) f α (t, x), α ∈ [0, M], and weighted Hermite basis functions H

We note that the coefficients encode information about how much the distribution function deviates from the equilibrium distribution function, which is given by the Maxwellian in Eq. (2) . The more coefficients are used, the better the deviations can be represented in the weighted Hermite basis. It is the goal of a moment model to use relatively few variables in comparison to a standard discretization of the velocity space while still allowing for high model accuracy. It can be attributed to the fact that the coefficients are a more clever choice of variables than normal point values of the distribution function.

The macroscopic variables as defined in Eq. (3) result in the following additional constraints that can be easily applied

Substituting expansion (6) into the kinetic equation (1) a system of equations can be derived by either matching coefficients of the basis functions [6] or by multiplying with a set of test functions (also called testing) [10] . The result is an explicit set of PDEs in space and time for the unknown coefficients f α of the expansion and the additional macroscopic variables. This set of PDEs is called the moment model.

The whole set of variables reads

. Using this definition the moment system can be given by

with system matrix A (w M ) ∈ R (M +1)×(M +1) depending on the variables. The vector S (w M ) ∈ R M +1 results from the corresponding discretization of the right-hand side collision operator [12] . Following the derivation in [6] , the so-called Grad model's system matrix A reads

where all other entries in the system matrix A are set to zero.

The right-hand side S is given by

Unfortunately, Grad's system (9) is not hyperbolic, as shown in [3, 4] . Numerical simulations using this model might suffer from instabilities and break down due to nonphysical values [3, 12] . New hyperbolic models have recently been developed that overcome these problems and make an efficient solution of Eq.

(1) possible. In this paper, we will use the Quadrature-Based Moment Equations (QBME), first developed in [10] . The model was derived using a different framework in [6] and was extended to the multi-dimensional case in [11] . It requires only a small modification of the Grad model in (10) to ensure hyperbolicity while still preserving the conservation laws and the most important physical properties of the model. For more information we refer to [9] .

In comparison to Grad's model, the QBME model adds regularization terms in the last and the second to last equation

Due to the additional terms, the two changed equations can no longer be written in conservative form. However, the effect of the non-conservative form on the solution was investigated in [3, 12] and it was found that stable and accurate non-equilibrium solutions can be computed despite the non-conservative form. In addition, we expect no problems for the equilibrium case → 0, as the nonconservative terms then vanish.

The expansion (6) is performed around local equilibrium (ρ, u, θ), which leads to a non-linear model. A simpler, linear model can be derived when considering the expansion around a Gaussian distribution function instead. The expansion then reads

where the weighted Hermite basis functions H α are defined as

where He α is the standard Hermite polynomial of degree α and the last factor is chosen for normalization of the basis functions.

Similar to the non-linear model, the following constraints hold for the linear model

The linear moment model is then derived in the same way as the non-linear model, by testing Eq. (1), and can be written as

using a constant system matrix A ∈ R (M +1)×(M +1) given by

The right-hand side vector S (w M ) ∈ R M +1 is given by

using the ansatz from (15) and can be computed analytically beforehand. We omit the details of the derivation here for conciseness.

In the case of the simple BGK collision operator, a first-order time splitting scheme can be easily implemented. In a splitting scheme the computation of a single time step with step size Δt = t n+1 − t n is formally split into two separate steps that are performed each after the other as follows

so that the first step solves only the hyperbolic transport part of the PDE and the second step solves only the relaxation with the right-hand side collision term. The second step can be solved exactly for the simple BGK model, both in the linear and in the non-linear case. We refer to [9] for more details. However, a splitting scheme highly depends on the specific form of the right-hand side operator and can become difficult for different (and more realistic) collision operators. Furthermore, a splitting scheme is not easily extendable to higher-order accuracy, which is a significant disadvantage if high-order solutions are necessary. In this paper we will thus use the splitting scheme only as a reference method to compare it to our new projective integration scheme for moment models. For this first application we will use a standard first-order splitting (see [12] for more details) to compare it with the first-order projective integration method.

Projective Integration (PI) is an asymptotic preserving time stepping scheme consisting of an inner integrator and an extrapolation step [14] . The first order PI scheme uses the standard forward Euler method as inner integrator, but higher order schemes using Runge-Kutta methods have been derived, see [13] . For the description of the scheme we use that the semi-discrete version of the model equation after spatial discretization is given by

where the term D x is the result of the spatial discretization and the second term represents the collision operator. Note that this form does not rely on any special form of the collision operator and it can also be used for models other than the BGK model. As inner integrator we use the explicit forward Euler scheme with time step size δt for K + 1 steps

After the K + 1 inner steps a discrete derivative using the last two values is obtained and used in an outer step to compute the value at the new time step 

The dynamic behavior of the semi-discrete system in Eq. (22) is governed by the eigenvalue spectrum of the right-hand side function D t (w M ). In [14] a spectral analysis was performed for the similar DVM model and a spectral gap could be shown analytically. Here, we do something similar numerically for the moment model. For the linear moment model from Eq. (18) using M = 5 and the UPRICE spatial discretization for 400 spatial discretization points, Fig. 1 shows the spectrum of the right-hand side operator. The spectrum shows a clear spectral gap, ideally suited for projective integration. The spectral gap increases with smaller , such that the time step constraint for a standard explicit Euler method becomes more and more severe. However, projective integration can be used to overcome that time step constraint. Following the stability criterion of the projective integration scheme [14] and the eigenvalue spectrum of the moment model, we can choose the inner time step size as δt = . Furthermore, we choose K = 2 so that three inner time steps are performed before the extrapolation step.

In this paper projective integration is used to speed up simulations of moment models close to equilibrium, where the stiffness of the model equation would normally require an extremely small time step size. Here we want to derive an estimate of the speedup factor when using projective integration in comparison to a standard time stepping method. A standard time stepping method needs to resolve the fast eigenvalues, the time step size is thus Δt ∼ . This method then needs 1 Δt ∼ 1 time steps to compute the solution over a unit time interval.

The projective integration method only does K + 1 small inner time steps and then extrapolates over the remainder of a standard CFL time step size Δt ∼ Δx. The projective integration method then needs K+1 Δt ∼ K+1 Δx time steps to compute the solution over a unit time interval.

Neglecting the computational overhead of the extrapolation step, the speedup S can be computed by the ratio of the respective number of time steps that need to be performed. For the projective integration method, we get S = time steps standard method time steps projective integration = Δx

In the test cases below, we used Δx = 0.001, K = 2 and up to = 10 −5 . In this case, the speedup can thus be estimated as S = 100 3 . This is mainly due to the larger time step of the projective integration method, which uses an outer time step size Δt ≈ 0.001 versus the standard method, which uses the finer Δt = 10 −5 .

We note that the speedup will be significantly higher for smaller values of or coarser spatial discretization, i.e. larger Δx. This will be possible when using higher-order spatial discretization.

For the numerical tests in this paper we consider a 1D shock tube test case, which is a standard benchmark problem in rarefied gases, see [3, 12] . The shock tube features a strong shock wave propagating forward. Close to the shock, the solution will be in non-equilibrium if the relaxation time is large. However, for small relaxation time, the solution will quickly relax to the equilibrium Maxwellian and in the limit it can be derived easily by the well known Euler equations. It is important for any model, to correctly predict the limit of vanishing relaxation time, as large parts of any simulation will typically feature equilibrium conditions. This is exactly the region in which the kinetic Boltzmann-BGK equation becomes stiff and is difficult to solve. We are thus interested in a speedup of moment models for simulations close to equilibrium like in this test case.

At t = 0, the gas is in exact equilibrium, whereas the density, velocity, and temperature are given by

modeling a jump in density at the discontinuity at x = 0. The computational domain is [−2, 2]. The simulations run until t END = 0.3 and Δt = 0.0001 corresponds to a CFL number of approximately 0.45 on a spatial grid discretized with 4000 cells as used in [12] . The first order, path conservative UPRICE method from [12] is used for the spatial discretization, see [17, 18] . The scheme needs the eigenvalues of the system matrix, but not the eigenvectors. It is therefore computationally more expensive for the non-linear model from Sect. 2.1, as the eigenvalues change throughout the simulation. The linear model from Sect. 2.2 has a constant system matrix and thus only requires one eigenvalue computation in the beginning of the simulation. Note that higherorder spatial discretization methods are available (for examples see [9] , but we will focus on the first order scheme for the first application of the new projective integration method.

According to the spectral analysis in Sect. 5, we will use δt = for the inner time step size use K = 2 for the number of inner time steps.

First we verify that the projective integration method results in the same solution as the reference method, which uses the splitting scheme. For both schemes we use the same HSM model (18) with M = 10 on a grid using 4000 cells and outer time step size Δt = 0.0001, so that the only cause of difference could be the time stepping method. Figure 2 shows the results for both = 10 −2 and = 10 −5 . There is no visual difference between the reference method and the projective integration method for both values of . With the help of the projective integration method we can thus get a reliable solution for smaller and smaller values of . In Fig. 3 , we compare the different solutions and can see that the limit approaches the shock structure of a standard Euler equations model. This has been shown for the discrete velocity model in [14] . It is important to note that the runtime of the new projective integration scheme for the HSM moment model does not depend on and the scheme is thus asymptotic preserving in terms of computational cost. A standard explicit forward Euler scheme would lead to an increasing number of time steps for smaller values of to ensure stability. 

The projective integration method can readily be extended towards using the non-linear QBME model in (9), as it requires no special treatment of the righthand side collision operator. In this case the non-linearity of the system matrix is not a problem, as the system is close to equilibrium. We choose the standard five moment case M = 4 [3, 9] and use the same settings as for the linear model from before to obtain stability from the previous linear stability analysis. The solutions in Fig. 4 show that the QBME model gives the same solution as the HSM model for this range of the parameter despite the simplicity of the linear model. This is due to the small value of that leads to solutions very close to equilibrium with no significant rarefaction effects. Due to that, the linear model yields accurate solutions and can be used equivalently for the more complex, non-linear model. We note that this is not the case for large values of , but this is not the regime of interest for the scope of this paper. For simulations further in the kinetic regime, we refer to [12] .

We showed the first application of an explicit, asymptotic preserving scheme for moment models based on projective integration. A linear stability analysis showed that the projective integration scheme is ideally suited for the eigenvalue structure due to the clear spectral gap. The spectrum is essentially the same as in the case of a discrete velocity model and the same parameters can be used. The projective integration scheme was validated in comparison to a splitting scheme that is specifically tailored to the type of the collision operator. However, the projective integration scheme is non-intrusive and does not depend on the implementation of the right-hand side collision operator. It can thus be used in further tests and applications using a more advanced collision operator. We have consistently extended the projective integration scheme to the non-linear moment model and achieved fast and accurate solutions. This allows for further applications of the scheme to full high-order simulations using the non-linear model.

Locally refined discrete velocity grids for deterministic rarefied flow simulations

A model for collision processes in gases. 1. Small amplitude processes in charged and neutral one-component systems

Globally hyperbolic regularization of Grad's moment system in one dimensional space

On hyperbolicity of 13-moment system

Numerical methods for kinetic equations

Model reduction of kinetic equations by operator projection

On the kinetic theory of rarefied gases

Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review. Rivista di Matematica della Università di Parma

Derivation and numerical solution of hyperbolic moment equations for rarefied gas flows. Dissertation

A framework for hyperbolic approximation of kinetic equations using quadrature-based projection methods

Hyperbolic moment equations using quadraturebased projection methods

Numerical study of partially conservative moment equations in kinetic theory

A high-order asymptotic-preserving scheme for kinetic equations using projective integration

Projective Integration for nonlinear BGK kinetic equations

Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics

Implicit-explicit runge-kutta schemes and applications to hyperbolic systems with relaxation

Upwind-biased FORCE schemes with applications to free-surface shallow flows

A finite volume upwind-biased centred scheme for hyperbolic systems of conservation laws: application to shallow water equations

Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory. Interaction of Mechanics and Mathematics

Solution of the boltzmann equation in stiff regime

Solution to the Boltzmann kinetic equation for high-speed flows

Conservative evaluation of boltzmann collision integral in discrete ordinates approximation

Modeling nonequilibrium gas flow based on moment equations