A simplified discrete unified gas kinetic scheme for incompressible flow
Mingliang Zhong, Sen Zou, Dongxin Pan, Congshan Zhuo, Chengwen Zhong
AA simplified discrete unified gas kinetic scheme for incompressibleflow
Mingliang ZHONG, ∗ Sen ZOU, † Dongxin Pan, ‡ Congshan Zhuo, § and Chengwen Zhong ¶ School of Aeronautics, Northwestern Polytechnical University,Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China. (Dated: July 16, 2020)
Abstract
The discrete unified gas kinetic scheme (DUGKS) is a new finite volume (FV) scheme for con-tinuum and rarefied flows which combines the benefits of both Lattice Boltzmann Method (LBM)and unified gas kinetic scheme (UGKS). By reconstruction of gas distribution function using par-ticle velocity characteristic line, flux contains more detailed information of fluid flow and moreconcrete physical nature. In this work, a simplified DUGKS is proposed with reconstruction stageon a whole time step instead of half time step in original DUGKS. Using temporal/spatial integralBoltzmann Bhatnagar-Gross-Krook (BGK) equation, the transformed distribution function withinclusion of collision effect is constructed. The macro and mesoscopic fluxes of the cell on nexttime step is predicted by reconstruction of transformed distribution function at interfaces alongparticle velocity characteristic lines. According to the conservation law, the macroscopic variablesof the cell on next time step can be updated through its macroscopic flux. Equilibrium distribu-tion function on next time step can also be updated. Gas distribution function is updated by FVscheme through its predicted mesoscopic flux in a time step. Compared with the original DUGKS,the computational process of the proposed method is more concise because of the omission of halftime step flux calculation. Numerical time step is only limited by the Courant-Friedrichs-Lewy(CFL) condition and relatively good stability has been preserved. Several test cases, including theCouette flow, lid-driven cavity flow, laminar flows over a flat plate, a circular cylinder, and anairfoil, as well as micro cavity flow cases are conducted to validate present scheme. The numericalsimulation results agree well with the references’ results.
PACS numbers: 47.45.Ab, 47.11.-j, 47.11.Df, 47.15.-x ∗ Zhong, M. L.: [email protected] † Zou, S.: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] J u l [email protected] § [email protected] ¶ Corresponding author: [email protected] . INTRODUCTION Lattice Boltzmann method (LBM) was developed from the lattice gas automation model.Its implementation can be divided into two steps: streaming and collision. Since the discretevelocity models are coupled with the computational mesh, the distribution function can beprecisely streaming from one node of mesh to its neighbor node in one time-step [1]. Themacroscopic flow variables can then be calculated using the distribution function. LBMrequires a small amount of coding effort while the computational and memory cost for largescale problems are significant [2].The standard LBM is constructed based on the Cartesian mesh, which leads to a short-coming in adaption to extensive application. For sufficient resolution, mesh points haveto be distributed in computational field with uniform spatial scale, which usually producesenormous mesh cells and huge computational cost. In practical applications, we have todispose of geometry with complex boundaries. Arbitrary mesh nodes distribution cannot bewell fitted by uniform lattices. These two issues reduce computational efficiency and limitits application in more complex flow cases.In the past few decades, there have been many attempts to apply LBM to arbitrarygeometrical computational domains. X. He, et al. proposed a mesh based mesoscopicinteraction for LBM [3]. The advantages of LBM are retained and simulation results agreewell with experimental data. J. Yao, et al. introduced an adaptive mesh refinement (AMR)method which refines meshes by constructing a linked-list of geometries for mesh levelsrefinement. Linked-list mesh data and LBM calculation are combined [4]. X. Guo, et al.applied AMR to immersed boundary LBM and bubble function in interpolation [5]. Thesemethods extend LBM to more universal cases. However, numerical time step is still restrictedby particle relaxation time. O. Aursjø, et al. proposed local lattice Boltzmann algorithm forincluding a general mass source term that results in Galilean invariant continuum equations[6]. D. Wang, et al. used LBM to simulate viscoelastic drops [7].Recently, finite volume LBM (FVLBM) is developed to apply LBM to hybrid meshes andovercome the time step restriction. In this method, the body-fitted mesh and hybrid meshcan be used. S. Succi, et al. were the first to propose a finite volume (FV) scheme for discreteBoltzmann equations [8]. G. Peng, et al. developed the FV scheme based on LBM which wasfirst applied to the unstructured mesh [9]. M. Stiebler, et al. developed an upwind discrete3cheme for the FVLBM [10]. W. Li, et al. proposed a grid-transparent FVLBM, whichshows low dissipation and high accuracy in simulation of viscous flows [11]. Y. Wang, et al.studied the performance of a FVLBM scheme which is discretized on a unstructured meshand simulated the steady and unsteady flows at relatively high Reynolds number. Resultsshow lower computational cost and good agreement with previous benchmark data [12, 13].L. Chen proposed a unified and preserved Dirichlet boundary treatment [14].Generally, in original FVLBM, numerical time step is still limited by relaxation time. Ona large time step, the computational process is no longer stable and non-physical oscilla-tion occurs. Due to temporal integration of collision term decoupled with advection term,the lack of real particle collision mechanism in flux evaluation restricts marching time stepto a same order of magnitude with collision time [15]. Discrete unified gas kinetic scheme(DUGKS) is a new kind of FV scheme for discrete velocity method (DVM), which com-bines the benefits of both lattice Boltzmann equation (LBE) and unified gas kinetic scheme(UGKS) methods [16]. In DUGKS, the evolution of the flux on the cell interface is simpli-fied by transformation of distribution function coupling with collision term. Half time stepadvection is used to replace the original one, which provides DUGKS semi-implicit property.Different from particle-based method [17, 18], multi-scale property in DUGKS is embodiedin temporal/spatial reconstruction in flux evaluation. Using particle velocity characteristicline, particle transport and collision processes are accurately traced. Gas nature is also de-scribed with high fidelity. In the recent study, DUGKS has been widely extended. DUGKShas been applied to rarefied gas flow, X. Zhao, et al. proposed a reduced order modeling-based DUGKS, which apply the reduced velocity space to the conventional DUGKS [19].Recently DUGKS has been extended to binary gas mixtures of Maxwell molecules, and Y.Zhang, et al. extended DUGKS to gas mixture flows based on the McCormack model [20].Many works have been done to develop efficient, accurate and robust methods based onDUGKS. P. Wang, et al developed a coupled DUGKS [21]. By implementation of kineticboundary condition and solving velocity and temperature independently, convection flowsfrom laminar to turbulent flows are accurately simulated. C. Wu, et al. retained particle ac-celeration term in original Boltzmann equation and considered force term in non-equilibriumgas distribution function [22]. Good accuracy and efficiency are validated in force-drivenflow cases. C. Shu, et al. constructed a third-order DUGKS [23]. Using Runge-Kutta tem-poral marching and high-order spatial reconstruction, more accurate results can be obtained4ith lower-resolution mesh. C. Zhang, et al. introduced a finite volume discretization of theanisotropic gray radiation equation based on DUGKS [24]. Radiative transfer in anisotropicscattering media was accurately simulated. J. Chen, et al. developed a conserved DUGKSfor multi-scale problems [25]. In the unstructured particle velocity space, macroscopic vari-ables are updated by moments of mesoscopic gas distribution function flux. Conservationproperty is greatly improved. Based on DUGKS solver, Z. Yang, et al. applied phase-fieldmethod to simulation of two-phase flows [26]. Several cases in which interfaces undergo se-vere deformation are efficiently simulated and many delicate details are accurately captured.D. Pan, et al. proposed an implicit DUGKS based on Lower-Upper Symmetric Gauss Seidel(LU-SGS) iteration [27]. In all flow regimes, computational efficiency can be improved byone or two orders of magnitude in comparison with explicit method.In present work, a kinetic flux reconstruction strategy is introduced in modeling flow evo-lution. Macroscopic flux is applied to flow evolution for improving macroscopic conservation[28]. Gas distribution function at next time step is applied to reconstruct particle motionin a mesh cell [29]. Different from original DUGKS, collision term is not implicitly includedin evolutionary process. Instead, the Bhatnagar-Gross-Krook (BGK) operator is integratedin time and space using Maxwellian and gas distribution function at next time step. Gasdistribution function is predicted by temporal and spatial reconstruction on particle velocitycharacteristic lines. Maxwellian is predicted by macroscopic variables which are updated bymoments of mesoscopic gas distribution function flux. To extend present method to unstruc-tured meshes, accurate and robust interpolation method is crucial. In previous gas kineticmodel, linear least squares regression (LLSR) was proved to be an efficient and accurate wayfor interpolation [30–32]. In present work, LLSR is applied in spatial interpolation. Thetime step in this flux reconstruction strategy is only constrained by Courant-Friedrichs-Lewy(CFL) condition which is more than ten times larger than relaxation time, which is moreefficient and robust than FVLBM. By using explicit discretization of collision term, presentmethod is simpler than original DUGKS.Our paper is organized as follows. Discretization of governing equation, reconstructionmethod and evaluation of both mesoscopic and macroscopic fluxes are presented in Sec. II.Boundary conditions we used are presented in Sec. III. In Sec. IV, several test cases, includingCouette flow, lid-driven cavity flow, and the flows over a flat plate, circular cylinder, andthe National Advisory Committee for Aeronautics (NACA) airfoil (NACA 0012 airfoil) are5onducted to validate present method. Micro cavity flows are carried out to verify the all flowregimes simulation ability of present method. Some remarks and discussions are concludedin Sec. V.
II. SIMPLIFIED DISCRETE UNIFIED GAS KINETIC SCHEME
The starting point of the DUGKS is the Boltzmann equation with BGK collision model.The Boltzmann BGK equation reads ∂f∂t + ξ · ∇ f = Ω = − f − f eq τ , (1)where f = f ( x , ξ , t ) is the gas distribution function for particles moving with velocity ξ at position x on time t , τ is the relaxation time, and f eq is the Maxwellian distributionfunction, f eq = ρ (2 πRT ) D/ exp (cid:32) − | ξ − u | RT (cid:33) , (2)where ρ is the density, R is the gas constant, T is the temperature, D is the spatial dimension, u is the flow velocity. W ( x , t ) = ρρ u ρE = (cid:90) ψ ( ξ ) f ( x , ξ , t ) d ξ , (3)where ρE = ρu / ρe is the total energy, e is the internal energy per unit mass, and ψ = (1 , ξ , ξ / T is the collision invariant. Eq. (1) can be discretized with a semi-implicitform in which flux and collision terms are predicted with kinetic reconstruction [33].By integrating Eq. (1) on the control volume V j from time t n to t n +1 = t n + ∆ t , discreteequation for solving gas distribution function f ( x , ξ , t ) at every time step can be written as f jn +1 − f jn + ∆ tV j F n +1 j,meso = ∆ t (cid:32) − f jn +1 − f jeq, ( n +1) τ j (cid:33) , (4)where F n +1 j,meso is the flux of distribution function f n +1j across every cell interfaces of the controlvolume V j at time t n +1 . Flux term F n +1 j,meso is applied to simplify the evaluation of distributionfunction which is one of the differences between our method and original DUGKS.6y simple derivation, update scheme for gas distribution function on next time step f jn +1 is easily written as (cid:18)
1+ ∆ tτ j (cid:19) f jn +1 = f jn − ∆ tV j F n +1 j,meso + ∆ tτ j f jeq, ( n +1) , (5) f jn +1 = (cid:20) f jn − ∆ tV j F n +1 j,meso + ∆ tτ j f jeq, ( n +1) (cid:21)(cid:30)(cid:18)
1+ ∆ tτ j (cid:19) . (6)The flux F n +1 j,meso and equilibrium distribution function f jeq, ( n +1) of the cell at time t n +1 areunknown. For better description of particle motion, Boltzmann BGK model equation isdiscretized and integrated on particle velocity characteristic line within one time step at thecell interface x b . Within a time step, trapezoidal rule is applied to integrate Eq. (1) f ( x b , ξ , t n + ∆ t ) − f ( x b − ξ ∆ t, ξ , t n ) = ∆ t x b , ξ , t n + ∆ t ) + Ω ( x b − ξ ∆ t, ξ , t n )] . (7)Then we move the time t n +1 term to the left hand side as well as moving the term on time t n to the right hand side. Two new transformed distribution functions are obtained¯ f j ( x b , ξ , t n +1 ) = ¯ f + j ( x b − ξ ∆ t, ξ , t n ) , (8)in which, ¯ f j ( x b , ξ , t n +1 ) = f j ( x b , ξ , t n +1 ) − ∆ t x b , ξ , t n + ∆ t )= 2 τ + ∆ t τ f j ( x b , ξ , t n +1 ) − ∆ t τ f eq,n +1 j ( x b , ξ , t n +1 ) , (9)¯ f + j ( x b − ξ ∆ t, ξ , t n ) = f j ( x b − ξ ∆ t, ξ , t n ) + ∆ t x b − ξ ∆ t, ξ , t n )= 2 τ − ∆ t τ f j ( x b − ξ ∆ t, ξ , t n ) + ∆ t τ f eq,nj ( x b − ξ ∆ t, ξ , t n ) . (10)¯ f j at next time step can be spatially reconstructed by ¯ f + j at current time step. With theTaylor expansion around the cell center, for smooth flow regime, ¯ f + j ( x b − ξ ∆ t, ξ , t n ) can beapproximated as¯ f + j ( x b − ξ ∆ t, ξ , t n ) = ¯ f + j ( x j , ξ , t n ) + ( x b − ξ ∆ t − x j ) · σ j . (11)To extend present method to more complex configuration, spatial derivative vector σ j isevaluated by LLSR which is more efficient on unstructured mesh. The fitting formula takesfollowing form ¯ f + j ( x , ξ , t n ) = ¯ f + ( x j , ξ , t n ) + ( x − x j ) · σ j . (12)7or 2D cases, with mesh cells adjacent to cell j , the linear least squares regression equationis constructed as N f (cid:88) k =1 (cid:16) x (1) k − x (1) j (cid:17) (cid:16) x (1) k − x (1) j (cid:17) (cid:16) x (2) k − x (2) j (cid:17)(cid:16) x (1) k − x (1) j (cid:17) (cid:16) x (2) k − x (2) j (cid:17) (cid:16) x (2) k − x (2) j (cid:17) σ j = N f (cid:88) k =1 (cid:16) x (1) k − x (1) j (cid:17) (cid:0) ¯ f + k − ¯ f + j (cid:1)(cid:16) x (2) k − x (2) j (cid:17) (cid:0) ¯ f + k − ¯ f + j (cid:1) , (13)where k represent cells around cell j and N f is the number of them. x (1) and x (2) are twocomponents of space vector x .The distribution function f j ( x b , ξ , t n ) at current time step is adopted to calculate the fluxof macro variables by following formula F nj,m acr o = (cid:73) ∂V j (cid:90) ψ ( ξ · n ) f j ( x b , ξ , t n ) d ξ dl, (14)where dl is circular integration variable around the interface of ∂V j .Macroscopic equations can be derived from moment of mesoscopic Boltzmann equation [34].By integration of Eq. (4), update scheme of macro variables is written as W jn +1 − W jn + ∆ tV j F nj,maro = . (15)After that, equilibrium distribution f jeq, ( n +1) on t n +1 is updated by macro variables W jn +1 .According to Eq. (9), distribution function f ( x b , ξ , t n +1 ) at the center of cell interface iswritten as. f ( x b , ξ , t n +1 ) = (cid:18) ¯ f j ( x b , ξ , t n +1 ) + ∆ t τ f eq,n +1 j ( x b , ξ , t n +1 ) (cid:19) · τ τ + ∆ t . (16)With distribution function f ( x b , ξ , t n +1 ), the flux F n +1 j,meso at the cell interface can be obtained. F n +1 j,meso = (cid:73) ∂V j ( ξ · n ) f ( x b , ξ , t n +1 ) dl. (17)Finally, gas distribution function at the cell center on t n +1 can be updated by Eq. (6) sinceall the unknown terms are obtained.In summary, the update rule of distribution function f from t n to t n + 1 in our simplifieddiscrete unified gas kinetic scheme (SDUGKS) is the following:1. F nj,macro is calculated by Eq. (14).2. Using Eq. (10) to obtain the ¯ f + at the cell center.8. With Eq. (11), we can obtain ¯ f + j ( x b − ξ ∆ t, ξ , t n ) by LLSR.4. According to Eq. (8), ¯ f j ( x b , ξ , t n +1 ) is obtained.5. Then the macro variables at cell interface on next time step can be updated using¯ f j ( x b , ξ , t n +1 ) by Eq. (3) based on conservation constraint.6. After the macro variables obtained, the equilibrium distribution function on the mid-point of interface at t n +1 can be updated by Eq. (2)7. According to Eq. (16), using the ¯ f j ( x b , ξ , t n +1 ) and f jeq,n +1 to get distribution function f ( x b , ξ , t n +1 ) at the center of cell interface.8. Then the F n +1 j,meso can be obtained by Eq. (17).9. Then update the macro variables at center of the cell on next time step by Eq. (15).10. Update the equilibrium distribution function f eq ( x j , ξ , t n +1 ) at the cell center byEq. (2).11. All the unknown values in Eq. (6) are obtained, we can update the distribution function f n +1 j at the cell center on next time step. III. BOUNDARY CONDITIONS
Boundary condition is decided by direction of particle velocity. Boundary outward normalvector is N bcf and tangential vector is A bcf . Sign of dot product between ξ and N bcf isapplied to distinguish particle information. If ξ · N bcf <
0, the implementation is as following.
A. Nonslip wall boundary condition ρ ( x bcf , t n ) = ρ ( x in , t n ) , u ( x bcf , t n ) = u wall ,f ( x bcf , ξ , t n ) = f eq ( x bcf , ξ , t n ) + ( f ( x in , ξ , t n ) − f eq ( x in , ξ , t n )) . (18)9onslip wall boundary condition is implemented based on non-equilibrium extrapolationmethod [35]. In above formulas bcf denotes interface on boundary and in denotes its neigh-bor cell. Maxwellian on boundary is computed by macroscopic variables of solid wall. B. Inlet flow boundary condition ρ ( x bcf , t n ) = ρ ∞ , u ( x bcf , t n ) = u ∞ ,f ( x bcf , ξ , t n ) = f eq ( x bcf , ξ , t n ) + ( f ( x in , ξ , t n ) − f eq ( x in , ξ , t n )) . (19)where ∞ represents freestream. C. Outlet flow ρ ( x bcf , t n ) = ρ ( x in , t n ) , u ( x bcf , t n ) = u ( x in , t n ) ,f ( x bcf , ξ , t n ) = f eq ( x bcf , ξ , t n ) + ( f ( x in , ξ , t n ) − f eq ( x in , ξ , t n )) . (20) D. Symmetric boundary condition
For consistency and conservation on symmetric boundary, ghost cells symmetric to innercells about boundaries are applied. Density in ghost cells should be equal to inner cells,which reads ρ ( x ghost , t n ) = ρ ( x in , t n ) . (21)Velocity vector in ghost cells is presented as mirror-symmetry, which is shown in Fig. 1.Relation of macroscopic velocity between ghost cells and inner cells is written as u ghost · N bcf = − u in · N bcf , u ghost · A bcf = u in · A bcf . (22)Macroscopic variables on boundaries are interpolated using inner cells and ghost cells aroundboundaries. Gas distribution function on symmetric boundaries is reconstructed by LLSR.10IG. 1: Velocity vector relation of symmetric boundary condition. E. Periodic boundary condition
As gas flows away from computational field, it will flow into the same field from theopposite side simultaneously. As shown in Fig. 2, cells A and B are ghost cells correspondingto inner cells A (cid:48) and B (cid:48) . The relation between them reads ρ ( x A , t n ) = ρ ( x A (cid:48) , t n ) , u ( x A , t n ) = u ( x A (cid:48) , t n ) ,f ( x A , ξ , t n ) = f ( x A (cid:48) , ξ , t n ) . (23) ρ ( x B , t n ) = ρ ( x B (cid:48) , t n ) , u ( x B , t n ) = u ( x B (cid:48) , t n ) ,f ( x B , ξ , t n ) = f ( x B (cid:48) , ξ , t n ) . (24)Mesoscopic and macroscopic variables on periodic boundaries are equal to the values on theinner interface of the corresponding periodic boundary.If ξ · N bcf ≥
0, variables on boundaries should be interpolated from inner flow field.In this way correct physical information is well preserved in reconstruction on the bound-ary. Consequently, evolution on boundaries is compatible with flow field, which increasesrobustness in computation. 11IG. 2: Implementation of periodic boundary condition.
F. Diffuse boundary condition
In rarefied flow simulation, diffuse-scattering rule is implemented on solid wall [36]. Gasdistribution function is computed by f ( x bcf , ξ , t ) = f eq ( x bcf , ξ , t ) (25)Density on wall is computed by [16] ρ ( x bcf , t ) = − (cid:88) ξ · N bcf > ( ξ · N bcf ) f eq ( x bcf , ξ , t ) − (cid:88) ξ · N bcf < ( ξ · N bcf ) f ( x bcf , ξ , t ) (26) IV. NUMERICAL RESULTS
In this section, several cases are conducted to validate the proposed SDUGKS. Thefirst Couette flow case verifies its spatial second-order accuracy. Lid-driven cavity flow iscomputed to validate present method in viscous flow simulation. On unstructured meshes,flows around NACA 0012 airfoil are accurately simulated. In unsteady cases, flows aroundcircular cylinder are simulated and many flow details are captured. The same accuracy withprevious work is proved with much larger time step and robustness is also pretty well. Micro12IG. 3: Meshes of Couette flow for spatial accuracy test. (a) 4 ×
8, (b) 4 ×
16, (c) 4 × ×
64, (e) 4 ×
1. Couette flow
Steady Couette flow case is computed for spatial accuracy validation. In this case, thegeometry details are introduced as follow. Unstructured meshes with different resolutionsof 4 ×
8, 4 ×
16, 4 ×
32, 4 ×
64, 4 ×
128 are used in this case, which are shown in Fig. 3.The height of computational domain is 1 .
0, the width of mesh is depending on the numberof cells it used.Nonslip wall boundary condition is implemented on the top and bottom walls, while thetop wall moves with constant velocity u w = 1 .
0, and the bottom wall is a static wall with u w = 0 .
0. The left and right boundaries are prescribed as periodic boundary condition.The analytical solution of velocity distribution in Couette flow field is shown as follow: u = u w y,v = 0 . . (27)Velocity profiles are plotted in Fig. 4 for comparison with analytical solution. In Fig. 5,second-order accuracy is proved by L2-norm of velocity errors varying with mesh resolution.
2. Laminar flow over a flat plate
As shown in Fig. ?? , In this case, computational domain is set to be x ∈ [ − , y ∈ [0 , x ∈ [ − ,
0] in set as symmetric condition, and in x ∈ [0 , u w is set to 0 . Re ∞ = 10 , ρ = 1 . ν = U ∞ L/Re ∞ , where L = 100 is the length of the flat plate.The boundary layer velocity distributions at position x = 5, x = 10, x = 20 and compar-isons with Blasius solution are shown in Fig. 7. The skin friction coefficients C f distributionalong the flat plate is shown in Fig. 8. Drag coefficient C d = 4 . × − fits the Blasiussolution well, which is 4 . × − . 14 a) (b) FIG. 6: Mesh for flat plate. (a) Mesh and computational domain (b) Mesh near leadingedge of flat plate. (a) (b)
FIG. 7: Comparison results for velocity profiles of flat plate flow. (a) U x and (b) U y .FIG. 8: Comparison result for skin friction coefficients of flat plate.15 . Lid-driven cavity flow Incompressible lid-driven cavity flow is a benchmark case for validation of numericalscheme in viscous flow simulation. Fluid in a square cavity is driven by a top moving wallwith other three static walls keeping. The top wall is moving with a fix speed u w . Thisincompressible case is only featured by Reynolds number. In our simulation, in order toguarantee a nearly incompressible flow, the driven velocity is set to be u w = 0 .
1, density ρ is set to be ρ = 1 .
0, the length of the cavity is taken to be the number of cell in x direction.The Reynolds number is defined as Re = u w Lν . Time step in this case is set to be ten timesof collision time τ . Convergence is decided by residual of the horizontal velocity U andvertical velocity V . The convergence criterion is 10 − which is being decided in every 1000iteration steps. (cid:15) = (cid:114)(cid:80) i (cid:104)(cid:0) U n +1000 i − U ni (cid:1) + (cid:0) V n +1000 i − V ni (cid:1) (cid:105)(cid:113)(cid:80) i (cid:0) ( U ni ) + ( V ni ) (cid:1) (28)The computational domain is discretized into 48 ×
48 cells. For comparison, The Ghia’sbenchmark solutions are included [37], in which 128 ×
128 mesh is adopted. It shows inFig. 9 that SDUGKS fits the benchmark data well. The streamline of the lid-driven cavityflow is shown in Fig. 10.
4. Flows around circular cylinder
The flow past a single stationary circular cylinder is a good test case to demonstratepresent SDUGKS in capturing boundary layer, detached structures and wake flows. Com-putational domain is set as a 75 d × d square zone with d -diameter cylinder at the center.The geometry and mesh of computational domain are shown in Fig. 12. In order to improvethe resolution of the tail vortex area and obtain more details of the flow field, the grid isrefined. The number of the mesh cells in Re ≤
47 cases is 17658, and the number of thecells in Re ≥
60 cases is 40029. The minimum mesh spacing for both of the mesh is 0 . d .On circular cylinder surface, nonslip wall boundary condition is implemented.In order to approximate the incompressible flow, a constant velocity u ∞ is specified at far-16 a) (b) FIG. 9: Comparison results of lid-driven cavity flow for u-velocity along the centralvertical line and v-velocity along the central horizontal line. (a) Re = 400 and (b) Re = 1000. (a) (b) FIG. 10: Streamline of the lid-driven cavity flow at (a) Re = 400 and (b) Re = 1000.field boundaries. u ∞ is set to be 0 . v ∞ is set to be 0 .
0, density ρ = 1 .
0. The Reynoldsnumber is defined as Re = ρU ∞ dµ , where Re = 10 , , , , , , ,
100 respectively.17 a) (b)
FIG. 11: (a) Pressure coefficient distributions and (b) skin friction coefficient distributionsfor flow around a NACA0012 airfoil at Ma = 0 .
8, Re = 500 and α = 10 ◦ .The drag coefficient and lift coefficient are defined as C d = F d . ρu ∞ d ,C l = F l . ρu ∞ d , (29)where F d is drag force, F l is lift force.The evolution of C d and C l at different Reynolds numbers is shown in Fig. 15. From theFig. 15, we can get that flow is steady when Re ≤
40, and C d can converge to a constant,and C l can also converge to a constant near zero. And the flow is unsteady when Re > C d and C l show sinusoidal fluctuation. In order to estimate the critical Reynolds number Re cr for the flow, we test the flow at Re = 45, and 47. C l shows slowly convergence while Re = 45, and for Re = 47, C l will develop into sinusoidal fluctuation, so 45 < Re cr < L is defined as the distance between the base point of thecylinder and the stagnation point on the center axis downstream of the cylinder, which isshown in Fig. 14.To validate temporal property of present method, time-evolution parameters are also18 a) (b)(c) (d) FIG. 12: Hybrid meshes for the simulation of laminar flow around circular cylinder: (a)Far field for Re = 10 , ,
40, and (b) mesh near wall of circular cylinder for Re = 10 , , Re = 60 , , Re = 60 , , C d and length of separate vortex L forflow around circular cylinder at Re = 10 , , Re Tritton [38] Park, et al. [39] Wang, et al. [12] Li, et al. [11] He, et al. [40] Present¯ C d ¯ C d L/R ¯ C d L/R ¯ C d L/R ¯ C d L/R ¯ C d L/R
10 2.926 2.78 0.476 2.88 0.505 3.003 0.6649 3.170 0.474 3.015 0.49320 2.103 2.01 1.814 2.072 1.866 2.118 2.0376 2.152 1.842 2.115 1.84240 1.605 1.51 4.502 1.545 4.609 1.568 4.7027 1.499 4.490 1.557 4.490 a) (b) FIG. 13: Evolution process of the lift coefficient for flows around circular cylinder at(a)Re=45, and (b) Re=47.FIG. 14: Definition of the separate vortex length L .compared to benchmark results. The Strouhal number ( St ) is a dimensionless number forunsteady temporal variation, which is denoted by St = f dU ∞ , (30)where f is the frequency of vortex shedding, d is the diameter of cylinder. The present St fits20 a) (b)(c) (d) FIG. 15: The evolutions at different Reynolds numbers of (a) (b) drag coefficient and (c)(d) lift coefficient.well with the experimental fitting curve for the parallel vortex shedding model in referencepaper [41], and the experimental fitting curve of St is defined as St = − . /Re +0 . . × − Re . The comparison of C d , St , and L/R are shown in Tab. I (for steady flows)and Tab. II (for unsteady flows).Finally, mean flow pressure distribution along circular cylinder surface is extracted andcompared with previous work. Pressure coefficient is defined as C p = p − p ∞ . ρ ∞ U ∞ , (31) C p distribution around the cylinder surface and its comparison with the reference data [39,42] is shown in Fig. 17. And it shows that all the results are in good agreement with the21 a) (b)(c) (d) FIG. 16: The streamlines at different Reynolds numbers: (a) Re = 10, (b) Re = 20,(c) Re = 40, (d) Re = 60.TABLE II: Comparison of mean drag coefficient ¯ C d , length of separate vortex L , andStrouhal number St for flows around circular cylinder at Re = 60 , , Re Tritton [38] Park, et al. [39] Wang, et al. [12] Present¯ C d ¯ C d L/R St ¯ C d L/R St ¯ C d L/R St
60 1.398 1.39 4.132 0.1353 1.422 4.155 0.1375 1.4228 4.2008 0.137780 1.316 1.35 3.312 0.1528 1.379 3.306 0.1550 1.375 3.314 0.1523100 1.271 1.33 2.782 0.1646 1.358 2.796 0.1670 1.349 2.832 0.1653 references’ data. 22 a) (b)
FIG. 17: Comparison results of mean pressure coefficients for flow around circular cylinderat (a) Re = 10 , ,
40 and (b) Re = 60 , , Re = 100.
5. Flows around NACA 0012 airfoil
To extend present SDUGKS to extensive applications, flows around a NACA 0012 airfoilare simulated on hybrid meshes. The geometry and meshes of computational domain areshown in Fig. ?? . Simulations were conducted with the angles of attack AOA = 0 ◦ and AOA = 8 ◦ . The total cells of the mesh is 38910, and minimum mesh spacing near thewall of airfoil is 4 . × − . The geometry and mesh of computational domain are shownin Fig. ?? . In airfoil case, meshes around wall and wake are refined to capture shear flowstructures. In farfield, rather coarse mesh is distributed for less computational cost. Differentresolutions of mesh are smoothly jointed by quadrilateral/triangle hybrid mesh. Rather23 a) (b) FIG. 19: Hybrid mesh for the simulation of laminar flow around a NACA 0012 airfoil: (a)Far field, and (b) mesh near wall of airfoil.good accuracy is obtained with rather few mesh cells. Boundary layer and detached trailingvortex are captured. On hybrid meshes, present method is accurate, robust and efficient incomputational fluid dynamics (CFD) research and fluids engineering.For the initial conditions, The u ∞ is set to be 0 .
1. The chord length of the airfoil isset to be L = 1 .
0, Reynolds number Re is set to be 500. The density of fluid is set as ρ = 1 .
0, and the velocity components are prescribed as u = U ∞ cos ( AOA · π/ , v = U ∞ sin ( AOA · π/ AOA = 0 ◦ , the drag and lift coefficient have been calculated,these calculation results are consistent with GILBM [43] and CFL3D, which are shown inTab. III. The pressure coefficient, pressure coefficient contour, and streamline around NACA0012 airfoil are shown in Fig. 21. The velocity profiles at different cross sections are shownin Fig. 20, which fits the CFL3D calculation results well in every cross section.Similar to AOA = 0 ◦ condition, we simulate AOA = 8 ◦ case. The drag and lift coefficientand its comparison with Fluent are shown in Tab. IV. The pressure coefficient C p , pressurecoefficient contour, and streamline around NACA 0012 airfoil are shown in Fig. 22. Thevelocity profiles at different cross section are shown in Fig. 23.
6. Micro cavity flow
The present SDUGKS is applied to micro cavity flow compared with simulations basedon Direct Simulation Monte Carlo (DSMC) Method. In the simulation, the Knudsen num-bers are set as Kn = 0 . , , ,
8. The velocity of top wall u w = 0 .
1. The velocity space is24 a)Cross sections (b) x = 0(c) x = 0 .
25 (d) x = 0 . x = 0 .
75 (f) x = 1 . FIG. 20: Comparison results for velocity profiles of laminar flow around a NACA 0012airfoil at
AOA = 0 ◦ at various cross sections x .25ABLE III: Drag and lift coefficient at AOA = 0 ◦ compared with GILBM and CFL3D Resolution (on airfoil) ∆ x min C d C l GILBM16705 (173) 4 . E −
03 0.1682 1 . E − . E −
04 0.1736 1 . E − . E −
03 0.1672 1 . E − . E −
04 0.1725 1 . E − . E −
04 0.1741 5 . E − . E − . E − . E − (a) (b) FIG. 21: (a) The pressure coefficient distribution over NACA 0012 airfoil at
AOA = 0 ◦ .(b) Pressure coefficient contour and streamline around NACA 0012 airfoil at AOA = 0 ◦ .discretized into 100 ×
100 nodes, and distribute uniformly in [ − , × [ − , ×
40 uniform mesh cells. The length of cavity is L = 1 . a) (b) FIG. 22: (a) The pressure coefficient distribution over NACA 0012 airfoil at
AOA = 8 ◦ .(b) Pressure coefficient contour and streamline around NACA 0012 airfoil at AOA = 8 ◦ .TABLE IV: Drag and lift coefficient at AOA = 8 ◦ compared with Fluent Resolution (on airfoil) ∆ x min C d C l Fluent38910(325) 4 . E − . E − The velocity profiles across the cavity center for different Knudsen numbers are shown inFig. 24, which shows that the present method fits the DSMC results well in all cases. Presentmethod is verified to be an accurate solver for all flow regimes. Its multi-scale property alsomakes it a reliable tool for all Knudsen number flows study.
V. CONCLUSION
In the present study, we propose a SDUGKS. In the spatial interpolation aspect, anaccurate and robust method is expected. And in the previous study, LLSR was proved tobe an efficient and accurate method for interpolation in reconstruction step. So, in this27 a)Cross sections (b) x = 0(c) x = 0 .
25 (d) x = 0 . x = 0 .
75 (f) x = 1 . FIG. 23: Comparison results for velocity profiles of laminar flow around a NACA 0012airfoil at
AOA = 8 ◦ at various cross sections x .28 a) (b)(c) (d) FIG. 24: Comparison results of micro cavity flow case for u-velocity along the centralvertical line and v-velocity along the central horizontal line at (a) Kn = 0 .
1, and(b) Kn = 1, and (c) Kn = 2, and (d) Kn = 8.work, LLSR is still adopted for spatial interpolation on unstructured mesh. In the proposedscheme, the flux at half time step in the original DUGKS is omitted. The mesoscopicflux of the cell on next time step is predicted by reconstruction of transformed distributionfunction at interfaces along particle velocity characteristic lines. Macroscopic variables ofthe cell on next time step are updated by the flux of gas distribution function at interfacesof the cell on current time step according to the conservation law. Conservation is improvedand larger time step can be used. The particle nature is retained by reconstruction basedon particle velocity characteristic line. The SDUGKS is validated by several test cases.29he Couette flow is simulated using unstructured mesh to verify that our new scheme hassecond order spatial accuracy. The lid-driven cavity flow, laminar flow over the flat plate,flows over a circular cylinder and a NACA 0012 airfoil are carried out to demonstrate theability of our proposed method to simulate complex flows. Micro cavity flow shows thatpresent method is applicable to rarefied flow simulations. By introducing prediction of gasdistribution function of interfaces on next time step in flux evaluation, particle motion istraced by kinetic temporal/spatial reconstruction. Physical evolution process is replaced bynumerical marching in a mesh cell within a whole time step. As a result, the stability ispreserved with the same accuracy as original DUGKS. A framework rather easy to learnand less complicated to code will also extend its application. ACKNOWLEDGEMENTS
This work is supported by National Numerical Wind tunnel project, the National NaturalScience Foundation of China (No. 11902264), the Natural Science Basic Research Planin Shaanxi Province of China (Program No. 2019JQ-315), and the 111 Project of China(B17037).
DATA AVAILABILITY
The data that support the findings of this study are available from the correspondingauthor upon reasonable request. 30
1] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of FluidMechanics 30 (1) (1998) 329–364.[2] Z. Guo, C. Shu, Lattice Boltzmann method and its applications in engineering, World Scien-tific, 2013.[3] X. He, L. S. Luo, M. Dembo, Some progress in lattice Boltzmann method. Part I. Nonuniformmesh grids, Journal of Computational Physics 129 (2) (1996) 357–363.[4] J. Yao, C. Zhong, K. Tang, An adaptive-gridding lattice Boltzmann method with linked-listdata structure for two-dimensional viscous flows, Progress in Computational Fluid DynamicsAn International Journal 17 (5) (2017) 267.[5] X. Guo, J. Y. Yao, C. Zhong, J. C. Cao, A hybrid adaptive-gridding immersed-boundarylattice Boltzmann method for viscous flow simulations, Applied Mathematics & Computation267 (2015) 529–553.[6] O. Aursjø, E. Jettestuen, J. L. Vinningland, A. Hiorth, On the inclusion of mass source termsin a single-relaxation-time lattice Boltzmann method, Physics of Fluids 30 (5) (2018) 057104.[7] D. Wang, D. Tan, N. Phan-Thien, A lattice Boltzmann method for simulating viscoelasticdrops, Physics of Fluids 31 (7) (2019) 073101.[8] F. Nannelli, S. Succi, The lattice Boltzmann equation on irregular lattices, Journal of Statis-tical Physics 68 (3-4) (1992) 401–407.[9] H. Xi, G. Peng, S.-H. Chou, Finite-volume lattice Boltzmann method, Physical Review E59 (5) (1999) 6202–6205.[10] M. Stiebler, J. T¨olke, M. Krafczyk, An upwind discretization scheme for the finite volumelattice Boltzmann method, Computers & Fluids 35 (2006) 814–819.[11] W. Li, L.-S. Luo, Finite volume lattice Boltzmann method for nearly incompressible flowson arbitrary unstructured meshes, Communications in Computational Physics 20 (02) (2016)301–324.[12] Y. Wang, C. Zhong, J. Cao, C. Zhuo, S. Liu, A simplified finite volume lattice Boltzmannmethod for simulations of fluid flows from laminar to turbulent regime, Part I: Numericalframework and its application to laminar flow simulation, Computers & Mathematics withApplications 79 (5) (2020) 1590–1618.
13] Y. Wang, C. Zhong, J. Cao, C. Zhuo, S. Liu, A simplified finite volume lattice Boltzmannmethod for simulations of fluid flows from laminar to turbulent regime, Part II: Extensiontowards turbulent flow simulation, Computers & Mathematics with Applications 79 (8) (2020)2133–2152.[14] L. Chen, L. Schaefer, A unified and preserved dirichlet boundary treatment for the cell-centered finite volume discrete Boltzmann method, Physics of Fluids 27 (2) (2015) 027104.[15] L. Zhu, Z. Guo, K. Xu, Discrete unified gas kinetic scheme on unstructured meshes, Computers& Fluids 127 (2016) 211–225.[16] Z. Guo, K. Xu, R. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows:Low-speed isothermal case, Physical Review E 88 (3) (2013) 1–11.[17] J. Zhang, B. John, M. Pfeiffer, F. Fei, D. Wen, Particle-based hybrid and multiscale methodsfor nonequilibrium gas flows, Advances in Aerodynamics 1 (1) (2019) 12.[18] M. Fang, Z.-H. Li, Z.-H. Li, J. Liang, Y.-H. Zhang, DSMC modeling of rarefied ionization re-actions and applications to hypervelocity spacecraft reentry flows, Advances in Aerodynamics2 (1) (2020) 1–25.[19] X. Zhao, C. Wu, Z. Chen, L. Yang, C. Shu, Reduced order modeling-based discrete unifiedgas kinetic scheme for rarefied gas flows, Physics of Fluids 32 (2020) 067108.[20] Y. Zhang, L. Zhu, P. Wang, Z. Guo, Discrete unified gas kinetic scheme for flows of binarygas mixture based on the McCormack model, Physics of Fluids 31 (1).[21] P. Wang, S. Tao, Z. Guo, A coupled discrete unified gas-kinetic scheme for Boussinesq flows,Computers & Fluids 120 (2015) 70–81.[22] C. Wu, B. Shi, Z. Chai, P. Wang, Discrete unified gas kinetic scheme with a force termfor incompressible fluid flows., Computers & Mathematics with Applications 71 (12) (2016)2608–2629.[23] C. Wu, B. Shi, C. Shu, Z. Chen, Third-order discrete unified gas kinetic scheme for continuumand rarefied flows: Low-speed isothermal case, Physical Review E 97 (2) (2018) 023306.[24] X. Song, C. Zhang, X. Zhou, Z. Guo, Discrete unified gas kinetic scheme for multiscaleanisotropic radiative heat transfer, Advances in Aerodynamics 2 (1) (2020) 1–15.[25] J. Chen, S. Liu, Y. Wang, C. Zhong, A conserved discrete unified gas-kinetic scheme withunstructured discrete velocity space, Physical Review E 100 (4) (2019) 043305.
26] Z. Yang, C. Zhong, C. Zhuo, Phase-field method based on discrete unified gas-kinetic schemefor large-density-ratio two-phase flows, Physical Review E 99 (4) (2019) 043302.[27] D. Pan, C. Zhong, C. Zhuo, An implicit discrete unified gas-kinetic scheme for simulationsof steady flow in all flow regimes, Communications in Computational Physics 25 (5) (2019)1469–1495.[28] L. Zhang, Z. Chen, L. Yang, C. Shu, Double distribution function-based discrete gas kineticscheme for viscous incompressible and compressible flows, Journal of Computational Physics(2020) 109428.[29] L. Yang, C. Shu, W. Yang, J. Wu, An improved three-dimensional implicit discrete velocitymethod on unstructured meshes for all Knudsen number flows, Journal of ComputationalPhysics 396 (2019) 738–760.[30] S. Lenz, M. Geier, M. Krafczyk, An explicit gas kinetic scheme algorithm on non-uniformcartesian meshes for GPGPU architectures, Computers & Fluids 186 (2019) 58–73.[31] J. Li, C. Zhong, D. Pan, C. Zhuo, A gas-kinetic scheme coupled with SST model for turbulentflows, Computers & Mathematics with Applications 78 (4) (2019) 1227–1242.[32] D. Pan, C. Zhong, L. Ji, C. Zhuo, A gas-kinetic scheme for the simulation of turbulent flows onunstructured meshes, International Journal for Numerical Methods in Fluids 82 (11) (2016)748–769.[33] R. Yuan, C. Zhong, A conservative implicit scheme for steady state solutions of diatomic gasflow in all flow regimes, Computer Physics Communications 247 (2020) 106972.[34] L. Wu, X.-J. Gu, On the accuracy of macroscopic equations for linearized rarefied gas flows,Advances in Aerodynamics 2 (1) (2020) 1–32.[35] Z.-L. Guo, C.-G. Zheng, B.-C. Shi, Non-equilibrium extrapolation method for velocity andpressure boundary conditions in the lattice Boltzmann method, Chinese Physics B 11 (4)(2002) 366–374.[36] J. Meng, Y. Zhang, Diffuse reflection boundary condition for high-order lattice Boltzmannmodels with streaming-collision mechanism, Journal of Computational Physics 258 (2014)601–612.[37] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, Journal of Computational Physics 48 (3) (1982)387–411.
38] D. J. Tritton, Experiments on the flow past a circular cylinder at low Reynolds numbers,Journal of Fluid Mechanics 6 (1959) 547–567.[39] J. Park, K. Kwon, H. Choi, Numerical solutions of flow past a circular cylinder at Reynoldsnumbers up to 160, KSME International Journal 12 (6) (1998) 1200–1205.[40] X. He, G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: Flow arounda circular cylinder, Journal of Computational Physics 134 (1997) 306–315.[41] C. H. K. Williamson, Oblique and parallel modes of the vortex shedding in the wake of acircular cylinder at low Reynolds numbers, Journal of Fluid Mechanics 206 (-1) (1989) 579–627.[42] B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, Journal ofFluid Mechanics 98 (4) (1980) 819–855.[43] T. Imamura, K. Suzuki, T. Nakamura, M. Yoshida, Flow simulation around an airfoil by latticeBoltzmann method on generalized coordinates, AIAA Journal 43 (9) (2005) 1968–1973.38] D. J. Tritton, Experiments on the flow past a circular cylinder at low Reynolds numbers,Journal of Fluid Mechanics 6 (1959) 547–567.[39] J. Park, K. Kwon, H. Choi, Numerical solutions of flow past a circular cylinder at Reynoldsnumbers up to 160, KSME International Journal 12 (6) (1998) 1200–1205.[40] X. He, G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: Flow arounda circular cylinder, Journal of Computational Physics 134 (1997) 306–315.[41] C. H. K. Williamson, Oblique and parallel modes of the vortex shedding in the wake of acircular cylinder at low Reynolds numbers, Journal of Fluid Mechanics 206 (-1) (1989) 579–627.[42] B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, Journal ofFluid Mechanics 98 (4) (1980) 819–855.[43] T. Imamura, K. Suzuki, T. Nakamura, M. Yoshida, Flow simulation around an airfoil by latticeBoltzmann method on generalized coordinates, AIAA Journal 43 (9) (2005) 1968–1973.