Unified Gas-kinetic Wave-Particle Methods II: Multiscale Simulation on Unstructured Mesh

TL;DR

This paper introduces a unified gas-kinetic wave-particle (UGKWP) method on unstructured meshes for efficient multiscale flow simulation, seamlessly bridging continuum and rarefied regimes with significant computational savings.

Contribution

The paper develops a novel UGKWP method that adaptively combines wave and particle descriptions, enabling accurate multiscale flow simulation on unstructured meshes without restrictive time step or mesh size constraints.

Findings

Achieves UGKS-level accuracy across all flow regimes.

Reduces computational cost and memory usage by several orders of magnitude.

Successfully validated through multiple numerical test cases.

Abstract

In this paper, we present a unified gas-kinetic wave-particle (UGKWP) method on unstructured mesh for multiscale simulation of continuum and rarefied flow. Inheriting from the multicale transport in the unified gas-kinetic scheme (UGKS), the integral solution of kinetic model equation is employed in the construction of UGKWP method to model the flow physics in the cell size and time step scales. A novel wave-particle adaptive formulation is introduced in the UGKWP method to describe the flow dynamics in each control volume. The local gas evolution is constructed through the dynamical interaction of the deterministic hydrodynamic wave and the stochastic kinetic particle. Within the resolution of cell size and time step, the decomposition, interaction, and evolution of the hydrodynamic wave and the kinetic particle depend on the ratio of the time step to the local particle collision time. In the rarefied flow regime, the flow physics is mainly recovered by the discrete particles and the UGKWP method performs as a stochastic particle method. In the continuum flow regime, the flow behavior is solely followed by macroscopic variable evolution and the UGKWP method becomes a gas-kinetic hydrodynamic flow solver for the viscous and heat-conducting Navier--Stokes solutions. In different flow regimes, many numerical test cases are computed to validate the UGKWP method on unstructured mesh. The UGKWP method can get the same UGKS solutions in all Knudsen regimes without the requirement of the time step and mesh size being less than than the particle collision time and mean free path. With an automatic wave-particle decomposition, the UGKWP method becomes very efficient. For example, at Mach number 30 and Knudsen number 0.1, in comparison with UGKS several-order-of-magnitude reductions in computational cost and memory requirement have been achieved by UGKWP.

Tables

Table 1. Table 1: Computational cost for hypersonic flow at Ma = 5 Ma 5 {\rm Ma}=5 and Kn = 0.1 Kn 0.1 {\rm Kn}=0.1 around a circle cylinder.

	UGKS	UGKWP	ratio
CPU time	$45h$ $16min$	$2h$ $10min$	$20.9$
Memory	$4.9\mathrm{GB}$	$277MB$	$18.1$
Steps	$5000$111$5000$ steps of first-order calculation for a better initial state in UGKS computation. $+$ $35000$	$30000$ $+$ $10000$222$10000$ steps of averaging process in the UGKWP simulation.

Table 2. Table 2: Computational cost for hypersonic flow at Ma = 5 Ma 5 {\rm Ma}=5 and Kn = 1 Kn 1 {\rm Kn}=1 around a circle cylinder.

	UGKS	UGKWP	ratio
CPU time	$45h$ $16min$	$2h$ $42min$	$16.8$
Memory	$4.9\mathrm{GB}$	$310MB$	$16.2$
Steps	$5000$111$5000$ steps of first-order calculation for a better initial state in UGKS computation. $+$ $35000$	$30000$ $+$ $10000$222$10000$ steps of averaging process in the UGKWP simulation.

Equations

$${\mathbf{w}}_i^{n+1}={\mathbf{w}}_i^n-\frac{\Delta t}{{\Omega }_i}\sum _{j\in N(i)}{\mathbf{F}}_{ij}{S}_{ij},$$

$$f_i^{n+1}=f_i^n-\frac{1}{{\Omega }_i}\sum _{j\in N(i)}{\int }_0^{\Delta t}\mathbf{u}\cdot {\mathbf{n}}_{ij}{f}_{ij}(t)S_{ij}dt+{\int }_0^{\Delta t}\mathcal{J}(f,f)dt,$$

$${\mathbf{w}}_i=\int f_i\mathbf{\psi }(\mathbf{u})d\mathbf{u},$$

$${\mathbf{F}}_{ij}=\frac{1}{\Delta t}{\int }_0^{\Delta t}\int \mathbf{u}\cdot {\mathbf{n}}_{ij}{f}_{ij}(t)\mathbf{\psi }(\mathbf{u})d\mathbf{u}dt,$$

$$\int \mathcal{J}(f,f)\mathbf{\psi }(\mathbf{u})d\mathbf{u}=0.$$

$$\mathcal{J}(f,f)=\frac{g-f}{\tau },$$

$$g=\rho {\left(\frac{\lambda }{\pi }\right)}^{\frac{d}{2}}\exp \left[-\lambda (\mathbf{u}-\mathbf{V}{)}^2\right],$$

$$f({\mathbf{x}}_0,t)=\frac{1}{\tau }{\int }_{t_0}^{t}g({\mathbf{x}}^{\prime },t^{\prime })e^{-(t-t^{\prime })/\tau }d{t}^{\prime }+e^{-(t-t_0)/\tau }{f}_0({\mathbf{x}}_0-\mathbf{u}(t-t_0)),$$

$$g(\mathbf{x},t)=g_0+g_{t}t+g_{\mathbf{x}}\cdot \mathbf{x},$$

$$f_0(\mathbf{x})=f_0+f_{\mathbf{x}}\cdot \mathbf{x}.$$

$$f_{ij}(t)$$

$$=f_{ij}^{eq}(t)+f_{ij}^{fr}(t)$$

$$c_1$$

$$c_2$$

$$c_3$$

$$c_4$$

$$c_5$$

$${\mathbf{w}}_0=\int g_0\mathbf{\psi }(\mathbf{u})d\mathbf{u}=\int f_0\mathbf{\psi }(\mathbf{u})d\mathbf{u},$$

$$\int g_{\mathbf{x}}\mathbf{\psi }(\mathbf{u})d\mathbf{u}$$

$$\int g_t\mathbf{\psi }(\mathbf{u})d\mathbf{u}$$

$${\int }_0^{\Delta t}\mathbf{u}\cdot {\mathbf{n}}_{ij}{f}_{ij}(t)dt$$

$$+\mathbf{u}\cdot {\mathbf{n}}_{ij}\left(q_4{f}_0+q_5{f}_{\mathbf{x}}\cdot \mathbf{u}\right)$$

$$={\mathcal{F}}_{ij}^{eq}+{\mathcal{F}}_{ij}^{fr}$$

$$q_1$$

$$q_2$$

$$q_3$$

$$q_4$$

$$q_5$$

$$f({\mathbf{x}}_e,t_e)$$

$$=(1-e^{-t_e/\tau })g_p+e^{-t_e/\tau }{f}_p$$

$${\mathbf{x}}^{\prime }=\mathbf{u}\left(\frac{t_e{e}^{-t_e/\tau }}{1-e^{-t_e/\tau }}-\tau \right),t^{\prime }=\left(\frac{t_e}{1-e^{-t_e/\tau }}-\tau \right),$$

$$\mathcal{G}(t)=e^{-t/\tau },$$

$$t_f=\min (-\tau \ln (r_0),\Delta t)$$

