Eulerian discrete kinetic framework in comoving reference frame for hypersonic flows

Abstract Flow physics vary in different regimes across the full Mach number range, with our knowledge being particularly poor about the hypersonic regime. An Eulerian realization of the particles on demand method, a kinetic model formulated in the comoving reference frame, is proposed to simulate hypersonic compressible flows. The present model allows for flux evaluation in different reference frames, in this case rescaled and shifted by local macroscopic quantities, i.e. fluid speed and temperature. The resulting system of coupled hyperbolic equations is discretized in physical space with a finite volume scheme ensuring exact conservation properties. Regularization via Grad expansion is introduced to implement distribution function and flux transformation between different reference frames. It is shown that the proposed method possesses Galilean invariance at a Mach number up to $100$. Different benchmarks including both inviscid and viscous flows are reproduced with the Mach number up to $198$ and pressure ratio up to $10^5$. Finally, the new model is demonstrated to be capable of simulating hypersonic reactive flows, including one-dimensional and two-dimensional detonations. The developed methodology opens up possibilities for the simulation of the full range of compressible flows, without or with chemical reactions, from the subsonic to hypersonic regimes, leading to enhanced understanding of flow behaviours across the full Mach number range.


Introduction
A variety of applications in science and engineering fall within the compressible flow regime.Generally speaking, compressible flows can be categorized as pertaining to different regimes based on Mach numbers, i.e. subsonic, transonic, supersonic and hypersonic flows.Considerable advances have been achieved in numerical methods, such as shock capturing and non-oscillatory schemes (Harten 1983;Harten et al. 1987;Liu, Osher & Chan 1994), to solve the compressible Navier-Stokes-Fourier (NSF) equations along with the advent of new and more powerful computing technologies.A remaining challenge is to develop a method capable of simulating the full range of flows from subsonic to hypersonic regimes.In particular, hypersonic flows with Mach numbers greater than five can, in some cases, be quite challenging to model using classical approaches relying on the NSF equations (Tsien 1946;Bird 1970).As the Knudsen number increases, the local state of the gas gets further away from the local thermodynamic equilibrium and the notion of the gas as a continuum-equilibrium fluid becomes questionable, in turn limiting the range of applicability of the NSF.For instance shock profiles for Mach numbers above two as obtained from the NSF are known not to match experimental observations and direct simulation Monte Carlo results, see for instance Greenshields & Reese (2007).The shortcomings of the NSF in these regimes have prompted the development of a variety of modified balance equations.One preferred, and rather successful approach to building balance laws beyond NSF, has been to derive reduced models based on the kinetic theory of gases, i.e. the Boltzmann equation and its variant.Over the years this has led to a variety of models, such as the more successful Grads moments method (Grad 1949) and its many variants, see Struchtrup (2005) for an overview.Another class of methods based on the kinetic theory of gases directly solves a discrete version of the Boltzmann equation in particle velocity space with quadratures guaranteeing correct recovery of moments of interest for the target balance equations.The now very popular lattice Boltzmann method (LBM) falls in that category.
The LBM, introduced in the early 1990s, was initially developed with the incompressible limit of the Navier-Stokes equations in mind (McNamara & Zanetti 1988;Benzi, Succi & Vergassola 1992).The past two decades have witnessed considerable efforts to extend the LBM to the compressible flow regime.These efforts can be categorized into two main approaches: (a) higher-order lattices and (b) two-population solvers with standard lattices.In the former Frapolli, Chikatamarla & Karlin (2015), the use of a larger number of discrete velocities and higher-order quadratures allows for the recovery of higher-order moments of the distribution function, in principle, correctly recovering the energy balance equation.While these approaches were successfully used to model compressible flows the computation and communication overhead brought about by the larger stencils is rather large especially considering the limited gain in maximum achievable Mach numbers.Furthermore, it has been observed that larger stencils result in more restricted stability domains in terms of temperature.It is worth noting that a number of recent studies have proposed solutions to improve numerical properties by using higher-order Hermite roots-based lattices.This has led to the development of semi-Lagrangian high-order lattice schemes which have been used for low supersonic flows (Wilde et al. 2020).As a way to reduce computational load for extension of the models to three-dimensional (3-D) models, the authors have also proposed strategies to reduce the number of discrete velocities (Wilde et al. 2021).The second category of approaches relies on the classical first-neighbour stencil supplemented with a second set of distribution functions carrying the total energy (Saadat et al. 2021).This class of models has been used to simulate a variety of flows in the transonic and supersonic regimes; however, it remains restricted to low supersonic speeds.
In parallel with the double distribution function approach, a number of studies have also proposed hybrid models replacing the second distribution function with discrete, finite volumes or finite differences solvers for the energy balance equation, see for instance (Feng, Sagaut & Tao 2015;Renard et al. 2021).While these approaches have been successful in modelling transonic flows and in some cases low supersonic flows, extension to higher Mach numbers has been quite challenging.A point common to all approaches is that their operation range is limited by the reference frame (Frapolli, Chikatamarla & Karlin 2016b), which is the static frame with a reference temperature dictated by the Gauss-Hermite quadrature.Other classes of solvers such as the discrete Boltzmann method (DBM), relying on the same principle as LBM for the discretization of particles speed space suffer from the same limitations, i.e. restricted operation range (see Xu et al. 2012;Gan et al. 2018b;Ji, Lin & Luo 2021).
In recent years, discrete velocity methods formulated on an adaptive reference frame, i.e. the particles on demand (PonD) method and its variants, have shown promising results (Frapolli et al. 2016b;Dorschner, Bösch & Karlin 2018).By adapting the reference frame of the discrete solver to the local fluid speed and temperature, higher-order moments are kept under control and the operation range of the solver is extended to very large Mach numbers.One consequence of the change in reference frame is that the streaming operator, for a Lagrangian solver, is no longer on the lattice.To that end, most Lagrangian realizations rely on an additional interpolation step to reconstruct distribution functions on the grid.While successfully applied to a number of flow configurations, they have been observed to suffer from conservativity issues due to the interpolation step, especially in flows involving large gradients (Kallikounis, Dorschner & Karlin 2022;Sawant, Dorschner & Karlin 2022).A flux-conserving realization of the PonD following the discrete unified gas kinetic scheme (Guo, Xu & Wang 2013) was recently introduced by Kallikounis et al. (2022).Restoration of conservative properties along with a regularized frame-to-frame transformation was shown to noticeably improve results and allow for large Mach number simulations.
Here, building upon the idea of the PonD method, different from the discrete unified gas kinetic scheme realization we propose a fully Eulerian finite-volume PonD solver for high-Mach-number flows.The outline of the paper is as follows.Section 2 presents a full procedure of model construction.The continuous kinetic equations are discretized in the phase space, which give adjustable Prandtl number and specific heat ratio.In addition, temporal and spatial discretization in an adaptive reference frame are presented and an overall algorithm is shown to clarify how to implement the new model.In § 3, a variety of simulations from simple low Mach one-dimensional (1-D) cases to more challenging 1-D and two-dimensional (2-D) configurations involving large variations in macro quantities are carried out to validate the new model.Furthermore, § 4 extends the proposed model to describe reactive flows where chemical reactions are coupled with fluid flows.Detonations are simulated with the developed model.In addition, a summary is provided in § 5.

