An Ultra-Sparse Approximation of Kinetic Solutions to Spatially Homogeneous Flows of Non-continuum Gas

We consider a compact approximation of the kinetic velocity distribution function by a sum of isotropic Gaussian densities in the problem of spatially homogeneous relaxation. Derivatives of the macroscopic parameters of the approximating Gaussians are obtained as solutions to a linear least squares problem derived from the Boltzmann equation with full collision integral. Our model performs well for flows obtained by mixing upstream and downstream conditions of normal shock wave with Mach number 3. The model was applied to explore the process of approaching equilibrium in a spatially homogeneous flow of gas. Convergence of solutions with respect to the model parameters is studied. ©2019 TheAuthors. Published by Elsevier B.V. This is an open access article under the CCBY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Non-continuum gas dynamics arises in a wide range of applications including high speed high altitude flight, reentry from space, and design of miniature electromechanical devices and those that operate in near vacuum conditions. Experimental studies of such flows are complicated due to challenges in reproducing the flow conditions in the laboratory setting, and in measuring quantities of interest. Consequently, use of numerical simulations to predict the devices' responses to interactions with non-continuum flows has become an indispensable tool in the design process. However, efficient simulation tools that can handle full complexity of real non-continuum gases remain elusive.
Non-continuum gases are best described using the kinetic approach by introducing the molecular velocity distribution ξ represents the number of molecules that are contained in a box with volume d⃗ x around point ⃗ x, whose velocities are contained in a box with volume d ⃗ ξ around point ⃗ ξ . Evolution of the velocity distribution function is governed by the Boltzmann equation, see e.g., [1]. An advantage of kinetic description is that it allows for prediction of flows at macroscopic scale by accounting for elementary interactions between molecules. A major disadvantage of kinetic models is their being multidimensional and expensive to compute. It is generally agreed, however, that a complete knowledge of the molecular distribution function f (t, ⃗ x, ⃗ ξ ) is redundant for most applications. Indeed, one only needs a handful of quantities that can be obtained from f (t, ⃗ x, ⃗ ξ ) by integration. As a result, there is an increasing interest in the development of low dimensional models that have low computational costs and can describe the physics of non-continuum gases accurately in specific regimes. A significant progress has been achieved recently in the development of low-dimensional models for the simulation of non-continuum gases using the method of moments, in particular, the Grad expansion [2][3][4]. Method of moments (MoM) methods are successfully implemented to simulate microflows [5][6][7][8], and flows with high Mach numbers [9,10]. Approaches have been proposed to guarantee hyperbolicity of the resulting macroscopic equations [11][12][13][14][15][16] and to construct well-posed boundary conditions for the macroscopic unknowns [17,18]. Models for non-equilibrium and polyatomic gases are reported in [8,19].
Our model is inspired by the well-known result of Mott-Smith on the structure of shock waves [20], and is related to the work of Kolyshkin-Ender-Ender on the method of Maxwellian expansion [21]. In particular, we approximate the kinetic solution by the sum of a small number of isotropic Gaussian densities, also known as the Maxwellian densities. The unknowns of the approximation are the number densities, bulk velocities, and temperatures of the Maxwellian streams entering the sums. Evolution equations for the unknown macro-parameters are obtained by substituting the approximation into the Boltzmann equation and reducing it to a least squares problem. The full Boltzmann collision operator with hard spheres potential is used to advance the unknowns in time. Evaluation of the collision operator is performed using the method of [22]. Our goal is to evaluate the potency of the approximation to capture non-continuum effects in spatially homogeneous flows. Hyperbolic formulations of MoM methods can be generalized to make our model suitable for simulation of spatially dependent flows, but is outside the scope of this paper.
Substitution of the sum of Maxwellian densities into the Boltzmann equation yields an overdetermined problem. We seek a solution that minimizes the discrete ∥ · ∥ 2 norm of the equation residual evaluated at randomly selected points in the velocity space. It is observed that the accuracy of the solutions depends on the choice of the velocity points for evaluation of the residual. The best result is achieved with uniform sampling of velocities from a ball around the point of bulk velocity of the solution. Accuracy also depends on the number of approximating streams and the number of sample velocity points. The derivatives of the unknown macro-parameters obtained by solving the least square problems generally do not conform to the conservation laws. Small corrections are added to the macro-parameters after each time advancement to enforce the conservation of mass, momentum, and energy in the solutions. Accuracy of the model is estimated by comparing the relaxation of moments of orders up to six to solutions to the full Boltzmann equation [22,23], and solutions to the ellipsoidal-statistical Bhatnagar-Gross-Krook (ES-BGK) model [24]. Our solutions show good agreement for flows with small deviation from continuum, and reasonable agreement for strongly non-continuum flows. The critical advantage of our ultra-sparse model is the small number of unknowns. We observe that our solutions have similar features to those of the ES-BGK model. In some cases, our model produces better approximations to solutions of the Boltzmann equation than those obtained by the model equations.

