Bhatnagar-Gross-Krook Research Articles

Evaluation of the non-Newtonian lattice Boltzmann model coupled with off-grid bounce-back scheme: Wall shear stress distributions in Ostwald-de Waele fluids flow.

We present a comprehensive analysis of the non-Newtonian lattice Boltzmann method (LBM) when it is used to simulate the distribution of wall shear stress (WSS). We systematically identify sources of numerical errors associated with non-Newtonian rheological behavior of fluids in off-grid geometries. We implement the single relaxation time, Bhatnagar-Gross-Krook (BGK), and multiple relaxation time (MRT) collision operators and investigate flow in a two-dimensional channel aligned with lattice directions and off-grid Hagen-Poiseuille flow of Ostwald-de Waele (power-law) fluids. As for boundary conditions, we implement constant body force-driven and pressure-driven flows. These two boundary conditions have different numerical challenges, which include numerical stability, accuracy, mass conservation, and compressibility effects, which are inherent in the LBM method. Our results indicate that MRT, when the relaxation times are adequately tuned in the non-Newtonian case, significantly improves the WSS distribution accuracy and the numerical stability of the LBM. MRT also enhances the stability and accuracy for non-Newtonian fluids compared with the Newtonian case, meaning that it is questionable if a BGK collision operator is appropriate to use in a non-Newtonian case with off-grid boundaries. When analyzing the non-Newtonian LBM in the context of staircase walls and interpolated bounce-back (IBB) walls, a MRT collision operator with the appropriate choice of tunable relaxation times makes it possible to achieve numerically accurate results without a significant increase in grid resolution for matching to the analytical solution of WSS distributions. In analyzing the non-Newtonian flows, we show that the viscosity dependency of bounce-back walls in the BGK-LBM deviates from the results obtained under Newtonian assumptions. The power-law index further influences these discrepancies, and errors caused by the viscosity dependency of the bounce-back boundary conditions can be effectively mitigated by implementing the MRT procedure. Results show that non-Newtonian fluids, in contrast with the Newtonian assumption, encounter a greater mass imbalance when flowing through a periodic system with IBB walls. MRT can address this challenge, as it allows for independent adjustments of physical relaxation times and enhances mass conservation in the case of non-Newtonian fluids. In pressure-driven non-Newtonian flows, there is a significant impact of bulk viscosity. This aspect is often overlooked in Newtonian simulations but can significantly impact fluid adapting to rapid changes in local effective viscosity. One of our main conclusions is that the MRT collision operator with tuned relaxation times can effectively resolve numerical problems caused by non-Newtonian rheological properties and off-grid geometries. We also provide practical guidelines for selecting the most suitable simulation approach.

Energy loss of a heavy fermion in a collisional QED plasma

We compute the energy loss of heavy fermions moving in a plasma, taking into account the modification of the photon collective modes induced by collisions using a Bhatnagar-Gross–Krook collisional kernel. We include contributions from both hard and soft scatterings of the heavy fermion using a collisionally modified hard-thermal-loop resummed propagator. Using this method, one does not need to introduce a separation scale between hard- and soft-momentum exchanges. To place our calculation in context, we review other theoretical approaches to computing the collisional energy loss of fermions and discuss the systematics and results obtained in each approach compared to using a resummed propagator for both hard and soft momentum exchanges. Our final results indicate that self-consistently including the effect of collisions in the self-energies of the resummed propagator results in an increased energy loss compared to using collisionless hard-thermal-loop propagators. The effect becomes larger as the magnitude of the coupling constant and the velocity of the fermion increase. Published by the American Physical Society 2024

Vortex-dipole impingement with convex and concave boundaries