Model description
In this part we will focus on model construction, starting from the Boltzmann equation with the Bhatnagar-Gross-Krook approximation for the collision operator (Bhatnagar, Gross & Krook 1954) where where D is the dimension of the physical space and R designates the gas constant.For the sake of convenience, R = 1 is used in the rest of this paper.
To take into account the effect of internal degrees of freedom allowing for variable specific heat ratios γ , following Rykov's original work (Rykov 1976) we supplement (2.1) with a balance equation for a second distribution function g carrying the excess internal energy stemming from non-translational degrees of freedom: with Ω g = (g eq (t, r, v) − g(t, r, v))/τ .Following the definition of the second distribution function, the local equilibrium state g eq can be defined as where C v indicates the specific heat at constant volume.As a result the total energy can be obtained as (2.5) While allowing for variable specific heat ratios, given that all moments relax at the same rate, the model still is limited to unity Prandtl number, defined as where C p is specific heat under constant pressure, μ the dynamic viscosity and κ the thermal conductivity.
To overcome this restriction, the relaxation to the equilibrium is split into two steps via the introduction of an intermediary state called the quasiequilibrium state (Gorban & Karlin 1994;Ansumali et al. 2007;Frapolli, Chikatamarla & Karlin 2014, 2016a;Hosseini & Karlin 2023).The collision terms are rewritten as where τ 1 and τ 2 are the two relaxation parameters related to the viscosity and thermal conductivity, and f * (t, r, v) and g * (t, r, v) are quasiequilibria.
In fact, the quasiequilibrium state is defined as the minimizer of the Boltzmann entropy function subject to a set of locally conserved fields and quasiconserved slow fields.Following Ansumali et al. (2007), in order to recover the NSF equations with variable Prandtl number, the quasiequilibrium distribution functions are required to satisfy the conservation of mass, momentum and total energy, (2.9) |v| 2 f eq + g eq dv. (2.10) In addition, the quasiequilibrium distribution functions are designed to conserve the third-order centred kinetic moments defined as where Qαβγ is the centred heat flux tensor and qα the centred heat flux.The constraints for the two quasiconserved quantities are Here, Q * αβγ and q * α are constructed by substituting f with f * in (2.11) and (2.12), respectively.These two constraints are the keys to realize an adjustable Prandtl number.In this way, Prandtl number Pr can be controlled by adjusting τ 1 and τ 2 independently as (2.15) Using the new collision terms (2.7) and (2.8), the kinetic equations (2.1) and (2.3) are partial differential equations (PDEs) in time, physical space and phase space -the space of the particles' velocities.In the upcoming sections we will introduce methods to discretize these equations while retaining conservativity and capturing the hydrodynamic regime, Euler and NSF, even under extreme conditions, i.e. large Mach number and temperature variations.

Frame-invariance of the Boltzmann equation
Before moving on into the details of the proposed discrete solver a brief discussion on the concept of frame is in order.A frame λ(u λ , T λ ), as intended here, is parameterized by a reference velocity, u λ , and a reference temperature, T λ .Particle velocities can therefore be written under this transformation of phase-space as It is straightforward to show, simply through the invariance of the corresponding moments system, the equivalence of Boltzmann equations on different frames, i.e. (2.17) where we have used the short-hand notation D λ t = ∂ t + v λ • ∇.Note that for the equivalence to stand, in each system the reference frame is fixed.

Frame-invariance of discrete velocity Boltzmann equations
Discretization of the particles' velocity space is usually achieved through finite-order quadratures, such as the Gauss-Hermite quadrature (Shan & He 1998), guaranteeing that moments of the continuous distribution function are matched by the discrete distribution function up to the order of the quadrature.Defining moments of the discrete distribution functions as we can introduce the set Φ of functions φ such that to define the set of moments supported by the lattice quadrature.Using (2.19) and the discussion in the previous section on frame-invariance of moments of the continuous distribution functions we can readily establish the frame-invariance of the set of moments defined by Φ of the discrete distribution functions.Moving one step forward one can also readily demonstrate the frame-invariance of the moments balance equations: (2.20) where φ λ defines the set of moments such that ( √ T λ c i + u λ )Φ λ ∈ Φ.For higher-order moments the equivalence is lost due to the limited order of the quadrature.Using these moments equations and further a multiscale analysis, see Kallikounis & Karlin (2024) for details, one can establish the order of quadrature needed to correctly recover the target balance equations.

