# Evaluation of the 1-D hyperbolic quadrature method of moments for non-equilibrium flows

Frédérique Laurent, Rodney O. Fox

<figure><img src="image_1.png" /></figure>

## To cite this version:

Frédérique Laurent, Rodney O. Fox. Evaluation of the 1-D hyperbolic quadrature method of moments for non-equilibrium flows. New Trends in Compex Flows, ESAIM: Proceedings, 2024. hal-04098742

HAL Id: hal-04098742

https://hal.science/hal-04098742

Submitted on 16 May 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# EVALUATION OF THE 1-D HYPERBOLIC QUADRATURE METHOD OF MOMENTS FOR NON-EQUILIBRIUM FLOWS

FRÉDÉRIQUE LAURENT<fnref n="1" /> AND RODNEY O. FOX<fnref n="2" />

**Abstract.** When considering moment methods for the resolution of the free-transport term of the 1-D kinetic equation, the hyperbolic quadrature method of moments (HyQMOM) closure introduced in [12] leads to a globally hyperbolic system of conservation equations. Here, the HLL scheme used for its resolution is first proved to be realizable, i.e., allows computed moments to remain moments of a positive measure, exploiting a new property shown for the closure. Then, the accuracy of HyQMOM is evaluated for two test cases from the literature. The first is related to rarefied gas dynamics, where the internal structure of stationary shock waves is simulated for non-equilibrium cases. The second is related to the description of a population of inertial particles, where two populations cross each other without interacting.

**Résumé.** En considérant les méthodes de moments pour la résolution du terme de transport libre de l'équation cinétique 1-D, la fermeture de la méthode de quadrature des moments hyperbolique (HyQMOM) introduite dans [12] conduit à un système globalement hyperbolique d'équations de conservation. Dans cet article, on prouve d'abord que le schéma HLL utilisé pour sa résolution est réalisable, c'est à dire qu'il permet aux moments calculés de rester des moments d'une mesure positive, en exploitant une nouvelle propriété montrée pour la fermeture. Ensuite, la précision de HyQMOM est évaluée pour deux cas tests de la littérature. Le premier est lié à la dynamique des gaz raréfiés, où la structure interne des ondes de choc stationnaires est simulée pour les cas hors équilibre. Le second est lié à la description d'une population de particules inertielles, où deux populations se croisent sans interagir.

## INTRODUCTION

Several applications can be described by a kinetic-type equation for a velocity distribution function (VDF). We can cite, for example, a gas with the Boltzmann equation, or a population of particles evolving in a gas flow with the generalized population balance equation [20]. In some cases, only a few moments of the velocity distribution, of order less than 2 or even 1, are sufficient to describe the system. This is the case for a gas near to the Maxellian equilibrium, i.e., when its Knudsen number is small enough. It is also the case for a population of particles with a small enough Stokes number [1, 7], the particle Stokes number being the ratio between the relaxation time of the particle velocity to the gas velocity and a characteristic time scale of the gas phase. Outside of these cases, e.g., moderately rarefied gases or a moderate-Stokes-number population of particles, moment methods can still be interesting compared to direct-simulation Monte-Carlo [3], which can be

<footnote n="1">¹ Laboratoire EM2C & Fédération de Mathématiques de CentraleSupélec, CNRS, CentraleSupélec, Université Paris-Saclay, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France (e-mail: frederique.laurent@centralesupelec.fr).</footnote>

<footnote n="2">² Department of Chemical and Biological Engineering, Iowa State University, 618 Bissell Road, Ames, IA 50011-1098, USA - Tocqueville-Fulbright Distinguished Chair, CentraleSupélec, Université Paris-Saclay, France (e-mail: rofox@iastate.edu).</footnote>
expensive. However, higher-order moments need to be used to describe the non-equilibrium effects for gas [21,22] or particle trajectory crossing for populations of non-interacting particles [8,30].

In general, high-order moment methods have unclosed terms. Thus, considering the free-transport term of the 1-D Boltzmann equation or of the 1-D generalized population balance equation, the equation for the $k^{th}$-order moment contains an unclosed spatial flux depending on the $(k+1)^{th}$-order moment. Broadly speaking, the closure of the latter can be accomplished (as proposed by Grad [15]) using a perturbative solution for the VDF valid near the equilibrium distribution, or using a non-perturbative reconstruction of the VDF such as entropy maximization [19,24].

Regardless of the method used to derive it, at least four properties of the obtained closure are desirable: (1) the moment set augmented by the predicted moments is realizable (i.e., they are moments on $\mathbb{R}$ of a positive measure), (2) the closure is well defined for any realizable moment set and its computation is affordable, (3) the moment system derived from the kinetic equation is globally hyperbolic, and (4) the equilibrium state is captured. The first property is automatically verified as soon as the closure is defined through the reconstruction of a positive VDF from the moments, but not the others. Regarding methods in the literature, these properties are rarely all satisfied. Grad's method satisfies them only for moments close to the ones of the Maxwellian distribution. Entropy maximization has good mathematical properties [19], but is costly to compute and it is not defined for some realizable moment sets along the so-called "Junk line" [18]. However, an affordable closure close to it was developed for the case of moment sets including moments up to fourth order [21]. Another example is the Gaussian-EQMOM closure [4], which leads to an hyperbolic system. However, both of these examples suffer from the same definition problem as the entropy maximization along the Junk line. Finally, the quadrature method of moment (QMOM) [9,23], used in the context of particle populations, is only weakly hyperbolic [5,17].

Recently, the hyperbolic quadrature method of moments (HyQMOM) [12] was designed to satisfy all the previously cited properties, reformulated and extended to arbitrary even-order moments on $\mathbb{R}$ using the method of [11]. The closure for the $(2n+1)^{th}$-order moment is directly computed from the moments up to order $2n$, without any VDF reconstruction, but ensuring the realizability of the augmented moment vector. This closure is possible for any realizable moment vector and an efficient algorithm allows to compute it. Moreover, the resulting moment system is globally hyperbolic (this was shown for $n \le 9$ and postulated for higher values of $n$) and captures the equilibrium state. However, in the original paper, a HLL scheme was used to solve the moment equations, but its realizability was not proved, i.e., the capability of the scheme to always lead to realizable moments when starting from realizable moments.

To date, HyQMOM has only been tested for the resolution of the 1-D free-transport term. If the behavior of the simulated moments was good, especially when increasing the considered moment set, the absence of source terms made this test case very difficult. Here, more realistic configurations are simulated, taken from the literature. In the context of rarefied gas dynamics, the internal structure of stationary shock waves is simulated for non-equilibrium cases [21,22]. Furthermore, in the context of particle populations, the 1-D version of the crossing-jets test case described in [28] is also considered, where two populations of inertial particles evolve in a gas with a significant strain rate and cross each other due to their inertia.

The remainder of this paper is then organized as follows. First, the HyQMOM closure is reviewed in the context of the 1-D free-transport equation. An additional property is also included, allowing then to prove in a second section, the realizability of the HLL scheme, which is used for the transport part of the equations. Finally, in the two following sections, the behavior of HyQMOM is studied for the two considered test cases.

# 1. REVIEW OF HYQMOM FOR THE 1-D TRANSPORT EQUATION

We consider the following kinetic equation for the velocity density function (VDF) $f(v; t, x)$, in 1-D, without source terms:

$$ \partial_t f + v \partial_x f = 0. \qquad (1) $$
Some source terms will be added in the result sections, depending on the application. In this section, after reviewing realizability based on orthogonal polynomials theory, the corresponding moment system is introduced and the HyQMOM closure is detailed, as well as its properties.

## 1.1. Moments and realizability definition

In this work, the moments of $f$ are considered:

$$M_k(t, x) := \int_{\mathbb{R}} v^k f(v; t, x) dv, \quad k \in \mathbb{N}. \qquad (2)$$

Since it will be useful for the results part, let us also introduce the central and standardized moments. As soon as $f$ is not zero (and then $M_0 > 0$), the central moments are