$${\mathbf{x}}_p={\mathbf{x}}_p^n+\mathbf{u}{t}_f,$$

$${\mathbf{W}}_i^{fr}=\sum _{k\in P(i)}{\mathbf{\varphi }}_k$$

$${\mathbf{W}}_i^{fr}=-\sum _{j\in N(i)}{S}_{ij}{\int }_0^{\Delta t}\int \mathbf{u}\cdot {\mathbf{n}}_{ij}{f}_{ij}^{fr}(t)\mathbf{\psi }(\mathbf{u})d\mathbf{u}dt=-\sum _{j\in N(i)}{S}_{ij}\int {\mathcal{F}}_{ij}^{fr}\mathbf{\psi }(\mathbf{u})d\mathbf{u}.$$

$${\mathbf{F}}_{ij}^{eq}={\int }_0^{\Delta t}\int \mathbf{u}\cdot {\mathbf{n}}_{ij}{f}_{ij}^{eq}(t)\mathbf{\psi }(\mathbf{u})d\mathbf{u}dt=\int {\mathcal{F}}_{ij}^{eq}\mathbf{\psi }(\mathbf{u})d\mathbf{u}.$$

Full text

multiscale transport, unified gas-kinetic scheme, wave-particle formulation, multi-efficiency

I Introduction

There are mainly two kinds of numerical methods for non-equilibrium gas flow simulations, i.e., the stochastic particle method and the deterministic method. The stochastic method employs discrete particles to simulate the statistical behavior of molecular gas dynamics. Since the Lagrangian-type stochastic particle method can be regarded as the best adaptive strategy in velocity space discretization, it requires low computational memory and gains high efficiency in rarefied flow computations, especially for high-speed flows in multidimensional cases. Due to the particle implementation, the stochastic method is very robust and shows great advantages on numerical modeling for complex flow physics, such as gas mixture and chemical reaction. However, by using a finite number of simulation particles to recover the gas distribution function, it suffers from the statistical noises especially for low-speed flow with small temperature variation. Moreover, in the continuum flow regime at small Knudsen number, the intensive particle collisions make the computational cost very high by explicitly following particle’s interaction physics. On the other hand, the deterministic method employs discrete distribution function to solve the kinetic equations. It can give accurate solutions without the statistical noises and is able to achieve high efficiency with equation-based numerical acceleration techniques. In order to capture non-equilibrium distribution, the whole velocity space has to be discretized with a high resolution, which leads to huge memory consumption and computational cost for high speed rarefied flow in three dimensional cases. In addition, for both the stochastic particle and the deterministic methods, once the gas evolution process is splitted into the collisionless trasport and instant collision, a numerical dissipation being proportional to the time step will be introduced implicitly. Therefore, the mesh size and the time step in these schemes have to be less than the mean free path and the particle collision time in order to properly control the numerical error and reduce its contamination to the physical dissipative effect. Otherwise, only the Euler limit can be recovered in these schemes once the physical dissipation is overwhelmingly taken by the numerical one in the continuum regime, such as the simulation of boundary layer or the cavity flow at high Reynolds numbers. This would severely constrain the applications of these kinetic methods to the continuum flow simulations.

The direct simulation Monte Carlo (DSMC) method bird1994book is the most successful particle simulation method in the engineering applications of rarefied flow oran1998direct . In the past decades, great efforts have been made to further improve the DSMC method on the aspects of accuracy and efficiency. In order to reduce the statistical error, the information preservation (IP) method fan2001statistical ; shen2006rarefied ; sun2002direct and the low-variance deviational simulation Monte Carlo (LVDSMC) method baker2005variance ; homolle2007low have been developed for low-speed microflows. By matching the moments of macroscopic equations, the moment-guided Monte Carlo method degond2011moment is able to reduce the fluctuations as well. In order to address the stiffness problem of the collision term in the continuum regimes, several asymptotic preserving (AP) Monte Carlo methods pareschi2001time ; pareschi2000asymptotic ; ren2014asymptotic ; dimarco2011exponential have been developed so that the Euler solution can be obtained in the hydrodynamic limit without the requirement on the time step and mesh size as that in the traditional DSMC method. In order to reduce the numerical diffusion error resulting from free molecular transport process, a low diffusion particle method burt2008low has been constructed for inviscid flow simulations. For further improving the computational efficiency, the stochastic particle methods based on kinetic model equations, such as the Bhatnagar–Gross–Krook (BGK), the ellipsoidal statistical BGK (ES-BGK) models nanbu1981simulation ; gallis2000application ; macrossan2001particle ; tumuklu2016particle ; fei2018particle , and the Fokker-Planck (FP) model jenny2010solution ; gorji2011fokker ; gorji2015fokker , have been constructed to simulate monatomic and diatomic gas flows. By taking account of the collision into the transport process, the particle methods of Jenny et al. jenny2010solution ; gorji2011fokker ; gorji2015fokker and Fei et al. fei2018particle are capable for flow simulation without the restriction on the time step and the mesh size to the kinetic level.

On the other hand, many studies have been carried out on the discrete velocity methods (DVM) to solve the Boltzmann equation chu1965kinetic ; aristov2012direct ; tcheremissine2005direct ; li2004gkua ; li2019 ; xu2010unified ; guo2013discrete . In order to reduce the computational cost and increase the numerical efficiency, many acceleration methods have been applied, including implicit algorithms yang1995rarefied ; mieussens2000discrete ; mieussens2000implicit ; zhu2016implicit ; zhu2017unified ; zhu2019implicit , high-order/low-order (HOLO) methods taitano2014moment ; chacon2017multiscale , memory reduction techniques chen2017memory ; yang2018memory , fast evaluation of the full Boltzmann collision term mouhot2006fast ; wu2013deterministic , and adaptive refinement method chen2012amr . Asymptotic preserving schemes filbet2010class ; dimarco2013imex are also developed to release the stiffness of the collision term in the small Knudsen number cases. However, for most AP schemes only the Euler solution in the hydrodynamic limit is recovered. By employing the integral solution of kinetic model equation, the unified gas-kinetic scheme (UGKS) xu2010unified ; xu2015direct and discrete UGKS (DUGKS) guo2013discrete are constructed with a true multiscale transport process which couples particles’ free streaming and collision. Similar to the gas-kinetic scheme (GKS) in the hydrodynamic regime xu2001gks , the NS solutions can be preserved by the UGKS without the cell resolution to the kinetic scale. The multiscale transport modeling in UGKS is general and can be applied in many transport problems, such as aerospace application jiang , radiative transfer sun2015gray ; sun2015asymptotic , phonon transport guo2016phonon ; luo2017discrete , plasma liu2017plasma , and granular flow liu2019granular . Since the time step and the mesh size in UGKS are not limited to be smaller than the particle collision time and the mean free path, the UGKS can achieve higher efficiency in the near continuum flow regime. However, as a deterministic method based on the discrete velocity distribution function, the computational cost and memory requirement of the UGKS are very high for multi-dimensional hypersonic rarefied flow due to the massive number of discrete points in the particle velocity space.