Particles on demand
To eliminate errors in higher-order moments of the distribution function not supported by the lattice quadrature, in the PonD method (Dorschner et al. 2018), the reference frame is chosen via the local macroscopic quantities, i.e. fluid speed and temperature (2.21a,b) In doing so equilibrium distribution functions are accurately evaluated and only depend on the local density (2.22) where W i are weights of the Gauss-Hermite quadrature (Frapolli et al. 2015).Based on this, the discrete quasiequilibrium distribution functions under constraints (2.9)-(2.10)and (2.13)-(2.14)have the following expression for Pr < 1: (2.24) where I is the unit tensor, ':' is the Frobenius inner tensorial product and ⊗ denotes tensor product.The choice of the reference frame based on local velocity and temperature would, at first sight, lead to a space-and time-dependent reference frame.
Different from approaches such as Kauf (2011) and Köllermeier (2013) where an adaptive reference-based Boltzmann equation is solved globally in the domain bringing in non-commutativity of the discrete velocities with the space-and time-derivatives and resulting in additional terms in the Boltzmann equation involving derivatives of the reference frame velocity and temperature, in PonD at every point in space and time, r and t, the Boltzmann equation is solved in the fixed reference frame of that point, here denoted as λ for the sake of readability, where M λ λ(r,t) is the reference frame transformation operator discussed in the next section.

Distribution function transformation: Grad's projection
Given the invariance of moments supported by the lattice quadrature, a straightforward approach to the transformation of distribution functions between different reference frames λ = (u λ , T λ ) and λ = (u λ , T λ ) would be to match moments, (2.26) Another approach, proposed in Zipunova et al. (2021) and further extended in Kallikounis et al. (2022) is the regularized reconstruction of distribution functions.In the regularized approach the distribution function in the new reference frame λ is sought in the form of a finite-order Grad expansion, where a λ (n) and H (n)  i (c i ) are tensors of rank n denoting the Hermite coefficient and polynomial of the corresponding order in the λ reference frame, N is the order of the Grad expansion and T L the lattice temperature.Solving (2.26) one can derive the Grad coefficients in the reference frame λ as a function of moment in λ; for instance, for the first three coefficients, (2.28) Table 1.Lattice temperature T L , roots of Hermite polynomials c iα and weights W iα for D2Q16. and Equations (2.26) to (2.31) allow us to go transform distribution function from reference frame λ to λ : where M λ λ is a short notation for transformation operation.To capture fundamental flow physics properties at the NSF level, up to third-order expansion is necessary for f i , while for the internal distribution function g i , the second-order expansion is sufficient.Furthermore, we use a fourth-order quadrature for the discrete velocities in this study, i.e. the D2Q16 model in 2-D (Ansumali, Karlin & Öttinger 2003;Kallikounis & Karlin 2024).The necessary characteristics of the D2Q16 lattice are listed in table 1 and the Hermite coefficients are provided in Appendix B.

Finite volume formulation
Different from the commonly used 'propagation + collision' scheme in most LB solvers, the present work we make use of a fully Eulerian finite volume discretization, therefore targeting the integral form of the local discrete Boltzmann PDE system of (2.25) in the reference frame λ: where n is the outwards unit vector normal to the surface S and δV is an infinitesimal control volume defining the grid and δS the corresponding surface.For a given unit volume bound by discrete surfaces σ one obtains the space-discretized system of equations as (2.34) where {F λ i , G λ i }(σ ) are the fluxes at interfaces σ surrounding the control volume and making up δS.The cell-averaged distribution functions appearing here are defined as As a result, the local macroscopic properties obtained as moments of the distribution function, are also cell-averaged values, i.e. (2.38)

Time discretization: first-order Euler
While the time derivative term can be discretized using any of the available higher-order approaches, e.g.Runge-Kutta schemes, etc., here for the sake of simplicity, time stepping is carried out using a fully explicit first-order Euler scheme, resulting in the following discrete equations: (2.40)

Collision operator
In the context of the present first order in-time realization the collision term is evaluated explicitly and on the cell-averaged reference frame of the control volume, i.e. λ: Y. Ji, S.A. Hosseini, B. Dorschner, K.H. Luo and I.V. Karlin (2.42) Note that though not the main focus of the present contribution, the Crank-Nicolson second-order approximation akin to the redefined fluxes in the LBM can also be implemented with the present FVM realization.

Streaming operator: flux reconstruction
Next comes the issue of reconstruction of fluxes at interfaces, and with it two challenges: one that is common to all PDEs i.e. reconstruction of fluxes at interfaces from volume-averaged fields and one specific to the present model, the choice of the reference frame λ in which fluxes are to be reconstructed.
To have exact conservation at cell interfaces, out-flux from a given cell must match exactly the in-flux from its neighbour, which is only guaranteed by setting λ to that of their interface.For the first issue, as for any other PDE, pointwise values at cell interfaces are computed on the basis of a pointwise field reconstructed through polynomial interpolation from the corresponding cell-averaged field.This is illustrated for the case of a 1-D system in figure 1. There, F(t, x − δx/2) and F(t, x + δx/2) are the left-and right-hand interface fluxes at x − δx/2 and x + δx/2, respectively.Assuming a linear shape function and using a centred scheme, (2.45) For the first-order upwind scheme, here the left-hand interface flux F λ(t,x−δx/2) i (t, x − δx/2) uses the similar recipe as Once the fluxes are computed at the interface reference frame, they can be put back into (2.40) using (2.49) In order to obtain a good compromise between numerical dissipation and accuracy, a monotonic upstream-centred scheme for conservation laws (MUSCL) is adopted in the following simulations with contact discontinuity.The flux and the interface reference frame are constructed by a combination of first-order upwind scheme (low order) and second-order central scheme (high order) with a limiter function φ, and the minmod limiter is used as the limiter function (Roe 1986) (2.52) where x represents the ratio of successive gradients.