The vortex-dipole interactions with convex and concave boundaries in a two-dimensional domain are analyzed using the Bhatnagar–Gross–Krook (BGK) collision-based lattice Boltzmann method (LBM). The formation and detachment of the boundary layers at the flat wall creates new dipoles of non-identical vortices, following cyclic trajectories and exhibiting logarithmic variations in the production of maximum vorticity with Reynolds number (Re). The vortex-dipole interactions with concave boundaries produce a series of secondary dipoles, whose relative strength linearly decreases for the given Re in subsequent vortex-releasing events. Oblique interactions with cavity corners cause secondary dipoles to undergo head-on collisions at the domain center, influencing partner exchange and orthogonal propagation. Conversely, convex boundaries split the primary vortex-dipole upon impingement, intensifying the vorticity production and strain effects. Whereas, the larger ingestion of vorticity at the convex corners of the “T-shaped” cavity forms new dipoles, which deflect, rebound, and follow a parabolic trajectory for the vortex exchange. Normalized enstrophy, Ω(t), and palinstrophy, P(t), show distinct peaks during dipole interactions with boundary walls, influencing enhanced kinetic energy, E(t), decay. Moreover, the evolution of E(t) and Ω(t) satisfies the relation valid for no-slip boundaries. The collision behavior, vorticity production, and vortex rebound are functions of Re. The convex boundaries modify the scaling results of maximum Ω(t) and P(t) to higher exponent values than the dipole interactions with flat and concave boundaries.

Perfectly matched layers for the Boltzmann equation: Stability and sensitivity analysis

We study the stability and sensitivity of an absorbing layer for the Boltzmann equation by examining the Bhatnagar–Gross–Krook (BGK) approximation and using the perfectly matched layer (PML) technique. To ensure stability, we discard some parameters in the model and calculate the total sensitivity indices of the remaining parameters using the ANOVA expansion of multivariate functions. We conduct extensive numerical experiments on two test cases to study stability and compute the total sensitivity indices, which allow us to identify the essential parameters of the model.

Relativistic BGK hydrodynamics

Bhatnagar–Gross–Krook (BGK) collision kernel is employed in the Boltzmann equation to formulate relativistic dissipative hydrodynamics. In this formulation, we find that there remains freedom of choosing a matching condition that affects the scalar transport in the system. We also propose a new collision kernel which, unlike BGK collision kernel, is valid in the limit of zero chemical potential and derive relativistic first-order dissipative hydrodynamics using it. We study the effects of this new formulation on the coefficient of bulk viscosity.

Exact Results for Non-Newtonian Transport Properties in Sheared Granular Suspensions: Inelastic Maxwell Models and BGK-Type Kinetic Model.

The Boltzmann kinetic equation for dilute granular suspensions under simple (or uniform) shear flow (USF) is considered to determine the non-Newtonian transport properties of the system. In contrast to previous attempts based on a coarse-grained description, our suspension model accounts for the real collisions between grains and particles of the surrounding molecular gas. The latter is modeled as a bath (or thermostat) of elastic hard spheres at a given temperature. Two independent but complementary approaches are followed to reach exact expressions for the rheological properties. First, the Boltzmann equation for the so-called inelastic Maxwell models (IMM) is considered. The fact that the collision rate of IMM is independent of the relative velocity of the colliding spheres allows us to exactly compute the collisional moments of the Boltzmann operator without the knowledge of the distribution function. Thanks to this property, the transport properties of the sheared granular suspension can be exactly determined. As a second approach, a Bhatnagar-Gross-Krook (BGK)-type kinetic model adapted to granular suspensions is solved to compute the velocity moments and the velocity distribution function of the system. The theoretical results (which are given in terms of the coefficient of restitution, the reduced shear rate, the reduced background temperature, and the diameter and mass ratios) show, in general, a good agreement with the approximate analytical results derived for inelastic hard spheres (IHS) by means of Grad's moment method and with computer simulations performed in the Brownian limiting case (m/mg→∞, where mg and m are the masses of the particles of the molecular and granular gases, respectively). In addition, as expected, the IMM and BGK results show that the temperature and non-Newtonian viscosity exhibit an S shape in a plane of stress-strain rate (discontinuous shear thickening, DST). The DST effect becomes more pronounced as the mass ratio m/mg increases.

Evaluation of stochastic particle Bhatnagar–Gross–Krook methods with a focus on velocity distribution function

For precise application of Bhatnagar–Gross–Krook (BGK) methods, assessing its accuracy in non-equilibrium flows is necessary. Generally, this assessment relies on macroscopic parameters, which are moments of the velocity distribution function (VDF). However, in non-equilibrium flows, the significance of each moment diminishes as the VDF deviates from the Maxwellian VDF. This study investigates the VDF in non-equilibrium flows. Two Prandtl-corrected BGK methods, the ellipsoidal statistical BGK and Shakhov BGK (SBGK), are compared with the direct simulation Monte Carlo method. To observe the VDF while excluding the effects of convection, the homogeneous relaxation of the initial non-equilibrium state is analyzed. The VDF in Couette flow and normal shock waves, where collision and convection coexist, is then examined. When comparing the accuracy of the BGK methods using higher-order moments, inconsistencies are observed. However, when comparing the VDFs, the SBGK method reproduces the non-equilibrium VDF more accurately. The results demonstrate the importance of the VDF in the evaluation of non-equilibrium flows.