Since the UGKS and the particle method have different, but complementary advantages and disadvantages, a particle implementation of the UGKS is preferred to design a more powerful tool for simulating multiscale transport. In this paper, we will present in details the construction of the unified gas-kinetic wave-particle (UGKWP) method on unstructured mesh. The direct particle implementation of the multiscale UGKS is the unified gas-kinetic particle (UGKP) method li2018ugkp , and the UGKWP method is a further improvement by incorporating waves into UGKP in order to reduce the computational cost and memory requirement liu2018ugkwp . In the UGKP method, the multiscale transport process in UGKS is recovered in the particles’ evolution process, where the collision effect is taken into account so that this particle method can present accurate solutions in different flow regimes without the DSMC requirement for the mesh size and time step. With the particle implementation, the capability of the UGKS to solve high-speed rarefied flow problems has been much enhanced. Since the particles’ interaction makes the distribution function to approach to the equilibrium state, some of the simulation particles can be replaced by the evolution of analytical distribution function after some collisions. Therefore, a wave-particle formulation is introduced in the construction of the UGKWP method, where these equilibrium particles will be expressed and computed in a deterministic way instead of by discrete simulation particles. As a result, for the near continuum flow, only a few particles are required and most computation can be handled analytically, so the computational efficiency could be greatly increased and the statistical noises from discrete particles could be efficiently reduced. The UGKWP method can adaptively become the particle simulation method in highly rarefied flow regimes and the hydrodynamic flow solver, same as the gas-kinetic scheme xu2001gks , in the continuum regimes. It should be pointed out that different from the hybrid methods wadsworth1990one ; sun2004hybrid ; degond2010multiscale which are based on the domain decomposition and solver hybridization, the UGKWP describes the physical state by an adaptive wave-particle decomposition in each cell with a unified treatment in the whole computational domain. Specifically, in the UGKWP method the physical state in a finite volume cell is separated into the hydrodynamic waves and discrete particles, which are expressed and transported through an analytic distribution for waves and the particles for the non-equilibrium part, respectively. According to the numerical resolution and local flow physics, the evolutions of the hydrodynamic waves and discrete particles are coupled dynamically. Therefore, with the multiscale transport process and the wave-particle adaptive property, the UGKWP method not only recovers the multiscale nature of UGKS, but also achieves multi-efficiency in all Knudsen number regimes. More specifically, the UGKWP method performs as a stochastic particle method in highly rarefied flows and becomes a hydrodynamic fluid solver for viscous flow simulations in the continuum limit. In the current paper, the UGKWP method is constructed and validated on the unstructured mesh. Numerical test cases, including the Sod shock tube, cavity flow, laminar boundary layer, and high-speed cylinder flow, will be computed across different flow regimes. For the hypersonic rarefied flow, such as Mach number $30$, the UGKWP has several-order-of-magnitude computational cost and memory reductions in comparison with UGKS.

The remainder of the paper is organized as follows. Section II presents the construction of the UGKWP method on the unstructured mesh. In Section III, numerical test cases are carried out to validate the present UGKWP method. Discussion and conclusion will be given in the last section.

II Numerical Method

In this section, we will introduce the unified gas-kinetic wave-particle (UGKWP) methods in details. Since the unified gas-kinetic particle (UGKP) method is a direct particle implementation of the unified gas-kinetic scheme (UGKS), and the UGKWP method is an enhanced UGKP method by employing the adaptive wave-particle decomposition, the UGKP method and UGKS will be introduced first before presenting the UGKWP method.

II.1 Unified gas-kinetic scheme (UGKS)

In the framework of finite volume method, the UGKS considers the conservations in the discretized space and time for both the macroscopic flow variables and microscopic gas distribution function. Specifically, for a discrete finite volume cell $i$ and discretized time step $\Delta t=t^{n+1}-t^{n}$, the governing equations are

[TABLE]

and

[TABLE]

where $\boldsymbol{w}$ is the conservative flow variables, i.e., the densities of mass, momentum and energy $(\rho,\rho\boldsymbol{V},\rho E)$, and $f$ is the gas distribution function. $N(i)$ denotes the set of the interface-adjacent neighboring cells of cell $i$, and cell $j$ is one of the neighbors. The interface between cells $i$ and $j$ is represented by the subscript $ij$. Hence, $S_{ij}$ and $\boldsymbol{n}_{ij}$ are referred to as the area and normal vector of the interface $ij$, $\boldsymbol{F}_{ij}$ denotes the macroscopic fluxes across the interface, and $f_{ij}(t)$ is the time-dependent distribution function on the interface. In addition, $\Omega_{i}$ is the volume of cell $i$, $\boldsymbol{u}$ denotes the microscopic velocity, and $\mathcal{J}(f,f)$ is the collision term. The connection between the macroscopic and microscopic governing equations (1) and (2) is that all the macroscopic variables can be obtained from the moments of distribution function, such as conservative variables and macroscopic fluxes

[TABLE]

and

[TABLE]

where $\boldsymbol{\psi}(\boldsymbol{u})=(1,\boldsymbol{u},\frac{1}{2}|\boldsymbol{u}|^{2})$. The collision term satisfies the compatibility condition

[TABLE]

In this paper, the BGK relaxation model bhatnagar1954model is considered for the collision term, i.e.,

[TABLE]

where the relaxation time or the mean collision time $\tau$ is computed by the ratio of dynamic viscosity to pressure, i.e., $\tau=\mu/p$. The equilibrium state $g$ is the Maxwellian distribution

[TABLE]

where $d$ is degree of freedoms, and $\lambda$ is related to the temperature $T$ by $\lambda=m_{0}/2k_{B}T$. Here $m_{0}$ and $k_{B}$ are the molecular mass and the Boltzmann constant.

It should be noted that Eq. (1) and Eq. (2) are the fundamental physical laws on the scale of mesh size and time step, which describe the general conservations of the macroscopic flow variables and microscopic gas distribution function. For finite volume method, the evolution of flow physics relies mainly on the construction of the flux function at the cell interfaces. In UGKS, the time-dependent distribution function $f_{ij}(t)$ for the flux function is constructed from the integral solution of kinetic model equation,

[TABLE]

where $\boldsymbol{x}_{0}$ is the point for the evaluation of local gas distribution function evolution, such as the center of a cell interface for flux evaluation. $f_{0}(\boldsymbol{x})$ is the initial distribution function around $\boldsymbol{x}_{0}$ at the beginning of each step $t_{0}=0$, and $g(\boldsymbol{x},t)$ is the equilibrium state distributed around $\boldsymbol{x}_{0}$ and $t_{0}$.

Specifically, for second order accuracy we have

[TABLE]

and

[TABLE]

The time dependent distribution function at the cell interface can be constructed

[TABLE]

with the coefficients

[TABLE]

where $f_{ij}^{eq}(t)$ and $f_{ij}^{fr}(t)$ are the terms related to the equilibrium state $g(\boldsymbol{x},t)$ and initial distribution function $f_{0}(\boldsymbol{x})$, respectively. The initial gas distribution function $f_{0}(\boldsymbol{x})$ in Eq. (10) is obtained from the spatial reconstruction of the distribution function at $t^{n}$. The equilibrium state $g_{0}$ is computed from the compatibility condition

[TABLE]

and the spatial and temporal derivatives of the equilibriums state can be obtained by

[TABLE]

where $\boldsymbol{w}_{\boldsymbol{x}}$ is the spatial derivatives of the conservative variables obtained from reconstruction.

Eq. (8) and Eq. (11) give a transition from the initial distribution function to the equilibrium state with the increment of time, which couples the particles’ free transport and collision processes. With an accumulating effect of the particle collision, the modeling scale changes from the kinetic scale to hydrodynamic scale. For a discretized space and time, the integral solution will adapt the physical solution on the scale of mesh size and time step according to the relation between the numerical resolution and flow physics, such as the ratio of the time step to the mean collision time $\Delta t/\tau$. Specifically, the integrated flux over a time step gives

