Cell Size Effect on Computational Fluid Dynamics: The Limitation Principle for Flow Simulation
LLimitation Principle for Computational Fluid Dynamics
Chang Liu a , Guangzhao Zhou d , Wei Shyy b , Kun Xu a,b,c, ∗ a Department of Mathematics, Hong Kong University of Science and Technology, Kowloon, Hong Kong, China b Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology,Kowloon, Hong Kong, China c Shenzhen Research Institute, Hong Kong University of Science and Technology, Shenzhen 518057, China d College of Engineering, Peking University, Beijing 100871, China
Abstract
For theoretical gas dynamics, the flow regimes are classified according to the physical Knud-sen number Kn s which is defined as the ratio of particle mean free path over the characteristiclength scale. The Navier-Stokes and Boltzmann equations are the corresponding modelingeqautions in the continuum Kn s ≤ − and rarefied Kn s ∼ Kn c , defined as the ratio of particle mean freepath over the cell size, will effect the solution as well. The final CFD result is controlled by anumerical Knudsen number Kn n , which is a function of Kn s and Kn c . The limitation principleis about their connections among Kn s , Kn c , and Kn n , and states the best result a numericalscheme can provide. This principle is the uncertainty theory for the capturing of a physicalsolution with a limited mesh resolution, which means the numerical flow regime changes withdifferent numerical resolution. In order to give the best result, a numerical scheme must bea multiscale method which correctly identifies flow physics in the corresponding cell Knudsennumber Kn c , such as capturing the hydrodynamic scale wave propagation in the coarse meshcase, to the kinetic scale particle transport in the fine mesh case. The unified gas-kinetic scheme(UGKS) is such a multiscale method which connects the NS to Boltzmann solutions seamlesslythrough the continuum, near-continuum, and the non-equilibrium solutions with a variationof cell Knudsen number. Even from a multiscale method, the CFD solution may not be thesame as the physical solution once the mesh resolution isn’t enough to resolve the physicalflow structure. In this paper, the mesh size effect on the representation of physical solutionis explicitly analyzed through examples. The current CFD practice, which targets on the cellconverged solution of a fixed governing equation, such as NS, is inadequate to recover a truly1 a r X i v : . [ phy s i c s . c o m p - ph ] M a y ultiple scale nature of gas dynamics. On the other hand, a scheme cannot target on theBoltzmann solution alone due to the limitation of computational resources. For an efficientCFD simulation, the representable flow physics should depend on the cell resolution under theconstraint of limitation principle. Keywords:
Cell Knudsen number, Grid-refinement, Multiscale modeling, Non-equilibriumflow
1. Introduction
The computational fluid dynamics (CFD) is more complicated than the theoretical one.For a computation, three Knudsen numbers, namely the physical, cell, and numerical ones,will effect the numerical solution and determine the quality of the scheme in the representationof a physical solution. According to physical Knudsen number, the flow regimes are classifiedinto rarefied ( Kn ≥ . < Kn ≤ . < Kn ≤ . Kn ≤ . ∗ Corresponding author
Email addresses: [email protected] (Chang Liu), [email protected] (Guangzhao Zhou), [email protected] (Wei Shyy), [email protected] (Kun Xu)
Preprint submitted to Elsevier May 28, 2018 umerical flux from the kinetic equation, and calculate the inviscid flux and viscous termsin a coupled way [3]. Comparing to the splitting approach, GKS is more physically reliable,especially under mesh refinement. As the cell size and the particle mean free path go to thesame order, there is no physical basis for the distinctive wave interaction from the macroscopicequations. However, GKS still represents NS solution due to the adoption of the Chapman-Enskog expansion for the initial gas distribution function at the beginning of each time step. Inreality, the flow physics modeling should be changed as well with the mesh refinement, insteadof solving purely the Euler, NS, or Boltzmann equations.For a flow field with high frequency mode or large gradients, the local Knudsen numbercan become large and non-equilibrium flow dynamics may emerge. The NS equations ceaseto appropriately describe such a phenomenon. Under such a circumstance, the converged NSsolution is not physically valid as the mesh size goes to the particle mean free path. With theincrease of mesh resolution, the CFD as a method for capturing physically valid solution mustbe associated with multiple scale modeling, such as solving different equations in different scales.The physical solution from a CFD algorithm should be synchronized with the cell resolution.Unfortunately, we don’t have such a traditional partial differential equation-based governingequation which could connect flow physics in different regimes smoothly. Numerically, as adirect modeling method, the unified gas kinetic scheme (UGKS) is constructed according tothe resolvable scale of a discretized scale, such as mesh size, and captures the solutions inall flow regimes through the variation of the cell Knudsen number [3]. However, even withperfect capturing of physical solution in different mesh size scale, there is still limitation due tothe mismatch between the physical flow structure thickness and the and the pre-defined meshresolution. In this paper, a limitation principle in CFD simulation will be proposed and it isrelated to the differences in the physical, cell, and numerical Knudsen numbers, to represent aphysical flow field in a discretized space.
2. The limitation principle in flow simulation
Theoretically, a physical flow field is a nature existence with different structure and dynamicsin different scales. The kinetic scale description follows the particle transport and collision, suchas the physical process in the Boltzmann equation, where the separate treatment of particletransport and collision determines its applicable scale within the particle mean free path. The3ydrodynamic equations, such as NS, are based on the macroscopic description and closedfluid elements with gigantic amount of particles inside are used in the modeling, even thoughthe scale of the fluid element has never been clearly defined [4]. Numerically, an algorithmin CFD becomes more complicated to represent a flow field in a discretized space with basicscales of cell size ∆ x and time step ∆ t , and to construct flow evolution in such a space. Thenumerical flow is the projection of the physical flow field on a discretized space. Both flowphysics and the cell resolution play dynamic roles for capturing multiple scale flow evolution.The limitation principle is about the limiting process in the representation of a physical flowfield in a discretized space.A multiscale method captures macroscopic flow behavior on a coarse mesh with the inclusionof a large amount of particles inside each control volume, and presents non-equilibrium flowdynamics on a fine mesh with the resolvable small kinetic scale physics. In a mesh refinementprocess, an idealized numerical method should provide a changeable dynamics in different scales.A multiscale method is the one to capture the corresponding flow physics on different meshresolution, such as the UGKS. Still, due to the possible mismatch between the cell resolutionand the physical flow structure thickness, there is a limitation principle about the best resultan idealized multiscale method can provide.The flow regime in CFD can be characterized by the numerical Knudsen number Kn n ,which is defined as the ratio of the molecular mean free path λ to the resolvable solution’slength scale L n , i.e., Kn n = λ/L n . The numerical Knudsen number depends on the physicalKnudsen number Kn s , and the cell Knudsen number Kn c . The relationship among Kn n , Kn s ,and Kn c gives the limitation principle about the real solution ( Kn s ) to be represented ( Kn n )numerically due to the cell resolution ( Kn c ).The physical Knudsen number is defined as the ratio of the molecular mean free path λ to the physical characteristic length scale L , i.e., Kn s = λ/L . The theoretical flow mechanicsstates that the NS equations are the modeling equations in the continuum flow regime at smallKnudsen number, with local thermodynamic equilibrium assumption. For the compressibleflow, the Knudsen number can be expressed as a function of Mach number Ma and the Reynoldsnumber Re by Kn s = (cid:114) γπ MaRe , (1)where γ is the specific heat ratio. Obviously Kn s reflects physical description of a flow field in4he characteristic length scale L . Physically, it is more reasonable to define the length scale of aflow structure, such as L s = ρ/ | (cid:53) ρ | based on density ρ , as the scale of a physical solution, i.e., Kn s = λ/L s . The NS modeling is applicable for flow structure with small Knudsen number,such as Kn s (cid:28) x , it is important to define a cell Knudsennumber: Kn c = λ ∆ x , (2)which is the ratio between the particle mean free path λ and the local cell size ∆ x . For aviscous flow, the particle mean free path, such as for a hard sphere molecule, is related to theviscosity coefficient [5], λ = (cid:114) πm kT µρ , (3)where m is the molecular mass, k is the Boltzmann constant, T is the temperature, µ is thedynamical viscosity coefficient, and ρ is the local density. Alternatively, the cell Knudsennumber can be expressed as Kn c = (cid:114) γπ MaRe c , (4)where Ma = | u | /a is now the local Mach number, Re c = | u | ∆ x/ν is the cell Reynolds number.Here u is the local flow velocity, a is the speed of sound, and ν is the kinematic viscositycoefficient ν = µ/ρ . The Kn c is related to the particle transport mechanism in the meshresolution, through the waves in the Riemann solution ( λ (cid:28) ∆ x ), or particle free transport inthe Boltzmann solver ( λ ∼ ∆ x ). With a mesh refinement, Kn c continuously increases. In thecurrent CFD practice, it isn’t surprising to see the mesh size approaching to the particle meanfree path in order to resolve the boundary layer solution in the hypersonic flow simulation. Foran ultimate CFD method, the dynamics in the algorithm should change with the cell resolution Kn c . When the mesh size ∆ x gets to the particle mean free path ∆ x (cid:39) λ , there are no intensiveparticle collisions around a cell interface to form the Riemann solution. The condition for theformation of distinct shock, contact, and rarefaction waves in the Riemann solver cannot besatisfied at all in the mesh refinement process. Therefore, to keep on using Riemann solver toget mesh converged solution is not physically founded, even for the NS solution.The CFD provides a physical evolution solution in a discretized space. The resolvable flowregimes are controlled by the numerical Knudsen number Kn n , instead of the physical and cell5nudsen numbers. The numerical Knudsen number is a function of physical and cell Knudsennumbers due to the representation of the flow physics in a cell size scale. As show in Fig. 1,for a well-resolved flow, the numerical characteristic length becomes the same as the physicalcharacteristic length. In such a case, the numerical Knudsen number goes to the physicalKnudsen number Kn n → Kn s . For a partially-solved flow, the numerical Knudsen number isin-between the physical Knudsen number and the cell Knudsen number. For an un-resolvedflow, the structure of a numerical solution is a purely numerical one with cell size thickness,and the characteristic length L n ∼ ∆ x . In such a case, the numerical Knudsen number goes tothe cell Knudsen number Kn n → Kn c . Therefore, the numerical Knudsen number is a functionof physical and cell Knudsen numbers by some soft minimum functions Kn n = softmin( Kn s , Kn c ) . (5)The solution obtained with the corresponding Kn n is the best result a best multiscale methodcan provide. The multiscale method is the one to capture flow dynamics in the mesh size scale.Here the multiscale method can smoothly connect the NS and Boltzmann solutions. With acontinuous variation of cell resolution, the multiscale method recovers the dynamics from therarefied particle transport and collision to the hydrodynamic wave propagation seamlessly.
3. Multiscale method: unified gas kinetic scheme
It is shown in Section 2 that the numerical solution depends on the numerical Knudsennumber instead of the physical and cell Knudsen numbers, but they are closely connected.Here we present a multiscale method, which could give the corresponding flow physics in thecell size scale, but this cell size may not be able to resolve the physical flow structure. The finalsolution from the multiscale method in any realistic computation will depend on the numericalKnudsen number, instead of the physical and cell Knudsen number.The multiscale method is different from the flow solvers targeting on the single scale NSand Boltzmann equations. Macroscopic equations such as the Euler and NS systems are validin hydrodynamic scale, so called continuum flow regime. The NS equations may fail to providea physically consistent solution in a highly non-equilibrium flow regime [4]. In order to capturethe physical solution under cell resolution, a multiscale method is needed. Following the directmodeling methodology [3, 6], the unified gas kinetic scheme (UGKS) has been constructed6ith the inclusion of cell size and time step effect, and becomes a multiscale method. In themodeling of UGKS, the coupling of the particle transport and collision is based on the localflow physics on the time step scale, which is crucial for the multiscale nature of UGKS. TheUGKS is constructed starting from the kinetic equation ∂f∂t + v · ∇ x f = Q, (6)where f ( x , t, v ) is the velocity distribution function and Q is collision term. For a discretephase space X × V = (cid:88) i,j Ω ij = (cid:88) i,j Ω x i × Ω v j , the evolution of cell averaged distribution function f ij = 1 | Ω ij | (cid:90) Ω ij f ( x , t, v ) d v d x , (7)is coupled with the evolution of cell averaged macroscopic conservative variables W i = 1 | Ω i | (cid:90) Ω i ρρ U ρE d x . (8)Note the above f ij and W i are cell averaged variables with a variable cell size. The governingequation for f ij is different from the Boltzmann equation, which is valid in the kinetic scaleonly, such as the scale less than the particle mean free path.The evolution equation for the cell-averaged velocity distribution function is f n +1 ij = f nij − | Ω i | (cid:90) t n +1 t n (cid:73) ∂ Ω i v · n f ∂ Ω i ( t, v j ) dsdt + β n Q n + (∆ t − β n ) Q n +1 , (9)and the evolution equation for the conservative variable is W n +1 i = W ni − | Ω i | (cid:90) t n +1 t n (cid:73) ∂ Ω i ψ v · n f ∂ Ω i ( t, v ) dsdt. (10)The above two equations are based on the direct modeling of conservation laws in a physicalscale of cell size and time step. In the following formulation, we choose Q n to be the Boltzmanncollision term and Q n +1 to be the Shakhov model. The explicit Boltzmann collision term iscalculated by the fast spectral method [7], and the evolution equation for the distributionfunction then becomes f n +1 ij = (cid:18) t − β n τ n +1 (cid:19) − (cid:34) f nij − | Ω i | (cid:90) t n +1 t n (cid:73) ∂ Ω i v · n f ∂ Ω i ( t, v j ) dsdt + β n Q ( f n , f n ) + ∆ t − β n τ n +1 f +( n +1) ij (cid:35) , (11) f + is defined as f + = g (cid:18) − Pr) c · q (cid:18) c mk B T − (cid:19) m pk B T (cid:19) , (12)with the Prandtl number Pr, the peculiar velocity c , and the heat flux q . The detail formulationof β n can be found in [6]. The local Maxwellian distribution g ( x , t, v ) is obtained from the localconservative flow variables. The numerical fluxes for the distribution function and conservativevariables are calculated from the integral solution f ∂ Ω i ( t, v ) of the kinetic equation Eq.(6) withShakhov collision term, f ∂ Ω i ( t, v ) = 1 τ (cid:90) t n + tt n f + ( x (cid:48) , t (cid:48) , v )e − ( t − t (cid:48) ) /τ dt (cid:48) + e − t/τ f ( x ∂ Ω i − v t, v ) , (13)where x (cid:48) = x ∂ Ω i − u ( t − t (cid:48) ) is the particle trajectory and f is the distribution function at time t n .Based on the integral solution Eq.(13), the UGKS couples two effects in the particle transportfor the numerical flux evaluation. The particle free transport and collision are connected basedon the ratio of local time step over the particle collision time, which is a function of local cellKnudsen number Kn c . The UGKS provides dynamics in different flow regime from the particlefree transport to the NS and Euler wave interactions [3]. The multiscale UGKS can onlyprovide a solution corresponding to the numerical Knudsen number. In different flow regimes,in order to resolve the physical solution the cell size can be varied significantly in comparisonwith the local particle mean free path. The possible discrepancy between the cell size and thedissipative flow structure thickness is the root of the limitation principle. Theoretically, theUGKS can always capture the physical solution in the mesh refinement, even in the highly non-equilibrium region. But, sometimes we don’t need such a refined mesh to resolve all dissipativewave structures, such as the shock structure. Therefore, the limitation principle still appliesto UGKS once there is discrepancy between the cell size resolution and physical structurethickness. As shown next in the numerical examples, the multiscale method provides the bestresult a numerical scheme can give efficiently. For a single scale flow solver, such as NS, themesh converged solution may not be the physical one at all. For the Boltzmann solver, eventhough it can provide accurate cell converged solution, but it always requires the cell size onthe scale of particle mean free path, which is too expensive to be used in many applications.8 . Numerical Experiments In this section, six numerical experiments will be used to demonstrate the multiscale solu-tions under different mesh resolutions. The examples include the viscous shock tube, densitysine wave propagation, 1-D shock interface interaction, 2-D shock interface interaction, shockbubble interaction, and force-driven Poiseuille flow. In the viscous shock tube problem, thehigh-order gas-kinetic scheme (GKS) is used for the NS solution. However, based on the cellKnudsen number distribution, it is realized that for a reasonable NS solution the mesh size caneasily go to the scale of the particle mean free path. The NS modeling becomes questionable incertain regions. For the density wave propagation and the shock-interface interaction problems,the NS solutions are compared with the multiscale results under different mesh resolution. Itcan be observed that when the cell Knudsen number is small, the NS solutions agree with themultiscale results, which means that the NS equations give reasonable solutions on large spaceand time scales. Relative large differences between the NS and multiscale solutions appear asthe cell Knudsen number increases in the mesh refinement process, and the NS equations don’tprovide accurate physical solution in the mesh refinement process. The shock-bubble interac-tion problem shows that for a fixed cell resolution, as the local Knudsen number increases, thesingle scale NS solutions deviate from the multiscale solutions. Even for small Knudsen num-ber, there is still deviation between the NS and multiscale solutions in the high order moment,such as the heat flux inside a shock layer. The study of force-driven Poiseuille flow concludesthat the conventional heat flux modeling, i.e. Fourier’s law, is incomplete at all for a flow underlarge external force field. The heat flux induced by the external force q Fx is on the same orderas the Fourier’s heat flux induced by the temperature gradient for Kn ≥ − and normalizedexternal force term F x ≥ − . With the decreasing of Knudsen number and external force,the force induced heat flux q Fx decays linearly with respect to Knudsen number and externalforce q Fx ∝ Kn · F x . This means that even in the continuum regime, the Navier-Stokes-Fourierequations are not adequate to provide a valid description of physical phenomena. These exam-ple show the multiple scale nature of gas dynamics, and the necessity of a multiscale method.At the same time, the final numerical solution is constrained by the limitation principle evenfrom a multiscale method. 9 .1. Viscous shock tube For any time-marching scheme, the time step ∆ t is determined by the CFL condition. Forthe Euler equations with the absence of dissipative terms, the time step is merely related tothe local convective wave speed, ∆ t ≤ σ (cid:18) | u | + a ∆ x (cid:19) − , (14)where σ is the maximum CFL number on the order of 1. For the NS equations, besides theabove time limitation, the viscosity coefficient will restrict the time step as well. With thedefinition of kinematic viscosity coefficient ν , the time step for the NS equations is determinedby [8]: ∆ t ≤ σ (cid:18) | u | + a ∆ x + C ν ∆ x (cid:19) − , (15)where C is a constant on the order of 1. The gas kinetic theory shows that the particle collisiontime τ is related to the dynamic viscosity coefficient µ (= νρ ) and local pressure p , τ = µp = γνa . (16)Then we get τ ∆ t ∼ γνσa (cid:18) | u | + a ∆ x + C ν ∆ x (cid:19) = γσ MaRe c (cid:18) C MaRe c + Ma + 1 (cid:19) . (17)With the definition cell Knudsen number in Eq. (2), Eq. (17) yields τ ∆ t ∼ Kn c + Ma + 1 C Kn c . (18)It indicates that for NS calculation τ / ∆ t is a quadratic function of the cell Knudsen number Kn c , and it is affected by the local Mach number Ma as well.In recent years, Daru and Tenaud [9] defined a viscous shock tube problem for high-speedNS solutions. In this test case, a diaphragm is vertically located in the middle of a square 2-Dshock tube with unit side length, separating the space into two parts with different initial states.Due to the symmetry, the flow in the lower half of the tube is simulated. The diaphragm isremoved instantly at time t = 0, resulting in a system of waves including a right-moving shockwith the Mach number 2 .
37 and their interactions. Non-slip conditions are used at all wallboundaries. Due to the boundary layer, incident shock, and reflecting shock wave interactions,10omplicated flow structures emerge with a large variation of wave structures. The flow fieldsat the Reynolds number 200 and 1000 and time t = 1 are used in the following investigation.The problem has been studied extensively and the results presented here are from a sim-plified high-order gas-kinetic scheme [10], which provides similar converged solutions as otherschemes on commonly agreed mesh size [11]. Though the GKS does not directly solve the NSequations, it is an accurate NS solver in the well-resolved region due to the use of Chapman-Enskog expansion for the initial gas distribution function reconstruction. Note that GKS is forNS solution and UGKS is a multiscale method for flow simulation in all regimes.Fig. 2(a) shows the variation of two dimensionless quantities λ/ ∆ x and τ / ∆ t in the wholecomputation domain at the Reynolds number 200. It is clear that the maximum τ / ∆ t on the1500 ×
750 mesh is around 2.5, indicating that the computational time step is already smallerthan the particle collision time. In many real engineering applications, it is common to see sucha mesh in the resolved hypersonic viscous flow computations.It should be noted that although a smaller kinematic viscosity ν leads to a larger ∆ t ac-cording to Eq. (15), in order to resolve the dissipative layer corresponding to smaller viscositycoefficient, a smaller cell size ∆ x is also needed. As a result, τ / ∆ t can be very large for higherReynolds number simulations. To make it clear, Fig. 2(b) plots the same curves as Fig. 2(a) butwith a Reynolds number 1000. The mesh for a converged flow field at Re = 200 is 500 × × τ / ∆ t is still large on a fine mesh for Re = 1000 case.For both cases, the physical Knudsen number for the whole problem is on the order of0 . ∼ .
01, because the flow has a Mach number on the order 1 and the Reynolds numberbetween 100 and 1000. However, with a fine mesh resolution the cell Knudsen number caneasily reach the order of 1. As shown in the figure, the ratio τ / ∆ t can go beyond 3. Thisclearly indicates the necessity to use a multiscale method for its solution at local high Kn c region.Eq. (4) shows that Kn c is inversely proportional to √ ρp ∆ x . Eq. (18) indicates that a largeMach number leads to a large τ / ∆ t . For the Re = 1000 case, the distributions of λ/ ∆ x and τ / ∆ t are presented in Fig. 3. And the distributions of density, pressure and Mach number areshown in Fig. 4. All of them are obtained with a 5000 × In the next three test cases, we compare the NS and multiscale solutions under different cellKnudsen numbers. First, we study the propagation of a density wave in argon gas. The initialcondition is set as ( ρ, U, p ) = (5 sin(30 x ) + sin( x ) + 10 , . , . . The computational domain is [0 , π ] with periodic boundary condition. The variable hardsphere (VHS) model is used with the viscosity temperature dependency index ω = 0 .
81. TheKnudsen number is set to be Kn = 5 × − . The velocity space is truncated from − Q , the deviation between the UGKSand NS solutions is calculated as Dev = (cid:107) Q NS − Q UGKS (cid:107) L (cid:107) Q UGKS (cid:107) L . (19)12ensity wave propagation Kn c τ / ∆ t Dev Fine mesh ∼ × − ∼ × − . × − Coarse mesh ∼ × − ∼ × − . × − Table 1: Parameters under different meshes, deviation between NS and multiscale solution for the density wavepropagation.
The initial density distribution is shown in Fig. 5, and the initial condition are numericallyresolved on both meshes. The density distribution as well as the cell Knudsen number at t = 4 π are shown in Fig. 6-7, and the deviations are listed in Table 1.On a relative coarse mesh, the NS solution agrees well with the UGKS solution. On such anumerical cell size scale, the cell Knudsen number is relatively small Kn c ∼ × − , and the cellaveraged velocity distribution is close to the local equilibrium. Under the limitation principle,the numerical solution is the one corresponding to cell Knuden number. When the grid pointsin physical space increase to 1000, the local non-equilibrium effect in the high frequency modeappears. Under such a resolution, the local flow physics is resolved by the mesh size, and thenumerical Knudsen number for the solution is the same as the physical Knudsen number. TheNS solution deviates from the physical one. The results show that the physical solution is moredissipative than the NS solution due to the inclusion of particle transport mechanism in such ascale. So, even with mesh converged solution, the NS solution still doesn’t capture the physicalsolution in the corresponding resolution. The UGKS can accurately recover the attenuation ofhigh frequency wave in the experiments [12]. In the shock interface interaction case, non-equilibrium effect beyond the NS modeling willemerge when the shock wave hits the interface with large gradients in the flow variables. TheNS solutions will be compared with the multiscale solutions under different cell resolution. Onthe relative coarse mesh case, the cell Knudsen number is small and the NS solver will have thesame results as the multiscale ones, because under limitation principle this is the best result amultiscale method can have under such a mesh resolution. When the cell size decreases andcell Knudsen number increases, the NS solutions deviate from the multiscale ones, which is dueto the limitation of NS modeling in capturing the non-equilibrium phenomena.For 1-D shock interface interaction, the working gas is argon modeled by VHS. As shown13hock interface interaction Kn c τ / ∆ t Dev Fine mesh ∼ ∼
50 3 . × − Coarse mesh ∼ . ∼ . × − Table 2: Parameters under different meshes, deviation between NS and multiscale solution for the shock interfaceinteraction. in Fig. 8, initially a normal shock wave with Ma = 3 . x = 0. The upstream initial condition is set as( ρ, u, T ) = (1 . , . , . x ≥ − , (2 . , . , . x < − . The x-axis is normalized by the shock upstream mean free path, which means the Knudsennumber is one. The NS solutions are compared with the multiscale solutions under a fine mesh∆ x = 0 . x = 2. For UGKS, the velocity domain is [ − ,
8] covered by100 velocity points. The solutions at t = 80 are shown in Fig. 9 when the interface movesto the origin and interacts with the shock wave. As shown in Table 2, a larger deviation isobserved under a fine mesh since the NS equations are not able to recover the physical structureof shock wave, which is supposedly resolved in such a fine mesh case. Under a coarse mesh, theshock wave structure is not resolved by the numerical resolution, and the physical shock waveis replaced by a numerical discontinuity, such as the shock capturing scheme. Therefore, boththe NS and the multiscale method give the same result under the limitation principle. The Richtmyer-Meshkov (RM) instability is caused by a shock wave passing though a contactinterface. Vortices will be generated during the passage of the shock wave and trigger interfaceinstability. In the following, the 2-D RM instability is numerically studied on different timescales. The interaction between shock and interface is studied on the time scale t ∼ τ , andthe development of the instability is studied on the time scale t ∼ τ . The NS solutions arecalculated by GKS while the UGKS provides the multiscale solutions. At the starting time, theinitial condition is shown in Fig. 10(a). The computational domain is [ − . , . × [ − . , . Kn c τ / ∆ t Dev Fine mesh ∼ ∼
50 1 . × − Coarse mesh ∼ . ∼ . × − Table 3: Parameters under different meshes, deviation between NS and multiscale solution for the Richtmyer-Meshkov instability.
Knudsen number is Kn = 5 × − . Initially, a shock wave with Ma = 3 . x = − .
15, and a contact interface is located at x = 0 as( ρ, u, T ) = (1 . , , . x < sin(2 πy + 1 . π ) , (5 . , , . x > sin(2 πy + 1 . π ) . The shock wave is pre-calculated and the fully developed shock structure is used as the initialcondition. For UGKS, the velocity domain is [ − , × [ − , ×
32 velocity grids areused. The NS and multiscale solutions are calculated under a fine mesh 240 ×
600 with Kn c ∼ ×
100 with Kn c ∼ .
1. The density distribution at t = 0 . ≈ τ along y = 0 is shown in Fig. 10(b-c). Under the coarse mesh, the NS solutions agree well with themultiscale solutions with Dev = 3 . × − , while large deviation can be observed under a finemesh with Dev = 1 . × − . The data is presented in Table 3. The UGKS solution is limitedby the limitation principle in the coarse mesh case.Next, we calculate the development of the interface instability for a longer time scale. Thecomputational domain is [ − . , . × [ − . , .
5] covered by 150 ×
100 cells in the physicaldomain. The Knudsen number is Kn = 1 . × − , and the Mach number of the shock wave is Ma = 1 .
3. The contact discontinuity located at x = 0 is( ρ, u, T ) = (1 . , − . , . x < sin(2 πy + 1 . π ) , (2 . , − . , . x > sin(2 πy + 1 . π ) . The solutions at t = 18 . ≈ . × τ are shown in Fig.11-12, where the density, cell Knudsennumber, vorticity magnitude, and streamline are plotted together with the NS solution (up)and the UGKS solution (down). On such large space and time scales, the NS solutions agreewell with the multiscale solutions. 15 .5. Shock Bubble interaction In this section, we study the process of a shock wave interacting with dense cold bubbleto show the capability of NS equations in describing the flow with large gradients. The initialcondition is shown in Fig. 13. A shock wave is situated with its density barycenter at x = − .
0, travelling in the positive x-direction into a flow field at rest with ( ρ, U, p ) = (1 . , , . x, y ) = (0 . , p = 0 .
5. Thecomputational domain is [ − , × [ − , − , × [ − ,
8] covered by 100 ×
100 velocity grid points.First, we set the Mach number of the shock wave to be Ma = 1 .
3, and set the density ofthe bubble to be ρ ( x, y ) = 1 + 0 . − x + y )) . (20)Two different Knudsen numbers are considered relative to the bubble radius, namely Kn =1 . × − and Kn = 0 .
3. For Kn = 1 . × − case, 250 ×
100 cells are used in the physicaldomain. The solutions of UGKS at t = 1 . Kn = 0 .
3, the physical domain is divided into 100 ×
40 mesh points. Solutions at t = 1 . Ma = 2 . ρ ( x, y ) = 1 + 1 . − x + y )) . (21)The solutions at small Knudsen number Kn = 1 . × − and at t = 1 . . The NSdensity and temperature profiles agree well with the multiscale solutions, which have a smallincrement in the deviation in comparison with Ma = 1 . Kn = 1 . × − case. For the16eat flux, large deviation can be clearly observed, especially for the x-directional heat flux with43% deviation. When the Knudsen number is increased to Kn = 0 .
3, relative large deviationin all flow variables can be observed in Fig. 18.For the shock bubble interaction at a relative large cell Knudsen number, the flow physicsis not in NS regime and the NS equations don’t provide accurate solutions. At a small cellKnudsen number, the low order moments from the NS, such as density and temperature, agreewith multiscale solution due to the cell averaged effect under the limitation principle. However,even at a small cell Knudsen number the high order moments, such as the heat flux, will notbe well predicted by NS equations at the region with large flow gradient. Therefore, even inthe near continuum flow regime, the multiscale numerical method can provide more accuratephysical solutions under the same cell resolution as NS solver.
In the following, through force-driven Poiseuille flow we study the capability of NS equationsin the description of non-equilibrium flow physics, even in the continuum flow regime. Thisexample clearly indicates the multiscale nature of the flow physics and the NS modeling is notfully complete for the physical solution.The Poiseuille flow is about the flow confined between two parallel isothermal plates. Aconstant body force is applied in the x-direction along the channel. Periodic boundary conditionis taken in the flow direction. The simulated gas is a hard sphere monatomic gas. The initialcondition is set as ρ = 1 . u x = 0, T = 1. In the following, the solutions at a variety ofKnudsen numbers and external force will be obtained. The Knudsen number is defined by theratio between the mean free path and the wall distance, Kn = λ/L y . The results from DSMCand NS have been previously presented [13].The first test has a Knudsen number Kn = 0 . F x = 1 . × − . Thephysical domain is divided into 51 cells in the y-direction, and velocity domain [ − , is dividedby 32 velocity points in x,z-direction and 64 velocity points in y-direction. For this set ofparameters, the time step is about 100 times smaller than the local collision time, therefore theBoltzmann collision term takes effect in Eq.(9). It can be observed in Fig. 19-21 that the UGKSsolutions agree well with DSMC results in the prediction of density, velocity, temperature,stress, and heat flux. The NS solution is calculated by GKS with the first order slip boundary17ondition u slip = αλ ∂u∂y , T slip = β γγ + 1 λ Pr ∂T∂y , with α = 1 .
11 and β = 1 .
13. In the flow resolved calculation, the non-equilibrium effect is notwell captured by NS equations at all, especially the stress and heat flux.Next, we reduce the Knudsen number and external force to F x = Kn = 5 × − and F x = Kn = 2 × − . For small Knudsen number and small acceleration, an asymptoticsolution from the BGK equation is derived by Aoki et al. [14]. The UGKS solutions arecompared with the second order asymptotic solutions and NS solutions. As shown in Fig. 22,asymptotic solution agrees well with DSMC solution especially for Kn = 2 × − , while theNS solutions always predict zero heat flux in the flow direction.Next, we study the Poiseuille flow under a relatively large external force. The calculationis performed with three set of parameters: Kn = 0 . , F x = 1, Kn = 5 × − , F x = 0 .
5, and Kn = 2 × − , F x = 0 .
2. The x-directional flux results are shown in Fig. 23. It is shownthat the heat flux induced by the external force is proportional to the first order of Knudsennumber, which is on the same order as the Fourier’s heat flux induced by the temperaturegradient. Therefore, with the external force field, the Fourier heat flux in the NS equationsis not complete in the description of energy transport, even in the continuum flow regime[15]. This test case indicates the importance of a multiscale method for the capturing ofintrinsically multiple scale flow dynamics. The limitation principle presented in this paper isfor the multiscale method. For the NS solver, the quality of the solution is far away from thelimiting solution constrained by the limitation principle.
5. Conclusion
Gas dynamics is associated with multiscale nature. The NS equations describe a flow physicsfor very small variation of flow variables, and this variation is measured by the physical Knudsennumber. At the same time, the NS uses the fluid element for modeling the dynamics, and thescale of the fluid element is not explicitly defined, which introduces inaccuracy in the predictionof physical solution even in the analytical level. For CFD computation, the situation becomesmore complicated in a discretized space. Here we have to have a clear picture about themodeling scale of the governing equation and the numerical scale, i.e., the mesh size and timestep, for the evaluation of the solution. With the variation of the mesh size, the coarse-grained18ow variables have a numerical dynamics with the corresponding numerical Knudsen number,and only a limited solution can be recovered by a numerical scheme. With the definitions ofcell Knudsen number Kn c and the physical solution Knudsen number Kn s , the representablenumerical solution in a mesh size scale is controlled by the numerical Knudsen number Kn n ,which is a function of Kn c and Kn s with the relationship Kn n = softmin( Kn s , Kn c ). This is thelimitation principle which states the best numerical solution can be obtained by a multiscaleCFD method, such as the UGKS. In order to understand such a limitation principle, in thispaper we first demonstrate the physical modeling deficiency in the NS equations in the meshrefinement process. As the mesh size goes to the particle mean free path in the high speedflow simulation, the wave interaction modeling, such as Riemann solver for the NS solution,doesn’t hold anymore. Instead, particle transport and collision will be the correct physicshere. Therefore, corresponding to different mesh size scale relative to the particle mean freepath, we need a consistent gas dynamics in a numerical scheme. The flow transport in UGKSdepends on the ratio of ∆ t/τ . With a variation of τ / ∆ t , the multiscale gas dynamics is coveredfrom the particle transport as τ (cid:39) ∆ t to the wave interaction as τ (cid:28) ∆ t , with a continuousvariation between them [4]. Even with a multiscale method, due to possible mismatch betweenthe mesh size and physical structure thickness, the final numerical solution is constrained bythe limitation principle with a newly defined the numerical Knudsen number instead of thephysical one. The numerical gas dynamic solution is related to the flow variable variation(physical Knudsen number) and mesh resolution (cell Knudsen number). In order to illustratesuch a limitation principle, many numerical examples are presented, such as shock interfaceinteraction and Richtmyer-Meshkov instability, where different solutions can be identified underdifferent mesh resolution. The force-driven Poiseuille flow clearly shows the multiscale natureof gas dynamics. Based on the UGKS solution, even in the continuum flow regime it showsthat the single scale NS equations are not adequate to present a physical solution. In a CFDsimulation, we need to consider both limitations from the governing equations and the cellsize resolution. The limitation principle is about the best result a best mutiscale method cansimulate in the mesh size scale. Any CFD simulation is associated with uncertainty principledue to the variation of mesh resolution. There is no faithful cell converged physical solution ifthe flow solver is targeting on the NS equations alone, even in the continuum flow regime.19 cknowledgement The current research was supported by Hong Kong research grant council (16206617,16207715,16211014),and National Science Foundation of China (91530319,11772281).
ReferencesReferences [1] E. Toro, Riemann solvers and numerical methods for fluid dynamics, Springer-Verlag, Berlin, 2009.[2] J. Li, Q. Li, K. Xu, Comparison of the generalized Riemann solver and the gas-kinetic scheme for inviscidcompressible flow simulations, Journal of Computational Physics 230 (12) (2011) 5080–5099.[3] K. Xu, Direct Modeling for Computational Fluid Dynamics: Construction and Application of UnifiedGas-kinetic Scheme, World Scientic, 2015.[4] K. Xu, C. Liu, A paradigm for modeling and computation of gas dynamics, Physics of Fluids 29 (2017)026101.[5] W. Vincenti, C. Kruger, Introduction to Physical Gas Dynamics, John Wiley & Sons, Inc., 1965.[6] C. Liu, K. Xu, Q. Sun, Q. Cai, A unified gas-kinetic scheme for continuum and rarefied flows IV: FullBoltzmann and model equations, Journal of Computational Physics 314 (2016) 305–340.[7] L. Wu, C. White, T. J. Scanlon, J. M. Reese, Y. Zhang, Deterministic numerical solutions of the Boltzmannequation using the fast spectral method, Journal of Computational Physics 250 (2013) 27–52.[8] J. Blazek, Computational Fluid Dynamics: Principles and Applications, 3rd Edition, Elsevier, 2015.[9] V. Daru, C. Tenaud, Evaluation of TVD high resolution schemes for unsteady viscous shocked flows,Computers & Fluids 30 (2001) 89–113.[10] G. Zhou, K. Xu, F. Liu, Simplification of the flux function for a high-order gas-kinetic evolution model,Journal of Computational Physics 339 (2017) 146–162.[11] G. Zhou, K. Xu, F. Liu, Grid-converged solution and physical dynamic process analysis of the viscousshock tube problem, submitted to Journal of Fluid Mechanics.[12] R.-J. Wang, K. Xu, The study of sound wave propagation in rarefied gases using unified gas-kinetic scheme,Acta Mechanica Sinica 28 (4) (2012) 1022–1029.[13] Y. Zheng, A. L. Garcia, B. J. Alder, Comparison of kinetic theory and hydrodynamics for poiseuille flow,Journal of Statistical Physics 109 (3-4) (2002) 495–505.[14] K. Aoki, S. Takata, T. Nakanishi, Poiseuille-type flow of a rarefied gas between two parallel plates drivenby a uniform external force, Physical Review E 65 (2) (2002) 026315.[15] T. Xiao, K. Xu, Q. Cai, T. Qian, An investigation of non-equilibrium heat transport in a gas system underexternal force field, International Journal of Heat and Mass Transfer 126 (2018) 362–379.
10 -8 -6 -4 -2 0 2 4 6 8 10 X Physical and Numerical Flow Field for a Well-Resolved Flow PhysicalFlow fieldNumericalFlow FieldPhysical and NumericalCharacteristic Length L x (a) -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 X Physical and Numerical Flow Field for a Partially-Resolved Flow PhysicalFlow fieldNumericalFlow FieldNumerical Characteristic scale L n Physical Characteristic Length L s x (b) -10 -8 -6 -4 -2 0 2 4 6 8 1000.10.20.30.40.50.60.70.80.91 Physical and Numerical Flow Field for an Under-Resolved Flow x Numerical Characteristic Length L n PhysicalFlow fieldNumericalFlow FieldPhysical CharacteristicLength L s (c) Figure 1: Diagram for physical flow field (solid line) and reconstructed numerical flow field (dotted line) for awell-resolved smooth flow (a), a partially-resolved flow (b), and an under-resolved flow with large gradient (c).
Number of cells in x direction
250 500 750 1000 1250 15000.511.522.5 max( λ / ∆ x)max( τ / ∆ t) (a) Number of cells in x direction λ / ∆ x)max( τ / ∆ t) (b)Figure 2: Maximum λ/ ∆ x and τ / ∆ t for different cell numbers at (a) Re = 200 and (b) Re = 1000. a) (b)Figure 3: Contour map of (a) λ/ ∆ x and (b) τ / ∆ t of the viscous shock tube problem at t = 1 with Re = 1000.5000 × t = 1 with Re = 1000. 5000 × a b Figure 5: The initial density distribution of the density wave propagation on coarse mesh (a) and on fine mesh(b). x D en s i t y Density distribution at t=4 ( x=2 /150)
NS (GKS)Multiscale (UGKS)
Dev=7.1514 10 -8 a x / x Loacl Knudsen number at t=4
NS (GKS)Multiscale(UGKS) b Figure 6: The density distribution at t = 4 π on coarse mesh (a), and the cell Knudsen number at t = 4 π oncoarse mesh (b). x D en s i t y Density distribution at t=4 ( x=2 /1000)
NS (GKS)Multiscale (UGKS) -6 a x / x Loacl Knudsen number at t=4
NS (GKS)Multiscale (UGKS) b Figure 7: The density distribution at t = 4 π on fine mesh (a), and the cell Knudsen number at t = 4 π on finemesh (b).
250 -200 -150 -100 -50 0 50 100 150 200 250 x y Initial density distribution
Figure 8: Initial condition of the shock-contact interaction -50 -40 -30 -20 -10 0 10 20 30 40 50 x D en s i t y Density distribution ( x=0.2 )
NS solution (GKS)Multiscale solution (UGKS)
Dev=3.24 10 -4 a -50 -40 -30 -20 -10 0 10 20 30 40 50 x D en s i t y Density distribution ( x=2 )
NS solution (GKS)Multiscale solution (UGKS)
Dev=1.47 10 -4 b Figure 9: Density of the shock-contact interaction at t = 80. The solutions of NS and multiscale method arecompared under a fine mesh (a) and a coarse mesh (b). -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 x D en s i t y Density distribution along y=0
NS SolutionMultiscale solution -0.12 -0.11 -0.1 -0.09 -0.08 -0.0734567891011
Dev=1.6 10 -3 Kn c b -0.5 0 0.5 1 1.5 x D en s i t y Density distribution along y=0
NS solution (GKS)Multscale solution (UGKS) -0.18 -0.16 -0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 033.544.555.56
Dev=3.15 10 -5 Kn c c Figure 10: Left figure shows the initial condition of RM instability. Right figures show the density distributionalong y = 0 at t = 0 . ensity Distribution -0.5 0 0.5 100.050.10.150.20.250.30.350.40.450.5 -0.5 0 0.5 1-0.5-0.45-0.4-0.35-0.3-0.25-0.2-0.15-0.1-0.050 Navier-Stokes Solution(GKS)Multiscale Solution(UGKS) a Cell Knudsen Number -0.5 0 0.5 100.050.10.150.20.250.30.350.40.450.5-0.5 0 0.5 1-0.5-0.45-0.4-0.35-0.3-0.25-0.2-0.15-0.1-0.050 -3 Navier-Stokes Solution(GKS)Multiscale Solution(UGKS) b Figure 11: The density (a) and cell Knudsen (b) at t=18.6, with the NS solution (up) and the down multiscalesolution (down).
Vorticity Magnitude -0.5 0 0.5 100.050.10.150.20.250.30.350.40.450.5 -0.5 0 0.5 1-0.5-0.45-0.4-0.35-0.3-0.25-0.2-0.15-0.1-0.050
Navier-Stokes Solution(GKS)Multiscale Solution(UGKS) a -0.5 0 0.5 100.050.10.150.20.250.30.350.40.450.5 Streamline -0.5 0 0.5 1-0.5-0.45-0.4-0.35-0.3-0.25-0.2-0.15-0.1-0.050
Navier-Stokes Solution(GKS)Multiscale Solution(UGKS) b Figure 12: The vorticity (a) and streamline (b) at t=18.6, with the NS solution (up) and the multiscale solution(down). igure 13: Initial condition for the shock bubble interaction process. Density
02 0 23 -14 10.5-2 0 0
Temputure
X-direction heat flux -3 Y-direction heat flux -4 Figure 14: Two-dimensional plots of the fields of density, temperature and heat fluxes at t = 1 .
28 b -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x D en s i t y Density distribution along x=0
NS (GKS)Multiscale (UGKS)
Dev=1.14 10 -8 c -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x T e m pe r a t u r e Temperature distribution along x=0
NS (GKS)Multiscale (UGKS)
Dev=8.92 10 -9 d -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -1-0.500.511.522.533.5 X - d i r e c t i on hea t f l u x -4 q x distribution along x=0 NS (GKS)Multiscale (UGKS)
Dev=0.0534 e -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -10123456 Y - d i r e c t i on hea t f l u x -5 q y distribution along y=-0.4 NS (GKS)Multiscale (UGKS)
Dev=0.0026 f Figure 15: Results of the shock bubble interaction with Ma = 1 . Kn = 1 . × − at t = 1 .
3. Top twofigures shows the cell Knudsen number (a) and total heat flux (b). (c)-(f) show the comparison between NSsolution profiles (solid line) and UGKS solution profiles (dash line).
29 b -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x D en s i t y Density distribution along y=0
NS (GKS)Multiscale (UGKS)
Dev=8.59 10 -5 c -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x T e m pe r a t u r e Temperature distribution along y=0
NS (GKS)Multiscale (UGKS)
Dev=4.16 10 -5 d -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -0.0100.010.020.030.040.050.06 X - d i r e c t i on hea t f l u x q x distribution along y=0 NS (GKS)Multiscale (UGKS)
Dev=0.0239 e -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -2024681012 Y - d i r e c t i on hea t f l u x -3 q y distribution along y=-0.4 NS (GKS)Multiscale (UGKS)
Dev=0.0327 f Figure 16: Results of the shock bubble interaction with Ma = 1 . Kn = 0 . t = 1 .
3. Top two figuresshows the cell Knudsen number (a) and total heat flux (b). (c)-(f) show the comparison between NS solutionprofiles (solid line) and UGKS solution profiles (dash line).
30 b -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x D en s i t y Density distribution along y=0
NS (GKS)Multscale (UGKS)
Dev=3.45 10 -7 c -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x T e m pe r a t u r e Temperature distribution along y=0
NS (GKS)Multscale (UGKS)
Dev=1.81 10 -7 d -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -0.500.511.522.53 X - d i r e c t i on hea t f l u x -3 q x distribution along y=0 NS (GKS)Multscale (UGKS)
Dev=0.43 e -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -10123456 Y - d i r e c t i on hea t f l u x -4 q y distribution along y=-0.4 NS (GKS)Multscale (UGKS)
Dev=0.0068 f Figure 17: Results of the shock bubble interaction with Ma = 2 and Kn = 1 . × − at t = 1 .
0. Top twofigures shows the cell Knudsen number (a) and total heat flux (b). (c)-(f) show the comparison between NSsolution profiles (solid line) and UGKS solution profiles (dash line).
31 b -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x D en s i t y Density distribution along y=0
NS (GKS)Multiscale (UGKS)
Dev=2.17 10 -4 c -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x T e m pe r a t u r e Temperature distribution along y=0
NS (GKS)Multiscale (UGKS)
Dev=5.79 10 -4 d -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -0.0500.050.10.150.20.250.30.350.4 X - d i r e c t i on hea t f l u x q x distribution along y=0 NS (GKS)Multiscale (UGKS)
Dev=0.1247 e -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 x -0.0200.020.040.060.080.10.120.14 Y - d i r e c t i on hea t f l u x q y distribution along y=-0.5 NS (GKS)Multiscale (UGKS)
Dev=0.0642 f Figure 18: Results of the shock bubble interaction with Ma = 2 and Kn = 0 . t = 1 .
0. Top two figures showsthe cell Knudsen number (a) and total heat flux (b). (c)-(f) show the comparison between NS solution profiles(solid line) and UGKS solution profiles (dash line). Y Density Distribution
DSMC SolutionUGKS SolutionNS Solution a Y U Velocity Distribution
DSMC SolutionUGKS SolutionNS Solution b Y T Temperature Distribution
DSMC SolutionUGKS SolutionNS Solution c Y p Pressure Distribution
DSMC SolutionUGKS SolutionNS Solution d Figure 19: The (a) density, (b) x velocity, (c) temperature, (d) pressure distribution for the force-driven Poiseuilleflow at steady state with Knudsen number Kn = 0 . F x = 0 . Y p xx Stress distribution p xx DSMC SolutionUGKS SolutionNS Solution a Y p yy Stress Distribution p yy DSMC SolutionUGKS SolutionNS Solution b Y -0.15-0.1-0.0500.050.10.15 p xy Stress Distribution p xy DSMC SolutionUGKS SolutionNS Solution c Y p zz Stress Distribution p zz DSMC SolutionUGKS SolutionNS Solution d Figure 20: The stress distribution for the force-driven Poiseuille flow at steady state: (a) p xx , (b) p yy , (c) p xy ,(d) p zz . The solid lines are the UGKS solutions, the dashed lines are the NS solutions and symbols stand forDSMC solutions. Y -0.1-0.08-0.06-0.04-0.0200.020.040.06 Q x X-directional heat flux
DSMC SolutionUGKS SolutionNS Solution a y -0.015-0.01-0.00500.0050.010.015 Q y Y-directional heat flux
DSMC SolutionUGKS SolutionNS Solution b Figure 21: The heat flux distribution for the force-driven Poiseuille flow at steady state: (a) x-directional heatflux Q x , (b) y-directional heat flux Q y . The solid lines are the UGKS solutions, the dashed lines are the NSsolutions and symbols stand for DSMC solutions. Y -1.5-1-0.500.511.522.53 Q x -3 X-directional heat flux (Kn=0.05)
Asymptotic Solution O(Kn )UGKS SolutionNS Solution O(Kn) a Y -2-1012345678 Q x -4 X-directional heat flux (Kn=0.02)
Asymptotic Solution O(Kn )UGKS SolutionNS Solution O(Kn) b Figure 22: X-directional heat flux with (a) Knudsen number Kn = 0 .
05 and acceleration F x = 0 .
05, (b)Knudsen number Kn = 0 .
02 and acceleration F x = 0 .
02. The UGKS solutions (solid lines) are compared withthe asymptotic solutions (symbols) and NS solutions (dashed lines). Y -0.100.10.20.30.40.50.6 Q x X-directional heat flux (Kn=0.1, F x =1) UGKS SolutionNS Solution a Y -0.0200.020.040.060.080.10.120.140.16 Q x X-directional heat flux (Kn=0.05, F x =0.5) UGKS SolutionNS Solution b Y -0.00500.0050.010.0150.020.0250.030.035 Q x X-directional flux (Kn=0.02, F x =0.2) UGKS SolutionNS Solution c Figure 23: X-directional heat flux with relative large external force (a) Knudsen number Kn = 0 . F x = 1 .
0, (b) Knudsen number Kn = 0 .
05 and acceleration F x = 0 .
5, (c) Knudsen number Kn = 0 .
02 andacceleration F x = 0 .
2. The solid lines represent UGKS solutions and dashed lines stand for the NS solutions.2. The solid lines represent UGKS solutions and dashed lines stand for the NS solutions.