Effect of a Wake-Field on the Dissociation of Quarkonium in Collisional Quark–Gluon Plasma

We have studied wake effects on the dissociation of heavy quarkonia states J/ψ and Y by introducing an in-medium modification to the inter-quark potential. The wakes in the quark–gluon plasma were modeled using linear response theory using a dynamic dielectric function obtained from kinetic theory (Boltzmann equation) with a Bhatnagar–Gross–Krook (BGK) collision term. The in-medium modified potential was used to investigate the dissociation character depending on various parameters such as the velocity of quarkonium moving through the medium and the collision frequency. We have also calculated critical values of the dissociation temperature. Modifications of the dissociation energy due to wake-field effects were found.

Rheology of a dilute binary mixture of inertial suspension under simple shear flow

Abstract The rheology of a dilute binary mixture of inertial suspension under simple shear flow is analyzed in the context of the Boltzmann kinetic equation. The effect of the surrounding viscous gas on the solid particles is accounted for by means of a deterministic viscous drag force plus a stochastic Langevin-like term defined in terms of the environmental temperature Tenv. Grad’s moment method is employed to determine the temperature ratio and the pressure tensor in terms of the coefficients of restitution, concentration, the masses and diameters of the components of the mixture, and the environmental temperature. Analytical results are compared against event-driven Langevin simulations for mixtures of hard spheres with the same mass density m1/m2 = (σ(1)/σ(2))3, mi and σ(1) being the mass and diameter, respectively, of the species i. It is confirmed that the theoretical predictions agree with simulations of various size ratios σ(1)/σ(2) and for elastic and inelastic collisions in a wide range of parameter space. It is remarkable that the temperature ratio T1/T2 and the viscosity ratio η1/η2 (ηi being the partial contribution of the species i to the total shear viscosity η = η1 + η2) discontinuously change at a certain shear rate as the size ratio increases; this feature (which is expected to occur in the thermodynamic limit) cannot be completely captured by simulations due to the small system size. In addition, a Bhatnagar–Gross–Krook (BGK)-type kinetic model adapted to mixtures of inelastic hard spheres is exactly solved when Tenv is much smaller than the kinetic temperature T. A comparison between the velocity distribution functions obtained from Grad’s method, the BGK model, and simulations is carried out.

Combination of Physics-Informed Neural Networks and Single-Relaxation-Time Lattice Boltzmann Method for Solving Inverse Problems in Fluid Mechanics

Physics-Informed Neural Networks (PINNs) improve the efficiency of data utilization by combining physical principles with neural network algorithms and thus ensure that their predictions are consistent and stable with the physical laws. PINNs open up a new approach to address inverse problems in fluid mechanics. Based on the single-relaxation-time lattice Boltzmann method (SRT-LBM) with the Bhatnagar–Gross–Krook (BGK) collision operator, the PINN-SRT-LBM model is proposed in this paper for solving the inverse problem in fluid mechanics. The PINN-SRT-LBM model consists of three components. The first component involves a deep neural network that predicts equilibrium control equations in different discrete velocity directions within the SRT-LBM. The second component employs another deep neural network to predict non-equilibrium control equations, enabling the inference of the fluid’s non-equilibrium characteristics. The third component, a physics-informed function, translates the outputs of the first two networks into physical information. By minimizing the residuals of the physical partial differential equations (PDEs), the physics-informed function infers relevant macroscopic quantities of the flow. The model evolves two sub-models that are applicable to different dimensions, named the PINN-SRT-LBM-I and PINN-SRT-LBM-II models according to the construction of the physics-informed function. The innovation of this work is the introduction of SRT-LBM and discrete velocity models as physical drivers into a neural network through the interpretation function. Therefore, the PINN-SRT-LBM allows a given neural network to handle inverse problems of various dimensions and focus on problem-specific solving. Our experimental results confirm the accurate prediction by this model of flow information at different Reynolds numbers within the computational domain. Relying on the PINN-SRT-LBM models, inverse problems in fluid mechanics can be solved efficiently.