Flux conservation properties
As discussed in the introduction, one motivation for the development of the present realization is to overcome conservation issues encountered with semi-Lagrangian schemes.To illustrate conservation properties of the present model, we consider a cell centred at x in a 1-D system and the interface with its left-neighbouring cell centred at x − δx, located at x − δx/2.For the system to be conservative, the change in the cell-averaged value of a given moment must match the change in cell-averaged value of its left neighbour, (2.55) meaning invariance with respect to the reference frame transformation operator is the necessary and sufficient condition for exact conservation.Following the discussion on invariance of the discrete Boltzmann equation in § 2.1.2this means that as long as the flux of a given moment is supported by the quadrature of the lattice, conservation is ensured.To illustrate, consider the case of the invariants of the collision operator, i.e. mass, momentum and energy.To ensure conservation, the lattice must support moments

Overall algorithm
We summarize this part with a flow chart in figure 2. The evolution scheme of the present model with the finite volume formulation is implemented as (i) compute the interfacial reference frame λ by (2.44); (ii) transform the distribution functions in the stencil from the corresponding cell-averaged frames to the interfacial frame λ; (iii) calculate the interfacial distribution function f λ i (t); (iv) calculate the interfacial flux F λ i (t) by (2.43); (v) obtain the flux F λ i (t) evaluated on the reference frame λ from F λ i (t) by (2.49); (vi) update the distribution functions through (2.40) and compute the comoving reference frame λ at each cell.

Results and discussion
In this section, the model is validated through several 1-D and 2-D benchmarks, probing model properties such as conservation, Galilean invariance of dissipation rates, showcasing its suitability for extreme conditions involving large Mach numbers and temperature variations.For the remainder of this section, unless otherwise indicated, the specific heat ratio is γ = 1.4,and the Prandtl number is Pr = 1 due to the reason that these benchmarks are Euler-level problems where Pr is irrelevant.The Courant-Friedrichs-Lewy (CFL) number is set as CFL = max|u|(δt/δr) = 0.2.

Numerical properties: conservation of mass
To illustrate the necessity of adopting a flux-conserving space-discretization scheme for the present discrete velocity Boltzmann solver we first consider the simple 1-D case of the Sod shock tube with initial conditions (ρ, p, u x ) = (0.5, 4, 0), 0 ≤ x ≤ 0.5, (0.5, 0.5, 0), 0.5 < x ≤ 1.0. (3.1) The configuration is studied both with the finite-volume realization introduced here and the original interpolation-based semi-Lagrangian PonD model Dorschner et al. (2018).
Given distribution functions: Given reference frames: Compute interfacial reference frame: Compute interfacial distribution function: Calculate macro quantities and renew reference frame: ρ, u, T, λ(t, x) Calculate interfacial flux: The resolution of the computational domain in both simulations is L x = 500δx.Figure 3 shows the evolution of total mass in the system over time for the two schemes.The total mass for the semi-Lagrangian solver is not conserved, leading in turn to incorrect density levels as shown in figure 4. Instead, the finite volume scheme makes the total mass conserved and restores correct density profiles, in agreement with the reference solution.

Assessment of physical properties
Next, to demonstrate that the model correctly recovers Euler and NSF-level dynamics we look into the dispersion and dissipation rates of hydrodynamic eigenmodes, i.e. shear, normal and entropic.

Dispersion of normal mode: speed of sound
As a first step and to validate the dispersion properties of normal modes, we look at the temperature-dependence of the speed of sound.To that end a freely travelling pressure front is simulated here.A quasi-1-D domain L x × 1 (with L x = 800) is separated into two parts with a pressure difference of p = 10 −4 and a uniform initial temperature T 0 and velocity u 0 = 0.The speed of sound is computed by tracking the shock front over time and compared with the theoretical value c s = √ γ T. The simulations are performed with two different specific heat ratios, i.e. γ = 1.4 and 1.8, in a wide range of temperatures.
From figure 5 it is obvious that the present model can correctly capture the sound speed for various temperatures.

Dissipation: shear, normal and entropic modes
Next, we look into the dissipation rates of the three hydrodynamic modes, i.e. shear, normal and entropic.From the Chapman-Enskog analysis, the kinematic shear viscosity ν and bulk viscosity ν B in the present model are related to the relaxation coefficient τ 1 and τ 2 as and the thermal diffusivity is defined as First, a plane shear wave is simulated to measure the kinematic shear viscosity.The wave corresponds to a small sinusoidal perturbation imposed to the initial velocity field.The initial conditions of the flows are where the initial density and temperature are (ρ 0 , T 0 ) = (1.0,1.0), the perturbation amplitude A = 0.001 and u 0 is derived from the Mach number as Ma = u 0 / √ γ T 0 .To simulate the cases with high Mach numbers, a fine resolution is required and the CFL number is set to 0.01.The simulation domain is L x × 1 (with L x = 1600).The shear viscosity is measured by monitoring the time evolution of the maximum velocity and fitting an exponential function to it, i.e.
The obtained results are shown in figure 6.For different advection Mach numbers up to 100, the measured viscosity is in excellent agreement with the imposed values.The bulk viscosity is investigated from the decay rate of sound waves.For that, a small perturbation is added to the density field (Dellar 2001;Hosseini, Darabiha & Thévenin 2020), which places the study in the linear regime, excluding effects such as nonlinear steepening.For a discussion on nonlinear acoustics, readers are referred to studies such as Buick et al. (2000).In the linear regime, via linearization of the Navier-Stokes and continuity equations, it is readily shown that wave dynamics are governed by the linear lossy wave equation, see Lamb (1924).The flow is initialized as where ρ 0 = 1.0 and the perturbation is ρ = ρ − ρ 0 with initial amplitude A = 0.001.The initial temperature is T 0 = 1.0.The decay of energy E(t) = u 2 x + u 2 y − u 2 0 + c 2 s ρ 2 over time is supposed to fit an exponential function depending upon an effective viscosity ν e which is a combination of kinematic shear and bulk viscosity: defined by Dellar (2001) as (3.9) The resolution is identical with the shear mode.The measured effective viscosity is extracted by a simple least square fit.In figure 7, the measured effective viscosities are compared with the intended values for varying Mach number.Clearly, the present model is able to correctly capture the bulk viscosity.
To assess the thermal diffusivity α, a different type of perturbation is introduced into the initial state of the system, where the initial density and temperature are (ρ 0 , T 0 ) = (1.0,1.0), the perturbation amplitude A = 0.001 and u 0 is derived from Mach number.The thermal diffusivity is