An ultra-sparse approximation in kinetic description of gas
We consider an approximation of the velocity distribution function f (t, ⃗ x, ⃗ ξ ) by a sum of isotropic Gaussian densities.
In the simplest case of a spatially homogeneous flow, the velocity distribution function is sought in the form where ) .
In this approximation, the velocity distribution function f (t, ⃗ ξ ) is given by 5k functions n i (t), ζ w i (t)), i = 1, . . . , k. Approximation (1) was introduced in [20] to obtain a closed-from solution to a normal shock wave problem. In this paper, molecular interactions are described using the Boltzmann collision operator. Hard spheres potential is used, but we note that our method is suitable for any potential.
Substituting (1) into (2) and using the chain rule, we obtain the following functional equation Hereφ(t) = d dt φ(t), whose components φ l (t) are the macro-parameters of the isotropic Gaussians ordered in the increasing index: The components of Λ l (t, ⃗ ξ ) are partial derivatives of the isotropic Gaussian densities with respect to the macroscopic variables n i (t), T i (t), ζ u i (t), ζ v i (t) and ζ w i (t) listed in their respective orders: ) .
To determine values ofφ l (t), the left and right sides of (3) are evaluated at selected points ⃗ ξ j , j = 1, . . . , r. Values r > 5k are used to guarantee well-posedness of the linear least squares problem Here Λ j,l = Λ l (t, ⃗ ξ j ) and the components Q j are the values of the Boltzmann collision operator, Q j = Q [f ](t, ⃗ ξ j ), j = 1, . . . , r, computed using the approach of [22,23,31]. The number of samples r and the algorithms to select ⃗ ξ j are parameters of the model and are discussed in the next section.
Once the values ofφ l (t) are computed, the first order Euler scheme is used to advance φ(t) to the next time step. The Euler scheme is chosen because the values ofφ l (t) are obtained through a stochastic procedure and are not expected to be smooth. Numerical experiments suggest that solutions can be made smooth by reducing the time step. Thus high order time integration schemes may be employed to produce favorable outcomes for sufficiently small time steps. For large time steps, about 1/50 of the mean free time, our solutions exhibit random oscillations but remain stable.
Computational costs of our approach are dominated by the costs of evaluating the collision operator Q [f ](t, ⃗ ξ ) at points ⃗ ξ j , j = 1, . . . , r. In this paper, nodal-DG discretizations of [22] are used to evaluate the Boltzmann collision operator.
Evaluation of the collision operator at a single point of the velocity space requires O(n 5 ) operations where n is the number of nodal velocity points in one velocity dimension. The collision operator is evaluated at r sample points. Therefore the computational cost of evaluating Q j is O(rn 5 ) operations. We note that approximation (1) does not by itself introduce discrete velocity points. Rather, a discrete mesh is introduced by the method of [22] to compute values of the collision The convolution form of the nodal-DG discretization [31] is used to evaluate Q [f ](t, ⃗ ξ ) at an arbitrary velocity point. In order to gauge the capacity of approximation (1) in the capture of non-continuum phenomena, we employ direct evaluation of the collision operator to transfer all available physical information to the model.