Stable LBM schemes for acoustic scaling simulations under high Reynolds to Mach ratio: Introduction to the DM-TS operator

Hydroacoustic analyses using direct methods are challenging in the field of computational physics. To date, hybrid methods are an alternative approach adopted for performing hydroacoustic analysis, although it demonstrates clear limitations in analyzing nonlinear turbulent noise. The lattice Boltzmann method (LBM) provides a solution to address the limitations of the hybrid method and is frequently applied in aeroacoustics. However, existing methods, such as the Bhatnagar–Gross–Krook (BGK) operator, suffer from instabilities that manifest in an oscillatory fashion in underwater conditions. Hydroacoustic situations are characterized by high Reynolds number (Re) to Mach number (Ma) ratios, and these conditions correspond to cases where stability problems occur. In this study, a direct-method-based two-step (DM-TS) LBM collision operator is introduced to solve the stability problem. The stability, accuracy, and consistency of the scheme are investigated and compared with the BGK-LBM. The DM-TS operator is developed to provide the same accuracy order as the BGK scheme while addressing the stability problem. Vortex shedding from a two-dimensional cylinder is simulated using the DM-TS operator to prove that the proposed method could be used to solve hydroacoustic problems. The results of the analysis indicate that dynamic and acoustic pressures can be obtained simultaneously using the DM–TS operator.

A phase-field-based multiple-relaxation-time lattice Boltzmann method for incompressible multiphase flows with density and viscosity contrasts

In this work, we present a unified framework of phase-field-based multiple-relaxation-time lattice Boltzmann (MRT-LB) method for incompressible multiphase flows with density and viscosity contrasts where a block-lower-triangular relaxation matrix is introduced. The present framework includes the classic MRT-LB model and central-moments-based LB model (CLBM, a particular version of the MRT-LB model). The governing equations of incompressible multiphase flows can be reproduced precisely through the direct Taylor expansion method at the second-order of expansion parameters. In this work, the pressure p is re-estimated, and we prove that the present pressure expression is suitable for classic MRT model, CLBM and various lattice structures with an accuracy of O(Δt2+ΔtMa2). In addition, the model construction and analysis are all carried out in the velocity space. As a result, the present model has a similar structure to the traditional Bhatnagar-Gross-Krook (BGK) model. Therefore the process of their analysis and implementation is also similar. We then test the present model by five classic physical problems. The results show that the present model is found to be accurate and has better numerical stability over the BGK model, and has the ability to simulate multiphase flow problems with large density ratio (1000).

Lattice Boltzmann solidification modeling of forced convection internal flows applied to Gen-IV nuclear reactor coolants

Generation-IV nuclear reactors involve novel coolants such as molten salts and lead. One of the issues that must be assessed with these coolants is the risk of solidification during reactor operation owing to their high melting temperature. However, the numerical simulation of solidification problems is computationally demanding. In this work, we assess the applicability of Lattice Boltzmann Method (LBM) to nuclear coolants solidification problems under laminar internal forced flow conditions. The double distribution Total Enthalpy-Partially Saturated Method is implemented in the open-source code OpenLB and compared against more standard Finite-Volume Computational Fluid Dynamics (FV-CFD) results, which uses the enthalpy-porosity Eulerian Multiphase method, in commercial code Star-CCM+. The test cases comprise the solidification of two fluids with significantly different Prandtl numbers, molten salt and lead, in a 3D pipe geometry in laminar flow conditions with an initial Reynolds number of 100. The solid interface advances from the cooled wall towards the bulk of the pipe reducing the area of flow passage as the transient evolves. Results of LBM velocity fields, solid fraction profiles, temperature fields, and pressure fields are in good agreement with FV-CFD results for both coolants. Furthermore, LBM stability and accuracy limits are independently explored for molten salts and lead using Two-Relaxation Time (TRT) and Bhatnagar-Gross-Krook (BGK) collision operators. Simulation time between FV-CFD and LBM is also contrasted for both coolants. In this work, LBM solvers use explicit time schemes with a regular grid, whereas FV-CFD solutions are able to include implicit time schemes and near-wall mesh refinement. For molten salt flows (high Prandtl number) FV-CFD allows for considerably larger timesteps than LBM, which is limited by the maximum lattice Mach number limitation and stability considerations. For lead flows (low Prandtl number), the implicit scheme in FV-CFD is limited in terms of accuracy for large timesteps due to the limiting Fourier-number. On the other hand, LBM does not present notable stability issues for lead and the lattice Mach number limit is not as relevant as for the molten salt case for this Reynolds number, allowing for faster simulation times than FV-CFD.