Dissipation: spectral behaviour
In this part, we go beyond dissipation rate for well-resolved features and look into the spectral dissipation of the solver.For the three modes introduced in the previous section, we run the dissipation rate tests for different perturbation wavelengths.The imposed values of dissipation are chosen as ν = 10 −2 , ν e = 10 −2 and α = 10 −2 for the shear, normal and entropic modes, respectively, with Ma = 1 and Pr = 1 .In figure 9, we plot the ratio θ of the measured value to the imposed value defined as for varying wavenumbers k x δ x .Here ψ is the tested physical parameter.The model recovers the correct dissipation rates in the limit of vanishing wavenumbers while introducing hyperviscosity for larger wavenumbers.The considerable positive hyperviscosity introduced by the numerical discretization at larger wavenumbers can be attributed to the activation of the MUSCL flux limiter switching to a first-order upwind reconstruction.

Thermal-viscous coupling: thermal Couette flow
Couette flow with a temperature gradient is simulated to probe the viscous heat dissipation and the Prandtl number in this part.The initial configuration is a viscous fluid flow between two infinite parallel flat plates, and the physical quantity of the flow is ρ = 1, T = 1 and u x = u y = 0.The bottom plate is fixed, and the top plate moves in the x-direction at a speed U.The temperatures for the bottom and top plates are fixed at T 0 and T 1 , respectively.The temperature distribution T along the y-direction satisfies the analytical solution: number is Ma = 0.1.The isothermal no-slip boundary conditions are adopted at both ends.With the variations of Prandtl numbers and Eckert numbers, the simulation results fit the analytical solutions very well in figure 10.When the velocity of the top plate is large enough, the compressibility of the fluid becomes appreciable and the velocity is no longer linearly distributed along the y-direction.Under the conditions of μ/μ 1 = T/T 1 and of adiabatic bottom plate condition (Xu 2001), there is an analytical solution for the horizontal velocity and temperature (Liepmann & Roshko 1957): where τ w and μ 1 are the shear stress and dynamic viscosity on the top plate, respectively.For this case with high Mach number Ma = 3.0 and Pr = 2/3, the isothermal no-slip boundary condition is adopted at the top boundary and the adiabatic boundary is used at the bottom boundary.The same number of 200 grids is used.In figure 11, a comparison of the density, velocity and temperature with reference from the literature (Liepmann & Roshko 1957) demonstrates a very good match.

Numerical simulations: 1-D cases
To further validate the proposed model and showcase its performances, a set of 1-D benchmarks are simulated.First, two classical shock tube problems are modelled to test the capability of the model to capture the shock and expansion waves.Then, the Shu-Osher problem with Ma = 3.0 is simulated to probe the ability of the present model to resolve discontinuities with fine structures.Furthermore, another two simulations with very strong discontinuities (i.e.strong shock problem and double rarefaction problem) are performed to demonstrate the performance of the shifted discrete velocities.For the above 1-D simulations, the left-hand boundary of the domain is set as an inflow and the right-hand flow is an outflow.

Sod shock tube
The Sod shock tube problem is a classical benchmark to verify the performance of the model for compressible flows (Sod 1978).The initial conditions are described by (3.16) The resolution of the computational domain is L x = 1000.Figure 12 indicates the results of the present model and the reference solutions at time t = 0.2.It can be observed that the present results are in excellent agreement with the reference solution (Sod 1978).

Lax shock tube
Similar to the Sod shock tube problem, the second test case is the Lax shock tube problem with a discontinuity in the velocity field (Lax 1954), with the following conditions: (ρ, p, u x ) = (0.445, 3.528, 0.698), 0 ≤ x ≤ 0.5, (0.5, 0.571, 0), 0.5 < x ≤ 1.0. (3.17) The same resolution is adopted with the Sod shock tube problem.Simulation results for the density, pressure and velocity at time t = 0.1 are plotted and compared with the reference solution in figure 13.The two sets of results match well with each other despite minor deviations at the shock front.

Shock-density wave interaction
We continue with the Shu-Osher problem which places sinusoidal density fluctuations upstream of a moving shock front.The Mach number is 3.0 and the flow is initialized as (3.857, 10.333, 2.629), 0 ≤ x ≤ 1.0, (1 + 0.2 sin(5(x − 5)), 1, 0), 1.0 < x ≤ 10.0. (3.18) Figure 14 displays the computed density, pressure and velocity profiles at t = 1.8 with a resolution of 5000 points.Clearly, the simulations are also in accordance with the reference solution.The characteristic structures, such as the shock front and fluctuations are correctly captured.

Strong shock tube
With respect to the former tests, a much stronger discontinuity is imposed in the flow to assess the robustness and accuracy of the present model with self-adjusted discrete velocities.We consider the strong shock problem with a ratio of 10 5 between the pressure of the left-and right-hand side, and the initial conditions are Exact solution of this problem includes a strong contact discontinuity and a left rarefaction wave, in which the shock wave has a shock Mach number Ma = 198.Figure 15 exhibits simulation results of the present model at t = 0.012 with a resolution of L x = 5000.The results compare well with the reference solution.