Conservation of mass, momentum, and energy
As is noted earlier, values ofφ(t) obtained by solving the least squares problem (4) generally do not guarantee conservation of mass, momentum, and energy in the solution. When the vector of moments φ(t) is updated using the Euler time stepping, small perturbations of the conservation laws are introduced. These perturbations accumulate and can make the solution non-physical in the period of three mean free times. To avoid accumulation of errors, a simple procedure is applied at each time step in which values of density, bulk velocity, and temperature of approximating streams are corrected so as to enforce conservation laws in the solution. We will briefly describe the procedure.
At the beginning of the time step, values of the solution density, n, bulk velocity ⃗ ζ = (ζ u , ζ v , ζ w ), and temperature T are computed analytically from the parameters of approximating Maxwellian streams and stored. After the time step, the value of density is computed again and the perturbation of density ∆n is determined. To enforce conservation of mass, values of the individual steams densities n i (t), i = 1, . . . , k are modified by adding ∆n(t)(n i (t)/n). Next, the bulk velocity of the updated solution is computed and the perturbation of bulk velocity, ∆ ⃗ ζ , is found. To enforce conservation of momentum, a correction is introduced in the bulk velocity moments of individual streams, ⃗ ζ i +∆ ⃗ ζ . Similarly, perturbations of the temperature ∆T are computed next and the temperature parameters of approximating streams are corrected, T i + ∆T . In this process, adjustments are introduced in parameters of approximating Maxwellians. As a result, positivity of solution is guaranteed in this procedure.

Algorithms for selection of sample velocity points ⃗ ξ j
It remains to describe how the sample points ⃗ ξ j are selected. Three approaches are considered for randomly generating sample velocities ⃗ v j . In the first approach, density of the distribution of the sample velocities is given by the solution f (t, ⃗ v) itself up to a normalization. This approach corresponds to randomly drawing molecules from the solution gas stream and using their velocities as the sample points ⃗ ξ j . In the second approach, the velocity values ⃗ ξ j are drawn randomly from the steady state solution which has the form of an isotropic Gaussian and the density, bulk velocity, and temperature of the solution f (t, ⃗ ξ ). In the third approach, points ⃗ ξ j are uniformly distributed in a ball with the center at the bulk velocity of the solution f (t, ⃗ ξ ) and the radius ρ √ T , where T is the temperature of the solution and ρ is a scaling parameter. Values of ρ = 2.58 and 3.20 were used in simulations. Steady state solutions are at most 0.5% of their height outside of the radius ρ √ T in the case of ρ = 2.58, and at most 0.005% of their height in the case of ρ = 3.20.
We note that the velocity sampling techniques used in this paper are not optimal. They were introduced to obtain proof of concept simulations only. To increase the method's efficiency, one wishes to replace them with an importance sampling technique that would guarantee good approximations of residual (4) and also yield well conditioned matrices Λ when possible. A promising approach to importance sampling for Boltzmann equation can be found in [32,33]. Implementation of importance sampling will be the authors future work. Additionally, one could employ quasi-stochastic sampling, see e.g., [34], which can lead to faster convergence in the case of smooth solutions.
In Fig. 1 solutions to spatially homogeneous problem (2) that use least squares system (4) and the three ways or drawing sample velocities are plotted. The initial data in these simulations is the sum of two Maxwellian streams, that is, the initial data has the form of (1). The values of the dimensionless densities, bulk velocities, and temperatures of the two Maxwellians are We note that while it is possible to compute the moments analytically using the form of solution (1), the results presented in this paper are computed numerically using nodal-DG discretizations of the velocity space native to the method for evaluating the collision operator [22]. In the results presented in Fig. 1, uniformly spaced meshes with 33 3 velocity cells were used in integration. Accuracy in evaluating the initial data was five significant digits in the parameter of density and four significant digits in temperature. We note that the obtained values of the moments are not used by the method, but are computed from solutions for visualization purposes.

