Hostname: page-component-5db58dd55d-d6ndz Total loading time: 0 Render date: 2026-07-08T08:08:13.139Z Has data issue: false hasContentIssue false
  • English
  • Français

Revisit viscous shock tube at low Reynolds number

Published online by Cambridge University Press:  08 July 2026

Yue Zhang Open the ORCID record for Yue Zhang [Opens in a new window]  and
Kun Xu Open the ORCID record for Kun Xu [Opens in a new window]
Yue Zhang
Affiliation:
Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
Kun Xu*
Affiliation:
Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong Shenzhen Research Institute, Hong Kong University of Science and Technology, Shenzhen, PR China
*
Corresponding author: Kun Xu, makxu@ust.hk
Rights & Permissions [Opens in a new window]

Abstract

Content of image described in text.

The viscous shock tube is a canonical test case for assessing Navier–Stokes (NS) solvers in the continuum-flow regime, widely used to validate numerical accuracy and probe flow physics. It features a rich set of interacting structures – shock and rarefaction waves, contact discontinuities, boundary layers and their couplings – spanning multiple spatial and temporal scales. However, NS-based modelling, which presumes near-equilibrium behaviour, may fail to capture important non-equilibrium effects even in nominally continuum conditions. This study investigates the viscous shock tube at low Reynolds numbers and demonstrates the presence of non-equilibrium phenomena within the conventional continuum regime. To obtain physically consistent solutions across scales, we employ the unified gas-kinetic scheme (UGKS) and compare its results with NS solutions computed using the gas-kinetic scheme (GKS). Discrepancies between UGKS and GKS solutions reveal pronounced non-equilibrium effects in regions where shock waves interact with boundary layers. For continuum flows at high Mach and low Reynolds numbers, such multiscale non-equilibrium transport becomes important, underscoring the need for multiscale methods in analysis and prediction.

JFM classification

Compressible Flows: Compressible boundary layers Compressible Flows: Shock waves

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 1038 , 10 July 2026 , A44
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NCCreative Common License - ND
This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives licence (https://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided that no alterations are made and the original article is properly cited. The written permission of Cambridge University Press or the rights holder(s) must be obtained prior to any commercial use and/or adaptation of the article.
Copyright
© The Author(s), 2026. Published by Cambridge University Press

1. Introduction

The viscous shock tube is a standard test case for assessing the accuracy of numerical methods in continuum-flow regimes. This configuration features a reflected shock interacting with a developing boundary layer and is relevant to applications such as high-speed aerodynamic facilities that investigate chemical relaxation in high-temperature gas mixtures or generate high-enthalpy reservoir conditions. In this setting, the reflected shock encounters an incident boundary layer whose properties have been altered by end-wall conditions, producing a complex, unsteady coupling between the shock and the boundary layer. When the boundary-layer stagnation pressure falls below that downstream of the normal reflected shock, separation occurs and a characteristic lambda-shaped shock system forms. The seminal analysis of this interaction was provided by Mark (Reference Mark1957), who used a simplified model, followed by experiments that corroborated his predictions. Subsequent detailed experimental studies in shock tubes in Davies & Wilson (Reference Davies and Wilson1969) and Matsuo, Kawagoe & Kage (Reference Matsuo, Kawagoe and Kage1974) refined and extended Mark’s original framework. However, at low Reynolds numbers, pronounced non-equilibrium effects arise that are not captured by the traditional Navier–Stokes equations.

Investigating viscous shock-tube problems at low Reynolds numbers (Re) is of paramount importance, driven by the growing demand for high-fidelity modelling in micro-scale gas dynamics and high-altitude flight regimes. In these environments, which are characterised by minute physical dimensions or extremely low gas densities, the flow physics deviates significantly from the classical continuum assumption. Prominent engineering applications include micro-thrusters in micro-electro-mechanical systems (MEMS) (Xu et al. Reference Xu, Zhang, Li, Liu, Ye and Shen2024), where pressure wave propagation is inherently viscous dominated, and the early-stage atmospheric re-entry of hypersonic vehicles (Li et al. Reference Li, Peng, Ma, Dang, Tang and Sun2019), where rarefaction effects are critical. Unlike high-Re flows, these low-Re conditions exhibit strong coupling between shock structures and thick boundary layers, with non-equilibrium effects and kinetic-scale interactions dictating the overall flow field. Consequently, the low-Re shock tube serves as an essential canonical model for evaluating the capability of multiscale numerical schemes to capture these non-classical features. Furthermore, its relatively simple physical set-up and distinct flow characteristics make it an ideal, easily standardised computational benchmark. Recent research (Tian & Wu Reference Tian and Wu2026) demonstrates that even subtle rarefaction effects can induce significant errors in surface stress and heat flux predictions when using the Navier–Stokes (NS) equations under turbulent conditions, thereby exposing the intrinsic limitations of the NS framework in non-equilibrium flows. Motivated by these challenges, this work investigates non-equilibrium phenomena within the viscous shock tube by comparing continuum NS solvers against a multiscale approach.

Multiscale gas flows arise frequently in aerospace engineering (Bird Reference Bird1994; Xu Reference Xu2021) and in MEMS (Senturia, Azuru & White Reference Senturia, Azuru and White1997; Wang et al. Reference Wang, Liu, Zhuo and Zhong2022). Across all Knudsen-number regimes, accurate prediction requires properly accounting for both free streaming and interparticle collisions. Because non-equilibrium effects demand additional degrees of freedom, they are naturally described within a kinetic formulation. The Boltzmann equation is the fundamental governing equation of rarefied gas dynamics. In principle, it captures multiscale physics over the full range of Knudsen numbers by resolving phenomena at the mean-free-path and collision time scales. Mainstream numerical approaches for non-equilibrium gas flows can be broadly classified as stochastic methods – such as the direct simulation Monte Carlo method (Bird Reference Bird1963, Reference Bird1998) – and deterministic methods – such as the discrete velocity method (DVM) (Chu Reference Chu1965; Yang & Huang Reference Yang and Huang1995). In recent years, a class of gas-kinetic methods has been developed, including the DVM-based unified gas-kinetic scheme (UGKS) (Xu & Huang Reference Xu and Huang2010; Huang, Xu & Yu Reference Huang, Xu and Yu2012), the particle-based unified gas-kinetic wave–particle (UGKWP) method (Zhu et al. Reference Zhu, Liu, Zhong and Xu2019; Liu, Zhu & Xu Reference Liu, Zhu and Xu2020), the discrete unified gas-kinetic scheme (DUGKS) (Guo, Xu & Wang Reference Guo, Xu and Wang2013) and the discrete UGKWP method (Yang et al. Reference Yang, Li, Shu, Liu, Liu and Wu2023). By coupling particle collisions and free transport within the evolution process, these methods relax the stringent mesh and time-step requirements for obtaining accurate solutions across regimes, especially in collision-dominated flows. Recently, Guo et al. (Reference Guo, Xu and Zhu2026a , Reference Guo, Zhu and Xub ) presented a unified perspective on kinetic schemes, such as UGKS, UGKWP and DUGKS, and clarified how their governing equations differ from the classical NS formulation. In this paper, we examine these differences through the viscous shock-tube problem by comparing a continuum solver (GKS) with a multiscale solver (UGKS).

The paper is organised as follows. Section 2 presents the finite-volume formulation for conservation laws and the Bhatnagar–Gross–Krook (BGK) model. The BGK-based gas-kinetic scheme (GKS) and UGKS will be presented. Section 3 presents the verification and validation of the GKS and UGKS in the cavity flow and one-dimensional shock-tube problem. Section 4 describes the viscous shock-tube configuration, including the geometry, boundary conditions and initial conditions. The GKS and UGKS results will be presented to highlight the non-equilibrium effects in the viscous shock-tube problem. Finally, § 5 summarises the main conclusions.

2. Methodology

Traditional computational fluid dynamics methods numerically solve the NS equations to determine the flow field. However, these equations lose their validity in certain regimes, such as in rarefied gas flows. To bridge this gap, direct modelling-based numerical methods are employed.

The finite-volume method is used in this paper. The domain is divided into cells, for each cell upper Omega Subscript i $\varOmega _i$ , the governing equation of cell-averaged conservative values bold italic upper W Subscript i Baseline equals left parenthesis rho Subscript i Baseline comma rho Subscript i Baseline bold italic upper V Subscript i Baseline comma rho Subscript i Baseline upper E Subscript i Baseline right parenthesis Superscript upper T $\boldsymbol{W}_i=(\rho _i,\rho _i \boldsymbol{V}_i,\rho _i E_i)^T$ , i.e. density, momentum and energy, is given by

(2.1) upper W Subscript i Superscript n plus 1 Baseline equals upper W Subscript i Superscript n Baseline minus StartFraction 1 Over StartAbsoluteValue upper Omega Subscript i Baseline EndAbsoluteValue EndFraction integral Subscript t Superscript n Baseline Superscript t Superscript n plus 1 Baseline Baseline sigma summation Underscript j element of upper N left parenthesis i right parenthesis Endscripts upper S Subscript italic ij Baseline bold italic upper F Subscript italic ij Baseline normal d t comma \begin{equation} W_i^{n+1}=W_i^n -\frac {1}{|\varOmega _i|}\int _{t^n}^{t^{n+1}}\sum _{j\in N(i)} S_{\textit{ij}} \boldsymbol{F}_{\textit{ij}} {\mathrm{\,d}} t, \end{equation}

where upper W Subscript i Superscript n $W_i^n$ and upper W Subscript i Superscript n plus 1 $W_i^{n+1}$ are cell-averaged conservative values at times t Superscript n $t^n$ and t Superscript n plus 1 $t^{n+1}$ , StartAbsoluteValue upper Omega Subscript i Baseline EndAbsoluteValue $|\varOmega _i|$ denotes the volume of cell i $i$ , upper N left parenthesis i right parenthesis $N(i)$ is the set of all interface-adjacent neighbouring cells of cell i $i$ and j $j$ is one of the neighbouring cells of i $i$ . The interface between cells i $i$ and j $j$ is labelled i j $ij$ and has area upper S Subscript italic ij $S_{\textit{ij}}$ . The flux bold italic upper F Subscript italic ij $\boldsymbol{F}_{\textit{ij}}$ is the macroscopic flux through interface i j $ij$ . It is important to note that (2.1) represents the conservation law of macroscopic conservative quantities at the discrete level, which is a fundamental physical law valid at all temporal and spatial scales. Thus, the key point is to accurately describe the macroscopic flux bold italic upper F Subscript italic ij $\boldsymbol{F}_{\textit{ij}}$ through interface i j $ij$ . According to the gas-kinetic theory, the macroscopic flux can be determined by taking moments of the microscopic flux at the cell interface i j $ij$ through

(2.2) StartLayout 1st Row bold italic upper F Subscript italic ij Baseline equals integral bold italic v bold dot bold italic n Subscript italic ij Baseline f Subscript italic ij Baseline upper Psi normal d bold italic v normal d bold italic xi comma EndLayout \begin{align} \boldsymbol{F}_{\textit{ij}}=\int \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}_{\textit{ij}} {f}_{\textit{ij}}\varPsi {\mathrm{\,d}} \boldsymbol{v} {\mathrm{\,d}} \boldsymbol{\xi }, \end{align}