Double rarefaction problem
In addition, the severe double rarefaction problem is simulated here (Hu & Khoo 2004;Hu, Adams & Shu 2013).The initial conditions are The solution of this case consists of two symmetric rarefaction waves and the centre region is close to a vacuum.The simulation adopts 5000 points and the results at t = 0.1 are plotted in figure 16.A good agreement of the present model with the reference solution is observed.Successful simulations of the rigorous problem demonstrate that the present model is robust and accurate enough to investigate compressible flows.

The 2-D cases
In this subsection, we further validate the model with a variety of 2-D problems, including the 2-D Riemann problem, double Mach reflection, flow over step and the shock-vortex interaction.

The 2-D Riemann configurations
The 2-D Riemann configurations consist of abundant geometric waves (Lax & Liu 1998;Kurganov & Tadmor 2002).The initial configuration is constructed by dividing a square domain into four quadrants with different densities, pressures and velocities.In detail, up to 19 different admissible configurations were extensively studied in the literature.In this study, four challenging configurations are solved with the following initial conditions in table 2. The four quadrants are distributed according to figure 17.Simulations are conducted on 1000 × 1000 grid nodes and zero-gradient boundary conditions are utilized Table 2.The initial conditions (ρ, p, u x , u y ) for the 2-D Riemann problems.
at each side.For each configuration the simulated result is illustrated by density contour plots in figure 18(a) and the reference solution of Lax & Liu (1998) is demonstrated in figure 18(b).
The typical complex phenomenology of the 2-D Riemann problems including the planar elementary waves including rarefaction waves and shock waves, as well as slip lines presented by Zhang & Zheng (1990) and Schulz-Rinne, Collins & Glaz (1993) are well reproduced in these contours.Generally, the simulated results agree well with the reference solutions.Specifically, for configuration (1), the intersection of shocks at interfaces creates a triple shock structure, which is well captured by the present model without spurious numerical oscillations.From configuration (2), the two shocks formed at interfaces with first quadrant propagate to the upper right corner and a pair of vortices is formed inside the third quadrant.All the pattern favourably compares with previous studies of Lax & Liu (1998) and Kurganov & Tadmor (2002).The characteristic vortex turning clockwise and a shape of a four-bladed propeller in configuration (3) is also well recovered.The slip lines of configuration (4) bisect the flow into a left-and right-hand region and join in a vortex near the centre.For configurations (1) and ( 2) where the initial states are symmetric along the diagonal line, the macroscopic fields should be symmetric all the time.From the simulation results, the symmetry property of terminal states is well retained.In addition, the ripples observed in the previous studies for configurations (3) and ( 4) are also well recovered.More results for different initial configurations can be seen in supplementary figures available at https://doi.org/10.1017/jfm.2024.94.These simulations further validate the new model with adaptive discrete velocities in the field of compressible flows.

Double Mach reflection
Next, we consider the problem of double Mach reflection which has been extensively studied in the literature (Woodward & Colella 1984;Ben-Dor 2007;Shirsat, Nayak & Patil 2022).The configuration consists of a Mach 10 shock wave colliding with a reflecting wall at an angle of π/6 with the wave propagation direction.The computational domain is a

Flow over step
This case consists of a uniform Mach 3 flow imposed in a 2-D wind tunnel with a step.The presence of the step leads to the formation of a transient shock wave reflecting at the wall.The domain is a rectangle of size 3 × 1 with a step of height 0.2 located at x = 0.6.Initial flow conditions in the domain are (ρ, u x , u y , p) = (1.4,3, 0, 1). (3.22) At the left-and right-hand boundaries, Dirichlet conditions corresponding to the the initial gas state are imposed, while at walls the flow is subject to reflecting boundary conditions.The simulation is conducted with 300 × 100 cells.The evolution of the shock over time is illustrated in figure 21 via six equidistant snapshots between t = 0.8 and t = 4.

Shock-vortex interaction
Sound generation by the interaction between a compressible vortex and a shock wave are simulated here.It helps assess the robustness and accuracy of the present model for unsteady and supersonic flows.We follow Inoue & Hattori (1999) and separate the main field by the Mach number of the shock, Ma s , and the left-and right-hand initial states satisfy the Rankine-Hugoniot condition.The left-hand region is set as upstream with (ρ, T, u x , u y ) l = (1, 0.05, Ma v √ γ T, 0).An isentropic vortex with vortex Mach number Ma v is passed through and perturbs the shock, where (x v , y v ) = (6, 12) is the vortex centre and r v is the vortex radius.The reduced radius is defined as r The simulation domain size is L x × L y = 28 × 24 discretized by a 1400 × 1200 mesh and initially the shock is located at x s = 8.Physical parameters are set to Ma s = 1.2, Ma v = 0.25, Pr = 0.75 and the Reynolds number Re = ρ l c s,l r v /μ = 800.
As common practice, we use a normalized pressure p = (p − p r )/p r to assess the present results where p r is the pressure behind the shock wave.Figure 22 demonstrates the normalized pressure contours at t * = 6 where we define t * = tc s,l /r v as the non-dimensional time.As seen from figure 22, the deformation of shock and the compression and rarefaction region are observed, showing very little difference with the reference solution (Inoue & Hattori 1999;Saadat et al. 2021).
In addition, distributions of the normalized pressure p along a radial cut of the fixed angle 45 • for t * = 6, t * = 8 and t * = 10 are plotted in figure 23   direct numerical simulation results by Inoue & Hattori (1999).Good agreement can be observed between the present and reference results and the propagation radially from the vortex of the precursor and the second sound are also captured.