Development and Validation of an in-House Direct Simulation Monte Carlo Code for Supersonic Rarefied Gas Microflows

An in-house Direct Simulation Monte Carlo (DSMC) code has been developed using an object-oriented interface for the purpose of eventually modelling gas and plasma microflows in micromechanical systems, by coupling it with an existing Particle-In-Cell (PIC) solver. The code is validated with canonical flow problems such as the micro-Couette and micro-Poiseuille flows. Since the emphasis is on supersonic and highly compressible flows, we perform validation studies of our algorithm for high velocity and Mach number cases, as compared to known analytical or semi-analytical profiles via asymptotic Bhatnagar-Gross-Krook (BGK) and Chapman-Enskog theory. Reasonably good agreement is found for all the test cases, which indicates the reliability of the present algorithm for modelling supersonic rarefied gas flows. Our results have implications for future validation studies of DSMC codes, where concrete quantitative theoretical and experimental comparisons across different flow configurations are required.

Unified stochastic particle simulation of polyatomic gas flows using SPARTACUS

Numerically simulating multiscale nonequilibrium gas flows that involve internal energy exchange is a challenging problem in engineering applications due to the complexity of the physical processes. The unified stochastic particle (USP) method, based on the Bhatnagar-Gross-Krook (BGK) model, is a promising approach for efficiently simulating multiscale flows. In this paper, we develop a USP method for polyatomic gases that can simulate hypersonic gas flows involving internal energy exchange across all flow regimes. This method is integrated into the open-source SPARTACUS solver, which was previously developed for monatomic gases. The newly developed USP method reassigns the velocity and internal energies of a subset of computational particles based on a precisely designed local target distribution function. A conservation algorithm is applied to ensure that the momentum and energy are conserved. To validate the accuracy and efficiency of the solver, we perform three benchmark cases, including flow past a two-dimensional cylinder, a three-dimensional 70° blunted cone, and an axisymmetric double cone configuration. The results demonstrate that SPARTACUS achieves excellent agreement with the reference results of DSMC and experimental data. Additionally, SPARTACUS provides significant computational efficiency advantages over the DSMC solver SPARTA, particularly in simulating three-dimensional gas flows.

A novel linear stability analysis method for plane Couette flow considering rarefaction effects

Following the stability analysis method in classic fluid dynamics, a linear stability equation (LSE) suitable for rarefied flows is derived based on the Bhatnagar–Gross–Krook (BGK) equation. The global method and singular value decomposition method are used for modal and non-modal analysis, respectively. This approach is validated by results obtained from Navier–Stokes (NS) equations. The modal analysis shows that LSEs based on NS equations (NS-LSEs) begin to fail when the Knudsen number ($Kn$) increases past$\sim$0.01, regardless of whether a slip model is used. When$Kn\geq 0.01$, the growth rate of the least stable mode is generally underestimated by the NS-LSEs. Under a fixed wavenumber, the pattern (travelling or standing wave) of the least stable mode changes with$Kn$; when the mode presents the same pattern, the growth rate decreases almost linearly with increasing$Kn$; otherwise, rarefaction effects may not stabilize the flow. The characteristic lengths of the different modes are different, and the single-scale classic stability analysis method cannot predict multiple modes accurately, even when combined with a slip model and even for continuum flow. However, non-modal analysis shows that this error does not affect the transient growth because modes with small growth rates offer little contribution to the transient growth. In rarefied flow, as long as the Mach number ($Ma$) is large enough, transient growth will occur in some wavenumber ranges. The rarefaction effect plays a stabilizing role in transient growth. The NS-LSEs-based method always overestimates the maximum transient growth.

Non-equilibrium flow of van der Waals fluids in nano-channels