[TABLE]

with the coefficients

[TABLE]

Here $\mathcal{F}_{ij}^{eq}$ and $\mathcal{F}_{ij}^{fr}$ are referred to as the equilibrium flux and the free transport flux, respectively. We can see that when $\Delta t\gg\tau$, only the terms $\mathcal{F}_{ij}^{eq}$ with $q_{1}\approx\Delta t$ and $q_{3}\approx\Delta t^{2}/2$ are remained for equilibrium wave interaction; when $\Delta t\ll\tau$, only $\mathcal{F}_{ij}^{fr}$ with $q_{4}\approx\Delta t$ and $q_{5}\approx-\Delta t^{2}/2$ are left for non-equilibrium particle free transport. With adaptive variation of $\Delta t/\tau$ in different regions, UGKS is able to provide multiscale flow solutions on the numerical scales. In comparison with the kinetic method based on particles’ free transport mechanics, the mesh size and time step of the UGKS are not constrained by the particles’ mean free path and mean collision time. With the scale-adaptive flux function, the UGKS is an efficient deterministic method for multiscale flow simulations in all regimes.

The algorithms for one time step evolution of the UGKS from $t^{n}$ to $t^{n+1}$ can be summarized as follows.

Step 1

Reconstruct the microscopic gas distribution function $f^{n}$, and obtain the initial gas distribution function $f_{0}(\boldsymbol{x})$ in Eq. (10).

Step 2

Obtain the equilibrium state $g_{0}$ at cell interface from the initial distribution function $f_{0}$ by the compatibility condition (13).

Step 3

Reconstruct the macroscopic flow variables $\boldsymbol{w}^{n}$ and obtain the spatial and temporal derivatives of the equilibrium state $g_{\boldsymbol{x}}$, $g_{t}$ from Eq. (14) with the reconstructed $\boldsymbol{w}_{\boldsymbol{x}}$.

Step 4

Compute the microscopic and macroscopic fluxes across cell interfaces by Eq. (15) and Eq. (4).

Step 5

Update the conservative flow variables $\boldsymbol{w}^{n+1}$ and the microscopic gas distribution function $f^{n+1}$ by Eq. (1) and Eq. (2).

Detailed implementation and analysis of the UGKS can be found in Refs. xu2010unified, ; huang2012unified, ; xu2015direct, .

II.2 Unified gas-kinetic particle (UGKP) method

Since the particles’ tracking and interaction can be regarded as an optimal strategy for the grid point adaption in the particle velocity space, the stochastic particle methods obtain very high efficiency for simulation of high-speed rarefied flows in three dimensional case. Therefore, in this section the particle implementation of the UGKS with multiscale transport process will be carried out to construct the UGKP method.

The physical picture for particles’ evolution in a time step $\Delta t$ is illustrated in Fig. 1. It shows that one particle will keep free transport until it encounters other particle and gets collided, and then it will continue its moving and colliding process. Before its first collision, the particle’s trajectory and the characteristic line are the same. The particle retains its initial discrete velocity. Once collision happens, the particle velocity changes and we cannot get the exact location and velocity of the particle unless it is traced step by step for each individual collision. It should be emphasized that Fig. 1 describes the physics on the numerical scale instead of the kinetic scales of mean free path and collision time. The free transport time $t_{f}$ changes with the local physics, and multiple or a huge number of collisions are allowed within the time interval $(t_{f},\Delta t)$, which specifies different flow regimes. Similar to the UGKS, direct modeling of the flow physics on the time step $\Delta t$ scale is the key to construct a multiscale scheme.

Here we re-write the integral solution (11) along the characteristic line for the end point $(\boldsymbol{x}_{e},t_{e})$ as

[TABLE]

where

[TABLE]

and $t_{e}$ is related to the time step for a numerical scheme. The point $(\boldsymbol{x}^{\prime},t^{\prime})$ locates on the characteristic line moving from the mid point to the end point as the increasing of the ratio $t_{e}/\tau$. Eq. (17) describes that the discrete distribution function at time $t_{e}$ is a combination of the initial distribution function $f_{p}$ and the Taylor expansion of the equilibrium state $g_{p}$. In analogy to the discrete distribution function in the deterministic methods, it is straightforward to evolve the particle in one time step through Eq. (17) by changing the mass-weight of particle instead of its microscopic velocity. However, this treatment will lose one of the most important advantages of particle methods, i.e., the adaptive property in the velocity space. In Eq. (17), the probability of maintaining the initial distribution function through particles’ free transport is given by $e^{-t/\tau}$. Statistically, the free transport time for each particle can be determined. Therefore, the motion of all the particles before their first collision can be accurately tracked. Although the subsequent collision and motion are not exactly followed, in a systematic point of view, all the particles encountering collision in a local region will approach to an equilibrium distribution $g_{p}$ according to the kinetic model.

In the framework of finite volume method, we will construct the UGKP method by directly modeling the above physical processes on the scale of mesh size and time step xu2015direct . Specifically, the free transport process of each particle before its first collision within a time step will be accurately tracked and the effect of collision is to annihilate the particles and be recovered by re-sampling them from a specific Maxwellian equilibrium state. In order to give a clear description, we will use the terms of free transport process and collision process as shown in Fig. 1 to denote the stages before and after the particles’ first collision in one time step, while the whole evolution process in one time step is denoted as the transport process or the multiscale transport process. As a finite volume method, the computation of the UGKP method for a time step evolution would contain

Macro level

Compute the numerical fluxes to update the conservative flow variables, which includes

(Pa)

computing the fluxes contributed from the particles’ motion in the free transport process;

(Pb)

computing the fluxes contributed from the particles’ motion in the collision process.

Micro level

Evolve the gas distribution function, i.e., update the simulation particles, which includes

(Pc)

tracking all the particles’ motion in the free transport process;

(Pd)

re-sampling collisional particles in the collision process.

II.2.1 Free transport process

Eq. (17) gives the cumulative distribution for particles’ free transport time

[TABLE]

so the free transport time of a particle within a time step $\Delta t$ can be determined by

[TABLE]

where $r_{0}$ is a random number generated from a uniform distribution between $(0,1)$. Given with the free transport time $t_{f}$, the particle can be accurately tracked by

[TABLE]

where the micro velocity $\boldsymbol{u}$ of the particle remains unchanged during the free transport process. Different from the DSMC method with $t_{f}=\Delta t$, where all the particles free stream with a whole time step, the free transport time in the UGKP method is constrained due to particles’ collision.

During the free transport process, the contribution to the numerical fluxes of cell $i$ can be obtained by counting the particles across the cell interfaces

[TABLE]

where $P(i)$ is the set of the particles moving across the interfaces of the cell $i$ during the free transport process. The vector $\boldsymbol{\phi}_{k}=(m_{p},m_{p}\boldsymbol{u}_{k},\frac{1}{2}m_{p}|\boldsymbol{u}_{k}^{2}|)$ denotes the mass, momentum,and energy carried by the particle $k$. In comparison with the multiscale transport process given in the UGKS in Eq. (15), the free transport process (21) only recovers the fluxes contributed by the initial distribution function $f_{0}(\boldsymbol{x})$, i.e., $\mathcal{F}_{ij}^{fr}$ with the terms $q_{4}$ and $q_{5}$. For comparison, the counterpart of $\boldsymbol{W}_{i}^{fr}$ in the deterministic UGKS would be