$$C_k(t, x) := \frac{1}{M_0} \int_{\mathbb{R}} (v-u)^k f(v; t, x) dv, \quad k \in \mathbb{N}, \qquad (3)$$

where $u = M_1/M_0$ is the mean velocity and in such a way that $C_0 = 1$ and $C_1 = 0$. They are the moments of $v \mapsto \frac{1}{M_0}f(v+u; t, x)$. If, in addition, $f$ has a nonzero variance, i.e., $C_2 > 0$, the standardized moments can also be introduced:

$$S_k(t, x) := \frac{C_k}{C_2^{k/2}}, \quad k \in \mathbb{N}, \qquad (4)$$

in such a way that $S_0 = 1$, $S_1 = 0$ and $S_2 = 1$. They are the moments of $v \mapsto \frac{\sqrt{C_2}}{M_0}f(u + v\sqrt{C_2}; t, x)$.

An arbitrary vector $(M_0, M_1, \dots, M_N)^t$ is not necessarily a moment vector of a non-negative function. Here moments of Radon measures are considered. The notion of realizability we use here is given in the following definition.

**Definition 1.1.** $(M_0, M_1, \dots, M_N)^t$ is said to be a **realizable moment vector** if there exists a Radon measure $d\mu$ such that $(M_0, M_1, \dots, M_N)^t = \int_{\mathbb{R}} (1, v, \dots, v^N)^t d\mu(v)$.

The set of all realizable moment vectors is called **moment space**. If a vector belongs to the interior of the moments space, it is said to be a **strictly realizable moment vector**.

The problem of the existence of such a measure $d\mu$ supported on $\mathbb{R}$, which is also the realizability characterization of $(M_0, M_1, \dots, M_N)^t$, is called the Hamburger truncated one-dimensional moment problem.

## 1.2. Orthogonal polynomials and realizability characterization

The one-dimensional moment problem is linked to the theory of orthogonal polynomials (see, for example, [14, 25]). For $M = (M_0, \dots, M_{2n})^t$, one can define the linear functional $\langle \cdot \rangle$ on the space $\mathbb{R}[X]_{2n}$ of the real polynomial function of degree smaller than $2n$ by

$$\langle X^k \rangle = M_k, \quad \text{for } k \in \{0, 1, \dots, 2n\}. \qquad (5)$$

This also allows to define the following bilinear functional on $\mathbb{R}[X]_n$: $(p, q) \mapsto \langle pq \rangle$. The strict realizability of the moments is equivalent to the definite positivity of this bilinear functional (see for example [25]), whose matrix in the basis of monomial polynomials is a Hankel matrix, in such a way that the strict realizability is characterized by the positivity of some Hankel determinants. Moreover, for a given strictly realizable moment vector, there are several possible corresponding measures. For a realizable moment vector at the boundary of moment space, some of these Hankel determinants are zero and the corresponding measure is unique, with a support having at most $N/2$ points.
For a strictly realizable moment vector $\mathbf{M}$, (.,.) is thus a scalar product. We can then define the family of monic orthogonal polynomials $Q_k$, deg($Q_k$) = $k$, for this scalar product: $\langle Q_k Q_m \rangle = \langle Q_k^2 \rangle \delta_{km}$. This family satisfies the following three-term recurrence relation [14]:

$$Q_{k+1} = (X - a_k)Q_k - b_kQ_{k-1} \qquad (6)$$

with $Q_{-1} = 0$ and $Q_0 = 1$. The coefficients $a_k$ and $b_k$ are given by the formulas:

$$a_k = \frac{\langle XQ_k^2 \rangle}{\langle Q_k^2 \rangle}, \quad b_k = \frac{\langle Q_k^2 \rangle}{\langle Q_{k-1}^2 \rangle}. \qquad (7)$$

Moreover, there is a one-to-one relation between $(b_0, a_0, \dots, b_{n-1}, a_{n-1}, b_n)$ and $\mathbf{M}$, when setting $b_0 = M_0$, and the Chebyshev algorithm [6, 14, 29] allows to compute them from the moments. This algorithm just uses (6) and the orthogonality of the $Q_k$ [29] to compute recursively the quantities $\langle X^i Q_j \rangle$.

It can be seen from (7) that, if $\mathbf{M}$ is strictly realizable, then $b_k > 0$ for $k = 0, \dots, n$. In fact, the following Theorem can be shown [25].

**Theorem 1.2.** Let $a_k$ and $b_k$ be some reals for $k = 0, \dots, n$. Let $Q_k$ be monic polynomials with deg($Q_k$) = $k$, for $k = 0, \dots, n$, defined by (6) with $Q_{-1} = 0$ and $Q_0 = 1$. Then, there exists a linear functional $\langle \cdot \rangle$ on $\mathbb{R}[X]_{2n}$ defining a scalar product $(p, q) \mapsto \langle pq \rangle$ on $\mathbb{R}[X]_n$ for which $(Q_k)_{k=0,\dots,n}$ is orthogonal if $b_k > 0$ for $k = 0, \dots, n$. Moreover, in this case, $\mathbf{M} = (\langle X^0 \rangle, \dots, \langle X^{2n} \rangle)^t$ is a strictly realizable moment vector.

This means that there is an equivalence between the strict realizability of $\mathbf{M}$ and the positivity of the $b_k$ computed from the moments by the Chebychev algorithm.

### 1.3. Moment equations and the HyQMOM closure

The 1-D moment system for the moment vector $\mathbf{M} := (M_0, M_1, \dots, M_{2n})^t$ with $n \in \mathbb{N}$ has the form

$$\partial_t \mathbf{M} + \partial_x \mathbf{F}(\mathbf{M}) = 0 \qquad (8)$$

with flux vector $\mathbf{F}(\mathbf{M}) := (M_1, M_2, \dots, M_{2n+1})^t$. Given a realizable moment vector $\mathbf{M}$, HyQMOM provides a closure, denoted $\overline{M}_{2n+1}(\mathbf{M})$ [12]. In general, HyQMOM does not use a reconstructed VDF to define the closure. However, the generalized quadrature moment of methods (GQMOM) can be used with the moments $(\mathbf{M}, \overline{M}_{2n+1}(\mathbf{M}))$ to find an $N$-node quadrature approximation for integrals involving the VDF where $N \ge n+1$ is arbitrary [10]. This provides a closure for any source term in integral form.

Instead of writing directly $\overline{M}_{2n+1}(\mathbf{M})$ from $\mathbf{M}$, the coefficients of the three-term recurrence relation of orthogonal polynomials are used. Then, choosing $\overline{M}_{2n+1}$ is equivalent to choosing $a_n$ and the Chebychev algorithm still makes the link between them.

However, for HyQMOM, the global hyperbolicity of the system in (8) is sought. Thus, let us consider the Jacobian matrix of the system $\mathbf{J}(\mathbf{M}) = \frac{DF}{DM}(\mathbf{M})$ and its characteristic polynomial $P_{2n+1}$. The HyQMOM closure was designed in such a way that the system is globally hyperbolic, i.e., $\mathbf{J}$ is diagonalizable with real eigenvalues. For that, a closure was sought to obtain a characteristic polynomial of the form $P_{2n+1} = Q_n[(X - \alpha_n)Q_n - \beta_n Q_{n-1}]$. In [12], the following theorem was then proved for $n \le 9$, but the result is conjectured to still be true for larger values of $n$.

**Theorem 1.3 (HyQMOM Closure - Hyperbolicity).** Let us assume that $\mathbf{M}$ is strictly realizable.

For all $n \in \{1, 2, \dots, 9\}$, the characteristic polynomial $P_{2n+1}$ of the Jacobian matrix of system (8) can be written as