where bold italic n Subscript italic ij $\boldsymbol{n}_{\textit{ij}}$ is the normal vector of interface i j $ij$ , f Subscript italic ij ${f}_{\textit{ij}}$ is the distribution function at the cell interface i j $ij$ and upper Psi equals left parenthesis 1 comma bold italic v comma 1 divided by 2 left parenthesis bold italic v squared plus bold italic xi squared right parenthesis right parenthesis Superscript upper T $\varPsi =(1,\boldsymbol{v},1/2(\boldsymbol{v}^2+\boldsymbol{\xi }^2))^T$ represents the collision invariants. In bold italic xi squared equals xi 1 squared plus xi 2 squared plus bold midline horizontal ellipsis plus xi Subscript upper K Superscript 2 $\boldsymbol{\xi }^2=\xi _1^2+\xi _2^2+\boldsymbol{\cdots} +\xi _K^2$ and normal d bold italic xi equals normal d xi 1 normal d xi 2 bold midline horizontal ellipsis normal d xi Subscript upper K Baseline ${\mathrm{\,d}}\boldsymbol{\xi }={\mathrm{\,d}}\xi _1{\mathrm{\,d}}\xi _2\boldsymbol{\cdots} {\mathrm{\,d}}\xi _K$ , upper K $K$ is the number of internal degrees of freedom. The parameter bold italic v $\boldsymbol{v}$ is the microscopic particle translation velocity, and bold italic xi $\boldsymbol{\xi }$ is the internal variable for the rotational and vibrational degrees of freedom. The conservative values bold italic upper W $\boldsymbol{W}$ are obtained by

(2.3) StartLayout 1st Row bold italic upper W equals integral f upper Psi normal d bold italic v normal d bold italic xi period EndLayout \begin{align} \boldsymbol{W}=\int {f}\varPsi {\mathrm{\,d}} \boldsymbol{v} {\mathrm{\,d}} \boldsymbol{\xi }. \end{align}