[TABLE]

So far, we have carried out the processes (Pa) and (Pc) by Eq. (22) and Eq. (21).

II.2.2 Collision process: macroscopic fluxes

In the collision process, the particles get collided and keep on moving and colliding process. During this process, once the particles move across the cell interfaces, they will contribute to the macroscopic fluxes as well. However, since we are not developing a full particle tracking method, the motion of simulation particles in the collision process will not be followed explicitly. So the macroscopic fluxes cannot be directly obtained from the discrete particles as in the free transport process. Fortunately, these fluxes have been already given in the UGKS in Eq. (15), i.e., $\mathcal{F}_{ij}^{eq}$ with the terms $q_{1}$, $q_{2}$ and $q_{3}$.

Hence, the macroscopic fluxes of the collision process across the cell interface $ij$ can be computed from the reconstructed macroscopic flow variables by

[TABLE]

In UGKS, $g_{0}$ is obtained from the reconstructed initial distribution function $f_{0}$ by the compatibility condition (13). However, in the particle method there is no explicit gas distribution function, so the equilibrium state $g_{0}$ at cell interface is computed from the reconstructed macroscopic flow variables, i.e.,

[TABLE]

where $g_{l}$ and $g_{r}$ are the equilibrium state on the left and right sides of cell interface, which are determined by the interpolated macroscopic flow variables $\boldsymbol{w}_{l}$ and $\boldsymbol{w}_{r}$. Same as that in the GKS and UGKS, $g_{t}$ and $g_{\boldsymbol{x}}$ can be obtained by Eq. (14), and then the equilibrium fluxes $\boldsymbol{F}_{ij}^{eq}$ can be analytically computed.

At this moment, the multiscale fluxes (15) in the UGKS have been fully recovered by free transport fluxes (22) and collisional fluxes (24). The macroscopic flow variables can be updated by the conservation laws

[TABLE]

So far, we have dealt with the process (Pb) and finished the update on the macroscopic level.

II.2.3 Collision process: microscopic particles

In the free transport process, the detailed motion of all the particles in the time interval $(0,t_{f})$ have been tracked. For the collisionless particles with $t_{f}=\Delta t$, the update of the microscopic state including particles’ velocity and location has been finished for the current step. While for the collisional particles with $t_{f}<\Delta t$, each of them will suffer at least one collision in the time interval $(t_{f},\Delta t)$, and the collective effect of the collisions is to force all the collisional particles in the local region to follow a specific equilibrium distribution $g_{p}$. According to the conservation, from the updated macroscopic flow variables and the streamed collisionless particles, we can easily obtain the the collisional particles at the end of the time step

[TABLE]

where $\boldsymbol{w}_{i}^{p}$ is the conservative flow variables carried by the collisionless particles in the currently investigated cell $i$ after their solely free transport process. Therefore, in the collision process $(t_{f},\Delta t)$, these collisional particles will be deleted first due to particles’ collision, and then re-sampled from their corresponding macroscopic flow variable $\boldsymbol{w}_{i}^{h}$ at each end of time step. A unique Maxwellian distribution (7) can be determined from $\boldsymbol{w}_{i}^{h}$. According to the macroscopic velocity and temperature, the collisional particles can be re-sampled in the cell $i$ to recover the distribution function. So far, we have carried out the process (Pd), and both the macroscopic flow variables and the microscopic particles have been updated.

II.2.4 Summary and discussions

In order to give a visual demonstration for the evolution of simulation particles, a series of diagrams are drawn in Fig. 2 to illustrate the composition of particles on different evolution stages within one time step. The explanation of Fig. 2 and the summary of the UGKP method will be given as follows.

Step 1

Give an initial state with macroscopic flow variables and microscopic particles as shown in Fig. 2(a), where the simulation particles could be an initial equilibrium distribution for the start of flow simulations, or a non-equilibrium distribution evolved from the previous step in the time evolution process.

Step 2

Free transport process, which includes

• •

Obtain the free transport time $t_{f}$ for each particle by Eq. (20) and classify the particles into collisionless particles (white circles in Fig. 2(b)) and collisional particles (solid circles in Fig. 2(b));

• •

Move all the particles for a free transport time $t_{f}$ in Eq. (21);

• •

Cumulate the free transport fluxes $\boldsymbol{W}_{i}^{fr}$ in Eq. (22) by counting the particles which move across the cell interfaces;

• •

Tally the collisionless particles with $t_{f}=\Delta t$ after streaming all the particles and calculate $\boldsymbol{w}_{i}^{p}$ of these freely transported particles in each cell $i$ (denoted by the particles on the bottom part of Fig. 2(c)), and delete the collisional particles.

Step 3

Collision process, which includes

• •

Reconstruct macroscopic flow variables $\boldsymbol{w}^{n}$ to obtain the conservative flow variables $\boldsymbol{w}_{l}$ and $\boldsymbol{w}_{r}$ on the left and right sides of cell interface;

• •

Obtain the equilibrium state $g_{0}$ at cell interface from Eq. (25), and compute the derivatives $g_{\boldsymbol{x}}$ and $g_{t}$ from the reconstructed spatial derivative of macroscopic flow variables $\boldsymbol{w}_{\boldsymbol{x}}$ by Eq. (14);

• •

Compute the collisional fluxes $\boldsymbol{F}_{ij}^{eq}$ in Eq. (24), i.e., the terms with $q_{1}$, $q_{2}$ and $q_{3}$ in Eq. (15).

• •

Update the conservative variables $\boldsymbol{w}_{i}^{n+1}$ by Eq. (26), and obtain the macroscopic variables $\boldsymbol{w}_{i}^{h}$ by Eq. (27) for the updated collisional particles (the grey area shown in Fig. 2(c)).

• •

Re-sample the collisional particles from $\boldsymbol{w}_{i}^{h}$ to finish the update process of the microscopic particles as illustrated in Fig. 2(d).

Step 4

Determine the computation to the next time step.

• •

If the solution is convergent for steady flows or the pre-described evolution time is achieved for unsteady flows, stop the program;

• •

Otherwise, go to Step 1 and continue the computation, where the updated state in Fig. 2(d) could be an initial state in Fig. 2(a) for next step evolution.

The UGKP method is a conservative finite volume method, where the simulation particles are employed to recover the underlying non-equilibrium distribution function. On the macroscopic level, the conservative variables are updated with the fluxes by conservation laws. On the microscopic level, all the particles are accurately tracked in the free transport process and the collisional particles are re-sampled from the updated equilibrium state in the collisional process. The maintenance of conservation laws is the key to the success of the current particle method. In addition, it should be noted that in the free transport process each particle moves over a free transport time $t_{f}$ instead of a whole time step $\Delta t$, and together with the subsequent collisional process a mult-scale transport process is constructed with recovering the multiscale nature of UGKS.

II.3 Unified gas-kinetic wave-particle (UGKWP) method