$$P_{2n+1} = Q_n R_{n+1}, \quad R_{n+1} = (X - \alpha_n) Q_n - \beta_n Q_{n-1} \qquad (9)$$
if and only if the coefficients $\alpha_n$ and $\beta_n$ and the closure for $M_{2n+1}$, defined through the coefficient $a_n$, are related to the recurrence coefficients $(a_k)_{k \in \{0, \dots, n-1\}}$ and $(b_k)_{k \in \{0, \dots, n\}}$ by

$$a_n = \alpha_n = \frac{1}{n} \sum_{k=0}^{n-1} a_k, \quad \beta_n = \frac{2n+1}{n} b_n.$$

Moreover, the $n+1$ roots of $R_{n+1}$ are real-valued and bound and separate the $n$ (real) roots of $Q_n$, in such a way that the system (8) is strictly hyperbolic.

This then allows system (8) to be globally hyperbolic, with good behavior of the eigenvalues of the Jacobian matrix of the flux when the moment vector tends to the boundary of the moment space, i.e., when $b_n$ tends to zero (see section 1.5).

Moreover, the roots of $Q_n$ and $R_{n+1}$ are the eigenvalues of the following two Jacobi matrices, $\mathbf{J}_n$ and $\mathbf{K}_{n+1}$ respectively:

$$\mathbf{J}_n = \begin{pmatrix} a_0 & \sqrt{b_1} & & & \\ \sqrt{b_1} & a_1 & \sqrt{b_2} & & \\ & \ddots & \ddots & \ddots & \\ & & \sqrt{b_{n-2}} & a_{n-2} & \sqrt{b_{n-1}} \\ & & & \sqrt{b_{n-1}} & a_{n-1} \end{pmatrix}, \quad \mathbf{K}_{n+1} = \begin{pmatrix} a_0 & \sqrt{b_1} & & & \\ \sqrt{b_1} & a_1 & \sqrt{b_2} & & \\ & \ddots & \ddots & \ddots & \\ & & \sqrt{b_{n-1}} & a_{n-1} & \sqrt{\frac{2n+1}{n} b_n} \\ & & & \sqrt{\frac{2n+1}{n} b_n} & a_n \end{pmatrix}. \quad (10)$$

Indeed, as seen in [14], the recurrence relation (6) can be written, for $P_k = Q_k / \sqrt{\prod_{i=1}^k b_i}$:

$$\sqrt{b_{k+1}} P_{k+1} = (X - a_k) P_k - \sqrt{b_k} P_{k-1}, \quad k = 0, \dots, n-1,$$

in such a way that, for $\mathbf{P} = (P_0, \dots, P_{n-1})^t$, one has $\mathbf{X}\mathbf{P} = \mathbf{J}_n\mathbf{P} + Q_n / \sqrt{\prod_{i=1}^{n-1} b_i} e_n$, with $e_n = (0, \dots, 0, 1)^t$. And the same kind of result is shown for $R_{n+1}$, given by (9), using $\beta_n$ instead of $b_n$ in the definition of $P_n$.

Let us denote $\lambda_{min}(\mathbf{M})$ and $\lambda_{max}(\mathbf{M})$ the minimal, respectively maximal, root of the characteristic polynomial $P_{2n+1}$ associated to the moment vector $\mathbf{M}$. They are, in fact, the minimal, respectively maximal, root of $R_{n+1}$ and then the minimal, respectively maximal, eigenvalues of $\mathbf{K}_{n+1}$. These eigenvalues of this symmetric tridiagonal matrix are computed by the implicit QL method, using a routine of EISPACK [13].

## 1.4. Additional property of the HyQMOM closure

Let us add a new result giving an interesting property of a measure that corresponds to the HyQMOM closure.

**Theorem 1.4 (Property of a representing measure).** Let us assume that the $\mathbf{M}$ is strictly realizable. Then, there exists a Radon measure $\mu$ of support included in $[\lambda_{min}(\mathbf{M}), \lambda_{max}(\mathbf{M})]$ such that $\mathbf{M} = \int_{\mathbb{R}} (1, v, \dots, v^{2n})^t d\mu(v)$, where $\lambda_{min}(\mathbf{M})$ and $\lambda_{max}(\mathbf{M})$ are the minimal and maximal roots of the characteristic polynomial $P_{2n+1}$ associated with $\mathbf{M}$.

*Proof.* The moment vector augmented by the HyQMOM closure $(\mathbf{M}, \overline{M}_{2n+1}(\mathbf{M}))$ is strictly realizable and a corresponding measure $\mu$ can be taken as the Gauss quadrature, i.e., $\mu$ is the sum of $n+1$ weighted Dirac delta functions such that $(\mathbf{M}, \overline{M}_{2n+1}(\mathbf{M})) = \int_{\mathbb{R}} (1, v, \dots, v^{2n}, v^{2n+1})^t d\mu(v)$. The abscissas of this quadrature, i.e., the points of the support of $\mu$, are then the zeros of $Q_{n+1} = (X - a_n)Q_n - b_nQ_{n-1}$ [14]. However,

$$Q_{n+1} - R_{n+1} = \frac{n+1}{n} b_n Q_{n-1}.$$
The first and last zeros of $R_{n+1}$ are $\lambda_{min}(\mathbf{M})$ and $\lambda_{max}(\mathbf{M})$. Since the zeros of $R_{n+1}$ separate the zeros of $Q_n$, which separate the zeros of $Q_{n-1}$ [14], the zeros of $Q_{n+1}$ are contained in $[\lambda_{min}(\mathbf{M}), \lambda_{max}(\mathbf{M})]$. Indeed, for $x > \lambda_{max}(\mathbf{M})$, $Q_{n-1}(x)$ is positive and then $Q_{n+1}(x) > R_{n+1}(x) > 0$. And for $x < \lambda_{min}(\mathbf{M})$, $Q_{n-1}(x)$ is of the same sign as $R_{n+1}(x)$ and then $Q_{n+1}(x)$ cannot be zero. $\square$

## 1.5. Link with the QMOM closure

In the QMOM closure [23], an even number of moments is considered $\mathbf{M} := (M_0, M_1, \dots, M_{2n-1})^t$ and the equations are closed by using the corresponding Gauss quadrature. When applied to (1), this means that a closure $\overline{M}_{2n}(\mathbf{M})$ is provided for the moment of order $2n$ in such a way that $(\mathbf{M}, \overline{M}_{2n}(\mathbf{M}))$ is at the boundary of moment space, i.e., its only corresponding measure has a finite support with at most $n$ points.

There is a one-to-one relation between $(b_0, a_0, \dots, b_{n-1}, a_{n-1})$ and $\mathbf{M}$, when setting $b_0 = M_0$. Then, choosing $\overline{M}_{2n}$ is equivalent to choosing $b_n$ and for QMOM, the choice is $b_n = 0$. Moreover, the characteristic polynomial $P_{2n}$ of the Jacobian matrix $\mathbf{J}$ of the corresponding system for the moments is $P_{2n} = Q_n^2$ [5, 17] and the system is only weakly hyperbolic.

We can then remark that for the HyQMOM closure, when $b_n$ tends to zero, the moment set corresponds to a unique measure, which is the Gaussian quadrature, thus imposing the values of higher-order moments. Moreover, the characteristic polynomial $P_{2n+1}$ tends to $(X - a_n)Q_n^2$ and the eigenvalues of the system then tend to the ones of the QMOM system (plus $a_n$), leading to a good limiting behavior for the system. This, and the ability to prove that the eigenvalues are real and distinct, were the main reasons for the choice of the form of the target characteristic polynomial.

## 2. REALIZABILITY OF THE HLL SCHEME

To solve the transport part of system (8), a first-order HLL scheme is used [16, 27]. Indeed, the structure of the solution of the Riemann problem is not known, but the eigenvalues can be easily computed from the coefficients $a_k$ and $b_k$, by computing the eigenvalues of the matrices defined by (10). Moreover, their minimal and maximal values $\lambda_{min}(\mathbf{M})$ and $\lambda_{max}(\mathbf{M})$ always belong to the eigenvalues of the second matrix. The aim of this section is to show that this scheme is realizable under some CFL-like conditions, i.e., that it ensures the realizability of the obtained moments.