Here, the BGK relaxation model is used to describe the process of distribution function evolution at the cell interface i j $ij$

(2.4) StartFraction partial differential f Over partial differential t EndFraction plus bold italic v bold dot bold nabla f equals StartFraction f Superscript plus Baseline minus f Over tau EndFraction comma \begin{equation} \frac {\partial f}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\!f =\frac {f^+-f}{\tau }, \end{equation}

where f equals f left parenthesis bold italic x comma t comma bold italic v comma bold italic xi right parenthesis $f=f(\boldsymbol{x},t,\boldsymbol{v},\boldsymbol{\xi })$ is the distribution function for gas molecules at physical space bold italic x $\boldsymbol{x}$ and time t $t$ . Here, tau $\tau$ is the particle collision time and f Superscript plus $f^+$ is the modified equilibrium distribution function. The modified equilibrium distribution function is given by the Shakhov model

(2.5) StartLayout 1st Row f Superscript plus Baseline equals g left bracket 1 plus left parenthesis 1 minus italic probability right parenthesis 2 lamda bold italic c bold dot bold italic q left parenthesis 2 lamda c squared minus 5 right parenthesis divided by left parenthesis 5 p right parenthesis right bracket equals g plus g Superscript plus Baseline comma EndLayout \begin{align} f^+=g\big [ 1+(1-{\textit{Pr}})2\lambda \boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{q}\big(2\lambda c^2-5\big)/(5p)\big]=g+g^+, \end{align}

where g $g$ is the Maxwellian distribution, italic probability ${\textit{Pr}}$ is the Prandtl number, bold italic c equals bold italic v minus bold italic upper V $\boldsymbol{c}=\boldsymbol{v}-\boldsymbol{V}$ is the peculiar velocity, bold italic upper V $\boldsymbol{V}$ is the macroscopic velocity, bold italic q $\boldsymbol{q}$ is the heat flux, lamda equals m divided by left parenthesis 2 k Subscript upper B Baseline upper T right parenthesis $\lambda =m/(2k_B T)$ , m $m$ is molecular mass, k Subscript upper B $k_B$ is the Boltzmann constant and upper T $T$ is the temperature. The Maxwellian distribution is

(2.6) StartLayout 1st Row g equals rho left parenthesis StartFraction lamda Over pi EndFraction right parenthesis Superscript StartFraction upper K plus upper D Over 2 EndFraction Baseline e Superscript minus lamda left parenthesis left parenthesis bold italic v minus bold italic upper V right parenthesis squared plus bold italic xi squared right parenthesis Baseline comma EndLayout \begin{align} g=\rho {\left (\frac {\lambda }{\pi }\right )}^{\frac {K+D}{2}}e^{-\lambda ((\boldsymbol{v}-\boldsymbol{V})^2+\boldsymbol{\xi }^2)}, \end{align}

where upper D $D$ is the spatial dimension. The heat flux is obtained through the distribution function

(2.7) StartLayout 1st Row bold italic q equals one half integral left parenthesis bold italic v minus bold italic upper V right parenthesis left parenthesis StartAbsoluteValue bold italic v minus bold italic upper V EndAbsoluteValue squared plus bold italic xi squared right parenthesis f normal d bold italic v normal d bold italic xi period EndLayout \begin{align} \boldsymbol{q}=\frac {1}{2}\int (\boldsymbol{v}-\boldsymbol{V})(|\boldsymbol{v}-\boldsymbol{V}|^2+\boldsymbol{\xi }^2)f{\mathrm{\,d}}\boldsymbol{v} {\mathrm{\,d}} \boldsymbol{\xi }. \end{align}

The collision term satisfies the compatibility condition

(2.8) StartLayout 1st Row integral left parenthesis f Superscript plus Baseline minus f right parenthesis upper Psi normal d bold italic v normal d bold italic xi equals 0 period EndLayout \begin{align} \int (f^+-f)\varPsi {\mathrm{\,d}} \boldsymbol{v}{\mathrm{\,d}}\boldsymbol{\xi }{=0}. \end{align}

To construct the numerical flux at the cell interface bold italic x 0 equals left parenthesis 0 comma 0 comma 0 right parenthesis Superscript upper T $\boldsymbol{x}_0=(0,0,0)^T$ , the integral solution along the characteristic line bold italic x prime equals bold italic x 0 minus bold italic v left parenthesis t minus t Superscript prime Baseline right parenthesis $\boldsymbol{x}^\prime =\boldsymbol{x}_0-\boldsymbol{v}(t-t^\prime )$ of the BGK equation (2.4) gives

(2.9) StartLayout 1st Row f left parenthesis bold italic x comma t comma bold italic xi right parenthesis equals StartFraction 1 Over tau EndFraction integral Subscript t 0 Superscript t Baseline f Superscript plus Baseline left parenthesis bold italic x prime comma t Superscript prime Baseline right parenthesis e Superscript minus left parenthesis t minus t Super Superscript prime Superscript right parenthesis divided by tau Baseline normal d t Superscript prime Baseline plus e Superscript negative t divided by tau Baseline f 0 left parenthesis bold italic x minus bold italic v t right parenthesis comma EndLayout \begin{align} f(\boldsymbol{x},t,\boldsymbol{\xi })=\frac {1}{\tau }\int _{t_0}^t f^+(\boldsymbol{x}^\prime ,t^\prime )e^{-(t-t^\prime )/\tau }{\mathrm{\,d}} t^\prime + e^{-t/\tau }f_{0}(\boldsymbol{x}-\boldsymbol{v} t), \end{align}

where f 0 left parenthesis bold italic x right parenthesis $f_{0}(\boldsymbol{x})$ is the initial gas distribution function at the beginning of each step t Superscript n $t^n$ , and f Superscript plus Baseline left parenthesis bold italic x comma t right parenthesis $f^+(\boldsymbol{x},t)$ is the effective equilibrium state distributed in space and time. The integral solution provides a multiscale model of the evolution from an initial non-equilibrium distribution f $f$ to an equilibrium distribution f Superscript plus $f^+$ via collisions. To achieve second-order accuracy, the initial distribution function f 0 left parenthesis bold italic x right parenthesis $f_{0}(\boldsymbol{x})$ is approximated as