Approaching continuum in the ultra-sparse model and the number of streams
An important parameter in the approximation (1) is the number k of Gaussian streams. Well known results from the neural network theory, see e.g. [35], state that as the number of approximating streams k increases and as the values of temperatures T i decrease, a very accurate approximation of the density solution can be obtained. In practice, however, we would like to limit k to a small value by sacrificing some accuracy. In simulations presented in the previous section, only two streams are needed to capture the initial data. However, as the solution starts to evolve, it is not clear if two streams are still sufficient for accurate modeling. This becomes apparent in solutions with stronger non-continuum. If one desires to introduce additional streams during the relaxation process, it is not clear how the macroscopic variables of the new streams are to be selected.
In this section we consider solutions to the problem of spatially homogeneous relaxation using k = 2 and k = 3 approximating streams. The initial data for the problem is the sum of two isotropic Gaussians. In the first problem, the dimensionless densities, bulk velocities, and temperatures of the two Gaussians are given by n 1 = 1, ⃗ ζ 1 = (1.225, 0, 0), T 1 = 0.2 and n 2 = 3, ⃗ ζ 2 = (0.408, 0, 0), T 2 = 0.733, respectively. The reference number density is N ref = 1.0E+20 and the reference temperature is T ref = 1500 K. These values correspond to downstream and upstream conditions of a normal shock wave with Mach number 3.0. In the second problem, the two Maxwellian streams are given by n 1 = 1.609,   The solutions are computed using k = 2 and k = 3 streams. In the case of k = 3, the two streams are initialized with the initial data and the third one is initialized with zero density, and the temperature and the bulk velocity are those of the steady state solution.
It can be observed that in both cases of k = 2 and k = 3 streams, the solutions using model (1)  Additional solutions of (1) with different number of approximating streams are presented in Fig. 3. These solutions are obtained using different numbers of sample points ⃗ ξ j . The sample sizes for solutions are 5kP, where P is the parameter specifying how many sample points are taken per degree of freedom. In simulations presented in Fig. 3, the numbers of samples are varied from P = 5 to P = 80. In the case of k = 3 streams, the solutions converge almost immediately and the plots corresponding to P = 10, 20, 40, and 80 all agree to at least three significant digits. Solutions for k = 2 behave differently. The case of P = 5 gives the best approximation to the Boltzmann solution. However, as the number Table 1 Convergence of selected moments in solutions corresponding to Mach 3.0 initial data. Values of moments f u,m (t), m = 2, 3, 6 are evaluated at t = 1.99 mean free times. The convergence order α is estimated using the formula α = ln(ε 1 /ε 2 )/ln(P 2 /P 1 ), were ε is the absolute value of the difference between a solution obtained for a particular value of P and the solution with the largest number of sample points, P = 80.  Convergence of solutions with respect to P is summarized in Table 1. Here, values of moments f 1,m (t), m = 2, 3, 6 at time t = 1.99 are compared for different values of P for the cases of two and three approximating Maxwellians. The rate of convergence is estimated using the formula α = ln(ε 1 /ε 2 )/ln(P 2 /P 1 ), were ε is the absolute value of the difference between a solution obtained for a particular value of P and the solution with the largest number of sample points. In Table 1, the highest number of sample points corresponds to P = 80. One can see that as the number of sample points increases, moments converge. The overall convergence rate appears to be about first order. However, asymptotic regime of convergence was not achieved for the used numbers of sample points. The oscillations in convergence rate can be caused by the sampling variance. However, they can also be caused by instabilities present in the model as will be seen in the case of Mach 6.5 initial data.
Limitations of the approximation (1) are observed in the case of stronger non-continuum. In Fig. 4-6 solutions corresponding to Mach 6.5 initial data are presented. In Fig. 4 convergence of solutions with respect to the number of samples is considered in the case of k = 2 approximating streams. It can be seen that in the case of P = 20 samples per degree of freedom, model (1) approximates solutions of the Boltzmann equation reasonably well for low order moments. However, the approximation is rather inaccurate for high order moments. Moreover, random oscillations can be observed in the plots of sixth moments for this solution. In fact, oscillations that are observed in the plots of the sixth moments are also present in low order moments, albeit with low amplitudes. These oscillations can be eliminated by dramatically reducing the time step at the expense of the model efficiency.
It can be seen in Fig. 4 that solutions using two approximating streams take a long time to converge with respect to the parameter P. The solutions appear to converge at about P = 160. Before that mark, the solutions exhibit differences that are visible on the graphs. Slow convergence in P is also observed for solutions using three approximating streams.
Not shown here, at the marks P = 80 and P = 120 the solutions still show significant differences. These results come in contrast to almost immediate convergence of the solutions in the case of Mach 3.0 initial data. An example of k = 3 Table 2 Convergence of selected moments in solutions corresponding to Mach 6.5 initial data. Values of moments f u,m (t), m = 2, 3, 6 are evaluated at t = 1.9 mean free times.  Convergence of moments of solutions at time t = 1.9 mean free times is considered in Table 2. As in the case of Mach 3.0 initial data, the solutions converge as value of P increases and the overall convergence appears to be first order. However, the estimated order of convergence fluctuates as P increases. It is possible that the number of velocity samples is still not sufficient to reduce the variance that is intrinsic to the sampling procedure. To reduce variance, other sampling techniques could be used including importance sampling and quasi-random sequences, see e.g., [34]. However, the simulations also suggest that there is an instability in the model that could be affecting convergence. A possible reason for these artifacts is that (1) cannot be used to accurately approximate solutions with strong noncontinuum. It is suspected that the process of approaching the continuum is complex and the solution may contain modes that cannot be approximated accurately by the sum of a small number of isotropic Gaussians. Instead, model (1) appears to produce solutions similar to those obtained in the ES-BGK model.