Let us introduce a spatial discretization in cells $[x_{j-1/2}, x_{j+1/2}]$ of size $\Delta x$ and a temporal discretization $t^p$. This finite-volume scheme provides an approximation $\mathbf{M}_j^p$ of the mean value of $\mathbf{M}$ in the $j^{th}$ cell and at time $t^p$, through the recurrence relation:

$$ \mathbf{M}_j^{p+1} = \mathbf{M}_j^p - \frac{\Delta t}{\Delta x} \left[ \mathcal{F}_{j+1/2}^p - \mathcal{F}_{j-1/2}^p \right]. \quad (11) $$

The HLL scheme uses an approximated Riemann solver, where only one intermediate state is considered. Thus, considering the Riemann problem between $\mathbf{M}_j^p$ and $\mathbf{M}_{j+1}^p$, we introduce the smallest propagation speed $S_L^{j+1/2}$ and the largest propagation speed $S_R^{j+1/2}$, given by

$$ S_L^{j+1/2} = \min\{\lambda_{min}(\mathbf{M}_j^p), \lambda_{min}(\mathbf{M}_{j+1}^p)\}, \quad S_R^{j+1/2} = \max\{\lambda_{max}(\mathbf{M}_j^p), \lambda_{max}(\mathbf{M}_{j+1}^p)\}. \quad (12) $$

The approximated solution of the Riemann problem is then