The Enskog–Vlasov equation provides a consistent description of the microscopic molecular interactions for real fluids based on the kinetic and mean-field theories. The fluid flows in nano-channels are investigated by the Bhatnagar–Gross–Krook (BGK) type Enskog–Vlasov model, which simplifies the complicated Enskog–Vlasov collision operator and enables large-scale engineering design simulations. The density distributions of real fluids are found to exhibit inhomogeneities across the nano-channel, particularly at large densities, as a direct consequence of the inhomogeneous force distributions caused by the real fluid effects including the fluid molecules' volume exclusion and the long-range molecular attraction. In contrast to the Navier–Stokes equation with the slip boundary condition, which fails to describe nano-scale flows due to the coexistence of confinement, non-equilibrium, and real fluid effects, the Enskog–Vlasov–BGK model is found to capture these effects accurately as confirmed by the corresponding molecular dynamics simulations for low and moderate fluid densities.

BGK Model for Two-Component Gases Near a Global Maxwellian

In this paper, we establish the existence of the unique global-in-time classical solutions to the two-component Bhatnagar–Gross–Krook (BGK) model suggested in [C. Klingenberg, M. Pirner, and G. Puppo, Kinet. Relat. Models, 10 (2017), pp. 445–465] when the initial data is a small perturbation of global equilibrium. For this, we carefully analyze the dissipative nature of the linearized two-component relaxation operator and observe that the partial dissipation from the intraspecies and the interspecies linearized relaxation operators are combined in a complementary manner to give rise to the desired dissipation estimate of the model. We also observe that the convergence rate of the distribution function increases as the momentum-energy interchange rate between the different components of the gas increases.

Seebeck and Nernst coefficients of a magnetized hot QCD medium with a number conserving kernel

We study the thermoelectric response of a hot and magnetized QCD medium created in the noncentral events at heavy-ion collider experiments. The collisional aspects of the medium have been embedded in the relativistic Boltzmann transport equation (RBTE) using the Bhatnagar-Gross-Krook (BGK) collision integral, which insures the particle number conservation, unlike the commonly used relaxation time approximation (RTA). We have incorporated the thermal medium effects in the guise of a quasiparticle model, where the interaction among the quarks and gluons is assimilated in the medium dependent masses of the quarks, which have been evaluated using imaginary-time formalism of thermal QCD with a background magnetic field. In the absence of $B$, the Seebeck coefficient for individual quark flavors gets slightly reduced in the BGK term in comparison to naive RTA, while it gets enhanced for the composite partonic medium. In the strong magnetic field ($B$), the BGK term enhances the Seebeck coefficient for the individual flavors, as well as that for the medium. The medium Seebeck coefficient rises with the strength of quark chemical potential ($\ensuremath{\mu}$) in the absence, as well as that in the strong $B$. We observe chirality dependence in the transport coefficients in the weak $B$ as the masses of chiral modes become nondegenerate. In the case of the $L$ modes, the BGK collision term causes a slight reduction in the Seebeck coefficient, while for $R$ modes both the collision integrals produce the same results. The Nernst coefficient gets reduced (enhanced) for $L$ ($R$) chiral modes in the BGK term.

Extension of the Shakhov Bhatnagar–Gross–Krook model for nonequilibrium gas flows

Bhatnagar–Gross–Krook (BGK) models are widely used to study rarefied gas dynamics. However, as simplified versions of the Boltzmann collision model, their performances are uncertain and need to be carefully investigated in highly nonequilibrium flows. In this study, several common BGK models, such as the ellipsoidal statistical BGK (ES-BGK) and Shakhov BGK (S-BGK) models, are theoretically analyzed using their moment equations. Then, numerical comparisons are performed between the Boltzmann collision model and BGK models based on various benchmarks, such as Fourier flow, Couette flow, and shock wave. The prediction performance of the ES-BGK model is better than that of the S-BGK model in Fourier flow, while prediction performance of the S-BGK model is better than that of the ES-BGK model in Couette flow and shock wave. However, with increasing Knudsen number or Mach number, the results of both ES-BGK and S-BGK deviate from the Boltzmann solutions. These phenomena are attributed to the incorrect governing equations of high-order moments of BGK models. To improve the performance of the current BGK models, the S-BGK model is extended by adding more high-order moments into the target distribution function of the original one. Our analytical and numerical results demonstrate that the extended S-BGK (S-BGK+) model provides the same relaxation coefficients as the Boltzmann collision model for the production terms of high-order moment equations. Compared with the other BGK models, the proposed S-BGK+ model exhibits better performance for various flow regimes.