In this section, the concept of the wave-particle will be introduced into the UGKP method for the further development of an efficient UGKWP method. In the UGKP method, the gas distribution function is fully represented by the simulation particles, and the collisional particles are deleted and re-sampled from a Maxwellian distribution in the collisional process. Theoretically, this portion of gas distribution function can be expressed in an analytic way. As shown in Fig. 2(c), the gas distribution function can be recovered by hydrodynamic waves and discrete kinetic particles, which correspond to the macroscopic variables $\boldsymbol{w}^{h}$ and $\boldsymbol{w}^{p}$. For the next time step evolution in the UGKP method, the re-sampled equilibrium particles will be re-classified into collisionless and the collisional particles according to the free transport time $t_{f}$ again. In the free transport process, both types of the particles will contribute to the free transport fluxes, but only the collisionless part particles will be remained as particles at the end of the time step to recover the non-equilibrium gas distribution function. Therefore, only the collisionless particles in the hydrodynamic waves should be re-sampled at the end of each time step and the contribution from these collisional particles, which are generated from the hydrodynamic wave previously, to the free transport fluxes in the next time step can be computed analytically.

From the cumulative distribution (19), we can easily obtain the expectation of the proportion of the collisionless particles in each cell, and the particles required to be sampled in the hydrodynamic waves $\boldsymbol{w}^{h}$ at the end of time step are

[TABLE]

The free transport fluxes contributed from the collisional particles of $(\boldsymbol{w}^{h}-\boldsymbol{w}^{hp})$ can be computed analytically by

[TABLE]

where $g_{0}^{h}$ is the Maxwellian distribution determined by $\boldsymbol{w}^{h}$ and $g_{\boldsymbol{x}}^{h}$ is the spatial derivative of the Maxwellian distribution, which can be obtained from the reconstruction of $\boldsymbol{w}^{h}$.

Therefore, in the UGKWP method the update process for the conservative variables will be

[TABLE]

In comparison with the update formula (26) in the UGKP method, the term of $\boldsymbol{F}_{ij}^{fr,wave}$ is the analytic part extracted from the particles’ free transport fluxes $\boldsymbol{W}_{i}^{fr}$. The combination of the last two terms in Eq.(30) for UGKWP is the same as the last term in Eq.(26) for UGKP.

A series of figures are drawn in Fig. 3 to illustrate the evolution of the UGKWP method. The algorithm of the UGKWP method can be summarized as follows.

Step 1

Get an initial state with the macroscopic flow variables $\boldsymbol{w}^{n}$ and the microscopic particles. The particles include the particles $\boldsymbol{w}^{p}$ evolved from the previous step and the collisionless particles sampled from the updated hydrodynamics waves $\boldsymbol{w}^{hp}$ as show in Fig. 3(d). For the first step, $\boldsymbol{w}^{p}=0$ and $\boldsymbol{w}^{h}=\boldsymbol{w}^{n=0}$ as shown in Fig. 3(a).

Step 2

Free transport process, which includes

• •

Classify the particles $\boldsymbol{w}^{p}$ into collisionless particles (white circles in Fig. 3(b)) and collisional particles (solid circles in Fig. 3(b)) according to the free transport time $t_{f}$ determined by Eq. (20).The free transport time of the re-sampled collisionless particles $\boldsymbol{w}^{hp}$ is always $t_{f}=\Delta t$ (white circles on the right top of Fig. 3(b));

• •

Stream all the particles over the free transport time $t_{f}$;

• •

Cumulate the free transport fluxes $\boldsymbol{W}_{i}^{fr}$ in Eq. (22) by counting the particles which move across cell interfaces;

• •

Tally the updated collisionless particles $\boldsymbol{w}^{p}$ and delete the streamed collisional particles;

• •

Compute the free transport fluxes $\boldsymbol{F}^{fr,wave}$ contributed from the un-sampled particles $(\boldsymbol{w}^{h}-\boldsymbol{w}^{hp})$ by Eq. (29).

Step 3

Collision process, which includes

• •

Compute the collisional fluxes $\boldsymbol{F}_{ij}^{eq}$ as same as in the UGKP method;

• •

Update the conservative variables $\boldsymbol{w}_{i}^{n+1}$ by the conservation laws (30), and obtain the updated hydrodynamic waves $\boldsymbol{w}_{i}^{h}$ by Eq. (27), as shown in Fig. 3(c);

• •

Sample the collisionless particles $\boldsymbol{w}_{i}^{hp}$ in Eq. (28) for next step evolution and finish the update process of the microscopic particles, see in Fig. 3(d).

Step 4

Determine the computation to the next time step.

• •

If the solution is convergent for steady flows or the pre-described evolution time is achieved for unsteady flows, stop the program;

• •

Otherwise, go to Step 1 and continue the computations.

In the UGKWP method, the wave-particle formulation is introduced to represent the non-equilibrium gas distribution function. The difference between the UGKP and UGKWP methods is that only the collisionless particles in the hydrodynamic waves are re-sampled at the end of time step; and for the next step evolution the free transport fluxes contributed from these un-sampled particles in the hydrodynamic waves are computed in a deterministic way. For near continuum flows where intense collisions are involved, i.e., $t_{f}\ll\Delta t$, the hydrodynamic waves will be dominant and only a few collionless particles are required to be sampled, which makes the current method very efficient. Therefore, in different flow regimes, the wave-particle decomposition will give an optimal formulation for the non-equilibrium gas distribution function, and achieve higher efficiency both in the continuum and rarefied flow regimes.

II.4 Unstructured mesh technique

II.4.1 Particle tracking

During the free transport process, trajectories of the simulation particles are fully tracked. For an arbitrary particle locating at point $O$ with microscopic velocity $\boldsymbol{u}$ in the triangular cell $\triangle{\rm ABC}$, its displacement in the free transport process would be $\boldsymbol{x}_{\rm OP}=\boldsymbol{u}t_{f}$. The displacement $\boldsymbol{x}_{\rm OP}$ may intersect with the faces satisfying $\boldsymbol{x}_{\rm OP}\cdot\boldsymbol{n}>0$, where $\boldsymbol{n}$ is the normal vector of cell interface. The intersection point, e.g., point $E$ on face ${\rm AC}$ satisfies

[TABLE]

where point $D$ is the centroid of face $\rm AC$. Similarly, the intersection point $F$ on face $\rm BC$ can be obtained as well. A minimum value

[TABLE]

can be used to determine the first intersection point of the trajectory and the cell interfaces haselbacher2007efficient . If $\alpha>1$, the particle is still inside the current cell, and the updated location will be $\boldsymbol{x}_{P}=\boldsymbol{x}_{O}+\boldsymbol{u}t_{f}$. If $\alpha\leq 1$, the particle will move out of the current cell. For this case, we will first move the particle to the intersection point $\boldsymbol{x}_{E}=\boldsymbol{x}_{O}+\boldsymbol{u}\alpha t_{f}$, and then track the particle in its neighboring cell using the same method for the remaining free transport time $(1-\alpha)t_{f}$.

II.4.2 Particle sampling

In the collision process, simulation particles will be re-sampled from a given Maxwellian distribution function to recover the gas distribution function on microscopic level. Given with a Maxwellian distribution determined by $(\rho_{s},\boldsymbol{V}_{s},\lambda_{s})$, the microscopic velocity for each particle to sample can be obtained from

[TABLE]