Reactive flows
In this section, the model is extended to reactive flows.It is well known that the computation of reactive flows is numerically challenging due to the complex interaction between fluid flow and chemical reactions.In particular, detonation, as a kind of violent reactive flow, involves a wave propagating with a supersonic speed sustained by chemical reaction, and is characterized by strong compressibility and non-equilibrium (Law 2010).
The intrinsic features of detonation lead to additional challenges for numerical simulation.
Compared with the kinetic methods, traditional continuum-based Euler or Navier-Stokes equations have limitations in describing the non-equilibrium characteristics (Meng et al. 2013).On the other hand, the DBM derived from the Boltzmann equation, is proven to be viable for 1-D, 2-D and 3-D detonation simulations up to Mach 5.422 (Yan et al. 2013;Lin & Luo 2019;Ji, Lin & Luo 2022).However, the DBM utilizes fixed discrete velocities whose values are empirically preset.Furthermore, the DBM adopts the matrix inversion method to calculate the discrete equilibrium distribution functions, which relies on a delicate compromise between sufficient isotropy of the discrete velocities and the invertibility of the transformation matrix, especially for the 3-D case (Gan et al. 2018a), thereby limiting the range of attainable fluid temperatures/speeds.We now demonstrate that the present model with a reaction term can overcome the limitations of the current DBM and enable supersonic reactive flow simulations.To that end, the discrete kinetic equations are coupled with a two-step reaction model to describe the dynamic process of detonation.The 1-D and 2-D detonations are simulated to demonstrate the robustness of the constructed model for high-Mach-number reactive flows.
4.1.Kinetic equations Detonation is a supersonic combustion wave sustained by the energy of a chemical reaction.The wave is initiated by shock compression accompanied by a large amount of heat release within a short time.Across the wave, the thermodynamic states increase sharply (Mader 1979;Lee 2008).To model the dynamic process of detonation, the reaction terms R f and R g are added to the right-hand sides of the kinetic equations (2.1) and (2.3).The reaction terms give the variation rate of the distribution functions due to the chemical reaction.Following the idea of the previous study (Ji et al. 2021(Ji et al. , 2022)), the reaction terms are derived from the variation of energy caused by reaction, which read where T * is the temperature after the chemical reaction and satisfies From (2.5), the change rate of the temperature due to the chemical reaction is obtained as where Q r indicates the chemical heat release per unit mass of reactants and λ r represents the mass fraction of products.To model the dynamics of detonation driven by chain-branching kinetics, a two-step reaction mechanism is considered here (Ng et al. 2005): where ξ r and λ r denote the reaction progress variable in the ignition and the heat release period, respectively, and k I and k R are the rate constants for the two processes.Generally, the activation energy required in the ignition stage is larger than the heat release stage ε I < ε R for typical hydrocarbon mixtures, because the strong chemical bonds of the reactants are broken in the ignition process (Ng et al. 2005).Hence, in the present study, typical values for the activation energy are adopted: 4.2.Detonation simulations In this part, 1-D detonation is simulated to validate the proposed model.The results are compared with the classical theory of Zel'dovich, von Neumann and Döring (ZND theory).Furthermore, to compare the present model with the previous DBM, 2-D detonations with higher Mach numbers are simulated, and the performances of the two models are evaluated.

The 1-D detonation
The classical ZND theory describes detonation as a steady 1-D solution with an inner structure that consists of an infinitesimally thin shock wave called the von Neumann spike and a zone of exothermic chemical reaction.At the von Neumann spike, the reactants are compressed to a high pressure and temperature to initiate reaction.A main reaction layer follows the shock, where the products are formed and chemical energy is released.While the ZND structure is rarely observed in practice, it can be used as a reference solution for numerical simulations of detonation.Therefore, a 1-D detonation is simulated and compared with the ZND solution in this part.
To demonstrate the performance of the proposed model for hypersonic reactive flows, a 1-D detonation simulation with parameters Q = 500, γ = 1.4,k I = 5 × 10 2 and k R = 10 3 is conducted.The initial configuration satisfies the Rankine-Hugoniot conditions: The Mach number for the detonation is 26.224.Under this condition, very fine temporal and spatial resolution is required.Therefore, the CFL number is set as 0.05 and L x = 5000 is employed in this case.Since the frame of reference is set on the detonation wave, the inflow and outflow boundary conditions are used in the simulation.and the computational domain is L x × L y = 1500 × 500.Initially, the planar wave is perturbed in the y direction with a sinusoidal perturbation where ϕ represents the physical quantity and L and R denote left-and right-hand values, respectively.The initial perturbation amplitude is chosen as A = 10.In the literature (Lin & Luo 2019), the DBM utilizes a discrete velocity model D2V16 where eight free parameters (v a , v b , v c , v d , η a , η b , η c , η d ) are included to ensure numerical stability.To conduct the 2-D detonation simulation, the parameters of the D2V16 model are set as (v a , v b , v c , v d , η a , η b , η c , η d ) = (4, 3.6, 2.2, 0.7, 0, 0, 0, 2.6).In fact, the second-order Runge-Kutta scheme for temporal discretization and the second-order non-oscillatory and non-free-parameter dissipation difference scheme for the space derivative are adopted in the literature (Lin & Luo 2019).For the sake of comparison, we perform simulations with the simplest first-order Euler scheme and the MUSCL scheme for the two models in the present study.In addition, the temporal and spatial resolutions t = 1 × 10 −6 and x = 2 × 10 −5 are adopted.Figure 25 depicts the contour of pressure at t = 0.1 from the D2V16 model by Lin & Luo (2019) and the present D2Q16 model.The cellular structures, including Mach stems and triple points, are observed in both results and are in good agreement with each other.In addition, the oscillation periods of the two results are both 8.6 × 10 −3 .Moreover, we perform a quantitative comparison by measuring the pressure along the horizontal centre axis.Figure 26 demonstrates the pressure distribution of two models.As is evident from the results, the shock location and the second pressure jump are captured very well, except for a slight difference in the amplitude.
Furthermore, we conduct a numerical experiment where the chemical heat release is increased from Q = 2 to Q = 40 with 2.164 ≤ Ma ≤ 7.539, to investigate the performance of the two models for higher-Mach-number detonations.Results show that (4, 3.6, 2.2, 0.7, 0, 0, 0, 2.6).On the other hand, the present model is capable of stably simulating 2-D detonations with Mach numbers in the full set-up range.Figure 27 shows the snapshots of the pressure at different times in the evolution of 2-D detonation with Ma = 7.539 by the proposed model.It can be found that the pressure distribution at t = 0.038 matches well with that at t = 0.0402 which indicates that an entire cycle of the detonation is 2.2 × 10 −3 .In the whole cycle, one cellular structure is observed where the triple points divide the shock front into several parts.Overall, the improvements achieved by the proposed model are remarkable.It represents a strong basis to further study more complex conditions.In fact, a higher Mach number could also be achieved, but finer resolution and more computational resources are required, which is beyond the scope of this paper.