(2.10) StartLayout 1st Row f 0 left parenthesis bold italic x right parenthesis equals StartLayout Enlarged left brace 1st Row 1st Column f Subscript upper G Superscript l Baseline plus bold italic x bold dot bold nabla f Subscript upper G Superscript l Baseline comma 2nd Column v Subscript n Baseline greater than 0 comma 2nd Row 1st Column f Subscript upper G Superscript r Baseline plus bold italic x bold dot bold nabla f Subscript upper G Superscript r Baseline comma 2nd Column v Subscript n Baseline less than 0 comma EndLayout EndLayout \begin{align} f_{0}(\boldsymbol{x}) = \begin{cases} f^l_G+\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{\nabla }\!f^l_G,&v_n\gt 0,\\[5pt] f^r_G+\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{\nabla }\!f^r_G,&v_n\lt 0,\\ \end{cases} \end{align}

where f Subscript upper G Superscript l $f^l_G$ and f Subscript upper G Superscript r $f^r_G$ are the reconstructed initial distribution functions at the left and right sides of the interface. The equilibrium state is approximated as

(2.11) StartLayout 1st Row f Superscript plus Baseline left parenthesis bold italic x comma t right parenthesis almost equals g 0 plus g 0 Superscript plus Baseline plus bold italic x bold dot bold nabla g 0 plus StartFraction partial differential g 0 Over partial differential t EndFraction t comma EndLayout \begin{align} f^+(\boldsymbol{x},t) \approx g_{0}+ g^+_{0} + \boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{\nabla }g_{0}+\frac {\partial g_{0}}{\partial t}t, \end{align}

where g 0 $g_{0}$ is the Maxwellian distribution function obtained from the conservative flow variables of colliding particles from both sides of the interface. For the gas-kinetic scheme, the initial distribution function f Subscript upper G $f_G$ is approximated by the Chapman–Enskog expansion

(2.12) StartLayout 1st Row f Subscript upper G Baseline equals g minus tau left parenthesis bold italic v bold dot bold nabla g plus StartFraction partial differential g Over partial differential t EndFraction right parenthesis period EndLayout \begin{align} f_G = g - \tau \! \left( \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }g + \frac {\partial g}{\partial t}\right)\!. \end{align}

In this limit, the GKS recovers the macroscopic NS solutions. To overcome the limitations of GKS in rarefied regions, the UGKS is implemented to enhance predictive accuracy. For the UGKS, the distribution function also needs to be updated with the following equation:

(2.13) StartLayout 1st Row f Subscript i comma k Superscript n plus 1 Baseline equals left parenthesis 1 plus StartFraction normal upper Delta t Over 2 tau EndFraction right parenthesis Superscript negative 1 Baseline left bracket f Subscript i comma k Superscript n Baseline plus StartFraction 1 Over upper Omega Subscript i Baseline EndFraction sigma summation Underscript j element of upper N left parenthesis i right parenthesis Endscripts upper S Subscript italic ij Baseline integral Subscript t Superscript n Baseline Superscript t Superscript n plus 1 Baseline Baseline bold italic v Subscript k Baseline bold dot bold italic n Subscript italic ij Baseline f Subscript i j comma k Baseline normal d t plus StartFraction normal upper Delta t Over 2 EndFraction left parenthesis StartFraction f Subscript i comma k Superscript plus comma n plus 1 Baseline Over tau EndFraction plus StartFraction f Subscript i comma k Superscript plus comma n Baseline minus f Subscript i comma k Superscript n Baseline Over tau EndFraction right parenthesis right bracket comma EndLayout \begin{align} f_{i,k}^{n+1} =\left (1+\frac {\Delta t}{2\tau }\right )^{-1}\left [f_{i,k}^{n} + \frac {1}{\varOmega _i}\sum _{j\in N(i)} S_{\textit{ij}}\int _{t^n}^{t^{n+1}} \boldsymbol{v}_k\boldsymbol{\cdot }\boldsymbol{n}_{\textit{ij}}f_{ij,k}{\mathrm{\,d}} t+\frac {\Delta t}{2} \left (\frac {f_{i,k}^{+,n+1}}{\tau }+\frac {f_{i,k}^{+,n}-f_{i,k}^{n}}{\tau }\right )\right ]\!, \end{align}

where the index k $k$ denotes the discrete points in the particle velocity space and the integral over the velocity space is approximated by the summation in the discrete velocity space. More details of the GKS and the UGKS are demonstrated in earlier work (Xu Reference Xu2001; Xu & Huang Reference Xu and Huang2010; Xu Reference Xu2015).

Figure 1.