where $U_{s}$, $V_{s}$ and $W_{s}$ are the components of $\boldsymbol{V}_{s}$, $r_{1}$, $r_{2}$, $r_{3}$ and $r_{4}$ are independent random numbers generated from the uniform distribution between the interval $(0,1)$. In the current study, a symmetric sampling process is employed. Specifically, from a group of $r_{1}$, $r_{2}$, $r_{3}$ and $r_{4}$, a pair of simulation particles with microscopic velocity $(u,v,w)$ and $(u^{\prime},v^{\prime},w^{\prime})$ are sampled, where the symmetric microscopic velocity is

[TABLE]

Given with a pre-set reference mass $m_{r}$ for the simulation particle, the number of particles to be sampled is determined by

[TABLE]

where $\Omega_{s}$ is the cell volume, and $m_{min}$ is the minimum mass to sample. In the sampling process, for the cases $N_{s}>0$ the mass weight actually sampled for each simulation particle is

[TABLE]

which guarantees the mass density $\rho_{s}$ in the volume $\Omega_{s}$ after sampling process.

Another way is to give a preference number $N_{r}$ for each cell, then the reference mass $m_{r}$ can be approximated by

[TABLE]

and then the number of particles to sample $N_{s}$ can be obtained by Eq. (35). By this way, the number of particles in each cell can be basically controlled around the given reference number $N_{r}$.

Besides the mass weight and the microscopic velocity, the location of each simulation particle is required as well in the sampling process. For an arbitrary triangular cell $\triangle{\rm ABC}$ as shown in Fig.5, a point inside can be denoted by $(\xi,\eta)$, which has ${\rm BF}=\xi{\rm BC}$, ${\rm AE}=\eta{\rm AB}$ and ${\rm EG}\parallel{\rm BC}$. The coordinates of the point $(\xi,\eta)$ in the global system are

[TABLE]

Assuming that the density is constant inside the cell, i.e., $\rho(\boldsymbol{x})=\rho_{s}$, the normalized cumulative distribution function up to the line $\xi=\xi_{0}$ is

[TABLE]

and along the line $\xi=\xi_{0}$ the cumulative distribution function up to the point $(\xi_{0},\eta_{0})$ is

[TABLE]

Therefore, generating two random numbers $r_{1}$ and $r_{2}$ from a standard uniform distribution, the particle location can be determined by Eq.(38) with $\xi=r_{1}$ and $\eta=\sqrt{r_{2}}$.

In the current study, although piecewise constant of density is assumed in a finite volume cell during the particle sampling process, we find the spatial accuracy is almost not reduced. This is due to the fact that only the portion $e^{-\Delta t/\tau}$ of the hydrodynamic waves are sampled, and the remaining part of hydrodynamic waves is computed analytically with second order accuracy. In the continuum regimes, although the hydrodynamic waves are dominant, the collisionless particles in the hydrodynamic waves are very few due to small value of $e^{-\Delta t/\tau}$; while in the rarefied flows, the hydrodynamic waves only take a small portion of the physical state due to mild collisions, so the particles in the hydrodynamic waves required to be re-sampled are not many as well. If the spatial derivative of density is considered in $\rho(\boldsymbol{x})$, the cumulative distribution can be derived as well, which would be more complicated than Eq. (39) and Eq. (40). Acceptance-rejection strategy bird1994book ; shen2006rarefied can be applied to handle the location sampling.

III Numerical examples

In this section, the UGKWP method will be tested in a wide range of multiscale flow problems. The performance of the method will be evaluated quantitatively.

III.1 Sod test case

On a two dimensional triangular mesh, the Sod shock tube problem has been computed at different Knudsen numbers to valid the current UGKWP method in the continuum and rarefied flows.

The initial condition is

[TABLE]

As shown in Fig. 6, the spatial discretization is carried out by an unstructured mesh with $100\times 5\times 2$ triangular cells. For the UGKWP computation, the pre-set reference mass $m_{r}$ for a simulation particle is $10^{-7}$; while for the UGKS simulation, $100\times 100$ velocity points are used to discretize the velocity space. The top and bottom boundaries are treated as symmetric planes. The results at the time $t=0.12$ for the cases at ${\rm Kn}=10^{-4}$, $10^{-3}$, $10^{-2}$, $0.1$, $1$ and $10$ in all flow regimes are presented.

The density, velocity and temperature obtained by the UGKS and the UGKWP method for different Knudsen number cases are plotted in Figs. 7–12, where the two dimensional flow field is projected to one dimensional data in the $x$ direction by taking average over the ten triangular cells along $y$ direction. In addition, in order to reduce the statistical noises, the unsteady flow solutions of the UGKWP method are averaged over $10$ times of computations. It can be seen that for all the cases in different flow regimes the UGKWP solutions agree well with the UGKS data. The capability of the UGKWP method for numerical simulations in continuum and rarefied flows are validated.

For the UGKS, once the discretization for the physical space and the velocity space is given, the computational costs for all Knudsen number cases will be same due to its unified treatment. The memory requirement and computational time in the UGKS simulations are $1.1$ GB and $15$ minutes. While for the UGKWP method, the overall CPU time of $10$ times of computations is about $65$ seconds for the cases with larger Knudsen numbers, and the memory cost is around $55$ MB. Moreover, for the case at ${\rm Kn}=10^{-4}$ in continuum flow, since the portion of hydrodynamic waves increases and much fewer discrete particles are needed to be sampled and tracked, the computational cost of the UGKWP method gets lower to $12$ seconds with $11$ MB. Generally speaking, in comparison with the UGKS, order-of-magnitude in efficiency increment and memory reduction can be achieved by the UGKWP method for the two dimensional Sod shock tube problem in the continuum and rarefied flows.

III.2 Cavity flow

The low-speed micro flow in a lid-driven cavity is computed at Knudsen numbers $0.1$, $1$ and $10$. The Knudsen number is defined as the ratio of molecular mean free path to the length of side wall. The argon gas with molecular mass $m_{0}=6.63\times 10^{-26}{\rm kg}$ is studied and the variable hard sphere model (VHS) is used for all three cases. The lid velocity is set to $50{\rm m/s}$. Isothermal boundary condition is applied with a fixed temperature $T_{w}=273{\rm K}$. The dynamic viscosity is computed by $\mu=\mu_{0}(T/T_{0})^{0.81}$.

The computational domain is discretized into $21\times 21\times 2$ triangular cells as shown in Fig. 16(a). For the UGKS computations, $100\times 100$ discrete velocity points are employed in the velocity space; and for the UGKWP method, we initially set the reference number of particles $N_{r}$ for each cell as $5000$. The numerical results are plotted in the Fig. 13 – Fig. 15, where the distributions of the velocity and temperature are compared between the UGKWP and UGKS solutions. Moreover, the velocity profiles along the central lines of the cavity are extracted by taking average of two neighboring triangular cells. From Fig. 16, it can be seen that satisfactory results are obtained for these three cases. For the low speed rarefied flow, we employ a large number of simulation particles and do many averaging process to reduce the statistical noises so that the high-order quantity, such as the temperature distribution, can be obtained. It takes about $5$ hours for the UGKWP method to obtain the current results. For the UGKS, with the acceleration techniques, such as implicit algorithm and multigrid method zhu2016implicit ; zhu2017unified , the convergent solution with no statistical noises can be obtained within $5$ minutes. Therefore, the deterministic method with acceleration techniques would still be a better choice for low-speed rarefied flow studies, which has much higher efficiency than the stochastic related method.

III.3 Laminar boundary layer