Figure 2 .
Figure 2. Flow chart for the proposed model: red lines represent transformation between different reference frames, blue lines stand for interpolation for interfacial values and black lines indicate other operations.

Figure 3 .Figure 4 .
Figure 3.Comparison of temporal history of the total mass with the finite volume (solid line) and the semi-Lagrangian (symbols) schemes.

Figure 5 .
Figure 5. Measurement of sound speed with different specific heat ratios.Blue and red colour denote the results with γ = 1.4 and γ = 1.8, respectively.Symbols denote the results of the present model, and solid lines represent theoretical solutions.

Figure 6 .
Figure 6.Measurement of the kinetic viscosity.Red and blue colour denote the results with ν = 1 × 10 −2 and ν = 5 × 10 −3 , respectively.Symbols denote the results of the present model, and solid lines represent the imposed viscosity.

Figure 7 .Figure 8 .Figure 8
Figure 7. Measurement of the effective viscosity.Red and blue colour denote the results with ν e = 5 × 10 −3 and ν e = 1 × 10 −3 , respectively.Symbols denote the results of the present model, and solid lines represent the imposed effective viscosity.

Figure 9 .
Figure9.Spectral dissipation analysis of the shear, normal and entropic modes.

Figure 10 .
Figure 10.Temperature ratio (T − T 0 )/(T 1 − T 0 ) in the Couette flow in the low Mach number case.(a) The Eckert number is fixed to 40, the Prandtl number takes the values 2.5, 1.0 and 0.72.(b) The Prandtl number is fixed to 0.5, the Eckert number takes the values 40, 20 and 4.Here lines are the reference solution and symbols are present solution.

Figure 11
Figure 11.(a) Velocity and (b) density and temperature distributions in high-speed Couette flow case with Pr = 2/3 and Ma = 3.0.Here lines are the reference solution (Liepmann & Roshko 1957) and symbols are the present solution.

Figure 12 .
Figure 12.Sod shock tube simulation results at t = 0.2: (a) density; (b) pressure; (c) velocity.Here the dashed lines are the reference solution (Sod 1978) and solid lines are the present solution.

Figure 13 .
Figure 13.Lax shock tube simulation results at t = 0.1: (a) density; (b) pressure; (c) velocity.Here the dashed lines are the reference solution (Lax 1954) and the solid lines are the present solution.

Figure 14 .
Figure 14.Shu-Osher problem simulation results at t = 1.8:(a) density; (b) pressure; (c) velocity.Here the dashed lines are solutions obtained using the numerical solver HyPar (see Debojyoti, John & Youngdae (2013) for more details on the code) and the solid lines are the present solution.

Figure 15 .
Figure 15.Strong shock simulation results at t = 0.012: (a) density; (b) pressure; (c) velocity.Here the dashed lines are the reference solution from an exact Riemann solver and the solid lines are from the present solution.

Figure 16 .Figure 17 .
Figure 16.Double rarefaction problem simulation results at t = 0.1: (a) density; (b) pressure; (c) velocity.Here the dashed lines are the reference solution(Hu et al. 2013) and the solid lines are the present solution.

Figure 20 .
Figure 20.Pressure profile for the double Mach reflection problem at y = 0.2, t = 0.25.Here the dashed line is the reference solution (Ben-Dor 2007; Shirsat et al. 2022) and the solid line is the present solution.
Figure 24 displays the comparison of the distribution of physical quantities around the front shock between the present results and the ZND solution.It can be clearly observed that the von Neumann spike is accurately captured by the developed model, and a good quantitative agreement is obtained between the present model and the ZND solution.This encouraging result shows the present model significantly extends the applicability of the previous DBM to higher Mach numbers.4.2.2.The 2-D detonationNow we proceed to a 2-D detonation simulation with Mach number 2.126 in the literatureLin & Luo (2019) to evaluate the robustness of the present model in two dimensions.The same reaction parameters as those in the literatureLin & Luo (2019) are used here: Q = 2, γ = 1.4,k I = 5 × 10 2 and k R = 10 4 .The initial configuration is (ρ, p, u x ) =(1.480, 3.054, −1.700), 0 ≤ x ≤ 0.027,

Figure 24 .
Figure 24.The 1-D detonation simulation results at t = 0.1: (a) density; (b) pressure; (c) velocity.Here the dashed lines are the reference solution (Law 2010) and the solid lines are the present solution.
Here v designates the particle velocity and r the position in space.Here, a single parameter, τ , controls the relaxation rate of the distribution function towards equilibrium state, known as local Maxwellian distribution parameterized by values of local temperature T, fluid velocity u and density ρ: and f eq (t, r, v) represent, respectively, the distribution function and the local equilibrium distribution function.