Two-step multi-resolution reconstruction-based compact gas-kinetic scheme on tetrahedral mesh
TTwo-step multi-resolution reconstruction-based compact gas-kineticscheme on tetrahedral mesh
Xing Ji a , Fengxiang Zhao b , Wei Shyy b , Kun Xu a,b,c, ∗ a Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay,Kowloon, Hong Kong b Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology,Clear Water Bay, Kowloon, Hong Kong c Shenzhen Research Institute, Hong Kong University of Science and Technology, Shenzhen, China
Abstract
In this paper, a third-order compact gas-kinetic scheme (GKS) on unstructured tetrahedralmesh is constructed for the compressible Euler and Navier-Stokes solutions. The time-dependent gas distribution function at a cell interface is used to calculate the fluxes for theupdating the cell-averaged flow variables and to evaluate the time accurate cell-averagedflow variables as well for evolving the cell-averaged gradients of flow variables. With theaccurate evolution model for both flow variables and their slopes, the quality of the schemedepends closely on the accuracy and reliability of the initial reconstruction of flow variables.The reconstruction scheme becomes more challenge on tetrahedral mesh, where the conven-tional second-order unlimited least-square reconstruction can make the scheme be linearlyunstable when using cell-averaged conservative variables alone with von Neumann neighbors.Benefiting from the evolved cell-averaged slopes, on tetrahedral mesh the GKS is linearlystable from a compact third-order smooth reconstruction with a large CFL number. In orderto further increase the robustness of the high-order compact GKS for capturing discontinu-ous solution, a new two-step multi-resolution weighted essentially non-oscillatory (WENO)reconstruction will be proposed. The novelty of the reconstruction includes the following.Firstly, it releases the stability issue from a second-order compact reconstruction throughthe introduction of a pre-reconstruction step. Secondly, in the third-order non-linear recon-struction, only one more large stencil is added beside those in the second-order one, whichsignificantly simplifies the high-order reconstruction. At the same time, the high-order wallboundary treatment is carefully designed by combining the constrained least-square tech-nique and the WENO procedure, where a quadratic element is adopted in the reconstructionfor the curved boundary. Numerical tests for both inviscid and viscous flow at low and highspeed are presented from the second-order and third-order compact GKS. The proposedthird-order scheme shows good robustness in high speed flow computation and favorablemesh adaptability in cases with complex geometry.
Keywords: compact gas-kinetic scheme, two-step reconstruction, multi-resolution WENO,two-stage time discretization, Navier-Stokes solution
Preprint submitted to Elsevier February 3, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] F e b . Introduction The simulation of compressible flow with complex geometry is important in the engi-neering applications. The use of unstructured mesh is especially favored from its geometricflexibility. On such a mesh even for a second-order finite volume method (FVM) it is noteasy to achieve the same performance as that in the structure mesh. Commonly, slope recon-struction schemes, such as cell-based Green-Gauss method [15] combined with different typesof limiters, are robust and widely used in the commercial or open-source software [36]. How-ever, these methods can easily deteriorate the spatial accuracy for skewed mesh and becomeover-dissipative for flow simulation with discontinuities. The least-square reconstructionwith cell-averaged variables and von Neumann neighbors only preserves a strictly second-order accuracy, but suffers from the linear instability on tetrahedral grid [5]. Therfore, thereconstruction stencil has to be further extended. In order to ensure linear stability, a two-step second-order weighted essentially non-oscillatory (WENO) method has been proposed[30] with the attempt of keeping a compact reconstruction even with an extended stencil.The high-order WENO-FVMs have been continuously developed and applied to large-scaleaeronautical simulations [1]. But, the compactness can be hardly kept in the high-orderFVMs. Even though the extended stencils in reconstruction can improve the robustnessof the schemes, difficulties still exist in the parallel programming and boundary treatment.The recently proposed multi-resolution reconstruction with only five equivalent sub-stencilsgreatly releases the above problems [37] even with the inclusion of neighbor-to-neighbor cells.The compact methods with the updates of multiple degrees of freedom (DOFs) for eachcell have been developed extensively in the past decades. Two main representatives are theDG [24] and the FR/CPR methods [6, 34], which hybridize the finite volume framework withthe finite element method or the finite difference method. These methods can achieve arbi-trary spatial order of accuracy with only the targeted cell as the reconstruction stencil, andyields great mesh adaptability and high scalability. Successful examples have been demon-strated in large eddy simulation (LES) [28] and RANS simulation [33] for subsonic flows.For the flow simulation with discontinuities, these methods have less robustness against thetraditional high-order FVMs. In addition, these methods have restricted explicit time stepsand high memory-consumption [14]. The P N P M [4] and reconstructed-DG (rDG) methods[14] were targeting to overcome the above weakness with the release of the compactness ofthe DG methods. Large time step and less memory requirement can be achieved in the rDGmethods in comparison with the same order DG ones.In recent years, a class of high-order compact GKS has been developed. The GKS isbased on a time accurate evolution model in the construction of the gas distribution functionat a cell interface [32]. The time-dependent solution provides not only the fluxes across a cellinterface, but also the corresponding flow variables. As a result, both the cell-averaged flowvariables and their slopes can be updated simultaneously through the divergence theorem. ∗ Corresponding author
Email addresses: [email protected] (Xing Ji), [email protected] (Fengxiang Zhao), [email protected] (Wei Shyy), [email protected] (Kun Xu)
2. Compact finite volume gas-kinetic scheme
The 3-D gas-kinetic BGK equation [2] is f t + u · ∇ f = g − fτ , (1)where f = f ( x , t, u , ξ ) is the gas distribution function, which is a function of space x , time t , particle velocity u , and internal variable ξ . g is the equilibrium state approached by f and τ is the collision time.The collision term satisfies the compatibility condition ˆ g − fτ ψψψ dΞ = 0 , (2)where ψψψ = (1 , u ,
12 ( u + ξ )) T , dΞ = d u d u d u d ξ ... d ξ K , K is the number of internal degreesof freedom, i.e. K = (5 − γ ) / ( γ −
1) in 3-d case, and γ is the specific heat ratio.In the continuum flow regime with the smoothness assumption, based on the Chapman-Enskog expansion of the BGK equation the gas distribution function can be expressed as[32], f = g − τ D u g + τ D u ( τ D u ) g − τ D u [ τ D u ( τ D u ) g ] + ..., where D u = ∂/∂t + u · ∇ . Different hydrodynamic equations can be derived by truncatingon different orders of τ . With the zeroth-order in truncated distribution function f = g ,the Euler equations can be recovered by multiplying ψψψ on Eq.(1) and integrating it over thephase space, W t + ∇ · F ( W ) = 0 . With the first-order truncation, i.e., f = g − τ ( u · ∇ g + g t ) , (3)4he N-S equations can be obtained, W t + ∇ · F ( W , ∇ W ) = 0 , with τ = µ/p and P r = 1.The conservative variables and their fluxes are the moments of the gas distributionfunction W ( x , t ) = ˆ ψψψf ( x , t, u , ξ )dΞ (4)and F ( x , t ) = ˆ u ψψψf ( x , t, u , ξ )dΞ . (5) Remark 1.
The cell-averaged conservative variables can be updated through the interfacefluxes under the finite volume framework. Besides the fluxes in Eq. (5) , the gas distributionfunction also provides the flow variables at the cell interface, such as that in Eq. (4) . It isthe key point for the possibility of constructing compact GKS. An obvious prerequisite is thetime accurate W ( x , t ) at a cell interface. In other words, it depends solely on the high-ordergas evolution model in the construction of high-order scheme.2.1. Compact gas-kinetic scheme on tetrahedral mesh For a tetrahedral cell Ω i in 3-D case, the boundary can be expressed as ∂ Ω i = N f (cid:91) p =1 Γ ip , where N f = 4 is the number of cell interfaces for cell Ω i .The update of the cell averaged conservative flow variables in a finite control volume ifrom t n to t n +1 can be expressed as W n +1 i | Ω i | = W ni | Ω i | − N f (cid:88) p =1 ˆ Γ ip ˆ t n +1 t n F ( x , t ) · n p d t d s, (6)with F ( x , t ) · n p = ˆ ψψψf ( x , t, u , ξ ) u · n p dΞ , (7)where W i is the cell averaged values over cell Ω i , | Ω i | is the volume of Ω i , F is the interfacefluxes, and n p = ( n , n , n ) T is the unit vector representing the outer normal direction ofΓ ip . The semi-discretized form of finite volume scheme can be written asd W i d t = L ( W i ) = − | Ω i | N f (cid:88) p =1 ˆ Γ ip F ( W ) · n p d s. (8)5or the interface fluxes F ip ( t ), the numerical quadrature can be adopted and Eq.(6) can berewritten as W n +1 i | Ω i | = W ni | Ω i | − N f (cid:88) p =1 | Γ ip | M (cid:88) k =1 ω k ˆ t n +1 t n F ( x p,k , t ) · n p d t. (9)Nowadays the curved mesh generation has been supported by popular commercial soft-ware. To be consistent with the spatial accuracy, the quadratic element is applied hereto describe the geometry. The controlling points for the quadratic triangle are shown inFig. 1(a). The iso-parametric transformation is used to evaluate the surface integral, whichis expressed as X ( ξ, η ) = (cid:88) l =0 x l φ l ( ξ, η ) , where x l is the location of the controlling points and φ l is the base function as follows [27] v = ( ξ + η − ξ + 2 η − , v = ξ (2 ξ − , v = η (2 η − ,v = − ξ ( ξ + η − , v = 4 ξη, v = − η ( ξ + η − . (10)The flux across Γ ip in Eq.(6) can be transferred to a standard isosceles right triangle ˜Γ ip F ip ( t ) = ˆ Γ ip F ( x , t ) · n p d s = ˆ ˜Γ ip F ( W ( X ( ξ, η ))) · n p (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, y, z ) ∂ ( ξ, η ) (cid:12)(cid:12)(cid:12)(cid:12) d ξ d η. To meet the requirement of a third-order spatial accuracy, the above equation can be ap-proximated through Gaussian quadrature as F ip ( t ) = 12 ∆ ξ ∆ η (cid:88) m =1 ˜ ω m F m ( t ) · ( n p ) m (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, y, z ) ∂ ( ξ, η ) (cid:12)(cid:12)(cid:12)(cid:12) m , where ∆ ξ = ∆ η = 1 and the local normal direction ( n p ) m = ( X ξ × X η ) / (cid:107) X ξ × X η (cid:107) . Thestandard Gaussian points are( ξ, η ) = ( 16 ∆ ξ,
16 ∆ η ) , ( ξ, η ) = ( 23 ∆ ξ,
16 ∆ η ) , ( ξ, η ) = ( 16 ∆ ξ,
23 ∆ η ) , with ˜ ω m = , m = 1 , ,
3. Compared with Eq.(9), we have ω m = 12 ˜ ω m (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, y, z ) ∂ ( ξ, η ) (cid:12)(cid:12)(cid:12)(cid:12) m , x p,k = x (( ξ, η )) m , n p,k = ( n p ) m , m = 1 , , . The quadratic element reduces to the linear element when every edge is a straight line.According to the coordinate transformation, the local coordinate for the cell interface6 (cid:2)(cid:3)(cid:4) (cid:5)(cid:6) (cid:1) (cid:2)(cid:3) (cid:4) (cid:5)(cid:6) (cid:7)(cid:8) (a) Quadratic triangluar surface (cid:1) (cid:2) (cid:3)(cid:4) (cid:4) (cid:3)(cid:2) (cid:1)(cid:5) (cid:6)(cid:7)(cid:8) (cid:9)(cid:10)(cid:11) (cid:12)(cid:13) (cid:8) (cid:9) (cid:10)(cid:11) (cid:12)(cid:13) (b) Quadratic tetrahedral volumeFigure 1: The controlling points and isoparametric transformation of the quadratic elements. Γ ip is expressed as ( (cid:101) x , (cid:101) x , (cid:101) x ) T = (0 , (cid:101) x , , (cid:101) x ) T , where ( (cid:101) x , , (cid:101) x ) T ∈ Γ ip , and the velocities inthe local coordinate are given by ˜ u = u n + u n + u n , ˜ u = − u n + u ( n + n n ) − u n n n , ˜ u = − u n − u n n n + u (1 − n n ) , n (cid:54) = − . (11)The macroscopic conservative flow variables in the local coordinate are expressed as (cid:102) W ( (cid:101) x , t ) = TW ( x , t ) , where T is the rotation matrix T = n n n − n n + n n − n n n − n − n n n − n n
00 0 0 0 1 , n (cid:54) = − . (12)Note that when n = −
1, Eq.(11) changes to ( ˜ u , ˜ u , ˜ u ) T = ( − u , − u , u ) T and the matrix(12) is replaced by a diagonal matrix Λ = diag(1 , − , − , , (cid:101) f ( (cid:101) x , t, (cid:101) u , ξ ) = f ( x , t, u , ξ ) and | d u | = | d (cid:101) u | , then the numerical fluxes can be transformed as F ( x , t ) = ˆ ψψψf ( x , t, u , ξ ) u · n p d u d ξ = ˆ ψψψ (cid:101) f ( (cid:101) x , t, (cid:101) u , ξ ) (cid:101) u d (cid:101) u d ξ. (13)In the computation, the fluxes are obtained firstly by taking moments of the gas distribution7unction in the local coordinates (cid:101) F ( (cid:101) x , t ) = ˆ (cid:101) ψψψ (cid:101) f ( (cid:101) x , t, (cid:101) u , ξ ) (cid:101) u d (cid:101) u d ξ, (14)where (cid:101) ψψψ = (1 , (cid:101) u ,
12 ( (cid:101) u + ξ )) T . According to Eq.(11), Eq.(13) and Eq.(14), the fluxes in theglobal coordinate can be expressed as a combination of the fluxes in the local coordinate F ( W ( x , t )) · n = T − (cid:101) F ( (cid:102) W ( (cid:101) x , t )) . (15) In order to construct the numerical fluxes at x = (0 , , T , the integral solution of BGKequation Eq.(1) is used f ( x , t, u , ξ ) = 1 τ ˆ t g ( x (cid:48) , t (cid:48) , u , ξ ) e − ( t − t (cid:48) ) /τ d t (cid:48) + e − t/τ f ( x − u t, u , ξ ) , (16)where x = x (cid:48) + u ( t − t (cid:48) ) is the particle trajectory. f is the initial gas distribution function, g is the corresponding equilibrium state in space and time. The integral solution basicallystates a physical process from the particle free transport in f in the kinetic scale to thehydrodynamic flow evolution in the integration of g term. The flow evolution at the cellinterface depends on the ratio of time step to the local particle collision time ∆ t/τ .To construct a time evolving gas distribution function at a cell interface, the followingnotations are introduced first a x i ≡ ( ∂g/∂x i ) /g = g x i /g, A ≡ ( ∂g/∂t ) /g = g t /g, where g is the equilibrium state. The variables ( a x i , A ), denoted by s , depend on particlevelocity in the form of [31] s = s j ψ j = s + s u + s u + s u + s
12 ( u + u + u + ξ ) . The initial gas distribution function in the solution (16) can be modeled as f = f l ( x , u ) H ( x ) + f r ( x , u )(1 − H ( x )) , where H ( x ) is the Heaviside function. Here f l and f r are the initial gas distributionfunctions on both sides of a cell interface, which have one to one correspondence with theinitially reconstructed macroscopic variables. The first-order Taylor expansion for the gasdistribution function in space around x = can be expressed as f k ( x ) = f kG ( ) + ∂f kG ∂x i ( ) x i = f kG ( ) + ∂f kG ∂x ( ) x + ∂f kG ∂x ( ) x + ∂f kG ∂x ( ) x , (17)8or k = l, r . According to Eq.(3), f kG has the form f kG ( ) = g k ( ) − τ ( u i g kx i ( ) + g kt ( )) , (18)where g k is the equilibrium state with the form of a Maxwell distribution. g k can be fullydetermined from the reconstructed macroscopic variables W l , W r at the left and right sidesof a cell interface ˆ ψψψg l dΞ = W l , ˆ ψψψg r dΞ = W r . (19)Substituting Eq.(17) and Eq.(18) into Eq.(16), the kinetic part for the integral solution canbe written as e − t/τ f k ( − u t, u , ξ ) = e − t/τ g k [1 − τ ( a kx i u i + A k ) − ta kx i u i ] , (20)where the coefficients a kx , ..., A k , k = l, r are defined according to the expansion of g k . Afterdetermining the kinetic part f , the equilibrium state g in the integral solution Eq.(16) canbe expanded in space and time as follows g ( x , t ) = g c ( ,
0) + ∂g c ∂x i ( , x i + ∂g c ∂t ( , t, (21)where g c is the Maxwellian equilibrium state located at an interface. Similarly, W c are themacroscopic flow variables for the determination of the equilibrium state g c ˆ ψψψg c dΞ = W c . (22)Substituting Eq.(21) into Eq.(16), the hydrodynamic part in the integral solution can bewritten as 1 τ ˆ t g ( x (cid:48) , t (cid:48) , u , ξ ) e − ( t − t (cid:48) ) /τ d t (cid:48) = C g c + C a cx i u i g c + C A c g c , (23)where the coefficients a cx i , A c are defined from the expansion of the equilibrium state g c . Thecoefficients C m , m = 1 , , C = 1 − e − t/τ , C = ( t + τ ) e − t/τ − τ, C = t − τ + τ e − t/τ . The coefficients in Eq.(20) and Eq.(23) can be determined by the spatial derivatives ofmacroscopic flow variables and the compatibility condition as follows (cid:104) a x (cid:105) = ∂ W ∂x = W x , (cid:104) a x (cid:105) = ∂ W ∂x = W x , (cid:104) a x (cid:105) = ∂ W ∂x = W x , (cid:104) A + a x u + a x u + a x u (cid:105) = 0 , (24)9here (cid:104) ... (cid:105) are the moments of a gas distribution function defined by (cid:104) ( ... ) (cid:105) = ˆ ψψψ ( ... ) g dΞ . (25)The details for the evaluation of each term from macroscopic variables can be found in [7].In smooth flow region, the collision time is determined by τ = µ/p , where µ is thedynamic viscosity coefficient and p is the pressure at the cell interface. In order to properlycapture the un-resolved shock structure, additional numerical dissipation is needed. Thephysical collision time τ in the exponential function part can be replaced by a numericalcollision time τ n . For the inviscid flow, the collision time τ n is modified as τ n = ε ∆ t + C | p l − p r p l + p r | ∆ t, where ε = 0 .
01 and C = 1. For the viscous flow, the collision time is related to the viscositycoefficient, τ n = µp + C | p l − p r p l + p r | ∆ t, where p l and p r denote the pressure on the left and right sides of the cell interface. Theinclusion of the pressure jump term is to increase the non-equilibrium transport mechanismin the flux function to mimic the physical process in the shock layer. Then substituteEq.(20) and Eq.(23) into Eq. (16) with τ and τ n , the final second-order time dependent gasdistribution function becomes f ( , t, u , ξ ) =(1 − e − t/τ n ) g c + [( t + τ ) e − t/τ n − τ ] a cx i u i g c + ( t − τ + τ e − t/τ n ) A c g c + e − t/τ n g l [1 − ( τ + t ) a lx i u i − τ A l ] H ( u )+ e − t/τ n g r [1 − ( τ + t ) a rx i u i − τ A r ](1 − H ( u )) . (26)For smooth flow, the time dependent solution in Eq. (26) can be simplified as [31], f ( , t, u , ξ ) = g c − τ ( a cx i u i + A c ) g c + A c g c t, (27)under the assumptions of g l,r = g c , a l,rx i = a cx i . The above gas-kinetic solver for smooth flowhas less numerical dissipation than that from the full GKS solver in Eq. (26). As shown in Eq. (26), a time evolution solution at a cell interface is provided by the gas-kinetic solver, which is distinguishable from the Riemann solver with a constant solution.Recall Eq.(4), the conservative variables at the Gaussian point x p,k can be updated by taking10oments ψψψ on the gas distribution function, W p,k ( t n +1 ) = ˆ ψψψf n ( x p,k , t n +1 , u , ξ )dΞ , k = 1 , ..., M. (28)Then the cell-averaged first-order derivatives within each element at t n +1 is given throughthe Gauss’s theorem, W n +1 x = ˆ V ∇ · ( W ( t n +1 ) , , V = 1∆ V ˆ ∂V (1 , , · n W ( t n +1 )d S = 1∆ V ˆ ∂V W ( t n +1 ) n d S = 1∆ V N f (cid:88) p =1 M (cid:88) k =1 ω p,k W n +1 p,k ( n ) p,k ∆ S p ,W n +1 y = ˆ V ∇ · (0 , W ( t n +1 ) , V = 1∆ V ˆ ∂V (0 , , · n W ( t n +1 )d S = 1∆ V ˆ ∂V W ( t n +1 ) n d S = 1∆ V N f (cid:88) p =1 M (cid:88) k =1 ω p,k W n +1 p,k ( n ) p,k ∆ S p ,W n +1 z = ˆ V ∇ · (0 , , W ( t n +1 ))d V = 1∆ V ˆ ∂V (0 , , · n W ( t n +1 )d S = 1∆ V ˆ ∂V W ( t n +1 ) n d S = 1∆ V N f (cid:88) p =1 M (cid:88) k =1 ω p,k W n +1 p,k ( n ) p,k ∆ S p , (29)where n p,k = (( n ) p,k , ( n ) p,k , ( n ) p,k ) is the outer unit normal direction at each Gaussianpoint x p,k .
3. Two-stage temporal discretization
The two-stage fourth-order (S2O4) temporal discretization is adopted here as that inthe previous compact GKSs [9, 35, 10]. Following the definition of Eq.(8), a fourth-ordertime-accurate solution for cell-averaged conservative flow variables W i are updated by W ∗ i = W ni + 12 ∆ t L ( W ni ) + 18 ∆ t ∂∂t L ( W ni ) , W n +1 i = W ni + ∆ t L ( W ni ) + 16 ∆ t (cid:0) ∂∂t L ( W ni ) + 2 ∂∂t L ( W ∗ i ) (cid:1) , (30)11here L ( W ni ) and ∂∂t L ( W ni ) are given by L ( W ni ) = − | Ω i | N f (cid:88) p =1 M (cid:88) k =1 ω p,k F ( x p,k , t n ) · n p,k ,∂∂t L ( W ni ) = − | Ω i | N f (cid:88) p =1 M (cid:88) k =1 ω p,k ∂ t F ( x p,k , t n ) · n p,k ,∂∂t L ( W ∗ i ) = − | Ω i | N f (cid:88) p =1 M (cid:88) k =1 ω p,k ∂ t F ( x p,k , t ∗ ) · n p,k . (31)The proof for the fourth-order accuracy in time is shown in [11].In order to obtain the numerical fluxes F p,k and their time derivatives ∂ t F p,k at t n and t ∗ = t n + ∆ t/
2, the time accurate solution in Eq.(26) can be approximated as a linearfunction of time. Let’s first introduce the following notation, F p,k ( W n , δ ) = ˆ t n + δt n F p,k ( W n , t )d t. For convenience, assume t n = 0, the flux in the time interval [ t n , t n + ∆ t ] is expanded in thelinear form F p,k ( W n , t ) = F np,k + t∂ t F np,k . The coefficients F np,k and ∂ t F np,k can be fully determined by F p,k ( W n , t n )∆ t + 12 ∂ t F p,k ( W n , t n )∆ t = F p,k ( W n , ∆ t ) , F p,k ( W n , t n )∆ t + 18 ∂ t F p,k ( W n , t n )∆ t = F p,k ( W n , ∆ t/ . By solving the linear system, we have F p,k ( W n , t n ) = (4 F p,k ( W n , ∆ t/ − F p,k ( W n , ∆ t )) / ∆ t,∂ t F p,k ( W n , t n ) = 4( F p,k ( W n , ∆ t ) − F p,k ( W n , ∆ t/ / ∆ t . (32)Finally, with Eq.(31) and (32), W n +1 i at t n +1 can be updated by Eq.(30).The time dependent gas distribution function at a cell interface is updated in a similarway, f ∗ = f n + 12 ∆ tf nt ,f n +1 = f n + ∆ tf ∗ t . (33)12n order to construct the first-order time derivative of the gas distribution function, thedistribution function in Eq.(26) is approximated by the linear function f ( t ) = f ( x p,k , t, u , ξ ) = f n + f nt ( t − t n ) . According to the gas-distribution function at t = 0 and ∆ tf n = f (0) ,f n + f nt ∆ t = f (∆ t ) , the coefficients f n , f nt can be determined by f n = f (0) ,f nt = ( f (∆ t ) − f (0)) / ∆ t. Thus, f ∗ and f n +1 are fully determined and the macroscopic flow variables at the cellinterface can be obtained by Eq. (4). Theoretically, a fourth-order temporal accuracy canbe achieved for the conservative flow variables on arbitrary mesh. The proof is given in [35].
4. Compact HWENO reconstruction
In this section, a compact HWENO-type reconstruction is designed to get the piecewisediscontinuous flow variables and their first-order derivatives at each Gaussian point on bothsides of a cell interface. The reconstruction procedure for an inner cell is given first, thenthe special treatment for the boundary cell is presented in subsequent section 4.4.
As a starting point of WENO reconstruction, a linear reconstruction will be presentedfirst. For a piecewise smooth function Q ( x ) over cell Ω , a polynomial P r ( x ) with degree r can be constructed to approximate Q ( x ) as follows P r ( x ) = Q ( x ) + O (∆ h r +1 ) , where ∆ h ∼ | Ω | is the equivalent cell size. In order to achieve a third-order accuracy andsatisfy conservative property, the following quadratic polynomial over cell Ω is obtained P ( x ) = Q + (cid:88) | k | =1 a k p k ( x ) , (34)13here Q is the cell averaged value of Q ( x ) over cell Ω , k = ( k , k , k ), | k | = k + k + k .The p k ( x ) are basis functions, which are given by p k ( x ) = x k x k x k − | Ω | ˚ Ω x k x k x k d V. (35)The controlling points for a quadratic tetrahedron are shown in Fig. 1(b). The iso-parametric transformation is used to evaluate the volume integral, which can be writtenas X ( ξ, η, ζ ) = (cid:88) i =0 x i φ i ( ξ, η, ζ ) , where x i is the location of the ith controlling point and φ i is the base function as follows[27], v = ( − ζ + η + ξ )( − ζ + 2 η + 2 ξ ) ,v = ξ ( − ξ ) , v = η ( − η ) , v = ζ ( − ζ ) ,v = − ξ ( − ζ + η + ξ ) , v = 4 ηξ,v = − η ( − ζ + η + ξ ) , v = − ζ ( − ζ + η + ξ ) ,v = 4 ζξ, v = 4 ζη. (36)Then, the integration of monomial in Eq. (35) becomes ´ Ω x k y k z k d x d y d z = ´ ˜Ω x k y k z k ( ξ, η, ζ ) (cid:12)(cid:12)(cid:12) ∂ ( x,y,z ) ∂ ( ξ,η,ζ ) (cid:12)(cid:12)(cid:12) d ξ d η d ζ. (37)It can be evaluated numerically as ˝ Ω x k y k z k d x d y d z = (cid:80) Mm =1 ω m x k y k z k ( ξ, η, ζ ) m (cid:12)(cid:12)(cid:12) ∂ ( x,y,z ) ∂ ( ξ,η,ζ ) (cid:12)(cid:12)(cid:12) m ∆ ξ ∆ η ∆ ζ, (38)where ω m is the quadrature weight at the Gaussian point ( ξ, η, ζ ) m and ∆ ξ = ∆ η = ∆ ζ = 1.A five-point Gaussian quadrature with fourth-order spatial accuracy is used with( ξ, η, ζ ) = ( 14 , ,
14 ) , ( ξ, η, ζ ) = ( 12 , ,
16 ) , ( ξ, η, ζ ) = ( 16 , ,
16 ) , ( ξ, η, ζ ) = ( 16 , ,
12 ) , ( ξ, η, ζ ) = ( 16 , ,
16 ) , with ω = − , ω m = , m = 2 , , P ( x )In order to achieve a third-order spatial accuracy, the quadratic polynomial P ( x ) on Ω is constructed on the compact stencil S including Ω and its all von Neumann neighbors,Ω m , m = 1 , ...,
4, where the averages of Q ( x ) and averaged derivatives of Q ( x ) over each cellare known.The following values on S are used to obtain P ( x ),14 cell averages Q for cell 0, 1, 2, 3, 4, • cell averages of the x -direction partial derivative Q x for cell 1, 2, 3, 4; • cell averages of the y -direction partial derivative Q x for cell 1, 2, 3, 4; • cell averages of the z -direction partial derivative Q x for cell 1, 2, 3, 4.The polynomial P ( x ) is required to exactly satisfy ˚ Ω m P ( x )d V = Q m | Ω m | , (39)where Q m is the cell averaged value over Ω m , m = 1 , ...,
4, with the following conditionsatisfied in a least-square sense ˚ Ω m ∂∂x P ( x )d V = ( Q x ) m | Ω m | , ˚ Ω m ∂∂x P ( x )d V = ( Q x ) m | Ω m | , ˚ Ω m ∂∂x P ( x )d V = ( Q x ) m | Ω m | , (40)where Q x i , i = 1 , , m in a global co-ordinate, respectively. On a regular mesh, the system has 16 independent equations. Theconstrained least-square method is used to solve the above linear system [12]. The abovereconstruction improves the linear stability of the scheme and reduces the numerical errors.The left and right states W l,r provided by the reconstructed P ( x ) yields a linearly stablethird-order compact GKS, as validated through numerical tests in Section 5. P ( x ) and P ( x )In order to deal with discontinuity, lower-order polynomials from the sub-stencils shouldbe determined. Following the multi-resolution reconstruction in [37], the first-order poly-nomial P ( x ) is determined from the same central stencil as the P ( x ) but only with thecell-averaged conservative variables • Q for cell 0, 1, 2, 3, 4.The polynomial P ( x ) is required to satisfy ˚ Ω m P ( x )d V = Q m | Ω m | , m = 1 , , , , (41)in a least-square sense. 15ote that the left and right states W l,r solely determined by the reconstructed P ( x )yield an unstable second-order GKS. The theoretical proof for such instability on the second-order Riemann solver-based-FVM can be found in [5]. The zeroth-order polynomial P ( x ) issimply determined by the cell-averaged conservative variables on the targeted cell Ω itself,i.e. P ( x ) = Q . The coefficient matrices for the above P j ( x ) , j = 0 , Define three polynomials p ( x ) = 1 γ , P ( x ) − (cid:88) (cid:96) =0 γ (cid:96), γ , p (cid:96) ( x ) ,p ( x ) = 1 γ , P ( x ) − γ , γ , P ( x ) ,p ( x ) = P ( x ) . (42)For a third-order reconstruction, the second-order polynomial P ( x ) can be rewritten as P ( x ) = γ , p + γ , p + γ , p (43)with arbitrary positive coefficients γ m,n satisfying γ , + γ , + γ , = 1 , γ , + γ , = 1.For a second-order reconstruction, the first-order polynomial P ( x ) can be rewritten as P ( x ) = γ , p + γ , p (44)with arbitrary positive coefficients γ m,n satisfying γ , + γ , = 1. The coefficients are chosenas γ , : γ , : γ , = 100 : 10 : 1, and γ , : γ , = 10 : 1 as suggested in [37].The smoothness indicators β j , j = 1 , β j = r j (cid:88) | α | =1 | Ω | | α |− ˚ Ω (cid:0) D α P j ( x ) (cid:1) d V, (45)where α is a multi-index and D is the derivative operator, r = 1, r = 2. The smoothnessindicators in Taylor series at ( x , y ) have the order β = O {| Ω | [1 + O ( | Ω | )] } = O ( | Ω | ) = O ( h ) ,β = O {| Ω | [1 + O ( | Ω | )] } = O ( | Ω | ) = O ( h ) . Assuming a suitable β , β = O {| Ω | [1 + O ( | Ω | )] } = O ( | Ω | ) = O ( h ) ,
16 global smoothness indicator σ similar to that in [37] can be defined σ rd = ( 12 ( | β − β | + | β − β | )) = O ( | Ω | ) = O ( h ) , and σ nd = | β − β | = O ( | Ω | ) = O ( h ) . Then, the corresponding non-linear weights are given by ω m,n = γ m,n (1 + ( σ(cid:15) + β m ) ) , ¯ ω m,n = ¯ ω m,n (cid:80) ω m,n = γ m,n + O ( h ) , (46)where m = 0 , , n = 2 and m = 0 , n = 1, and (cid:15) takes 10 − to avoid zero in thedenominator.Replacing γ m,n by the normalized non-linear weights ¯ ω m,n in Eq. (43) and Eq. (44), thefinal reconstructed polynomials are given by R rd ( x ) = ¯ ω , p + ¯ ω , p + ¯ ω , p (47)for a third-order spatial accuracy, and R nd ( x ) = ¯ ω , p + ¯ ω , p (48)for a second-order spatial accuracy.As a result, the non-linear reconstruction meets the requirement for a third-order accu-racy R ( x ) = P ( x ) + O ( h ). If any of these values yield negative density or pressure, thefirst-order reconstruction is used instead. The desired non-equilibrium states at Gaussianpoints can be obtained from the weighted polynomials Q l,rp,k = R l,r ( x p,k ) , ( Q l,rx i ) p,k = ∂R l,r ∂x i ( x p,k ) . (49) According to the definition in Eq. (45), the smooth indicator of the zeroth-order polyno-mial P ( x ) is always 0. So, a new smooth indicator for P ( x ) has to be defined and can begiven as a non-linear combination of the first-order biased sub-stencils as suggested in [37].One of the choices is P on S = { ¯ Q , ¯ Q , ¯ Q , ¯ Q } , P on S = { ¯ Q , ¯ Q , ¯ Q , ¯ Q } ,P on S = { ¯ Q , ¯ Q , ¯ Q , ¯ Q } , P on S = { ¯ Q , ¯ Q , ¯ Q , ¯ Q } . In this plan, a total of 16 × × × × P can be greater than those of P and P under irregular gird and the WENOwill fail to suppress oscillations.Inspired by the method in [30], a two-step reconstruction is designed as follows to main-tain the compact manner of the scheme: • Reconstruction Step 1: Construct the first-order polynomial P ( x ) in each cell byEq. (41). Compute the slopes, e.g., b , b , b for each component, and store them. Acoefficient matrix with dimension 3 × P ( x ) and another matrix withdimension 5 × • Reconstruction Step 2: Conduct the multi-resolution reconstruction for each cell, andthe smooth indicator β for P ( x ) is given as a non-linear combination of the smoothindicators of P j ( x ) from the neighbor cells, i.e., β ,j = | ∆ | ( b ,j + b ,j + b ,j ) ,σ st = ( 16 ( (cid:88) | β ,j − β ,k | )) ,ω stj = 1 + σ st (cid:15) + β j , ¯ ω stj = ω j (cid:80) ω j ,β = (cid:88) ¯ ω stj β ,j , (50)where j, k = 1 , , , j > k . Then, the second-order reconstruction is complete.For the third-order reconstruction, only one extra beta β in Eq. (45) is needed. A coefficientmatrix with dimension 9 ×
16 is stored for P ( x ).Through the reconstruction step 1, the sub-stencils are extended to neighboring cells ofneighbors. Compared with the first plan, the robustness and mesh adaptability is signifi-cantly improved. 16 × × × The strategy of the two-step and multi-resolution reconstruction is extended to theboundary condition treatment. The one-sided reconstruction without ghost cell is adoptedhere with special care on Dirichlet boundary condition, i.e. the non-slip adiabatic walland the non-slip isothermal one. For the non-slip adiabatic boundary, the velocities areconstrained. For the non-slip isothermal boundary condition, both the velocities and thetemperature are constrained. 18
Reconstruction Step 1:For ith (i=0,...,4) conservative variables: – If there is no constraint for all the boundary faces: ∗ If the neighboring cell number is no less than 3, construct the first-orderpolynomial P ( x ). If the coefficient matrix is found to be nearly singular,which suggests a poor mesh quality, set P ( x ) = P ( x ). ∗ If the neighboring cell number is less than 3, using the cell-averaged slopesas the slopes of the first-order polynomial instead. – If there exists at least one constraint for all the boundary faces on the targeted cell,the weighted constrained least square reconstruction involving all the neighborcells and boundary faces are conducted. The weights for those boundary facesthat do not have constraint are set to be zero. In this step, each constrainedboundary face has one constraint, which is located at the geometric center of theface. ∗ If the sum of the neighboring cell number and the constraint number is noless than 3, and the constraint number is no greater than 3, construct thefirst-order polynomial P ( x ) by using constrained least-square method. ∗ If the constraint number is greater than 3 (which is impossible for tetrahedronmesh), construct the first-order polynomial P ( x ) by using the least-squaremethod. ∗ If the sum of the neighboring cell number and the constraint number is lessthan 3, use the cell-averaged slopes as the slopes of the first-order polynomialinstead.For non-slip adiabatic wall, each component is reconstructed in the followingorder. ∗ Step 1. One-sided reconstruction for density. ∗ Step 2. One-sided constrained reconstruction for momentum ρ U = ρ U wall where the reconstructed density is used. ∗ Step 3. One-sided reconstruction for energy.For non-slip isothermal wall, each component is reconstructed in the followingorder. ∗ Step 1. One-sided reconstruction for density. ∗ Step 2. One-sided constrained reconstruction for momentum ρ U = ρ U wall where the reconstructed density is used. ∗ Step 3. One-sided constrained reconstruction for energy ρE = ρ U wall + ρT wall / ( r − • Reconstruction Step 2. 19or a second-order reconstruction, the WENO procedure in Eq. (48) is complete. Ifa third-order reconstruction is adopted, the second-order polynomial P ( x ) is needed.For the ith (i=0,...,4) conservative variables: – If there is no constraint for all the boundary faces: ∗ If the neighboring cell number is no less than 3 and the coefficient matrixis not singular, construct the second-order polynomial P ( x ), by constrain-ing the cell-averaged values. Otherwise, the first-order polynomial from the P ( x ) stencils will be constructed in a least-squares sense. ∗ For the smooth reconstruction, the first-order polynomial using the P ( x )stencils is reconstructed instead. – If there exists at least one constrained face on the targeted cell, the weighted con-strained least square reconstruction involving all the neighbor cells and boundaryfaces are conducted. The weights for those boundary faces that do not have con-straint are set to be zero. Each constrained triangular face has three constraints,which are located at the corresponding Gaussian points. ∗ If the sum of neighboring cell-averaged data and the constraint number is noless than 9, and the constraint number is no greater than 9, construct thefirst-order polynomial P ( x ) by using constrained least-square method. ∗ If the constraint number is greater than 9 (which never happens in the testsof this paper), construct the first-order polynomial using the P ( x ) stencilsby the weighted least-square instead. ∗ If the sum of neighboring cell-averaged data and the constraint number isless than 9, construct the first-order polynomial using the P ( x ) stencil bythe weighted least-square instead.The constrained quantities and reconstruction order for the non-slip wall bound-aries are the same as those in the reconstruction step 1.It should be emphasized that the above criteria are general for other types of mesh andhybrid mesh. After obtaining the inner state (assume as ˜ W r ) at a boundary Gaussian point,a ghost state (assume as ˜ W l ) can be assigned according to boundary condition under localcoordinates. There is possible discontinuity between ˜ W l and ˜ W r if the WENO reconstruc-tion is used. The ghost state setting at the solid wall boundary is given as follows (the tildeis omitted). • Slip wall. The conservative variables under local coordinate ( ρ, ρU , ρU , ρU , ρE ) l =( ρ, − ρU , ρU , ρU , ρE ) r . The normal derivatives ( ρU ) lx = ( ρU ) rx while the normalderivatives for other components are W ilx = − W irx , i = 0 , , ,
4. The tangentialderivatives ( ρU ) lx j = − ( ρU ) rx j , j = 2 , W ilx j = W irx j , i = 0 , , , , j = 2 , Non-slip adiabatic wall. The conservative variables under local coordinate are given as( ρ, ρU , ρU , ρU , ρE ) l = ( ρ, − ρU , − ρU , − ρU , ρE ) r . The derivatives for all momenta( ρU i ) lx j = ( ρU i ) rx j , i = 1 , , , j = 1 , ,
3, while the normal derivatives for othercomponents are W ilx j = − W irx j , i = 0 , , j = 1 , , • Non-slip isothermal wall. – Assume the same pressure p l = p r . The velocity is opposite U li = − U ri , i =1 , ,
3. The temperature is set as T l = 2 T − T r , where T = T wall . Then, ρ l = p l /RT l = p r /RT l = p r /R (2 T − T r ). – Use primitive variables ( ρ, U , p ) l to get ( ρ, U , ρE ) l . – From the chain rule, ∂U i = ∂ ( ρU i ) − ∂ρU i ρ i = 1 , ,
3. Denote Q = (cid:80) U i , ∂Q = (cid:80) ∂U i U i . Then, ∂ρE = ∂ρQ + ∂Qρ + γ − ∂p and ∂p = ( γ − ∂ρE − ∂ρQ − ∂Qρ ).From p = ρRT , ∂T = ∂p − R∂ρTRρ is obtained. And ∂U li , i = 1 , , , ∂p l , and ∂T l aredetermined. – The derivatives of the primitive variables for the ghost states are set as ∂U li = ∂U ri , i = 1 , , ∂T l = ∂T r . ∂p l = − ∂p r . – Then, obtain ∂ρ l by ∂ρ = ∂p − R∂T ρRT and ∂ ( ρU i ) l by ∂ ( ρU i ) = ∂ρU i + ∂U i ρ . – Finally, get ∂ ( ρE ) l by ∂ρE = ∂ρQ + ∂Qρ + γ − ∂p .A summary for the reconstruction procedure is shown in Fig. 2. Yes No
Boundary cell? Construct by Eq.(41) and storethe slopes
Yes No
Is there any constrain for all the cellinterfaces? Contrain the requiredvariables on thebarycenter of the cellinterface One sidedreconstruction for each componentand store the slopes.Reduce to zeroth-order if necessary Construct from thestored slopes on thetargeted cell by Eq.(45) Construct from thestored slopes onneighbor cells by Eq.(50)
NoYes
Third-orderreconstruction? Obtain the desiredquantities by Eq.(49)Determine the finalWENO polynomial byEq.(46) and Eq.(48) Compute by Eq.(45) Determine the finalWENO polynomial byEq.(46-47)
Yes No
Boundary cell?
YesNo
Is there any constrain for all the cellinterfaces? One sidedreconstruction for each componentand store the slopes.Reduce to first-orderif necessary Constrain therequired variables onthe Gaussian pointsof the cell interface
Reconstrctuion Step 1 Reconstrctuion Step 2
Construct byEq.(40-41).
Figure 2: The process of the compact two-step multi-resolution reconstruction.
The reconstructions for the non-equilibrium states have the same spatial order and canbe used to get the equilibrium state g c , g cx i directly by a suitable average of g l,r , g l,rx i . To be21onsistent with the construction of g c , we make an analogy of the kinetic-based weightingmethod for g cx i , which are given by ˆ ψψψg c dΞ = W c = ˆ u> ψψψg l dΞ + ˆ u< ψψψg r dΞ , ˆ ψψψg cx i dΞ = W cx i = ˆ u> ψψψg lx i dΞ + ˆ u< ψψψg rx i dΞ . (51)The data for this method has compact support. In programming, this procedure is includedinside the subroutine of the gas distribution function, since it is performed at the localcoordinate. Thus, it is also cache-friendly. This method has been validated in the non-compact WENO5-GKS [8]. In this way, all components of the microscopic slopes in Eq.(26)have been fully obtained. It is worth to remark that the above reconstruction procedure canbe directly implemented to arbitrary mesh.
5. Numerical examples
In this section, numerical tests will be presented to validate the proposed scheme. Thetime step is determined by ∆ t = C CF L
Min( ∆ r i | U i | + ( a s ) i , (∆ r i ) ν i ) , (52)where C CF L is the CFL number, and | U i | , ( a s ) i , and ν i = ( µ/ρ ) i are the magnitude ofvelocities, sound speed, and kinematic viscosity coefficient for cell i. The ∆ r i is taken as theapproximated inscribed sphere radius of a tetrahedron,∆ r i = 3 | Ω i | (cid:80) | Γ ip | . All reconstructions will be performed on the conservative variables. Quadratic elements anda
CF L = 1 are used if no specified. An algorithm flowchart of the compact GKS is given inFig. 3.
The advection of density perturbation is tested with the initial condition ρ ( x, y, z ) = 1 + 0 . π ( x + y + z )) , U ( x, y, z ) = (1 , , , p ( x, y, z ) = 1 , within a cubic domain [0 , × [0 , × [0 , × N cells are used. With the periodic boundary condition in all directions, the22 o Yes
Smooth flow? Determine the equalibrium state by Eq.(51) Macroflow variables Update via two-stagetime discretizationCell-averagedconservative variables andtheir first-orderderivatives Linearreconstructionwith constrainedboundarycondition Gas evolution modelObtain interfacevalues , YesNoSmooth flow?Determine the time-accurate distributionfunction byEq.(26) Determine the time-accurate distributionfunction byEq.(27) The time-accurate fluxes andconservative variables , ,S2O4 in Eq.(30) and (33)Gauss's theorm in Eq.(29)Determine the non-equalibriumstate , by Eq.(13) and (14)
Yes No
Boundary cell?
NoYes
Boundary cell?Two-step non-linear 3rd-orderHWENOreconstructionLinear 3rd-orderreconstruction Two-step non-linear 3rd-orderHWENOreconstructionwith constrainedboundaryconditionTimestep inEq.(52)
Figure 3: The brief algorithm of the compact GKS. analytic solution is ρ ( x, y, z, t ) = 1 + 0 . π ( x + y + z − t )) , U ( x, y, z ) = (1 , , , p ( x, y, z, t ) = 1 . The collision time τ = 0 is set since the flow is smooth and inviscid. The L , L and L ∞ errors and the corresponding orders with linear and non-linear Z-type weights at t = 2 forthe third-order compact GKS are given in Table 1 and Table 2. The result with non-linearZ-type weights for the second-order scheme is also presented in Table 3. Expected accuracyis achieved for all cases. Mesh number L error Order L error Order L ∞ error Order6 × × × × Table 1: Accuracy test for the 3D sin-wave propagation by the linear third-order compact reconstruction.CFL=1.0. (a) Shu-Osher problem
This is the Shu-Osher problem [25] with the initial condition( ρ, U, p ) = (cid:40) (3 . , . , . , < x ≤ , (1 + 0 . x ) , , , < x < . esh number L error Order L error Order L ∞ error Order6 × × × × Table 2: Accuracy test for the 3D sin-wave propagation by the third-order compact HWENO reconstructionwith d : d : d = 100 : 10 : 1. CFL=1.0. mesh number L error Order L error Order L ∞ error Order6 × × × × Table 3: Accuracy test for the 3D sin-wave propagation by the second-order WENO reconstruction with d : d = 10 : 1. CFL=1.0. The computational domain is [0 , t = 1 . /
200 and 1 / x D e n s i t y Reference2nd3rd x D e n s i t y Reference2nd3rd
Figure 4: Shu-Osher problem. Mesh number: 6 × × × (b) Blast wave problem D e n s i t y Reference2nd3rd x D e n s i t y Reference2nd3rd
Figure 5: Shu-Osher problem. Mesh number: 6 × × × The initial conditions for the blast wave problem [29] are given as follows( ρ, u, p ) = (1 , , , ≤ x < . , (1 , , . , . ≤ x < . , (1 , , , . ≤ x ≤ . In the computational domain, 6 × × × × × × t =0 .
038 are presented in Fig. 6 and Fig. 7. Both the second-order and third-order schemesshow good robustness for such a strong shock-shock interaction. x D e n s i t y Reference2nd3rd x D e n s i t y Reference2nd3rd
Figure 6: Blastwave problem. Mesh number: 6 × × × D e n s i t y Reference2nd3rd x D e n s i t y Reference2nd3rd
Figure 7: Blastwave problem. Mesh number: 6 × × × A 3-D cavity is bounded in a unit cube and is driven by a uniform translation of the topboundary. In this case, the flow is simulated with Mach number
M a = 0 .
15 and γ = 5 /
3. Allboundaries are isothermal and nonslip. The computational domain [ − . , . × [ − . , . × [ − . , .
5] is covered by a uniform mesh with 6 × × ×
32 points and a refined uniformmesh with 5 × × ×
40 points, as shown in Fig. 8. A CFL number of 0.5 is used. Theflow is initialized with ρ = 1, U = 0 . U = U = 0, and p = 1 /γ . Since the flow is nearlyincompressible and mesh is regular, the smooth reconstruction and the simplified solver inEq. (26) are adopted in the computations. Figure 8: Lid-driven cavity flow. Left: uniform mesh with near wall size h = 1 /
32. Right: non-uniformmesh with near wall size h = 1 / (a) Re=1,000 For the Reynolds number Re = 1 , U -velocities alongthe line x = 0 , z = 0, and V -velocities along the line y = 0 , z = 0, are shown in Fig. 9. Thevelocity profiles from the third-order scheme match very well with the benchmark data [23].The velocity magnitude contours and streamlines by the third-order scheme are shown inFig. 10. The cavity case demonstrates the high-order accuracy of the compact GKS. y U Reference2nd3rd x V Reference2nd3rd
Figure 9: Lid-driven cavity flow: Re=1,000. The velocities profiles compared with the reference data in[23]. (b) Re=3,200
The flow becomes transient when
Re > U velocity component along the line x = 0 , z = 0, and the V velocity com-ponent along the line y = 0 , z = 0 are presented in Fig. 11. Both results under the uniformmesh and the refined non-uniform mesh agree well with the experimental data. A betteragreement in the RMS U-velocity U rms can be observed in the refined mesh calculation. A low-speed viscous flow passing through a sphere is tested. The Reynolds number basedon the diameter of the sphere D = 1 is 118. In such case, a drag coefficient C D = 1 wasreported from the experiment in [26]. The far-field flow condition outside boundary of thedomain is set with the free stream condition( ρ, U, V, W, p ) ∞ = (1 , . , , , γ ) , with γ = 1 . M a ∞ = 0 . h ≈ . × − D , as shown in Fig. 12. Bothsecond-order and third-order schemes with non-linear reconstructions are tested. A cleanand symmetric velocity contour is observed from the third-order compact GKS, as shown inFig. 13. The pressure contour and the 3-D streamline are also presented in Fig. 14, where27 igure 10: Lid-driven cavity flow: Re=1,000. Top left: The velocity magnitude contours. Others: Thestreamlines on x = 0, y = 0, and z = 0 planes. Ux V y
1 0.5 0 0.5 1 0.4 0.2 0 0.2 0.4 1 0.500.51 0.4 0.200.20.4
Expcoarsefine U rms x V r m s y
1 0.5 0 0.5 1 0.4 0.2 0 0.2 0.4 1 0.500.51 0.4 0.200.20.4
Expcoarsefine
Figure 11: Lid-driven cavity flow: Re=3,200. The mean and RMS velocity profiles obtained by thethird-order compact GKS are compared with the experimental data in [19, 20]. C D , the separation angle θ , and the closest wake length L , as defined in [10]. Figure 12: Subsonic flow passing through a viscous sphere. Mesh number: 399,546.Figure 13: Subsonic flow passing through a viscous sphere. Ma=0.2535. Re=118. Left: The second-orderGKS. Right: The third-order GKS.
To validate the robustness of the current scheme for the high-speed viscous flow, asupersonic flow passing through a sphere with
M a = 1 . D = 1. The Prandtl number is P r = 1. The tetrahedral meshwith an upstream length of 5 and a downstream length of 40 is shown in Fig. 15. The firstmesh size at the wall has a thickness 2 . × − D . The numerical results obtained by the29 igure 14: Subsonic flow passing through a viscous sphere by the third-order compact GKS. Ma=0.2535.Re=118. Scheme Mesh number Cd θ L ClExperiment [26] – 1.0 151 1.07 –Current 2nd 399,546 1.027 126.9 1.00 1.5e-2Current 3rd 399,546 1.018 127.4 1.00 1.7e-3Implicit third-order DDG [3] 160,868 1.016 123.7 0.96 –Implicit fourth-order VFV [27] 458,915 1.014 – – 2.0e-5Implicit third-order AMR-VFV [17] 621,440 1.016 – – –
Table 4: Quantitative comparisons among different compact schemes for the viscous flow over a sphere.
Figure 15: Supersonic flow passing through a viscous sphere. Mesh number: 665,914.Figure 16: Supersonic flow passing through a viscous sphere by the third-order compact GKS. Ma=1.2.Re=300.
As a classic validation case for compressible external flow, the transonic flow over theONERA M6 wing is tested. Experimental data are reported in [21], where the flow is fullyturbulent. Same as the inviscid calculation in [13], an incoming Mach number Ma=0.8395and an angle of attack AOA=3.06 ◦ are used, which correspond to a rough prediction of theflow field under a very high Reynolds number. In the computation, the wing has a slipwall boundary condition, and the Riemann invariants are applied 10 times of the root chord31 cheme Grid Number Cd θ L Shock stand-offWENO6 [16] 909,072 1.282 126.9 1.61 0.69Current 3rd 665,914 1.274 126.3 1.64 0.66-0.82 (0.72)
Table 5: Quantitative comparisons between the current scheme and the reference solution for the supersonicviscous flow over a sphere. length away from the wing. Two sets of meshes are used to test the mesh sensitivity, as shownin Fig. 17. For each mesh, the results from the second and third-order GKS are presented.The surface pressure distributions and Mach number slices at different wing sections underMesh I for both schemes are shown in Fig. 18. The “Lambda” shock is resolved from bothschemes. Third-order scheme presents accurate solutions with high resolution in pressureand Mach contours in smooth region. Similar conclusions can be drawn from the resultsobtained from Mesh II, as shown in Fig. 20. Quantitatively comparisons on the pressuredistributions at six different locations on the wing are given in Fig. 19 and Fig. 21. Abetter agreement in the secondary shock position is obtained with Mesh II at the semi-spanlocations Y/B = 0.20, 0.44, and 0.65.
Figure 17: Mesh for the inviscid ONERA M6 wing. Left: Mesh I with 294,216 cells. Right: Mesh II with347,094 cells.
The inviscid supersonic flow passing through a complete aircraft model is computed.The computational mesh for a YF-17 (”Cobra”) fighter model is shown in Fig. 22 whichis provided at “https://cgns.github.io/CGNSFiles.html”. A free stream at a Mach numberMa=1 . .
25 are adopted as the initial conditions. The surfacepressure, Mach number distributions, and streamlines are presented in Fig. 23 for the GKSwith the second-order WENO reconstruction. Complicated shocks appear in the locationsincluding the nose, cockpit-canopy wing, horizontal stabilizer, and vertical stabilizer. Aslightly smoother solution is obtained by the compact GKS with the third-order HWENO32 igure 18: Transonic flow over an inviscid ONERA M6 wing under Mesh I. Ma=0.8935. AOA=3.06 ◦ . Left:The second-order GKS. Right: The third-order GKS. X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
Figure 19: Pressure distributions for wing section at different semi-span locations Y/B on the ONERAM6 wing under Mesh I. Ma=0.8935. AOA=3.06 ◦ . Top: Y/B=0.20, 0.44, 0.65 from left to right. Bottom:Y/B=0.80, 0.90, 0.95 from left to right.Figure 20: Transonic flow over an inviscid ONERA M6 wing under Mesh II. Ma=0.8935. AOA=3.06 ◦ .Left: The second-order GKS. Right: The third-order GKS. /C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
X/C C p Experiment2nd3rd
Figure 21: Pressure distributions for wing section at different semi-span locations Y/B on the ONERAM6 wing under Mesh II. Ma=0.8935. AOA=3.06 ◦ . Top: Y/B=0.20, 0.44, 0.65 from left to right. Bottom:Y/B=0.80, 0.90, 0.95 from left to right. reconstruction, as shown in Fig. 24. The maximum Mach number on the surface is 2.4for the second-order scheme and 2.28 for the third-order one. The current algorithm canhandle complicated geometry, such as the mesh skewness near the wing tips and the lack ofneighboring cell for the cell near boundary corners. The compact GKS demonstrates goodmesh adaptability in the computation. A space-shuttle-like blunt-body model is considered to test the robustness of the schemesfor the hypersonic inviscid flow. The initial condition has Ma=5 and AOA=0 ◦ . The surfacemesh is given in Fig. 25, where the controlling points of the quadratic elements are shown.The pressure distributions are shown in Fig. 26, where no significant differences are observedin the results from the second-order and the third-order GKS. The Mach number distributionand streamlines are also plotted in Fig. 27.
6. Conclusion
A 3rd-order compact GKS is developed for 3-D tetrahedral mesh. The main difficultyfrom the structured to tetrahedral mesh is related to the linear instability in unstructuredmesh. On a compact stencil with von Neumann neighbors only, even a second-order FVM intetrahedral mesh can become linearly unstable. The high-order method based on Riemann-solvers with a compact stencil is also associated with instability. Also, the traditional WENO34 YZ Figure 22: Supersonic flow passing through a YF-17 (”Cobra”) model. Ma=1.8. AOA=1.25. Mesh number:325,096.Figure 23: Supersonic flow passing through a YF-17 (”Cobra”) model by the second-order GKS. Ma=1.8.AOA=1 . ◦ .Figure 24: Supersonic flow passing through YF-17 (”Cobra”) model by the third-order GKS. Ma=1.8.AOA=1 . ◦ . igure 25: Hypersonic inviscid flow over a blunt body. Mesh number: 117,221.Figure 26: Pressure distributions around the blunt body. Ma=5.0. AOA=0 . ◦ . Left: the second-orderGKS. Right: the third-order GKS. igure 27: Mach number distribution, stream-lines, and surface density distributions around the bluntbody. Ma=5.0. AOA=0 . ◦ . Left: the second-order GKS. Right: the third-order GKS. strategy based on the above compact stencil fails in dealing with discontinuities. However,benefiting from the direct evolution of the cell-averaged first-order spatial derivatives inthe compact gas-kinetic scheme, the linear stability of the compact third-order GKS hasbeen validated on tetrahedral mesh through the smooth inviscid and viscous tests. Tofurther improve the mesh adaptability and robustness of the scheme, a new reconstructionbased on the two-step and multi-resolution WENO methods is proposed. At the sametime, a new second-order GKS can be naturally obtained as a byproduct. Both second andthird-order schemes keep the compactness. The reconstruction in this paper is carefullydesigned with the consideration of possible singularities from the mesh distortion or theboundary corner, and it becomes suitable for arbitrary mesh. The compact GKS also usesthe two-stage time discretization as a building block for temporal accuracy, which becomesefficient in comparison with the Runge-Kutta time stepping method for the same third-ordertemporal accuracy. Various numerical examples from low-speed smooth flow to hypersonicflow are tested. The compact GKS shows properties of robustness, high accuracy, and lowdissipation. Reliable mesh adaptability is also validated in the supersonic flow computationover a complete aircraft model. Moreover, a large explicit time step with a CFL of 1can be used for most test cases. The proposed compact GKS with the two-stage timediscretization and the two-step multi-resolution WENO reconstruction exhibits excellentnumerical performance among the current existing compact schemes on tetrahedral mesh.The compact GKS is currently extended to hybrid mesh with high aspect ratio for theboundary layer flow computation in supersonic and hypersonic viscous flow. Acknowledgments
The authors would like to thank Dr. Jun Zhu for helpful discussion, and be grateful toMr. Nianhua Wang, Dr. Yangyang Liu, and Dr. Liming Yang for providing computational37esh. The current research is supported by National Numerical Windtunnel project andNational Science Foundation of China (11772281, 91852114).
ReferencesReferences [1] Antonis F. Antoniadis, Panagiotis Tsoutsanis, and Dimitris Drikakis. Assessment of high-order finitevolume methods on unstructured meshes for rans solutions of aeronautical configurations.
Computers& Fluids , 146:86 – 104, 2017.[2] Prabhu Lal Bhatnagar, Eugene P Gross, and Max Krook. A model for collision processes in gases. i.Small amplitude processes in charged and neutral one-component systems.
Physical Review , 94(3):511,1954.[3] Jian Cheng, Xiaodong Liu, Tiegang Liu, and Hong Luo. A parallel, high-order direct discontinuousGalerkin method for the Navier-Stokes equations on 3d hybrid grids.
Communications in ComputationalPhysics , 21(5):1231–1257, 2017.[4] Michael Dumbser. Arbitrary high order PNPM schemes on unstructured meshes for the compressibleNavier–Stokes equations.
Computers & Fluids , 39(1):60–76, 2010.[5] F. Haider, J.-P. Croisille, and B. Courbet. Stability analysis of the cell centered finite-volume MUSCLmethod on unstructured grids.
Numerische Mathematik , 113(4):555–600, 2009.[6] Hung T Huynh. A flux reconstruction approach to high-order schemes including discontinuous Galerkinmethods. In , page 4079, 2007.[7] Xing Ji.
High-order non-compact and compact gas-kinetic schemes . PhD thesis, Hong Kong Univeristyof Science and Technology, 2019.[8] Xing Ji and Kun Xu. Performance enhancement for high-order gas-kinetic scheme based on WENO-adaptive-order reconstruction.
Communications in Computational Physics , 28(2):539–590, 2020.[9] Xing Ji, Fengxiang Zhao, Wei Shyy, and Kun Xu. A HWENO reconstruction based high-order compactgas-kinetic scheme on unstructured mesh.
Journal of Computational Physics , page 109367, 2020.[10] Xing Ji, Fengxiang Zhao, Wei Shyy, and Kun Xu. A three-dimensional compact high-order gas-kineticscheme on structured mesh. arXiv preprint arXiv:2009.02908 , 2020.[11] Jiequan Li and Zhifang Du. A two-stage fourth order time-accurate discretization for Lax–Wendrofftype flow solvers I. hyperbolic conservation laws.
SIAM Journal on Scientific Computing , 38(5):A3046–A3069, 2016.[12] Wanai Li.
Efficient implementation of high-order accurate numerical methods on unstructured grids .Berlin, Heidelberg: Springer, 2014.[13] Yangyang Liu, Liming Yang, Chang Shu, and Huangwei Zhang. Three-dimensional high-orderleast square-based finite difference-finite volume method on unstructured grids.
Physics of Fluids ,32(12):123604, 2020.[14] Hong Luo, Luqing Luo, Robert Nourgaliev, Vincent A Mousseau, and Nam Dinh. A reconstructeddiscontinuous Galerkin method for the compressible Navier–Stokes equations on arbitrary grids.
Journalof Computational Physics , 229(19):6961–6978, 2010.[15] Dimitri Mavriplis. Revisiting the least-squares procedure for gradient reconstruction on unstructuredmeshes. In , page 3986, 2003.[16] T. Nagata, T. Nonomura, S. Takahashi, Y. Mizuno, and K. Fukuda. Investigation on subsonic tosupersonic flow around a sphere at low reynolds number of between 50 and 300 by direct numericalsimulation.
Physics of Fluids , 28(5):056101, 2016.[17] Jianhua Pan, Qian Wang, Yusi Zhang, and Yuxin Ren. High-order compact finite volume methodson unstructured grids with adaptive mesh refinement for solving inviscid and viscous flows.
ChineseJournal of Aeronautics , 31(9):1829–1841, 2018.
18] Liang Pan, Kun Xu, Qibing Li, and Jiequan Li. An efficient and accurate two-stage fourth-ordergas-kinetic scheme for the Euler and Navier–Stokes equations.
Journal of Computational Physics ,326:197–221, 2016.[19] Ajay Prasad, Chin-Yuan Perng, and Jeffrey Koseff.
Some Observations on the Influence of LongitudinalVortices in a Lid-Driven Cavity Flow , page 3654. 1988.[20] Ajay K Prasad and Jeffrey R Koseff. Reynolds number and end-wall effects on a lid-driven cavity flow.
Physics of Fluids A: Fluid Dynamics , 1(2):208–218, 1989.[21] V Schmitt. Pressure distributions on the ONERA M6-wing at transonic mach numbers, experimentaldata base for computer program assessment.
AGARD AR-138 , 1979.[22] David C Seal, Yaman G¨u¸cl¨u, and Andrew J Christlieb. High-order multiderivative time integrators forhyperbolic conservation laws.
Journal of Scientific Computing , 60(1):101–140, 2014.[23] C Shu, L Wang, and YT Chew. Numerical computation of three-dimensional incompressible Navier–Stokes equations in primitive variable form by DQ method.
International Journal for Numerical Meth-ods in Fluids , 43(4):345–368, 2003.[24] Chi-Wang Shu. High order WENO and DG methods for time-dependent convection-dominated PDEs:A brief survey of several recent developments.
Journal of Computational Physics , 316:598 – 613, 2016.[25] Chi-Wang Shu and Stanley Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. In
Upwind and High-Resolution Schemes , pages 328–374. Springer, 1989.[26] Sadatoshi Taneda. Experimental investigation of the wakes behind cylinders and plates at low Reynoldsnumbers.
Journal of the Physical Society of Japan , 11(3):302–307, 1956.[27] Qian Wang.
Compact High-Order Finite Volume Method on Unstructured Grids . PhD thesis, TsinghuaUniversity, 6 2017.[28] ZJ Wang, Y Li, F Jia, GM Laskowski, J Kopriva, U Paliath, and R Bhaskaran. Towards industriallarge eddy simulation using the FR/CPR method.
Computers & Fluids , 156:579–589, 2017.[29] Paul Woodward and Phillip Colella. The numerical simulation of two-dimensional fluid flow with strongshocks.
Journal of computational physics , 54(1):115–173, 1984.[30] Yidong Xia, Xiaodong Liu, and Hong Luo. A finite volume method based on WENO reconstructionfor compressible flows on hybrid grids. In , page 0939, 2014.[31] Kun Xu. A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificialdissipation and Godunov method.
Journal of Computational Physics , 171(1):289–335, 2001.[32] Kun Xu.
Direct Modeling for Computational Fluid Dynamics: Construction and Application of UnifiedGas-Kinetic Schemes . World Scientific, 2014.[33] Xiaoquan Yang, Jian Cheng, Hong Luo, and Qijun Zhao. Robust implicit direct discontinuous Galerkinmethod for simulating the compressible turbulent flows.
AIAA Journal , 57(3):1113–1132, 2019.[34] Meilin Yu, Z.J. Wang, and Yen Liu. On the accuracy and efficiency of discontinuous Galerkin, spectraldifference and correction procedure via reconstruction methods.
Journal of Computational Physics ,259:70 – 95, 2014.[35] Fengxiang Zhao, Xing Ji, Wei Shyy, and Kun Xu. A compact high-order gas-kinetic scheme on un-structured mesh for acoustic and shock wave computations. arXiv preprint arXiv:2010.05717 , 2020.[36] Zhong Zhao, Lei He, and Xianyao He. Design of general CFD software PHengLEI (in chinese).
ComputerEngineering & Science , 42(2):210–219, 2020.[37] Jun Zhu and Chi-Wang Shu. A new type of third-order finite volume multi-resolution WENO schemeson tetrahedral meshes.
Journal of Computational Physics , 406:109212, 2020., 406:109212, 2020.