The laminar boundary layer over a flat plate is computed to validate the current multiscale method for viscous NS solutions in the continuum limit. The computational domain is $[-50,120]\times[0,50]$ as shown in Fig. 17. A non-uniform mesh with $120\times 50$ cells is employed. The free stream is monatomic gas flow at Reynolds number ${\rm Re}=10^{5}$ and Mach number ${\rm Ma}=0.3$ with constant viscosity. The Reynolds number and Mach number is defined with respect to the length of the flat plate $L=120L_{0}$. The reference variables $U_{0}$ and $t_{0}$ to non-dimensionalize the velocity and time are given by $U_{0}=\sqrt{2k_{B}T_{0}/m_{0}}$ and $t_{0}={L_{0}}/{U_{0}}$, where $T_{0}$ is the temperature in the free stream. The flow field at time $t=1000t_{0}$ is given as the convergent steady state solution in Fig. 18, where the distribution of the density and velocity around the leading edge is enlarged in the $y$ direction. Comparison between the UGKWP results and the Blasius solutions is given in the Fig. 19.

In the computation, the time step $\Delta t$ and particle collision time $\tau$ are $0.02t_{0}$ and $6.57\times 10^{-4}t_{0}$, respectively. Since the ratio $e^{-\Delta t/\tau}$ has a very small value of $6\times 10^{-14}$, the hydrodynamic wave is dominant and the particle contribution can be neglected. The computational time for $50000$ step simulations is $15$ minutes and the memory cost is $24$ MB. Under such condition, the present UGKWP method automatically becomes a hydrodynamic fluid solver, such as the GKS xu2001gks ; xu1994numerical . Due to the multiscale transport, the UGKWP method can recover NS solutions without the requirement of the mesh size and the time step being less than the mean free path and the particle collision time. Moreover, the computational cost is comparable to the hydrodynamic fluid solver in the continuum regime.

III.4 Cylinder flow

Hypersonic flow past a circle cylinder at ${\rm Ma}=5$ and $30$ are simulated to show the capability of the current method for high-speed rarefied flow simulations. The free stream is initialized with the monatomic gas flow of argon with an initial temperature $T_{\infty}=273{\rm K}$. The diameter $D$ of the cylinder is $1{\rm m}$ long. The solid boundaries are isothermal walls with a constant temperature $T_{w}=273{\rm K}$. The Knudsen number is defined with respect to the diameter of the circle cylinder.

For the free stream with a relatively low Mach number ${\rm Ma}=5$, the cases at the Knudsen numbers $0.1$ and $1$ are computed. The computational domain is discretized by $50\times 50\times 2$ triangular cells as shown in Fig. 20, which covers a region of $\pi(15D)^{2}$. Along the radial direction, the minimum heights of the triangles near the boundaries are $0.01{\rm m}$ and $0.03{\rm m}$ for ${\rm Kn}=0.1$ and $1$, respectively. The UGKS employs $100\times 100$ velocity points in the velocity space, and the initial reference number of particles $N_{r}$ for the UGKWP method is set as $400$. In comparison with the UGKS solutions, the flow fields computed by the UGKWP method are shown in Fig. 21 and Fig. 23. It can be seen that the UGKWP results agree well with those obtained from the UGKS computations. Detailed comparisons of the surface quantities, such as the pressure, shear stress, and heat flux, are given in Fig. 22 and Fig. 24. The computational cost is listed in Tables  1 and 2. The UGKS solutions are fully recovered by the UGKWP method on the unstructured meshes, but with one-order-of-magnitude lower in computational cost and memory consumption from UGKWP.

Furthermore, a very high speed flow at ${\rm Ma}=30$ is computed for the case ${\rm Kn}=0.1$ on the same unstructured mesh. Since the memory requirement of the discrete velocity points for the UGKS is unaffordable for such high Mach number computation, we only show the results of the UGKWP method in Fig. 25. In the computation, the memory cost of the UGKWP method is only $375$ MB. The advantage of the particle method with a nature adaptivity in the phase space through particles is well inherited by the UGKWP method for high speed rarefied flow computations.

IV Discussion and Conclusion

In this paper, we introduce the unified gas-kinetic wave-particle (UGKWP) method on unstructured mesh for flow simulation in all Knudsen regimes. Similar to the UGKS methodology, the direct modeling of the flow physics on numerical mesh size and time step is carried out to construct the multiscale algorithm. The early discrete velocity-based UGKS is further developed to the purely particle-based UGKP and wave-particle-based UGKWP methods. In the UGKP method, based on the integral solution of the kinetic model equation, the free transport and collision processes are well described for the evolution of simulation particles in a statistical point of view. Different from the DSMC method where simulation particles stream for a whole time step $\Delta t$ and then get possible collision, the free transport time of the simulation particles in both UGKP and UGKWP methods is obtained from the integral solution, and the free streaming convection is constrained due to particles’ interaction in different flow regimes. The collision process is handled by re-sampling simulation particles from a Maxwellian distribution according to the conservation laws. Due to the multiscale transport modeling in the UGKS methodology, the UGKWP method has no kinetic scale related time step and cell size limitations which are imposed in many other kinetic equation solvers and particle methods.

A novel wave-particle adaptive formulation is introduced to describe the microscopic gas distribution function. Specifically, the flow state in each cell contains the deterministic hydrodynamic waves and the stochastic simulation particles, and the proportion between the waves and particles evolves adaptively according to the local flow physics. Unified treatment can be carried out for all finite volume cells in the computational domain. In the continuum regimes, the hydrodynamic waves are dominant and the UGKWP method goes to a hydrodynamic fluid solver, such as GKS; while in the highly rarefied flow, the UGKWP method performs as same as the stochastic particle method. The wave-particle adaptivity makes the UGKWP method very efficient in different flow regimes by inheriting the advantages of the deterministic method and the stochastic method, and the advantages of kinetic particle transport and hydrodynamic continuum wave evolution.

Numerical test cases, including Sod problem at different Knudser numbers, low-speed micro cavity flow, laminar boundary layer for viscous NS solutions, and high-speed cylinder flow, are computed to validate the current method. It shows that the UGKWP method can recover the UGKS solutions in all flow regimes. For low-speed micro flow with a small temperature variation, the UGKS with acceleration techniques shows obvious advantages over the stochastic UGKWP method on the aspects of efficiency and accuracy because the deterministic UGKS does not suffer from the statistical noises. For the continuum flows at small Knudsen numbers and the rarefied gas flow at high Mach numbers, the UGKWP method can achieve much higher efficiency and requires lower memory cost. The unified treatment, multiscale property, high efficiency of the UGKWP method make it a very promising tool to study the multiscale problems in real engineering applications, such as the reentry of space vehicles and the high-speed near-space flights.

In the current study, we only consider the BGK model equation with a unit Prandtl number. It would not be difficult to apply other kinetic models in the current method to obtain more accurate results or to develop a more realistic method for real gas simulations. Moreover, the concepts of the wave-particle adaptive formulation and the direct modeling on the mesh size and time step scales could be used in other multiscale transport processes, such as plasma, granular flow, and radiation, in the construction of multiscale multi-efficiency methods.

V Acknowledgement

The work of Zhu and Zhong is supported by the National Natural Science Foundation of China (Grant No. 11472219), the 111 Project of China (B17037) as well as the ATCFD Project (2015-F-016). The research of Liu and Xu is supported Hong Kong research grant council (16206617) and National Natural Science Foundation of China (11772281, 91852114).