Approaching continuum in the ultra-sparse model
To illustrate the process of approaching the state of continuum in model (1), we consider time evolution of the approximating streams in solutions to the problem of spatially homogeneous relaxation with Mach 3.0 initial data. These solutions are schematically plotted in Fig. 7 as follows. One horizontal axis corresponds to the time variable measured in mean free paths. The other horizontal axis and the vertical axis represent the u and w variables of the velocity space, respectively. The component v of velocity is suppressed in the plots. However, dependencies of the presented streams in w are identical to their dependencies in v due to the symmetry in the initial data. Approximating Maxwellian streams are represented by tubes that stretch in the temporal direction. A cross section of the tube by a plane t = const has   Table 3 Macroscopic variables of the exact solution, and macroscopic variables of the dominant streams in solutions using two and three approximating streams. * In the case of three approximating streams, the relative error in density is calculated using the combined density of the dominating streams. proportional to the stream density and its color is determined by the stream temperature. In Fig. 7, colder streams appear as blue and hotter streams appear as dark red. It can be observed that during the process of relaxation, streams exchange molecules: in particular, warmer streams gain density as is indicated by the increasing radii of the tubes, while colder streams lose density as is indicated by the decreasing radii of the tubes. One can also notice that the bulk velocities of the streams change with time. The bulk velocities of the streams drift toward the bulk velocity of the steady state solution.
Similarly, the temperatures of the streams are approaching the temperature of the steady state solution. The final values of the macroscopic variables of the dominating streams are given in Table 3. It can be seen that the temperatures and bulk velocities of the dominating streams can deviate up to 3% from the temperature and the bulk velocity of the exact solution. The difference is more prominent in the case of three streams. At the same time, the macroscopic variables of the combined solution are accurate to ten digits.
Similar observations were made in [31] where the authors approximated solutions to the full Boltzmann equation by sums of isotropic Gaussians (1) using a likelihood maximization algorithm. The approximation algorithms identified two streams with the colder stream initially being depleted and the hotter stream gaining mass. This process continued for about two mean free times with the temperatures and bulk velocities of the streams drifting toward their steady state values. At the mark of eight mean free times, the approximation algorithms returned multiple copies of a single Maxwellian. Between two and eight mean free times, the approximation returned multiple streams with the bulk velocities and temperatures approaching the steady state. However, densities of the approximating streams were not monotonic in time. Again, the macroscopic variables of the identified Maxwellian streams missed their true values by 2%-5% in the case of the initial data and the steady state solution. At the same time, the L 1 norms of the approximation errors were within 1%.

Conclusion
We consider a compact approximation of the kinetic velocity distribution function by a sum of a small number of isotropic Gaussians densities. The approximation was applied to the solution of the problem of spatially homogeneous relaxation of non-continuum gas. Derivatives of the macro-parameters of the approximating streams are obtained as solutions to a linear least squares problem that is derived from the Boltzmann equation. The full Boltzmann collision operator was used in the model to advance the macro-parameters of the approximating streams in time. Our model has a very small number of degrees of freedom and therefore offers computational advantages. The obtained solutions yield good approximations to the solutions of the Boltzmann equation in the gas flows with small to medium deviations from the continuum. These solutions use 10 or 15 degrees of freedom only. In the cases of strong deviations from continuum, approximations are less accurate. Convergence of solutions with respect to time step, number of approximating streams, and size of the sample velocities used to generate the least squares problem is considered. High order moments of the solutions behave similar to those in solutions to the ES-BGK model.
An obvious limitation of the proposed model is its low accuracy in approximating strongly non-continuum flows. While the results are promising in the spatially homogeneous case, the model needs to be tested in spatially dependent flows. Generalizations of the model to spatially dependent cases are straightforward. An important question that remains unanswered is the choice of the number of approximating streams and initialization of the new streams during the time evolution. These issues will be addressed in a subsequent paper.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.