Cavity simulation using GKS and UGKS at upper R e equals 20 $Re=20$ . (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The upper U $U$ -velocity along the central vertical line and the upper V $V$ -velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Three panels show temperature contours and velocity profiles in a cavity simulation.

3. Scheme validation and verification

In this section, we evaluate the GKS and UGKS methods and validate the simulation results across a range of Mach and Reynolds numbers.

3.1. Cavity flow

In this subsection, the cavity flow problem, a classic benchmark for validating numerical methods, is employed to verify the accuracy of the GKS and UGKS. The computational set-up consists of a two-dimensional flow within a square domain containing a centrally located square obstacle. To comprehensively assess the methods, simulations are conducted at three Reynolds numbers: 20 comma 50 $20, 50$ and 100 $100$ . The working fluid is argon gas with a Prandtl number of italic probability equals 0.73 ${\textit{Pr}}=0.73$ . The isothermal wall temperature is fixed at upper T Subscript w Baseline equals 273 $T_w=273$ K, and the upper wall moves with a velocity of upper U Subscript w Baseline equals 50.0 $U_w=50.0$ m s−1. For both the GKS and UGKS, the physical domain is discretised using a 60 times 60 $60\times 60$ spatial mesh, with the first layer of near-wall cells refined to a height of 1 times 10 Superscript negative 3 $1\times 10^{-3}$ . Additionally, for the UGKS, the velocity space is discretised into 60 times 60 $60\times 60$ elements.

The results of the GKS and UGKS are shown in figures 1, 2 and 3. For all three Reynolds numbers, both GKS and UGKS exhibit similar velocity distributions. However, differences arise in temperature and heat flux distributions between the two methods. Notably, at a low Reynolds number (e.g. upper R e equals 20 $Re=20$ ), the heat flux predicted by UGKS deviates from Fourier’s law, exhibiting a counter-gradient phenomenon where heat flows from low-temperature to high-temperature regions. As the Reynolds number increases, the UGKS results gradually get back to the NS solutions.

Figure 2.

Cavity simulation using GKS and UGKS at upper R e equals 50 $Re=50$ . (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The upper U $U$ -velocity along the central vertical line and the upper V $V$ -velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Three panels show temperature contours and velocity profiles for cavity simulation using different methods.

Figure 3.

Cavity simulation using GKS and UGKS at upper R e equals 100 $Re=100$ . (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The upper U $U$ -velocity along the central vertical line and the upper V $V$ -velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Three contour and line graphs depict temperature and velocity simulations in cavity flows.

3.2. One-dimensional viscous shock tube

In order to validate the grid convergence solution of UGKS, the following one-dimensional shock-tube problem will be studied. In this study, the computational domain is defined as left bracket 0 comma 1 right bracket $[0,1]$ with all boundaries being adiabatic walls. The initial conditions are specified as

(3.1) StartLayout 1st Row left parenthesis rho comma p comma u comma v right parenthesis equals StartLayout Enlarged left brace 1st Row 1st Column left parenthesis 120.0 comma 120.0 divided by gamma comma 0.0 comma 0.0 right parenthesis comma 2nd Column if x less than 0.5 comma 2nd Row 1st Column left parenthesis 1.2 comma 1.2 divided by gamma comma 0.0 comma 0.0 right parenthesis comma 2nd Column if x greater than 0.5 comma EndLayout EndLayout \begin{align} (\rho ,p,u,v)=\begin{cases} (120.0, 120.0/\gamma , 0.0, 0.0), & \text{if } x\lt 0.5, \\[5pt] (1.2, 1.2/\gamma , 0.0, 0.0), & \text{if } x\gt 0.5, \\ \end{cases} \end{align}

where rho $\rho$ is the density, p $p$ is the pressure and left parenthesis u comma v right parenthesis $(u,v)$ are the velocities in the x $x$ and y $y$ directions, respectively. The non-dimensional gas equation of state is gamma p equals rho upper T $\gamma p = \rho T$ with a specific heat ratio of gamma equals 1.4 $\gamma =1.4$ . The Reynolds number is defined based on the domain length upper L equals 1 $L=1$ , reference density rho 0 equals 1.0 $\rho _0=1.0$ , reference velocity upper U 0 equals 1.0 $U_0=1.0$ and a constant viscosity, i.e. mu equals 1.0 divided by upper R e $\mu =1.0/{Re}$ . The mean free path is given by l equals mu divided by left parenthesis rho a right parenthesis equals 1 divided by left parenthesis rho a upper R e right parenthesis $l = \mu /(\rho a) = 1/(\rho a{Re})$ , where a equals StartRoot upper T EndRoot $a = \sqrt {T}$ denotes the speed of sound. The Prandtl number is set to italic probability equals 0.73 ${\textit{Pr}} = 0.73$ . The simulation time is t equals 1.0 $t=1.0$ , with the time step determined by a Courant-Friedrichs-Lewy number 0.5 $0.5$ .

The grid configuration for both the physical and velocity spaces is as follows. For two-dimensional simulations, the physical domain is discretised into 700 times 700 $700\times 700$ cells, and for UGKS, the velocity space left bracket negative 6 comma 6 right bracket squared $[-6, 6]^2$ is discretised into 80 times 80 $80 \times 80$ uniform grid points. For the one-dimensional reference case, which captures the evolution of key flow features such as the shock wave, contact discontinuity and rarefaction wave, and their interactions, the computational domain is divided into 700 cells along the x $x$ -axis, with the velocity space discretised into 80 elements. Only the left and right boundaries are considered in the one-dimensional case.

To validate the adequacy of the mesh resolution and ensure mesh independence, systematic studies have been conducted for both the physical and velocity-space discretisations. Figure 4 presents the density profiles at different times computed by UGKS with two levels of DVS resolutions: 80 elements in the range left bracket negative 6 comma 6 right bracket $[-6,6]$ , and 200 elements in the range left bracket negative 10 comma 10 right bracket $[-10,10]$ . The comparison confirms that the results are insensitive to further refinement of the velocity mesh, with good agreement between the fine and coarse discretisations.

Figure 4.

Discrete velocity-space (DVS) mesh independence test: comparison of density profiles using 80 elements in the range left bracket negative 6 comma 6 right bracket $[-6, 6]$ and 200 elements in the range left bracket negative 10 comma 10 right bracket $[-10, 10]$ for (a) upper R e equals 50 $Re=50$ , (b) upper R e equals 100 $Re=100$ and (c) upper R e equals 200 $Re=200$ .

Three line graphs compare density profiles using different mesh elements.

Figure 5 demonstrates the mesh independence solution from UGKS with 700 cells and 1400 cells in the one-dimensional physical domain. The test is performed at italic Re equals 200 $\textit{Re}=200$ . The results show that the converged physical solution can be obtained by UGKS using 700 mesh points.

Figure 5.

Physical-domain mesh independence test: comparison of density profiles using 700 and 1400 cells in the physical domain for upper R e equals 200 $Re=200$ .

Line graph comparing density profiles using different numbers of cells.

Therefore, unless otherwise stated, all one-dimensional simulations adopt 700 cells in the physical domain and 80 elements in the velocity space, while two-dimensional simulations use 700 times 700 $700 \times 700$ cells and 80 times 80 $80 \times 80$ grid points in the velocity space.

The one-dimensional case is considered to illustrate the difference in the main flow structure between GKS and UGKS for different Reynolds numbers. The evolution of density in the x minus t $x{-}t$ diagram is shown in figure 6.

Figure 6.

The x minus t $x{-}t$ diagram of the density by GKS and UGKS with different Reynolds numbers. The first row is the results of GKS, and the second row is the results of UGKS. From (a–f), the Reynolds numbers are 50, 100 and 200.

Six contour plots showing density variations with different Reynolds numbers.

A shock wave, followed by a contact discontinuity, propagates toward the low-pressure region (the right side), while a rarefaction wave propagates toward the high-pressure region (the left side). The shock Mach number is 2.37. The incident shock wave is relatively weak and reflects off the right wall at time t almost equals 0.2 $t\approx 0.2$ . Following this reflection, it interacts with the contact discontinuity. At later times, the contact discontinuity remains stationary near the right wall. Conversely, the reflected shock wave continues to propagate, beginning to interact with the incident rarefaction wave at time t equals 0.4 $t=0.4$ , and subsequently with the reflected rarefaction wave. During these interactions, the reflected shock wave maintains an approximately constant velocity.

A primary distinction between the two methods is that the GKS generally exhibits lower numerical dissipation. For instance, at the head of the rarefaction wave, the GKS may produce an overshoot (i.e. density values exceeding expected physical limits). This suggests that the GKS lacks sufficient dissipation to accurately capture physical diffusion within extreme flow structures, such as shocks and rarefaction waves, potentially leading to non-physical numerical artefacts. In contrast, the UGKS, as a multi-scale method, more comprehensively accounts for rarefaction effects and various diffusion processes. This effectively suppresses overshoots and provides a more accurate resolution of the shock and rarefaction structures. A more detailed comparison of these methods will be presented in the subsequent sections, where numerical solutions at various times are analysed to highlight specific differences in main flow structures, shock resolution and dissipation characteristics.

Figure 7.

The density profile of Reynolds numbers 50 and 100 by GKS and UGKS at different times.

Multiple line graphs depict density profiles at different times and Reynolds numbers.

The density profile at different times is shown in figure 7. At t equals 0.175 $t=0.175$ , the shock wave has not reached the right wall, so the density profile is not affected by it. The GKS and UGKS predict the same position of the shock wave and the contact discontinuity, and, for upper R e equals 50 $Re=50$ , the UGKS gives a thickness of the shock wave that is slightly thicker than the GKS results. At the time of t equals 0.4 $t=0.4$ , the shock wave has reflected off the right wall and interacts with the contact discontinuity. At Reynolds number 100, the GKS and UGKS results are close. At Reynolds number 50, the UGKS density is higher than the GKS results at the front of the shock wave. At the times of t equals 0.6 $t=0.6$ and t equals 1.0 $t=1.0$ , the shock wave interacts with the rarefaction wave. As shown in the figure enlargement, the density is higher in front of the shock wave and remains so at subsequent times. And little oscillation is observed in the shock structure at upper R e equals 100 $Re=100$ by GKS.

4. Viscous shock tube

Whereas the previous section analysed the evolution of the primary flow structures in the one-dimensional case, this section focuses on two-dimensional configurations to investigate the interaction between the main flow and the boundary layer. Transitioning to a two-dimensional framework enables the resolution of more complex flow phenomena, such as vortices. When the incident shock wave impinges on the right wall and reflects, the total pressure in the boundary layer ahead of the shock is lower than that in the region behind it. Consequently, the boundary layer cannot smoothly penetrate the reflected-shock region; instead, it separates from the surface, forming a separated flow region near the wall. The fluid passing through the oblique foot and tail shocks retains greater forward momentum than the fluid decelerated by the main reflected shock, resulting in the formation of a wall jet. The bifurcated foot and tail shocks merge with the main reflected shock at the triple point, forming a lamda $\lambda$ -shock wave. Finally, the entropy difference between the gas traversing the bifurcation zone and the gas processed by the single-shock system generates a distinct slip line or shear layer.

The density gradient of the viscous shock-tube problem with GKS and UGKS at the time of t equals 0.6 $t=0.6$ is shown in figure 8. At this time, the lamda $\lambda$ -shock wave was formed, and other lamda $\lambda$ -shock waves were formed by the interaction of the tail shock wave and the separated flow region. The second lamda $\lambda$ -shock wave was more obvious at the Reynolds number of 100. In addition, the triple point of the lamda $\lambda$ -shock predicted by UGKS is lower than the GKS results. The temperature profile along the line of y equals 0.095 $y=0.095$ at the time of t equals 0.6 $t=0.6$ is shown in figure 9(a). From the figure, the different positions of the tail shock wave by GKS and UGKS are observed. The position of the tail shock wave predicted by UGKS is behind the GKS results, which is caused by the different positions of the triple point of the lamda $\lambda$ -shock wave predicted by GKS and UGKS, and the difference is more obvious at the Reynolds number of 50, which means that, at lower Reynolds numbers, the rarefied effect is apparent and UGKS can resolve it. For the Reynolds number of 100, three shock waves caused by the two lamda $\lambda$ -shock waves are observed, while for the Reynolds number of 50, only two shock waves are observed. However, in the GKS results at Reynolds number 50, a flat temperature region is located at x equals 0.84 $x=0.84$ , indicating that the second lamda $\lambda$ -shock wave is weakly resolved; in the UGKS results, this flat region is resolved. The temperature contour of GKS and UGKS at the time of t equals 0.6 $t=0.6$ is shown in figure 10. The temperature contours show the low-temperature jet from the gas that encounters the oblique foot and tail shocks. The x $x$ -direction velocity profile along the line of x equals 0.9 $x=0.9$ at the time of t equals 0.6 $t=0.6$ is shown in figure 9(b). This figure shows the velocity profile of the low-temperature jet from the gas that encounters the oblique foot and tail shocks. The velocity predicted by GKS is higher than that by UGKS, consistent with the temperature profile, and the difference is more pronounced at a Reynolds number of 50.

Figure 8.

The contour of density gradient by GKS and UGKS at the time of t equals 0.6 $t=0.6$ . The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Four contour plots compare density gradients at different Reynolds numbers using two methods.

Figure 9.

The line profile at the time of t equals 0.6 $t=0.6$ . (a) The temperature profile along the line of y equals 0.095 $y=0.095$ , (b) the x $x$ -direction velocity profile along the line of x equals 0.9 $x=0.9$ .

Two line graphs compare temperature and velocity profiles using different methods and Reynolds numbers.

Figure 10.

The temperature contour of GKS and UGKS at the time of t equals 0.6 $t=0.6$ . The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Four contour plots compare temperature distributions for different Reynolds numbers and methods.

At the time of t equals 1.0 $t=1.0$ , shock waves move away from the wall, and more complicated flow structures are observed at the boundary layer. The density contours of GKS and UGKS at the time of t equals 1.0 $t=1.0$ are shown in figure 11. In addition, the density profile along the line of y equals 0.0 $y=0.0$ at the time of t equals 1.0 $t=1.0$ is shown in the figure 12. This profile is usually set as the reference profile to evaluate the accuracy of the numerical method. However, an unusual phenomenon is observed. The difference between GKS and UGKS is more obvious at the Reynolds number of 100, which is consistent with the density contour. The reason is that at a Reynolds number of 100, even though the flow condition is more continuous, the non-equilibrium effect is more apparent than at a Reynolds number of 50. At the Reynolds number of 50, the viscous effect is more apparent, but the flow is smoother; thus, the difference is not as obvious as at the Reynolds number of 100. As shown in the figure 13, the Knudsen-number ( Kn Subscript upper G l l Baseline equals l StartAbsoluteValue bold nabla rho EndAbsoluteValue divided by rho $\text{Kn}_{Gll}=l|\boldsymbol{\nabla }\rho |/\rho$ ) contours of UGKS with Reynolds numbers of 50 and 100 are shown, which indicates the non-equilibrium effect. The maximum Knudsen number is 0.077 for Reynolds number 50 and 0.053 for Reynolds number 100. The non-equilibrium effect appears in the region of the main shock wave, lamda $\lambda$ -shock wave and the boundary layer. Although the non-equilibrium effect of the main shock wave and lamda $\lambda$ -shock wave is not as obvious as the Reynolds number of 50, the non-equilibrium effect of the boundary layer is more apparent at the Reynolds number of 100. The main difference between GKS and UGKS at the Reynolds number of 100 is caused by the non-equilibrium effect of the boundary layer. This case shows that, in the continuum-flow regime, common computational fluid dynamics methods do not account for non-equilibrium effects, leading to deviations between numerical and physical results.

Figure 11.

The density contour of GKS and UGKS at the time of t equals 1.0 $t=1.0$ . The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Four contour plots compare density distributions for different Reynolds numbers and methods.

Figure 12.

The density profile along the line of y equals 0.0 $y=0.0$ at the time of t equals 1.0 $t=1.0$ .

A line graph showing density profiles for different Reynolds numbers.

Figure 13.

The Knudsen-number Kn Subscript upper G l l $\text{Kn}_{Gll}$ contour at the time of t equals 1.0 $t=1.0$ . (a) The results of Reynolds number 50. (b) The results of Reynolds number 100.

Heat maps showing Knudsen-number contours at different Reynolds numbers.

Figures 14 and 15 present the stress and heat flux tensors obtained using the UGKS at t equals 1.0 $t=1.0$ for Reynolds numbers of 50 and 100. Two formulations for these macroscopic quantities are compared: one derived directly from the moments of the velocity distribution function, and the other calculated using the Newtonian viscous stress model and Fourier’s law of heat conduction. The most significant deviations between the two methods occur at the shock waves, where strong non-equilibrium effects render the continuum hypothesis inadequate for accurately approximating stress and heat flux. Furthermore, at Re = 50, substantial discrepancies arise both near the wall and within the interaction region between the vortex and the main flow. In certain areas, the two methods even predict opposite signs for the stress and heat flux, a phenomenon particularly pronounced in the x $x$ -direction heat flux profiles. At Re = 100, however, the differences are primarily concentrated near the wall, which aligns with the discrepancies observed in the near-wall density distribution shown in figure 12.

Figure 14.

The sigma Subscript x y $\sigma _{xy}$ and q Subscript x $q_x$ and q Subscript y $q_y$ contours of UGKS at the time of t equals 1.0 $t=1.0$ at the Reynolds number of 50. (a) Calculated by the momentum of distribution function. (b) Calculated by the Newton stress tensor and Fourier’s law of heat conduction. (c) The ratio of (a) and (b).

Multiple contour plots depict stress and heat flux at a Reynolds number of 50.

Figure 15.

The sigma Subscript x y $\sigma _{xy}$ and q Subscript x $q_x$ and q Subscript y $q_y$ contours of UGKS at the time of t equals 1.0 $t=1.0$ at the Reynolds number of 100. (a) Calculated by the momentum of distribution function. (b) Calculated by the Newton stress tensor and Fourier’s law of heat conduction. (c) The ratio of (a) and (b).

Multiple contour plots depict stress and heat flux at a Reynolds number of 100.

5. Conclusion

In this study, we investigated the viscous shock-tube problem in the low-Reynolds-number regime, demonstrating that non-equilibrium effects can be significant even within the nominal continuum-flow regime. To this end, we compared a continuum solver, the GKS, with a multiscale solver, the UGKS, across both one- and two-dimensional configurations. The results indicate that non-equilibrium effects are most pronounced in the one-dimensional case at upper R e equals 50 $Re=50$ and during the early stages of the two-dimensional evolution, primarily driven by non-equilibrium behaviour within the shock structure. In contrast, during the later stages of the two-dimensional evolution, the discrepancies between the solvers become more pronounced at upper R e equals 100 $Re=100$ , where non-equilibrium effects associated with the developing boundary layer dominate. Overall, these findings reveal that the viscous shock-tube problem – commonly employed as a standard benchmark in the continuum regime – exhibits inherently multiscale characteristics. Consequently, non-equilibrium effects must be carefully accounted for in numerical simulations, even when the flow is traditionally classified as a continuum. This observation challenges conventional assumptions and underscores the necessity of multiscale kinetic approaches for accurate flow prediction.

Acknowledgements

We would like to thank X. Ji for the helpful discussion and valuable suggestions.

Funding

This work was supported by the National Key R&D Program of China (Grant No. 2022YFA1004500), the National Natural Science Foundation of China (92371107) and the Hong Kong Research Grant Council (16208324).

Declaration of interests

The authors report no conflict of interest.

References

Bird, G.A. 1963 Approach to translational equilibrium in a rigid sphere gas. Phys. Fluids 6, 1518–1519.10.1063/1.1710976CrossRefGoogle Scholar
Bird, G.A. 1994 Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford University Press.10.1093/oso/9780198561958.001.0001CrossRefGoogle Scholar
Bird, G.A. 1998 Recent advances and current challenges for DSMC. Comput. Math. Appl. 35 (1–2), 1–14.10.1016/S0898-1221(97)00254-XCrossRefGoogle Scholar
Chu, C.K. 1965 Kinetic-theoretic description of the formation of a shock wave. Phys. Fluids 8 (1), 12–22.10.1063/1.1761077CrossRefGoogle Scholar
Davies, L. & Wilson, J.L. 1969 Influence of reflected shock and boundary-layer interaction on shock-tube flows. Phys. Fluids 12 (5), I–37.10.1063/1.1692625CrossRefGoogle Scholar
Guo, Z., Xu, K. & Wang, R. 2013 Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case. Phys. Rev. E 88 (3), 033305.10.1103/PhysRevE.88.033305CrossRefGoogle ScholarPubMed
Guo, Z., Xu, K. & Zhu, Y. 2026 a A unified gas-kinetic framework from Boltzmann to Navier–Stokes scales. Adv. Aerodyn. 8 (1), 8.10.1186/s42774-025-00252-1CrossRefGoogle Scholar
Guo, Z., Zhu, Y. & Xu, K. 2026 b Kinetic representation of the unified gas-kinetic wave-particle method and beyond. Commun. Comput. Phys. 39 (5), 1512–1535.10.4208/cicp.OA-2024-0192CrossRefGoogle Scholar
Huang, J.-C., Xu, K. & Yu, P. 2012 A unified gas-kinetic scheme for continuum and rarefied flows II: multi-dimensional cases. Commun. Comput. Phys. 12 (3), 662–690.10.4208/cicp.030511.220911aCrossRefGoogle Scholar
Li, Z.-H., Peng, A.-P., Ma, Q., Dang, L.-N., Tang, X.-W. & Sun, X.-Z. 2019 Gas-kinetic unified algorithm for computable modeling of boltzmann equation and application to aerothermodynamics for falling disintegration of uncontrolled Tiangong No. 1 spacecraft. Adv. Aerodyn. 1 (1), 4.10.1186/s42774-019-0009-4CrossRefGoogle Scholar
Liu, C., Zhu, Y. & Xu, K. 2020 Unified gas-kinetic wave-particle methods I: continuum and rarefied gas flow. J. Comput. Phys. 401, 108977.10.1016/j.jcp.2019.108977CrossRefGoogle Scholar
Mark, H. 1957 The Interaction of a Reflected Shock Wave with the Boundary Layer in a Shock Tube. Cornell University.Google Scholar
Matsuo, K., Kawagoe, S. & Kage, K. 1974 The interaction of a reflected shock wave with the boundary layer in a shock tube. Bull. JSME 17 (110), 1039–1046.10.1299/jsme1958.17.1039CrossRefGoogle Scholar
Senturia, S.D., Azuru, N. & White, J. 1997 Simulating the behavior of MEMS devices: computational methods and needs. IEEE Comput. Sci. Engng 4 (1), 30–43.10.1109/99.590854CrossRefGoogle Scholar
Tian, S. & Wu, L. 2026 Small rarefaction, large consequences: limits of Navier Stokes turbulence simulations. arXiv preprint arXiv: 2602.08770.Google Scholar
Wang, Y., Liu, S., Zhuo, C. & Zhong, C. 2022 Investigation of nonlinear squeeze-film damping involving rarefied gas effect in micro-electro-mechanical systems. Comput. Math. Appl. 114, 188–209.10.1016/j.camwa.2022.03.045CrossRefGoogle Scholar
Xu, K. 2001 A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method. J. Comput. Phys. 171 (1), 289–335.10.1006/jcph.2001.6790CrossRefGoogle Scholar
Xu, K. 2015 Direct modeling for computational fluid dynamics. Acta Mechanica Sin. 31, 303–318.10.1007/s10409-015-0453-2CrossRefGoogle Scholar
Xu, K. 2021 A Unified Computational Fluid Dynamics Framework from Rarefied to Continuum Regimes. Cambridge University Press.10.1017/9781108877534CrossRefGoogle Scholar
Xu, K. & Huang, J.-C. 2010 A unified gas-kinetic scheme for continuum and rarefied flows. J. Comput. Phys. 229 (20), 7747–7764.10.1016/j.jcp.2010.06.032CrossRefGoogle Scholar
Xu, J., Zhang, J., Li, F., Liu, S., Ye, Y. & Shen, R. 2024 A review on solid propellant micro-thruster array based on MEMS technology. FirePhysChem 4 (2), 95–106.10.1016/j.fpc.2023.03.002CrossRefGoogle Scholar
Yang, J.-Y. & Huang, J.-C. 1995 Rarefied flow computations using nonlinear model Boltzmann equations. J. Comput. Phys. 120 (2), 323–339.10.1006/jcph.1995.1168CrossRefGoogle Scholar
Yang, L., Li, Z., Shu, C., Liu, Y., Liu, W. & Wu, J. 2023 Discrete unified gas-kinetic wave-particle method for flows in all flow regimes. Phys. Rev. E 108 (1), 015302.10.1103/PhysRevE.108.015302CrossRefGoogle ScholarPubMed
Zhu, Y., Liu, C., Zhong, C. & Xu, K. 2019 Unified gas-kinetic wave-particle methods. II. Multiscale simulation on unstructured mesh. Phys. Fluids 31 (6), 067105.10.1063/1.5097645CrossRefGoogle Scholar
Figure 0

Figure 1. Cavity simulation using GKS and UGKS at Re=20$Re=20$. (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The U$U$-velocity along the central vertical line and the V$V$-velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Figure 1

Figure 2. Cavity simulation using GKS and UGKS at Re=50$Re=50$. (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The U$U$-velocity along the central vertical line and the V$V$-velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Figure 2

Figure 3. Cavity simulation using GKS and UGKS at Re=100$Re=100$. (a) The temperature contour and heat flux using GKS. (b) The temperature contour and heat flux using UGKS. (c) The U$U$-velocity along the central vertical line and the V$V$-velocity along the central horizontal line; symbols: GKS, lines: UGKS.

Figure 3

Figure 4. Discrete velocity-space (DVS) mesh independence test: comparison of density profiles using 80 elements in the range [−6,6]$[-6, 6]$ and 200 elements in the range [−10,10]$[-10, 10]$ for (a) Re=50$Re=50$, (b) Re=100$Re=100$ and (c) Re=200$Re=200$.

Figure 4

Figure 5. Physical-domain mesh independence test: comparison of density profiles using 700 and 1400 cells in the physical domain for Re=200$Re=200$.

Figure 5

Figure 6. The x−t$x{-}t$ diagram of the density by GKS and UGKS with different Reynolds numbers. The first row is the results of GKS, and the second row is the results of UGKS. From (a–f), the Reynolds numbers are 50, 100 and 200.

Figure 6

Figure 7. The density profile of Reynolds numbers 50 and 100 by GKS and UGKS at different times.

Figure 7

Figure 8. The contour of density gradient by GKS and UGKS at the time of t=0.6$t=0.6$. The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Figure 8

Figure 9. The line profile at the time of t=0.6$t=0.6$. (a) The temperature profile along the line of y=0.095$y=0.095$, (b) the x$x$-direction velocity profile along the line of x=0.9$x=0.9$.

Figure 9

Figure 10. The temperature contour of GKS and UGKS at the time of t=0.6$t=0.6$. The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Figure 10

Figure 11. The density contour of GKS and UGKS at the time of t=1.0$t=1.0$. The first two figures are the results of Reynolds number 50, and the last two figures are the results of Reynolds number 100. For each Reynolds number, the first figure is the result of UGKS, and the second figure is the result of GKS.

Figure 11

Figure 12. The density profile along the line of y=0.0$y=0.0$ at the time of t=1.0$t=1.0$.

Figure 12

Figure 13. The Knudsen-number KnGll$\text{Kn}_{Gll}$ contour at the time of t=1.0$t=1.0$. (a) The results of Reynolds number 50. (b) The results of Reynolds number 100.

Figure 13

Figure 14. The σxy$\sigma _{xy}$ and qx$q_x$ and qy$q_y$ contours of UGKS at the time of t=1.0$t=1.0$ at the Reynolds number of 50. (a) Calculated by the momentum of distribution function. (b) Calculated by the Newton stress tensor and Fourier’s law of heat conduction. (c) The ratio of (a) and (b).

Figure 14

Figure 15. The σxy$\sigma _{xy}$ and qx$q_x$ and qy$q_y$ contours of UGKS at the time of t=1.0$t=1.0$ at the Reynolds number of 100. (a) Calculated by the momentum of distribution function. (b) Calculated by the Newton stress tensor and Fourier’s law of heat conduction. (c) The ratio of (a) and (b).