$$ \left\{ \begin{array}{ll} \mathbf{M}_j^p & \text{if } \frac{x}{t} \le S_L^{j+1/2}, \\ \mathbf{M}_{j+1/2}^* & \text{if } S_L^{j+1/2} \le \frac{x}{t} \le S_R^{j+1/2}, \\ \mathbf{M}_{j+1}^p & \text{if } \frac{x}{t} \ge S_R^{j+1/2}, \end{array} \right. \quad \text{with } \mathbf{M}_{j+1/2}^* = \frac{S_R^{j+1/2} \mathbf{M}_{j+1}^p - S_L^{j+1/2} \mathbf{M}_j^p - \mathbf{F}(\mathbf{M}_{j+1}^p) + \mathbf{F}(\mathbf{M}_j^p)}{S_R^{j+1/2} - S_L^{j+1/2}} \quad (13) $$
and the flux is given by

$$ \mathcal{F}_{j+1/2}^{p} = \begin{cases} \mathbf{F}(\mathbf{M}_{j}^{p}) & \text{if } 0 \le S_{L}^{j+1/2}, \\ \mathbf{F}_{j+1/2}^{*} & \text{if } S_{L}^{j+1/2} \le 0 \le S_{R}^{j+1/2}, \\ \mathbf{F}(\mathbf{M}_{j+1}^{p}) & \text{if } 0 \ge S_{R}^{j+1/2}. \end{cases} \quad (14) $$

There are three equivalent ways of writing $\mathbf{F}_{j+1/2}^*$, obtained by a space and time integration of (8):

$$ \mathbf{F}_{j+1/2}^* = \frac{1}{S_R^{j+1/2} - S_L^{j+1/2}} \left[ S_R^{j+1/2} \mathbf{F}(\mathbf{M}_j^p) - S_L^{j+1/2} \mathbf{F}(\mathbf{M}_{j+1}^p) + S_L^{j+1/2} S_R^{j+1/2} (\mathbf{M}_{j+1}^p - \mathbf{M}_j^p) \right] \quad (15) $$

$$ = \mathbf{F}(\mathbf{M}_j^p) + S_L^{j+1/2} (\mathbf{M}_{j+1/2}^* - \mathbf{M}_j^p) \quad (16) $$

$$ = \mathbf{F}(\mathbf{M}_{j+1}^p) + S_R^{j+1/2} (\mathbf{M}_{j+1/2}^* - \mathbf{M}_{j+1}^p). \quad (17) $$

For the computations, an equivalent flux formula can be used [28]:

$$ \mathcal{F}_{j+1/2}^{p} = \frac{1}{2} [\mathbf{F}(\mathbf{M}_{j+1}^{p}) + \mathbf{F}(\mathbf{M}_{j}^{p})] - \frac{1}{2} |S_{L}^{j+1/2}| (\mathbf{M}_{j+1/2}^{*} - \mathbf{M}_{j}^{p}) - \frac{1}{2} |S_{R}^{j+1/2}| (\mathbf{M}_{j+1}^{p} - \mathbf{M}_{j+1/2}^{*}), \quad (18) $$

The realizability of this scheme was shown in [28] for the case $n = 1$, where the closure corresponds to a Gaussian distribution. Here, a more general result is shown in a quite similar way.

**Theorem 2.1** (Realizability of HLL). The scheme defined by (11) with (14) for the HyQMOM flux $\mathbf{F}$ is realizable if $\Delta t \le \frac{\Delta x}{\lambda^+ - \lambda^-}$, where $\lambda^+ = \max\{0, (\lambda_{max}(\mathbf{M}_j^p))_j\}$ and $\lambda^- = \min\{0, (\lambda_{min}(\mathbf{M}_j^p))_j\}$.

*Proof.* Let us assume that the $\mathbf{M}_j^p$ are strictly realizable for any $j$ and let us prove that the $\mathbf{M}_j^{p+1}$ are then also realizable. In Theorem 1.4, it was shown that there exists a Radon measure $\mu_j$ of support included in $[\lambda_{min}(\mathbf{M}_j^p), \lambda_{max}(\mathbf{M}_j^p)]$ and corresponding to the moment vector augmented by the closure $(\mathbf{M}_j^p, \overline{M}_{2n+1}(\mathbf{M}))$.

First, let us prove that $\mathbf{M}_{j+1/2}^*$ is realizable. It can be written:

$$ \mathbf{M}_{j+1/2}^* = \frac{1}{S_R^{j+1/2} - S_L^{j+1/2}} \left[ \int_{\mathbb{R}} (S_R^{j+1/2} - v) (1, v, \dots, v^{2n})^t d\mu_{j+1}(v) + \int_{\mathbb{R}} (v - S_L^{j+1/2}) (1, v, \dots, v^{2n})^t d\mu_j(v) \right]. $$

The supports of $\mu_j$ and $\mu_{j+1}$ are included in $[S_L^{j+1/2}, S_R^{j+1/2}]$. Then, on their supports $S_R^{j+1/2} - v$ and $v - S_L^{j+1/2}$ are positive and $\mathbf{M}_{j+1/2}^*$ is realizable.

Next, let us then consider the different cases. If $0 \le S_R^{j+1/2}$ and $S_L^{j-1/2} \le 0$, then let us use the following values for the fluxes:

$$ \mathcal{F}_{j+1/2}^{p} = \mathbf{F}(\mathbf{M}_{j}^{p}) + (S_{L}^{j+1/2})^{-} (\mathbf{M}_{j+1/2}^{*} - \mathbf{M}_{j}^{p}), \quad \mathcal{F}_{j-1/2}^{p} = \mathbf{F}(\mathbf{M}_{j}^{p}) + (S_{R}^{j-1/2})^{+} (\mathbf{M}_{j-1/2}^{*} - \mathbf{M}_{j}^{p}) $$

where for any $y$, we define $y^+ = \max\{0, y\}$ and $y^- = \min\{0, y\}$. Then

$$ \mathbf{M}_{j}^{p+1} = \left[ 1 - \frac{\Delta t}{\Delta x} (S_{R}^{j-1/2})^{+} + \frac{\Delta t}{\Delta x} (S_{L}^{j+1/2})^{-} \right] \mathbf{M}_{j}^{p} - \frac{\Delta t}{\Delta x} (S_{L}^{j+1/2})^{-} \mathbf{M}_{j+1/2}^{*} + \frac{\Delta t}{\Delta x} (S_{R}^{j-1/2})^{+} \mathbf{M}_{j-1/2}^{*}, $$

and it is realizable as soon as $1 - \frac{\Delta t}{\Delta x} (S_R^{j-1/2})^{+} + \frac{\Delta t}{\Delta x} (S_L^{j+1/2})^{-} \ge 0$.

If $0 > S_R^{j+1/2}$ and $S_L^{j-1/2} > 0$, then $\mathcal{F}_{j+1/2}^p = \mathbf{F}(\mathbf{M}_{j+1}^p)$, $\mathcal{F}_{j-1/2}^p = \mathbf{F}(\mathbf{M}_{j-1}^p)$, so that

$$ \mathbf{M}_{j}^{p+1} = \mathbf{M}_{j}^{p} + \frac{\Delta t}{\Delta x} \int_{\mathbb{R}} (-v)(1, v, \dots, v^{2n})^{t} d\mu_{j+1}(v) + \frac{\Delta t}{\Delta x} \int_{\mathbb{R}} v(1, v, \dots, v^{2n})^{t} d\mu_{j-1}(v). $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Stationary shock. Normalized density (left) and dimensionless heat flux (bottom) for Ma = 2 (left), Ma = 4 (middle) and Ma = 8 (right): solution of the BGK equation (black solid line) and HyQMOM system with n = 2 (blue dash-dot line) and n = 3 (red solid line).</figcaption></figure>

Since $v < 0$ on the support of $\mu_{j+1}$ and $v > 0$ on the support of $\mu_{j-1}$, $M_j^{p+1}$ is realizable.

If $0 \le S_R^{j+1/2}$ and $S_L^{j-1/2} > 0$, then let us use the following values for the fluxes:

$$ \mathcal{F}_{j+1/2}^{p} = \mathbf{F}(\mathbf{M}_{j}^{p}) + (S_{L}^{j+1/2})^{-} (\mathbf{M}_{j+1/2}^{*} - \mathbf{M}_{j}^{p}), \quad \mathcal{F}_{j-1/2}^{p} = \mathbf{F}(\mathbf{M}_{j-1}^{p}). $$

Then

$$ \mathbf{M}_{j}^{p+1} = \left[ 1 + \frac{\Delta t}{\Delta x} (S_{L}^{j+1/2})^{-} \right] \mathbf{M}_{j}^{p} - \frac{\Delta t}{\Delta x} (S_{L}^{j+1/2})^{-} \mathbf{M}_{j+1/2}^{*} + \frac{\Delta t}{\Delta x} \int_{\mathbb{R}} v(1, v, \dots, v^{2n})^{t} d\mu_{j-1}(v). $$

Similarly, as for the previous cases, it can be shown that $M_j^{p+1}$ is realizable as soon as $1 + \frac{\Delta t}{\Delta x} (S_L^{j+1/2})^{-} \ge 0$.

Finally, if $0 > S_R^{j+1/2}$ and $S_L^{j-1/2} \le 0$, it can be proved in the same way that $M_j^{p+1}$ is realizable as soon as $1 - \frac{\Delta t}{\Delta x} (S_R^{j-1/2})^{+} \ge 0$.

Then, for all cases, the condition $1 - \frac{\Delta t}{\Delta x} \lambda^{+} + \frac{\Delta t}{\Delta x} \lambda^{-} \ge 0$ is sufficient, which concludes the proof. $\square$

### 3. STATIONARY SHOCK SIMULATION

In the context of rarefied gas dynamics, it is interesting to investigate the behavior of the HyQMOM closure by simulating the internal structure of stationary shock waves in non-equilibrium cases [21, 22]. Here, not considering the model error caused by the collision operator, a very simple and common collision model is chosen: the BGK or relaxation collision operator [2]. Thus, the following Riemann problem will be simulated, written here in dimensionless form:

$$ \partial_t f + v \partial_x f = \frac{f_M - f}{\text{Kn}}, \quad f(v; 0, x) = f^0(v; x) = f_L(v) \mathbf{1}_{\mathbb{R}_-}(x) + f_R(v) \mathbf{1}_{\mathbb{R}+}(x), \quad (19) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Stationary shock. Orbits of velocity moments ($S_3, S_4$) corresponding the transition and internal structure for stationary shock-wave solutions with shock Mach numbers of Ma = 2 (left), Ma = 4 (middle) and Ma = 8 (right): solution of the BGK equation (black solid line) and HyQMOM system with $n = 2$ (blue dash-dot line) and $n = 3$ (red solid line).</figcaption></figure>

where Kn is the Knudsen number, i.e., the dimensionless characteristic relaxation time to the equilibrium distribution $f_M$, which is the Gaussian distribution with the same moments of order 0, 1 and 2 as $f$. On the left, the initial condition $f_L$ is a Gaussian distribution of dimensionless density $\rho_L = 1$, mean velocity $u_L = \sqrt{3}Ma$ and variance $T_L = 1$, where Ma is then the Mach number on the left (the speed of sound is $\sqrt{3}T$ for the corresponding 1-D Euler equations). On the right, $f_R$ is also a Gaussian distribution, of dimensionless density $\rho_R$, mean velocity $u_R$ and variance $T_R$. These values are, satisfying the Rankine-Hugoniot relations for the 1-D Euler equations:

$$ \rho_R = \frac{2Ma}{Ma^2 + 1}, \quad u_R = \frac{\sqrt{3} Ma^2 + 1}{2 Ma}, \quad T_L = \frac{(3Ma^2 - 1)(Ma^2 + 1)}{4Ma^2}. \qquad (20) $$

This initial distribution represents the shock obtained for a very small Knudsen number, corresponding to the equilibrium distribution $f_M$. However, the Knudsen number is here chosen as $Kn = 0.05$, high enough so that the solution deviates from this distribution. A stationary solution is attained for a large enough time, that we will try to approximate with HyQMOM. Three cases are considered: Ma = 2, Ma = 4 and Ma = 8.

A reference solution for (19) is computed by discretizing both the physical space and velocity. A large enough support is considered for the velocity to obtain very good accuracy for the moments of $f$, up to order 4: the moments of the initial condition outside this interval are smaller than $10^{-6}$. A large enough spatial domain is also considered so that there is no influence of the boundary conditions: [-5, 5] for Ma = 2, [-7, 7] for Ma = 4, [-11, 11] for Ma = 8. Then, a MUSCL scheme is used for the flux computations, coupled with an explicit Euler method for the source term. A refined uniform discretization is considered for good accuracy with 2000 cells in the spatial direction and 200 for the velocity. The simulations are done until time $t = 15$, after which the solution no longer evolves.

The system for the corresponding moment vector $\mathbf{M} = (M_0, M_1, \dots, M_{2n})^t$ is

$$ \partial_t \mathbf{M} + \partial_x \mathbf{F}(\mathbf{M}) = \frac{\mathbf{M}^e - \mathbf{M}}{Kn}, \qquad (21) $$

where the flux $\mathbf{F}(\mathbf{M})$ is given by the HyQMOM closure, with $n \in \{2, 3, 4, 5\}$ and $\mathbf{M}^e = (M_0^e, M_1^e, \dots, M_{2n}^e)^t$ are the moments of the equilibrium Gaussian distribution $f_M$, in such a way that $M_k^e = M_k$ for $k = 0, 1, 2$. For its numerical resolution, the considered spatial domain is the same as for the reference solution and a refined discretization with 8000 cells is used, to ensure spatial convergence. The HLL flux is used, coupled with an implicit Euler method for the source term so that the CFL condition does not need to be modified to ensure
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Stationary shock. Normalized density (top) and dimensionless heat flux (bottom) for Ma = 2 (left), Ma = 4 (middle) and Ma = 8 (right): solution of the BGK equation (black solid line) and HyQMOM system with n = 4 (blue dash-dot line) and n = 5 (red solid line).</figcaption></figure>

the realizability. Indeed, the numerical scheme is then written

$$M_j^{p+1} = M_j^p - \frac{\Delta t}{\Delta x} [\mathcal{F}_{j+1/2}^p - \mathcal{F}_{j-1/2}^p] + \Delta t \frac{M_j^{ep+1} - M_j^{p+1}}{Kn}, \quad (22)$$

where $M_j^{ep+1}$ is the moment vector of the Gaussian distribution with the same moments of order 0 to 2 as $M_j^{p+1}$. Thus, we first compute the moment vector obtained by the transport $M_j^* = M_j^p - \frac{\Delta t}{\Delta x} [\mathcal{F}_{j+1/2}^p - \mathcal{F}_{j-1/2}^p]$. This gives the first three components of $M_j^{p+1}$, allowing to compute $M_j^{ep+1}$. Then, (22) can be written

$$M_j^{p+1} = M_j^* + \Delta t \frac{M_j^{ep+1} - M_j^{p+1}}{Kn}, \quad \Rightarrow \quad M_j^{p+1} = \left( M_j^* + \frac{\Delta t}{Kn} M_j^{ep+1} \right) / \left( 1 + \frac{\Delta t}{Kn} \right).$$

This allows to compute $M_j^{p+1}$ as a convex combination of the moment vectors $M_j^*$ and $M_j^{ep+1}$, thus implying its realizability.

The solution of HyQMOM with $n = 1$ is not shown, since the corresponding equations are the 1-D Euler equation and the initial solution is then the stationary solution of the system. In Figure 1, for the reference solution and HyQMOM simulations with $n = 2$ and $n = 3$, the normalized density $\frac{M_0 - \rho_L}{\rho_R - \rho_L}$ and the dimensionless heat flux $S_3$ are represented as functions of $x/\lambda$, where the dimensionless mean free path is $\lambda = \frac{16Kn}{5\sqrt{2\pi}}$, considering the same formula as in [21]. It can be seen that, for the considered Mach numbers, the HyQMOM closure with $n = 2$ is sufficient to accurately predict the density, but not the heat flux in the left part of the shock. The solution with this closure is however more accurate than the solution of the Navier-Stokes equations represented in [21], which was not able to accurately describe the density or the heat flux, not even in the right part of the shock. These equations, obtained from a Chapman-Enskog expansion in Knudsen, show their limitations in this non-equilibrium configuration, while moments methods allow to go further from the equilibrium case. The HyQMOM solution with $N = 2$ is also less accurate than the simulations with the interpolated maximum entropy closure of [21, 22], which uses the same number of moments, but the computations with HyQMOM
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Stationary shock. Orbits of velocity moments ($S_3, S_4$) corresponding to the transition and internal structure for stationary shock-wave solutions with shock Mach numbers of Ma = 2 (left), Ma = 4 (middle) and Ma = 8 (right): solution of the BGK equation (black solid line) and HyQMOM system with $n = 4$ (blue dash-dot line) and $n = 5$ (red solid line).</figcaption></figure>

are much easier to do, since it avoids the singularity of the flux near the Junk line [18] and the corresponding arising numerical difficulties [22]. The HyQMOM closure can be also used with a higher number of moments and HyQMOM results are improved by using $n = 3$, with very good accuracy in the case Ma = 8. In Figure 2, the orbits ($S_3, S_4$) are also represented for the same simulations. This describes the structure of shock waves and its non-equilibrium behaviour for the reference solution, that is tried to be reproduced with the HyQMOM simulations. It shows the accuracy of HyQMOM with $n = 3$, especially in the case Ma = 8.

In order to observe the convergence of HyQMOM when the number of moments increases, results for $n = 4$ and $n = 5$ are also presented: the normalized density and the dimensionless heat flux in Figure 3 and the orbits ($S_3, S_4$) in Figure 4. In particular, we can see the excellent agreement with the reference solution for HyQMOM with $n = 5$.

## 4. CROSSING-JETS SIMULATIONS

We reproduce here a 1-D version of the crossing-jets test case described in [28], in such a way that our spatial direction $x$ corresponds to the coordinate $y$ of [28] and our time $t$ to the coordinate $x$ of [28]. This is possible because in [28], all the systems evolve at a constant velocity in the $x$ direction and a steady state is attained on the considered domain. Crossing jets are challenging to treat with moment methods because at crossing points the distribution function is strongly bimodal even when away from such points it is nearly Gaussian.

Here, we consider two populations of particles, also called sprays, in a gas, subject to a Stokes drag force. These particles have all the same size, characterized by their Stokes number $St$, the ratio of the characteristic relaxation time to the gas velocity and the characteristic time of the gas. This number is proportional to the particle surface area. These sprays are initially located in the symmetric intervals $[-x_b, -x_a]$ and $[x_a, x_b]$, with $0 < x_a < x_b$ and have an initial Gaussian velocity distribution with variance $\sigma$ and zero mean. Let us remark that in [28], the injected particles have no velocity in the $y$ direction, which can be seen as the case $\sigma = 0$. Then, here, an agitation can be considered in this direction, as can be generated in turbulent flows.

The strain rate of the gaseous carrier field will generate particle trajectory crossing. Indeed, the gas velocity is $v_g(x) = -\epsilon x$, where $\epsilon$ is the dimensionless strain rate and the particles are inertial enough ($St > \frac{1}{4\epsilon}$) to induce crossing: they are first set in motion by the gas, and then slowed down, but not fast enough to prevent them from passing the point $x = 0$.

Then, the NDF $f(v; t, x)$ of the particles satisfies the following equation, in its dimensionless form:

$$ \partial_t f + v \partial_x f + \partial_v \left( \frac{v_g(x) - v}{\text{St}} f \right) = 0, \quad f(v; 0, x) = f^0(v; x), \qquad (23) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Number density ($M_0$) for the analytical solution of the crossing jets with $\sigma = 0$ (left) and $\sigma = 0.001$ (right).</figcaption></figure>

with the initial condition $f^0(v; x) = \lambda(x)g_\sigma(v)$, where

$$ \lambda(x) = \varphi\left(\frac{x - x_a}{x_b - x_a}\right) + \varphi\left(\frac{x + x_b}{x_b - x_a}\right), \quad \varphi(x) = 16x^2(1-x)^2 \mathbf{1}_{[0,1]}(x), \quad g_\sigma(v) = \frac{1}{\sqrt{2\pi\sigma}} \exp\left(-\frac{v^2}{2\sigma}\right). $$

In the numerical examples below, we take $x_a = 0.3$, $x_b = 0.7$, $\sigma = 0$ or $\sigma = 10^{-3}$, $\epsilon = 1$, and $\text{St} = 20$. In the case $\sigma = 0$, $g_\sigma$ corresponds in fact to a Dirac distribution.

## 4.1. Analytical solution

An analytical solution to (23) can be found by the method of characteristics. (See examples in Figure 5.) Then, let us define the characteristics $(X, U)(t; s, x, u)$, which are the solution of

$$ d_t(X, U)(t; s, x, u) = \left( U(t; s, x, u), \frac{v_g(X(t; s, x, u)) - U(t; s, x, u)}{\text{St}} \right), \quad (X, U)(s; s, x, u) = (x, u). \quad (24) $$

Denoting $\omega = \frac{\sqrt{4\epsilon\text{St}-1}}{2\text{St}}$ and $\mu(t) = \frac{\sin(\omega t)}{2\omega\text{St}} + \cos(\omega t)$, it is easy to see that

$$ \begin{pmatrix} X \\ U \end{pmatrix} (t; s, x, u) = e^{-\frac{t-s}{2\text{St}}} \begin{pmatrix} \mu(t-s) & \frac{\sin(\omega(t-s))}{\omega} \\ -\epsilon \frac{\sin(\omega(t-s))}{\omega\text{St}} & \mu(s-t) \end{pmatrix} \begin{pmatrix} x \\ u \end{pmatrix}. \quad (25) $$

Since $\mu(t)\mu(-t) = 1 - \epsilon \frac{\sin(\omega t)^2}{\omega^2\text{St}}$, the Jacobian $J(t; s, x, u)$ of the transformation $(x, u) \mapsto (X, U)(t; s, x, u)$ is

$$ J(t; s, x, u) = \exp\left(-\frac{t-s}{\text{St}}\right). $$

The solution of (23) is given by $f(v; t, x) = J(0; t, x, v)f^0(U(0; t, x, v); X(0; t, x, v))$, i.e.,

$$ f(v; t, x) = e^{\frac{t}{\text{St}}} g_{\sigma} \left( e^{\frac{t}{2\text{St}}} \left[ \frac{\epsilon \sin(\omega t)}{\omega \text{St}} x + \mu(t)v \right] \right) \lambda \left( e^{\frac{t}{2\text{St}}} \left[ \mu(-t)x - \frac{\sin(\omega t)}{\omega \text{St}} v \right] \right). \quad (26) $$

Moreover, for any function $\Phi(v, x)$, one has, with the analytical solution,

$$ \langle \Phi, f \rangle \equiv \iint \Phi(v, x) f(v; t, x) dv dx = \iint \Phi(U(t; 0, x, v); X(t; 0, x, v)) f^0(v; x) dv dx. $$
Then, with $\Phi(v, x) = v^k \frac{1}{\Delta x} \mathbf{1}_{[x_{j-1/2}, x_{j+1/2}]}(x)$, one obtains the mean value of the $k^{th}$-order moment on the cell $[x_{j-1/2}, x_{j+1/2}]$:

$$m_{k,j}(t) = \iint \frac{\mathbf{1}_{[x_{j-1/2}, x_{j+1/2}]}\left(e^{\frac{-t}{2\text{St}}} \left[\mu(t)x + \frac{\sin(\omega t)}{\omega\text{St}}v\right]\right)}{\Delta x} \left(e^{\frac{-t}{2\text{St}}} \left[\frac{-\epsilon \sin(\omega t)}{\omega\text{St}}x + \mu(-t)v\right]\right)^k \lambda(x)g_\sigma(v)dx dv. \quad (27)$$

Let us remark that the analytical solution given by (26) goes from a Gaussian distribution in $v$ when $\sin(\omega t) = 0$ ($|\mu(t)| = 1$) to a compact support distribution in $v$ when $\mu(t) = 0$. To accurately compute the moments $m_{k,j}$ given by (27), a Gauss-Hermite quadrature with 400 abscissas is then used for the velocity integral when $|\mu(t)| \ge 0.35$, whereas, when $|\mu(t)| < 0.35$, Gauss-Legendre quadratures are used on $[x_a, x_b]$ and $[-x_b, -x_a]$ for the spatial integral, with $2N+3$ abscissas on each interval. In each case, the integration over the other direction is analytical. In the case $\sigma = 0$, the Gauss-Hermite quadrature just need one abscissa.

In the case $\sigma = 0$, let us remark that $X(t; 0, x, 0) = \mu(t)e^{\frac{-t}{2\text{St}}}$ also gives the position of the particles initially at $x$. They all cancel when $\mu(t) = 0$: the trajectories of the particles cross at the same points. In fact, the analytical solution is a mono-kinetic distribution (i.e., a Dirac delta function in $v$) at each point and for any $t$, except when $\mu(t) = 0$, where the distribution is then concentrated at $x = 0$ with a continuous velocity distribution. This fact makes this case quite difficult to simulate using moment methods.

## 4.2. Results with HyQMOM

The equation for the $k^{th}$-order moment, $k \in \{0, \dots, 2n\}$, is

$$\partial_t M_k + \partial_x M_{k+1} = k \frac{v_g(x) M_{k-1} - M_k}{\text{St}}, \quad (28)$$

where $M_{2n+1}(\mathbf{M})$ is closed with HyQMOM. To solve this system, a Strang splitting [26] is used: the drag term is solved for half a time step, then the transport during a time step and finally the drag term for half a time step. Moreover, the drag operator can be solved analytically, assuming that the gas velocity is constant in each cell $j$, since it is a system of linear ODEs of the form: $d_t \mathbf{M} = A(x_j)\mathbf{M}$. The matrix exponential $\exp(A(x_j)t)$ could be used for that, but numerical error can then lead to unrealizable moments when the moment vector is initially close to the boundary of the moment space. That is why, instead, a quadrature is used, associated to the moment vector completed by the closure $(M_0, M_1, \dots, M_{2n}, M_{2n+1}(\mathbf{M}))^t$. Let us call $w_i^p$ and $v_i^p$, with $i \in \{1, \dots, n+1\}$, the weights and abscissas corresponding to this moment vector at time $t^p$ in cell $j$. Then, assuming $M_k(t) = \sum_{i=1}^N w_i^p(v_i^p)^k$, there is the following equivalence:

$$\left\{ \begin{array}{l} d_t M_k = k \frac{v_g(x_j) M_{k-1} - M_k}{\text{St}}, \\ M_k(t^p) = \sum_{i=1}^N w_i^p(v_i^p)^k, \quad k \in \{0, \dots, 2n\} \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} d_t v_i = \frac{v_g(x_j) - v_i}{\text{St}}, \\ v_i(t^p) = v_i^p, \quad i \in \{1, \dots, n+1\}. \end{array} \right.$$

Thus, after a time step, the weights are unchanged, $w_i^{p+1} = w_i^p$ and the abscissas are given by

$$v_i^{p+1} = v_g(x_j) + (v_i^p - v_g(x_j)) \exp\left(-\frac{\Delta t}{\text{St}}\right). \quad (29)$$

Since the moments at step $p+1$ are written $\sum_{i=1}^N w_i^{p+1}(v_i^{p+1})^k$, they are more likely numerically realizable.

Moreover, in most of the domain, there are no particles so that the moments are zero and in the case $\sigma = 0$, moments are at the boundary of the moment space in a part of the domain (i.e., $b_2$ is zero). To manage this, without introducing any small moments or an artificial velocity dispersion in all the domain, the HyQMOM closure algorithm is slightly modified, allowing very small modifications of the moments through a projection: considering a small quantity $\epsilon$, the moment vector is set to zero if $M_0 < \epsilon$ and when computing the $b_j$, the moment vector is considered at the boundary of the moment space if there exists $k \le n$ such that $\prod_{j=0}^k b_j < \epsilon$.
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Number density ($M_0$) for the crossing jets using HyQMOM with $\sigma = 0$ (left) and $\sigma = 0.001$ (right): $n = 2$ (first row), $n = 3$ (second row), $n = 4$ (third row).</figcaption></figure>

(Recall that on this boundary, the exact distribution consists of $k$ Dirac delta function.) In this last case, $b_j$ and $a_j$ are set to zero for $j \ge k$ and the moments are computed from the $a_j$ and $b_j$ through the reverse Chebychev algorithm. In practice $\epsilon = 10^{-16}$. This induces a small change of the moments of order greater than $2k$. The use of this procedure is necessary because the Jacobian of the transformation of the $a_j$ and $b_j$ to the moments is the product of the $b_j$, thus leading to important errors on the inverse transformation as soon as this product becomes too small. In practice, this procedure places an upper limit on the number of moments needed to describe the distribution, since, for a Gaussian distribution of variance $\sigma$, the values of the $b_k$ are $b_k = k\sigma$.

Results for $M_0$ with $\sigma = 0$ and $\sigma = 0.001$ for $n = 2, 3, 4$ are presented in Figure 6 for HyQMOM and in Figure 7 for QMOM. These can be compared to the analytical results in Figure 5. It can be seen that HyQMOM produces better results than QMOM. Recall that QMOM uses one less moment as compared to HyQMOM. For the case $n = 3$, at the first crossing point ($x = 0$) QMOM produces a nonzero weight with zero velocity that propagates towards the right, increasing in magnitude at each crossing. This behavior is not observed with HyQMOM, even for $n = 2$, where a possible reconstruction is a sum of three Dirac distribution, one at $v = 0$, due to symmetry at the crossing points. With QMOM, such behavior is observed whenever $n$ is an odd integer, the worst case being with $n = 1$ where all the weight lands on the double eigenvalues $v = 0$ at the first crossing point, producing a delta shock. In contrast, with $n = 1$ HyQMOM has three distinct eigenvalues so that no delta shock is produced. The results of this last case are not presented here since the closure then corresponds to a Gaussian reconstruction, like in [28] where this crossing-jet test case was introduced. Moreover, we can remark that, in the case $\sigma = 0$, a small part of the density stay at $x = 0$ after the third crossing for the HyQMOM simulation with $N = 4$. This is probably due to the accumulation of numerical errors in this hard test case.

Of the two presented, the case with $\sigma = 0.001$ is clearly more difficult to reproduce with moment methods. Nonetheless, the HyQMOM results are superior to QMOM, which produces distinct “packets of particles” moving with the (repeated) eigenvalues of the moment system. Similar, but less pronounced, behavior is observed with HyQMOM, which is closer to the analytical solution. We should note that increasing $n$ with
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Number density ($M_0$) for the crossing jets using QMOM with $\sigma = 0$ (left) and $\sigma = 0.001$ (right): $n = 2$ (first row), $n = 3$ (second row), $n = 4$ (third row).</figcaption></figure>

HyQMOM eventually leads to moments on the boundary of moment space. However, as noted earlier, this process occurs smoothly as the simulation progresses and the $b_j$ decrease due to fluid drag. In summary, HyQMOM provides a robust moment closure for 1-D crossing-jets simulations, and yields superior results as compared to QMOM.

## 5. CONCLUSIONS

In this work, HyQMOM and an associated realizable HLL scheme have been further tested for approximating solutions to 1-D kinetic equations. Extending prior work on the free-transport equation, which corresponds to the 1-D Boltzmann equation with infinite Knudsen number, we have shown that HyQMOM with a relatively small number of moments can accurately capture the first moments of the solution of the kinetic equation for the stationary shock problem. Interestingly, the HyQMOM solution for larger Mach numbers appears to converge relatively quickly in terms of the number of moments. In a second example, a highly non-equilibrium particle-trajectory-crossing case was considered. There, the particle dynamics are driven by the coupling to the gas phase, which depends on the Stokes number. Due to spatial transport, for sufficiently large Stokes numbers the VDF exhibits bimodal behavior at crossing points that is difficult to capture with moment methods. HyQMOM with a relatively small number of moments is shown to well capture the analytical solution. In contrast, QMOM generates delta shocks due to its weakly hyperbolic nature. Overall, the realizable HLL scheme combined with the HyQMOM closure provides a robust computational tool for solving 1-D kinetic equations. In future work, the treatment of nonlinear source terms using GQMOM [10] and multivariate versions of HyQMOM [11] will be investigated.

## REFERENCES

[1] S. Balachandar and J. K. Eaton. Turbulent dispersed multiphase flow. *Annu. Rev. Fluid Mech.*, 42:111–133, 2010.
[2] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. *Phys. Rev.*, 94:511–525, May 1954.

[3] G. A. Bird. Molecular gas dynamics and the direct simulation of gas flows. *Oxford Science Publications*, 42, 1994.

[4] C. Chalons, R. O. Fox, F. Laurent, M. Massot, and A. Vié. Multivariate Gaussian extended quadrature method of moments for turbulent disperse multiphase flow. *Multiscale Model. Simul.*, 15(4):1553–1583, 2017.

[5] C. Chalons, D. Kah, and M. Massot. Beyond pressureless gas dynamics: quadrature-based velocity moment models. *Commun. Math. Sci.*, 10(4):1241–1272, 2012.

[6] P. L. Chebyshev. Sur l'interpolation par la méthode des moindres carrés. *Mém. Acad. Impér. Sci. St. Petersburg*, 1(15):1–24, 1859. Also in œuvres I pp. 473–498.

[7] S. de Chaisemartin. *Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays*. PhD thesis, Ecole Centrale Paris, France, 2009. Available online at http://tel.archives-ouvertes.fr/tel-00443982/en/.

[8] Stephane de Chaisemartin, Lucie Fréret, Damien Kah, Frédérique Laurent, Rodney O. Fox, Julien Reveillon, and Marc Massot. Eulerian models for turbulent spray combustion with polydispersity and droplet crossing. *Comptes Rendus Mécanique*, 337(Numéro Thématique Spécial “Combustion for aerospace propulsion”):438–448, 2009.

[9] R. O. Fox. A quadrature-based third-order moment method for dilute gas-particle flow. *J. Comput. Phys.*, 227(12):6313–6350, 2008.

[10] R. O. Fox, F. Laurent, and A. Passalacqua. The generalized quadrature method of moments. *J. Aerosol Sci.*, 167:106096, 2023.

[11] R. O. Fox, F. Laurent, and A. Vié. Conditional hyperbolic quadrature method of moments for kinetic equations. *J. Comput. Phys.*, 365:269–293, 2018.

[12] Rodney O. Fox and Frédérique Laurent. Hyperbolic quadrature method of moments for the one-dimensional kinetic equation. *SIAM Journal on Applied Mathematics*, 82(2):750–771, 2022.

[13] Burton S. Garbow. Eispack — a package of matrix eigensystem routines. *Computer Physics Communications*, 7(4):179–184, 1974.

[14] W. Gautschi. *Orthogonal Polynomials: Computation and Approximation*. Oxford University Press, Oxford, UK, 2004.

[15] H. Grad. On the kinetic theory of rarefied gases. *Commun. Pure Appl. Math.*, 2(4):331–407, 1949.

[16] Amiram Harten, Peter D. Lax, and Bram van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Rev.*, 25(1):35–61, 1983.

[17] Q. Huang, S. Li, and W.-A. Yong. Stability analysis of quadrature-based moment methods for kinetic equations. *SIAM J. Appl. Math.*, 80(1):206–231, 2020.

[18] M. Junk. Domain of definition of Levermore's five-moment system. *J. Stat. Phys.*, 93(5/6):1143–1167, 1998.

[19] C. D. Levermore. Moment closure hierarchies for kinetic theories. *J. Stat. Phys.*, 83:1021–1065, 1996.

[20] D. L. Marchisio and R. O. Fox. *Computational Models for Polydisperse Particulate and Multiphase Systems*. Cambridge University Press, Cambridge, UK, 2013.

[21] J. McDonald and M. Torrilhon. Affordable robust moment closures for CFD based on the maximum-entropy hierarchy. *J. Comput. Phys.*, 251:500–523, 2013.

[22] James G. McDonald and Clinton P. T. Groth. Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution. *Continuum Mechanics and Thermodynamics*, 25(5):573–603, 2012.

[23] R. McGraw. Description of aerosol dynamics by the quadrature method of moments. *Aerosol Science and Technology*, 27:255–265, 1997.

[24] I. Müller and T. Ruggeri. *Rational Extended Thermodynamics*. Springer-Verlag, New York, 1998.

[25] Konrad Schmüden. *The Moment Problem*, volume 277 of *Graduate Texts in Mathematics*. Springer, Cham, 2017.

[26] G. Strang. On the construction and comparison of difference schemes. *SIAM Journal of Numerical Analysis*, 5(3):506–517, 1968.

[27] E. F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. Springer, New York, 1999.

[28] A. Vié, F. Doisneau, and M. Massot. On the anisotropic Gaussian velocity closure for inertial-particle laden flows. *Commun. Comput. Phys.*, 17(01):1–46, January 2015.

[29] J. C. Wheeler. Modified moments and Gaussian quadratures. *Rocky Mt. J. Math.*, 4:287–296, 1974.

[30] C. Yuan and R. O. Fox. Conditional quadrature method of moments for kinetic equations. *J. Comput. Phys.*, 230(22):8216–8246, 2011.