A three-dimensional compact high-order gas-kinetic scheme on structured mesh
AA three-dimensional compact high-order gas-kinetic scheme onstructured 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 is firstly proposed for three-dimensionalcomputation for the compressible Euler and Navier-Stokes solutions. As an extension of theprevious compact gas-kinetic scheme (GKS) on two-dimensional structured/unstructuredmesh [12, 15], the new scheme is based on three key ingredients: the time-accurate gas-kineticevolution solution, the Hermite weighted essentially non-oscillatory (HWENO) reconstruc-tion, and the two-stage temporal discretization. The scheme achieves its compactness dueto the time-dependent gas distribution function in GKS, which provides not only the fluxesbut also the time accurate flow variables in the next time level at a cell interface. As a result,the cell averaged first-order spatial derivatives of flow variables can be obtained naturallythrough the Gauss’s theorem. Then, a third-order compact reconstruction involving the cellaveraged values and their first-order spatial derivatives can be achieved. The trilinear inter-polation is used to treat possible non-coplanar elements on general hexahedral mesh. Theconstrained least-square technique is applied to improve the accuracy in the smooth case. Todeal with both smooth and discontinuous flows, a new HWENO reconstruction is designedin the current scheme by following the ideas in [48]. No identification of troubled cells isneeded in the current scheme. In contrast to the Riemann solver-based method, the compactscheme can achieve a third-order temporal accuracy with the two-stage two-derivative tem-poral discretization, instead of the three-stage Runge-Kutta method. Overall, the proposedscheme inherits the high accuracy and efficiency of the previous ones in two-dimensional case.The desired third-order accuracy can be obtained with curved boundary. The robustness ofthe scheme has been validated through many cases, including strong shocks in both inviscidand viscous flow computations. Quantitative comparisons for both smooth and discontinu-ous cases show that the current third-order scheme can give competitive results against thefifth-order non-compact GKS under the same mesh. A large CFL number around 0.5 canbe used in the present scheme.
Keywords: compact gas-kinetic scheme, Hermite WENO reconstruction, two-stage timediscretization, Navier-Stokes solution
Preprint submitted to Elsevier September 8, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p . Introduction Tremendous efforts have been paid on the development of high-order computational fluiddynamics (CFD) methods for the compressible Euler and Navier-Stokes (N-S) equations inpast decades. Representative methods include weighted essentially non-oscillatory (WENO)methods [28], discontinuous Galerkin (DG) methods [29], the flux reconstruction (FR) [10]or correction procedure via reconstruction (CPR) [43] methods, etc. WENO-type recon-struction has been widely applied on structured grids, which can keep good robustness withvery high order of accuracy in space [1]. However, it is not a trivial task to extend theWENO approach to unstructured mesh or non-uniform mesh with the same high-order ac-curacy, since the stencil is too large. Zhu et al. has proposed a class of new WENO schemesin the attempt for releasing the problems recently [49].On the other hand, high-order methods based on the compact stencil which only involvesthe target cell and its von Neumann neighbors are attractive because of their good meshadaptability, high scalability, which becomes a hot topic in CFD research [34]. Two mainrepresentatives of compact methods are the DG methods, which combine the finite volumeframework and the finite element framework and the FR/CPR methods, which hybridizethe finite difference and finite volume discretization originally. These methods can achievearbitrary spatial order of accuracy with only the targeted cell as the reconstruction stencil,which yields a great mesh adaptability. Numerical results have demonstrated their power inlarge eddy simulation (LES) [35] and RANS simulation [42] for smooth flow. However, in theflow simulations with strong discontinuities, these methods seem to lack robustness. Manytechniques have been used to limit the troubled internal degree of freedom [26]. At the sametime, they have more restricted explicit time step for linear stability against the traditionalhigh-order finite volume method. New compact schemes have been proposed to improvethe above shortcomings, such as the multi-moment constrained finite volume scheme [36], P N P M [41] scheme. Most of these methods use Riemann solvers or approximate Riemannsolvers for the flux evaluation, and the Runge-Kutta time-stepping methods are adopted forthe temporal accuracy.In recent years, a class of compact high-order gas-kinetic scheme (HGKS) has beendeveloped [22]. With the adoption of the two-stage time discretization [18, 7], and theHermite WENO (HWENO) reconstruction [25], a fourth-order compact gas-kinetic scheme(GKS) on two-dimensional structured mesh is constructed [12]. It has higher resolution thanthe conventional non-compact fourth-order GKS [23] while performs competitive robustnessas the second-order scheme. Most importantly, it allows a CFL number around 0.5 inthe computation. The success of the above scheme lies in the gas-kinetic framework. Incomparison with traditional Riemann solver based high-order CFD methods, it includesthe following distinguishable points: (i) GKS is based on an analytical integral solutionof the BGK equation, which can recover the N-S equations from the Chapman-Enskog ∗ 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 three-dimensional gas-kinetic BGK equation [2] can be written as f t + u · ∇ f = g − fτ , (1)where f is the gas distribution function, g is the corresponding equilibrium state, and τ isthe collision time. f = f ( x , t, u , ξ ), where x is location in physical space, t is time and u isparticle velocity in phase space.The collision term satisfies the following 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 degreeof freedom, i.e. K = (5 − γ ) / ( γ −
1) in three-dimensional case, and γ is the specific heatratio. 4ased on the Chapman-Enskog expansion for BGK equation [39], the gas distributionfunction in the continuum regime can be expanded as f = g − τ D u g + τ D u ( τ D u ) g − τ D u [ τ D u ( τ D u ) g ] + ..., where D u = ∂/∂t + u · ∇ . By truncating on different orders of τ , the corresponding macro-scopic equations can be derived. If the zeroth-order truncation is taken, i.e., f = g , theEuler equations can be recovered by multiplying ψψψ on Eq.(1) and integrating it over thephase space, W t + ∇ · F ( W ) = 0 . If the first-order truncated distribution function is applied, i.e., f = g − τ ( u · ∇ g + g t ) , (3)the N-S equations can be obtained, W t + ∇ · F ( W , ∇ W ) = 0 , with τ = µ/p and a fixed Prandtl number P r = 1.Taking moments of the time-dependent distribution function f , the conservative variablescan be obtained W ( x , t ) = ˆ ψψψf ( x , t, u , ξ )dΞ . (4)so as the macroscopic fluxes F ( W ( x , t )), F ( x , t ) = ˆ u ψψψf ( x , t, u , ξ )dΞ . (5) Remark 1.
It is well known that the cell-averaged conservative variables can be updatedthrough the interface fluxes under the finite volume framework. Beside the fluxes, Eq. (4) provides additional information, which can be updated at the next time level. It is the keypart in constructing the compact GKS, which will be introduced in detail in subsection 2.3.An obvious prerequisite is that the W ( x , t ) must be time-accurate, while the time-independentRiemann solution cannot provide such a time accurate evolution solution2.1. Finite volume scheme on general structured mesh For a polyhedron cell Ω i in 3-D case, the boundary can be expressed as ∂ Ω i = N f (cid:91) p =1 Γ ip , N f is the number of cell interfaces for cell Ω i , e.g., N f = 4 for tetrahedron and N f = 6for cuboid or general hexahedron.The increment of the cell averaged conservative flow variables in a finite control volumei in a time interval [ t n , 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 value over cell Ω i , | Ω i | is the volume of Ω i , F is the interfaceflux, 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)Numerical quadratures can be adopted to give a high-order spatial approximation for F ip ( t ) or (cid:101) F ip ( t ), where Eq.(6) can be rewritten 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) Remark 2.
The bilinear interpolation is used to describe a given quadrilateral interface Γ ip with coplanar or non-coplanar vertexes, X ( ξ, η ) = (cid:88) l =1 x l φ l ( ξ, η ) , where ( ξ, η ) ∈ [ − / , / , x l is the locations of the mth vertex and φ l is the base functionas follows φ = 14 (1 − ξ )(1 − η ) , φ = 14 (1 − ξ )(1 − η ) ,φ = 14 (1 − ξ )(1 + 2 η ) , φ = 14 (1 − ξ )(1 + 2 η ) , The flux across Γ ip in Eq. (6) can be rewritten as F ip ( t ) = ˆ Γ ip F ( x , t ) · n p d s = ˆ / − / ˆ / − / 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 approx- mated through Gaussian quadrature as F ip ( t ) = (cid:88) m,n =1 ω m,n F m,n ( t ) · ( n p ) m,n (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, y, z ) ∂ ( ξ, η ) (cid:12)(cid:12)(cid:12)(cid:12) m,n ∆ ξ ∆ η, where ∆ ξ = ∆ η = 1 , and the local normal direction ( n p ) m,n = ( X ξ × X η ) / (cid:107) X ξ × X η (cid:107) .The standard Gaussian points are ( ξ, η ) , = ( − √ ξ, − √ η ) , ( ξ, η ) , = ( − √ ξ, √ η ) , ( ξ, η ) , = ( 12 √ ξ, − √ η ) , ( ξ, η ) , = ( 12 √ ξ, √ η ) , with ω , = ω , = ω , = ω , = .Compared with Eq. (9) , we have ω k = ω m,n (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, y, z ) ∂ ( ξ, η ) (cid:12)(cid:12)(cid:12)(cid:12) m,n , x p,k = x (( ξ, η ) m,n ) , n p,k = ( n p ) m,n , with k = 2 m + n , m = 1 , , n = 1 , . According to the coordinate transformation, the local coordinate for the cell interfaceΓ 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) = − . (10)And the macroscopic conservative variables in the local coordinate are given 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) = − . (11)Note that when n = −
1, the Eq.(10) changes to ( ˜ u , ˜ u , ˜ u ) T = ( − u , − u , u ) T and thematrix (11) 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 ξ. (12)In the computation, the fluxes in the local coordinate are obtained first by taking momentsof the gas distribution function in the local coordinate (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 ξ, (13)where (cid:101) ψψψ = (1 , (cid:101) u ,
12 ( (cid:101) u + ξ )) T . According to Eq.(10), Eq.(12) and Eq.(13), 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 )) . (14) 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 , ξ ) , (15)where x = x (cid:48) + u ( t − t (cid:48) ) is the trajectory of particle. f is the initial gas distribution function, g is the corresponding equilibrium state. The integral solution basically states a physicalprocess from the particle free transport in f in the kinetic scale to the hydrodynamic flowevolution in the integral of g term. The flow evolution at the cell interface depends on theratio of time step to the local particle collision time ∆ t/τ .To construct a time evolution solution of gas distribution function at a cell interface, thefollowing notations 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 [38] s = s j ψ j = s + s u + s u + s u + s
12 ( u + u + u + ξ ) . For the kinetic part of the integral solution Eq.(15), the initial gas distribution function canbe constructed 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 distribution8unctions 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 = is 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 , (16)for k = l, r . According to Eq.(3), f kG has the form f kG ( ) = g k ( ) − τ ( u i g kx i ( ) + g kt ( )) , (17)where g k are the equilibrium states 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 as ˆ ψψψg l dΞ = W l , ˆ ψψψg r dΞ = W r . (18)Substituting Eq.(16) and Eq.(17) into Eq.(15), 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 ] , (19)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.(15) canbe expanded in space and time as follows g ( x , t ) = g c ( ,
0) + ∂g c ∂x i ( , x i + ∂g c ∂t ( , t, (20)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 . (21)Substituting Eq.(20) into Eq.(15), the hydrodynamic part for 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 , (22)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.(19) and Eq.(22) can be determined by the spatial derivatives of9acroscopic 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 , (23)where (cid:104) ... (cid:105) are the moments of a gas distribution function defined by (cid:104) ( ... ) (cid:105) = ˆ ψψψ ( ... ) g dΞ . (24)The details for calculation of each microscopic term from macroscopic quantities can refer[11]. Then the second-order time dependent gas distribution function is given as f ( , t, u , ξ ) =(1 − e − t/τ ) g c + [( t + τ ) e − t/τ − τ ] a cx i u i g c + ( t − τ + τ e − t/τ ) A c g c + e − t/τ g l [1 − ( τ + t ) a lx i u i − τ A l ] H ( u )+ e − t/τ g r [1 − ( τ + t ) a rx i u i − τ A r ](1 − H ( u )) . Finally in order to properly capture the un-resolved shock structure, additional numericaldissipation is needed. The physical collision time τ in the exponential function part can bereplaced by a numerical collision time τ n , which will be defined later, 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 )) . (25) Distinguished from the approximate Riemann solver with a constant state at a cell in-terface, the gas-kinetic scheme provides a time evolution solution. Recall Eq.(4), the con-servative variables at the Gaussian point x p,k can be updated by the moments with ψψψ , W p,k ( t n +1 ) = ˆ ψψψf n ( x p,k , t n +1 , u , ξ )dΞ , k = 1 , ..., M (26)According to the Gauss’s theorem, the cell-averaged first-order derivatives within each10lement at t n +1 can be given by 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 , (27)where n m,k = (( n ) m,k , ( n ) m,k , ( n ) m,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 which has been adopted in theprevious compact schemes on two-dimensional cases is implemented here [12, 15]. Followingthe definition of Eq.(8), a fourth-order time-accurate solution for cell-averaged conservativeflow variables W i is 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) , (28)where L ( W ni ) and ∂∂t L ( W ni ) are 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 . (29)11he proof for fourth-order accuracy can be found in [18].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 flux function, as shown in Eq.(25), can be approximated asa linear function of time within a time interval. 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 as thefollowing linear 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 . (30)Finally, with Eq.(29), and Eq.(30), W n +1 i at t n +1 can be updated by Eq.(28). Remark 3.
The use of S2O4 time integration in GKS is due to the following reasons.Firstly, the fourth-order two-derivative method with two-stage is unique and is efficient incomparison with the RK method. Secondly, the accuracy and robustness of S2O4 time in-tegration has been validated by numerical tests for both non-compact and compact schemes[24, 14, 12, 23, 13], and has been extended to compressible multi-component flow [21], hyper-sonic non-equilibrium multi-temperature flow [3], and direction simulation of compressiblehomogeneous turbulence [4]. The S2O4 method becomes a building block of a family of HGKSfor time integration.
Similar to the two-stage temporal discretization in the flux evaluation, the time depen-dent gas distribution function at a cell interface is updated as f ∗ = f n + 12 ∆ tf nt ,f n +1 = f n + ∆ tf ∗ t , (31)where a second-order evolution model is used for the update of gas distribution function onthe cell interface and for the evaluation of flow variables.12n order to construct the first-order time derivative of the gas distribution function, thedistribution function in Eq.(25) 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 at the cell interface for the evaluation of macroscopicflow variables. This temporal evolution for the interface value is similar to the one used inGRP solver [7] on Cartesian grid. There is no rigorous proof for its temporal accuracyon general mesh so far. However, numerical tests demonstrate that the proposed compactscheme achieves a third-order temporal accuracy.
4. Compact HWENO reconstruction
In this section, the discontinuous values and their first-order derivatives of flow variablesat each Gaussian point of a cell interface will be constructed by a newly designed compactHWENO-type reconstruction. Then, based on such an initial condition the time-dependentgas distribution function in Eq.25 can be fully determined.As a starting point of WENO reconstruction, a linear reconstruction will be presentedfirst. For a piecewise smooth function Q ( x ) (Q can be conservative or characteristic vari-ables) over cell Ω , a polynomial P r ( x ) with degrees 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 ) , (32)where Q is the cell averaged value of Q ( x ) over cell Ω , k = ( k , k , k ), | k | = k + k + k .13he 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. (33) Remark 4.
The trilinear interpolation is used to describe a given hexahedron Ω with copla-nar or non-coplanar vertexes, X ( ξ, η, ζ ) = (cid:88) i =1 x i φ i ( ξ, η, ζ ) , where ( ξ, η, ζ ) ∈ [ − / , / , x i is the locations of the ith vertex and φ i is the base functionas follows φ = 18 (1 − ξ )(1 − η )(1 − ζ ) ,φ = 18 (1 − ξ )(1 − η )(1 + 2 ζ ) ,φ = 18 (1 − ξ )(1 + 2 η )(1 − ζ ) ,φ = 18 (1 − ξ )(1 + 2 η )(1 + 2 ζ ) ,φ = 18 (1 + 2 ξ )(1 − η )(1 − ζ ) ,φ = 18 (1 + 2 ξ )(1 − η )(1 + 2 ζ ) ,φ = 18 (1 + 2 ξ )(1 + 2 η )(1 − ζ ) ,φ = 18 (1 + 2 ξ )(1 + 2 η )(1 + 2 ζ ) . (34) Then the integration of monomial in Eq. (33) can be given by ˝ Ω 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 ζ. (35) It can be evaluated numerically as ˝ Ω x k y k z k d x d y d z = (cid:80) Ml,m,n =1 ω l,m,n x k y k z k ( ξ l , η m , ζ n ) (cid:12)(cid:12)(cid:12) ∂ ( x,y,z ) ∂ ( ξ,η,ζ ) (cid:12)(cid:12)(cid:12) ( ξ l ,η m ,ζ n ) ∆ ξ ∆ η ∆ ζ, (36) where ω l,m,n is the quadrature weight for the Gaussian point ( ξ l , η m , ζ n ) . For the currentthird-order scheme, a × × Gaussian quadrature with fourth-order spatial accuracy isused with ω l,m,n = ( ) and ( ξ l , η m , ζ n ) = ( ± √ , ± √ , ± √ ) .4.1. Large stencil and sub-stencils In order to reconstruct the quadratic polynomial P ( x ) on Ω , the large stencil forreconstruction includes Ω and all its von-Neumann neighbors, Ω m , m = 1 , ..., igure 1: The compact stencils of cell Ω for HWENO-type reconstruction. Center: the large stencil. Others:the eight sub-stencils. averages of Q ( x ) and averaged derivatives of Q ( x ) over each cell are known.The following values are used to obtain P ( x ), • cell averages Q for cell 0, 1, 2, 3, 4, 5, 6, • cell averages of the x -direction partial derivative Q x for cell 1, 2, 3, 4, 5, 6, • cell averages of the y -direction partial derivative Q y for cell 1, 2, 3, 4, 5, 6, • cell averages of the z -direction partial derivative Q z for cell 1, 2, 3, 4, 5, 6.The polynomial P ( x, y ) is required to exactly satisfy ˚ Ω m P ( x, y )d V = Q m | Ω m | , (37)where Q m is the cell averaged value over Ω m , m = 1 , ...,
6. Then we require the followingcondition satisfied in a least-squares 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 | , (38)where Q x i , i = 1 , , m in a global co-ordinate, respectively. On a regular mesh, the system has 24 independent equations. Theconstrained least-square method is used to solve the above linear system [19].15 emark 5. The constraints introduced in Eq. (37) can improve the linear stability of thereconstruction and reduce the numerical errors. The information from the cell-averagedvalues takes up a high proportion using the technique. Under uniform mesh, the coefficientsin Eq. (32) for cell Ω i,j,k are given in a concise and elegant form a , , = Q i +1 ,j,k − Q i − ,j,k h , a , , = Q i,j +1 ,k − Q i,j − ,k h , a , , = Q i,j,k +1 − Q i,j,k − h ,a , , = Q i +1 ,j,k + Q i − ,j,k − Q i,j,k h , a , , = Q i,j +1 ,k + Q i,j +1 ,k − Q i,j,k h ,a , , = Q i,j,k +1 + Q i,j +1 ,k − − Q i,j,k h ,a , , = ( Q x ) i,j +1 ,k − ( Q x ) i,j − ,k + ( Q x ) i +1 ,j,k − ( Q x ) i − ,j,k h ,a , , = ( Q x ) i,j,k +1 − ( Q x ) i,j,k − + ( Q x ) i,j +1 ,k − ( Q x ) i,j − ,k h ,a , , = ( Q x ) i,j,k +1 − ( Q x ) i,j,k − + ( Q x ) i +1 ,j,k − ( Q x ) i − ,j,k h , (39) The information from the derivatives only shows in the cross terms. In addition, it willreduce to a third-order scheme involving the cell-averaged values only in the 1-D case.
In order to deal with discontinuity, eight sub-stencils S j , j = 1 , ..., P j ( x ), 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 } ,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 } , There is always one sub-stencil in smooth region with the appearance of discontinuity nearany one of the interfaces of the target cell. The method in [46] can be used to obtain P j ( x, y ), which avoids the singularity caused by mesh irregularity, and the linear polynomialis expressed as P j ( x ) = ¯ Q + (cid:88) | k | =1 a j,k p k ( x ) . (40)Note that the choice of the large and sub-stencils is not unique.16 .2. Define the values of linear weights P ( x ) is written as P ( x ) = d [ 1 d P ( x ) − (cid:88) j =1 d j d P j ( x )] + (cid:88) j =1 d j P j ( x ) , where the linear weights are chosen as γ = 0 . , γ j = 0 . , j = 1 , ..., The smoothness indicators β j , j = 0 , ..., β j = r j (cid:88) | α | =1 | Ω | | α |− ˚ Ω (cid:0) D α P j ( x ) (cid:1) d V, where α is a multi-index and D is the derivative operator, r = 2, r j = 1 , j = 1 , ...,
8. Thesmoothness indicators in Taylor series at ( x , y ) have the order β = O {| Ω | [1 + O ( | Ω | )] } = O ( | Ω | ) = O ( h ) ,β j = O {| Ω | [1 + O ( | Ω | )] } = O ( | Ω | ) = O ( h ) , j = 1 , ..., . By using a similar technique [48], a global smoothness indicator σ can be defined σ = ( 18 (cid:88) | β − β j | ) = O ( | Ω | ) = O ( h ) , then the corresponding non-linear weights are given by ω j = d j (1 + σ(cid:15) + β j ) , j = 0 , ..., ,δ j = ω j (cid:80) l =0 ω l = d j + O ( h ) , (41)where (cid:15) takes 10 − to avoid zero in the denominator.The final reconstruction polynomial for the approximation of Q ( x ) yields R ( x ) = δ [ 1 d P ( x ) − (cid:88) j =1 d j d P j ( x )] + (cid:88) j =1 δ j P j ( x ) . (42)As a result, the non-linear reconstruction achieves a third-order accuracy R ( x ) = Q ( x ) + O ( h ). If any of these values yield negative density or pressure, the first-order reconstruction17s used instead. So all the desired quantities at Gaussian points can be fully determined as Q l,rp,k = R l,r ( x p,k ) , ( Q l,rx i ) p,k = ∂R l,r ∂x i ( x p,k ) . The reconstructions for the non-equilibrium states have the uniform order and can beused to get the equilibrium state directly, such as g c , g cx i by a suitable average of g l,r , g l,rx i . Thesimplest way is to use the arithmetic average, but it is only applicable for smooth flow. Tobe consistent 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Ξ . (43)This method has been used in an early version of second-order GKS [37] and validated inthe non-compact WENO5-GKS [13]. In this way, all components of the microscopic slopesin Eq.(25) have been fully obtained.
5. Numerical tests
In this section, numerical tests will be presented to validate the compact high-order GKS.For the inviscid flow, the collision time τ n is defined by τ n = (cid:15) ∆ 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, µ isthe dynamic viscosity coefficient, and p is the pressure at the cell interface. In smooth flowregion, it reduces to τ n = τ = µ/p . The ratio of specific heats takes γ = 1 .
4. The inclusionof the pressure jump term is to enlarge the collision time in the discontinuous region, wherethe numerical cell size is not enough to resolve the shock structure. It increases the non-equilibrium transport mechanism in the flux function to mimic the physical process in theshock layer.All reconstructions will be performed on the characteristic variables. Ghost cells aremainly adopted in the current scheme for boundary treatment. After we obtain the innerstate at the boundary, a ghost state can be assigned according to boundary condition,18nd the corresponding gas distribution function in Eq.(25) can be determined. A high-order boundary reconstruction is only applied to the test of subsonic flow passing througha circular cylinder. Explorations on the construction of accurate and stable condition oncurved boundaries will continue. The time step is determined by∆ t = C CF L
Min( ∆ r i || U i || + ( a s ) i , (∆ r i ) ν i ) , 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∆ r i = | Ω i | Max | Γ ip | . The CFL number is set as 0 . (a) 3-D sinusoidal wave propagation The advection of density perturbation is tested, and the initial condition is given asfollows ρ ( 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, the analyticsolution 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 CF L = 0 . L , L and L ∞ errors and the corresponding orders with linear andnon-linear Z-type weights at t = 2 are given in Tab.1 and Tab.2. The expected accuracy isconfirmed. mesh number L error Order L error Order L ∞ error Order5 Table 1: Accuracy test for the 3D sin-wave propagation by the linear third-order compact reconstruction.
CF L = 0 . (b) Subsonic flow past a circular cylinder esh number L error Order L error Order L ∞ error Order5 Table 2: Accuracy test for the 3D sin-wave propagation by the non-linear third-order compact HWENOreconstruction with d = 0 . d i = 0 . , i = 1 , ..., CF L = 0 . This 2-D test has been widely used to test the spatial accuracy for a high-order schemewith curved wall boundary [17, 20, 33].A circular cylinder is put in the center of the computational domain with a radius of r = 0 .
5. The concentric computational domain is bounded by a circle r out = 20. Foursuccessively refined meshes with 16 × ×
4, 32 × ×
4, 64 × ×
4, and 128 × × ρ, U, V, W, p ) ∞ = (1 , . , , , γ ) , with γ = 1 .
4. The periodic boundary condition is given in the Z-direction. It describesa subsonic inviscid flow at
M a ∞ = 0 .
38 passing through a cylinder. Ideally, the flow isisentropic with S ( x, y, t ) = S ∞ . Thus, an entropy error, defined as (cid:15) s = S − S ∞ S ∞ = pp ∞ ( ρ ∞ ρ ) γ − , is used for measuring the error of the numerical solution. The simulation is initialized withthe free stream value. The error is recorded when the flow gets to a steady state. To achievea third-order accuracy on the cylinder wall, a one-side compact stencil with six cells 0 , , ..., • cell averages ¯ W for cell 0, 1, 2, 3, 4, 5, • cell averages of the x -direction partial derivative ¯ W x for cell 0, 1, 2, 3, 4, 5, • cell averages of the y -direction partial derivative ¯ W y for cell 0, 1, 2, 3, 4, 5, • cell averages of the z -direction partial derivative ¯ W z for cell 0, 1, 2, 3, 4, 5,20nd a quadratic polynomial can determined in a least-squares sense. Moreover, the curvedboundary modification proposed in [17] is adopted, by adjusting the normal directions on theboundary at Gaussian point. A third-order convergence rate is achieved through the abovetreatment, as shown in Tab.3. The numerical result with low-order boundary reconstructionhas a more visible wake than that with the high-order one, as shown in Fig.3. XY Z x y
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1 0.8 0.6 0.4 0.200.20.40.60.81
Figure 2: The mesh sample with 128 × × x y
2 1.5 1 0.5 0 0.5 1 1.5 2 2 1.5 1 0.500.511.52 x y
2 1.5 1 0.5 0 0.5 1 1.5 2 2 1.5 1 0.500.511.52
Figure 3: Subsonic flow passing through a circular cylinder: 20 equidistant Mach contours from 0 .
038 to0 .
76. Mesh 128 × ×
4. Left: the second-order one-sided boundary reconstruction. Right: the third-ordercompact boundary reconstruction. The curved boundary treatment is applied in both cases.
This test case is used to test the capability of the proposed method in resolving low-speedviscous flow. The Reynolds number based on the diameter of the sphere is 118. In suchcase, the drag coefficient C D = 1 according to the experimental work [31].21 econd-order Third-ordermesh number L error Order L error Order64 ×
16 1.15e-03 1.40e-4128 ×
32 5.02e-5 4.51 1.01e-5 3.79256 ×
64 8.75e-6 2.52 1.19e-6 3.21
Table 3: Accuracy test for the subsonic flow past a circular cylinder by different boundary reconstructions.The curved boundary treatment is applied in both cases.
The far-field boundary condition is set around the outside of the domain, which has afree stream condition ( ρ, U, V, W, p ) ∞ = (1 , . , , , γ ) , with γ = 1 . M a ∞ = 0 . D = 1 and the first grid off the wall isabout 2 . × − D . The height of the grid grows from the wall with a constant ratio andstop at 20 D in the radial direction. The computational streamlines are compared with theexperimental streamlines [31], as shown in Fig.6. The shape of the steady separation bubbleagrees well with each other. The quantitative results are given in Tab.4, including the dragcoefficient C D , the separation angle θ , and the closed wake length L . The drag coefficientare defined as C D = F D ρ ∞ U ∞ S , where S = πD . The definition of θ and L is given in Fig.5. With similar DOFs, thecurrent compact scheme has the closest drag coefficient to the experiment data. The closedwake length has a visible difference with the experiment, but very close to the DG’s result[5]. It is probably due to the compressible effect since the experiment is conducted in thelow-speed water tank. Scheme DOF Cd θ LExperiment [31] – 1.0 151 1.07Current 524,288 1.009 125.1 0.95Implicit third-order DDG [5] 1,608,680 1.016 123.7 0.96Implicit fourth-order VFV [32] 458,915 1.014 – –Implicit third-order AMR-VFV [16] 621,440 1.016 – –Fourth-order FR [30] – – 123.6 1.04
Table 4: Quantitative comparisons among different compact schemes for the viscous flow past a sphere. igure 4: The mesh for viscous flow passing through a sphere with 32 × × × θ and the closed wake length L . igure 6: Comparison of streamlines between the experiment [31] and computation. The implicit large eddy simulation (ILES) of a three-dimensional Taylor-Green vortex[6] is conducted to validate the new compact GKS for nearly incompressible viscous flow.The initial flow field is given by U = V sin( xL ) cos( yL ) cos( zL ) ,U = − V cos( xL ) sin( yL ) cos( zL ) ,U = 0 ,p = p + ρ V
16 (cos( 2 xL ) + cos( 2 yL ))(cos( 2 zL ) + 2) , within a periodic cubic box − πL ≤ x, y, z ≤ πL . The density distribution is given bykeeping a constant temperature. In the computation, L = 1 , V = 1 , ρ = 1, and the Machnumber takes M = V /a s = 0 .
1, where a s is the sound speed. The characteristic convectivetime t c = L/V . The specific heat ratio γ = 1 . P r = 1. TwoReynolds number Re = 280 and 1600 are studied here. The linear weights of reconstructionand the smooth flux function are adopted in this case. The equilibrium state is obtainedby the arithmetic average of the non-equilibrium states to further reduce the numericaldissipations. The CFL is set as 0.3. 24wo quantities are investigated in the current study as the flow evolves in time. Thefirst one is the volume-averaged kinetic energy E k = 1 ρ Ω ˆ Ω ρ U · U dΩ , where Ω is the volume of the computational domain. ε k = − d E k d t = − d E k d t , By using the data E k ( t i ) , t i ∈ [0 , t stop ], the dissipation rate of the kinetic energy is given bya second-order interpolation ε k ( t i ) = − E k ( t i +1 ) − E k ( t i − ) t i +1 − t i − . The numerical results are compared with the reference DNS data in [6] and the non-compactfifth-order GKS [13]. The iso-surfaces of Q criterion colored by Mach number at t = 5 and10 for both Reynolds number are shown in Fig.7 and Fig.8. With the time increment, thevortex structures become denser and smaller. The case with higher Reynolds number hasricher structures, which requires high resolution for a numerical scheme.For Re = 280, the quantitative result agrees nicely with the reference data under acoarse mesh 96 , as shown in Fig.9. The new third-order compact scheme even shows betterability against the traditional fifth-order GKS with the same mesh, seen as Fig.9(c). Thetime history of the normalized volume-averaged kinetic energy and dissipation rate under Re = 1600 with 128 and 196 mesh points are presented in Fig.10. A detailed zoom-inplot is given in Fig.11. The compact scheme is capable to capture the complicated vortexstructure as the non-compact one under the same mesh. The one-dimensional Sod test case is performed in three-dimensional simulation. Theinitial condition is given by( ρ, U, V, W, p ) = (cid:40) (1 , , , , , < x < . , (0 . , , , , . , . ≤ x < , where 100 × × , × [0 , × [0 , t = 0 .
2. Non-reflection boundary condition is adopted at theleft and right boundaries of the computational domain, and periodic boundary conditionis adopted at the rest of the boundaries. The 3-D plot of density distribution in Fig.12(a)shows the uniformity in the flow distributions along x direction. The density distributionat the center horizontal line is also extracted, as shown in Fig.12(b). For this mild case,the numerical result agrees well with the exact solution. However, the contact discontinuity25 a) Re=280, t=5 (b) Re=280, t=10Figure 7: Taylor-Green vortex: the iso-surfaces of Q criterion colored by Mach number at time t = 5, 10 forRe = 280. Cell number = 128 . The x-y plane is shown.(a) Re=1600, t=5 (b) Re=1600, t=10Figure 8: Taylor-Green vortex: the iso-surfaces of Q criterion colored by Mach number at time t = 5, 10 forRe = 1600. Cell number = 196 . The x-y plane is shown. k i n e t i c e n e r g y Refweno5 linear 128^396^3128^3 (a) E k t d k Refweno5 linear 128^396^3128^3 (b) ε k t d k Refweno5 linear 128^396^3128^3 (c) Local enrlargement of ε k Figure 9: Taylor-Green vortex: Re=280. The time history of kinetic energy and its dissipation rate. t k i n e t i c e n e r g y Refweno5 linear 196^3128^3196^3 (a) E k t d k Refweno5 linear 196^3128^3196^3 (b) ε k Figure 10: Taylor-Green vortex: Re=1600. The time history of kinetic energy and its dissipation rate. d k Refweno5 linear 196^3196^3 (a) E k t d k Refweno5 linear 196^3196^3 (b) ε k t d k
11 12 13 14 150.0040.0060.0080.01
Refweno5 linear 196^3196^3 (c) ε k Figure 11: Taylor-Green vortex: Re=1600. The local enlargement of the time history of kinetic energydissipation rate.
28s not as sharp as the traditional fifth-order WENO-GKS. It is consistent with Remark 5,which shows the scheme will recover to a third-order non-compact GKS under 1-D smoothcase. (a) x d e n s i t y ExactCurrent (b)Figure 12: Sod problem: 3-D view of density distribution and its extraction along the center x-horizontalline. The mesh size along the x-direction is 1 / As an extension of the Sod problem, the spherical explosion test problem is considered.The initial conditions are given by( ρ, U, V, W, p ) = (cid:40) (1 , , , , , < r < . , (0 . , , , , . , . ≤ r < , where r = (cid:112) ( x + y + z ). The solution contains a spherical shock wave and a contactsurface traveling away from the center and a spherical rarefaction wave moving towards theorigin (0 , , x, y, z ) ∈ [0 , × [0 , × [0 ,
1] and the uniform grid with dx = dy = dz = 1 /
100 is used. Thesymmetric boundary conditions are imposed on the planes x = 0, y = 0, and z = 0, whilethe outflow boundary conditions are imposed on the planes x = 1, y = 1, z = 1. The densityand pressure profiles along different radial directions at t = 0 .
25 are given in Fig.13. Thecompact scheme resolves the wave profiles crisply. Slightly overshoot can be observed at thefront of the rarefaction wave.
To validate the robustness of the current scheme with non-coplanar meshes, a high-speedinviscid flow pasting through a sphere is tested. The structured grids with six blocks are29 d e n s i t y Refy=0,z=0x=0,z=0x=0,y=0 r p r ess u r e Refy=0,z=0x=0,z=0x=0,y=0
Figure 13: Three-dimensional explosion problem with 100 cells. Left and middle: the density and pressureprofiles extracted along the y = 0 , z = 0, x = 0 , z = 0, and x = 0 , y = 0 lines. Right: 3-D view of pressuredistribution. used in the computation, as shown in Fig.14. The diameter of the sphere is D = 1 and thefirst grid off the wall is 1 × − D . The slip boundary condition is imposed on the surfaceof the sphere. The outer domain is around 3 . D in the radial direction. The supersonicinlet/outlet is adopted on the outside boundary, which is set according to the angle betweenthe outer-pointing normal vector of each boundary interface and the incoming velocity. Asupersonic flow with M a = 3 is tested first. The simulation starts with the free stream flowcondition. The computation is directly started. The pressure and Mach distributions atsteady state are shown in Fig.15. Then a hypersonic flow with
M a = 5 is tested. A primaryflow field calculated by the first-order kinetic method [37] is used as the initial field. Thenumerical results are shown in Fig.16. The shock is captured sharply and the carbunclephenomenon does not appear in both cases. The result is essentially axis-symmetric. Theasymmetric pattern can be observed at the leeward side of the sphere. A nearly vacuumstate forms in this region and the HWENO reconstruction will easily give a negative densityor pressure and then reduce to a low-order reconstruction.
Figure 14: The mesh for inviscid flow passing through a sphere with 32 × × × igure 15: Supersonic inviscid flow passing through a sphere: Ma=3. Left: the pressure and Mach distribu-tions for the x-z and y-z planes. Right: the pressure and Mach distributions for the y-z plane. The sphereis colored by the temperature. CFL=0.5. The maximum Mach number M a
Max ≈ . igure 16: Hypersonic inviscid flow passing through a sphere: Ma=5. Left: the pressure and Mach distribu-tions for the x-z and y-z planes. Right: the pressure and Mach distributions for the y-z plane. The sphereis colored by the temperature. CFL=0.5. The maximum Mach number M a
Max ≈ . .7. Compressible isentropic turbulence A decaying homogeneous isotropic compressible turbulence is computed within a squarebox defined as − π ≤ x, y, z ≤ π , and the periodic boundary conditions are used in alldirections [27]. Given spectrum with a specified root mean square U (cid:48) U (cid:48) = (cid:28) U · U (cid:29) / , a divergence-free random velocity field U is initialized, where (cid:28) ... (cid:29) is a volume averageover the whole computational domain. The specified spectrum for velocity is given by E ( k ) = A k exp( − k /k ) , where A is a constant to set initial kinetic energy, k is the wavenumber, k is the wavenumberat spectrum peaks. The initial volume averaged turbulent kinetic energy K and the initiallarge-eddy-turnover time τ are given by K = 3 A √ πk , τ = (cid:114) A (2 π ) / k − / . The Taylor micro-scale λ and corresponding Reynolds number Re λ and are given by λ = ( U (cid:48) ) (cid:28) ( ∂ U ) (cid:29) , Re λ = (cid:28) ρ (cid:29) U (cid:48) λ (cid:28) µ (cid:29) = (2 π ) / ρ µ (cid:112) A k / . The turbulence Mach number
M a t is defined as M a t = √ U (cid:48) (cid:28) c s (cid:29) = √ U (cid:48) √ γT . The dynamic viscosity is determined by the power law µ = µ (cid:0) TT (cid:1) . , where µ and T can be determined from Re λ and M a t with initialized U (cid:48) and ρ = 1.A fixed Re λ = 72 is investigated in the current work. The other parameters, i.e., A =1 . × − , k = 8, are chosen according to [4]. The random initial flow field evolves complexlocal structures, as shown in the direct numerical simulations by the conventional WENO-GKS [4]. When M a t = 0 .
5, the flow is initially transonic, since the maximum Mach numberin the flow filed is about three times of the initial turbulent Mach number. The puresmooth GKS solver and the WENO-AO reconstruction with linear weights can be usedunder such mild case to achieve a higher resolution. Uniform meshes with 64 and 128 cellsare used in the simulations. Initially the cell-averaged slopes are constructed automaticallyby setting the first explicit time step ∆ t = 0. The kinetic energy, root-mean-square of33ensity fluctuation, and the skew factor are calculated, which are given by K ( t ) = 12 (cid:28) ρ U · U (cid:29) , ρ rms ( t ) = (cid:112) (cid:28) ( ρ − ρ ) (cid:29) , S u ( t ) = (cid:88) i (cid:28) ( ∂ i u i ) (cid:29)(cid:28) ( ∂ i u i ) (cid:29) / . Since the cell-averaged W x i is stored in each cell, the ∂ i u i can be calculated convenientlythrough the chain rule. The time history of normalized kinetic energy K ( t ) /K , normalizedroot-mean-square of density fluctuation ρ rms ( t ) /M a t agree well with the reference data andthe traditional WENO-GKS with 64 cells, as shown in Fig.17, and Fig.18. For the higher-order moment, the skew factor needs finer meshes to resolve, as shown in Fig.19.When the Mach number gets higher, the flow becomes supersonic, and the strongershocklets are generated, followed by complex shock-vortex interactions. It is nontrivial forhigh-order methods to survive under M a t ≥
1. In addition, it becomes more challengingunder a coarse mesh, since the discontinuities become stronger due to the limitation of themesh resolution. Thus, a series of turbulent Mach numbers have been chosen to test therobustness of the current scheme with the mesh 64 . The full GKS solver and the non-linearHWENO reconstruction are used. Considering the large velocity jump in the initial field, amodified τ is taken as τ = µp + (cid:80) δQ ∆ t , where Q means all five primitive variables, operator δQ = | Q l − Q r || Q l | + | Q r | +1 e − . The statistical quantities with respect to different Mach numbers arepresented in Fig.20. With the increase of M a t , the kinetic energy gets dissipated morerapidly. The visualized results, i.e., the iso-surfaces of Q criterion and the selected surfaceslice of Mach number distribution at z = − π are plotted in Fig.21. The complex vortexesand widespread shocklets can be observed clearly. t/ τ K (t) / K Refweno5 linear 128^364^3128^3 t/ τ K (t) / K Refweno5 linear 128^364^3128^3
Figure 17: Compressible isotropic turbulence: K ( t ) /K . Re λ = 72, M a t = 0 . / τ ρ r m s / M a t Refweno5 linear 128^364^3128^3 t/ τ ρ r m s / M a t Refweno5 linear 128^364^3128^3
Figure 18: Compressible isotropic turbulence: ρ rms ( t ) /M a t . Re λ = 72, M a t = 0 . t/ τ S u Refweno5 linear 128^364^3128^3 t/ τ S u Refweno5 linear 128^364^3128^3
Figure 19: Compressible isotropic turbulence: S u ( t ). Re λ = 72, M a t = 0 . / τ ρ r m s / M a t ma t =0.5ma t =0.8ma t =1.0ma t =1.2 t/ τ K (t) / K ma t =0.5ma t =0.8ma t =1.0ma t =1.2 Figure 20: Compressible isotropic turbulence: comparison with different
M a t numbers by the new compactGKS. Mesh: 64 . CFL=0.5.Figure 21: Compressible homogeneous turbulence with M a t = 1 at time t/τ = 1. Left: iso-surfaces of Qcriterion colored by Mach number; right: the Mach number distribution with z = − π . Mesh: 64 . Conclusion In this paper, a compact high-order gas-kinetic scheme for three-dimensional flow sim-ulation is presented. The distinguishable feature of the scheme is that the high-order GKSevolution model at a cell interface provides not only the fluxes, but also the time accurateflow variables. As a result, based on the cell interface values the first-order spatial deriva-tives of flow variables inside each control volume can be directly obtained through Gauss’stheorem at the next time level. The way for the updates of gradients in GKS is differentfrom the weak formulation for the updates of similar degree of freedom in the compact DG-type methods. Therefore, equipped with the cell-averaged values and their gradients, a newHWENO reconstruction with less sub-stencils and all positive weights has been designed forthe initial data reconstruction in the scheme. At the same time, the multi-stage and multi-derivative technique is used as a time marching strategy in the scheme, which subsequentlyleads to a high efficiency in comparison with the traditional Runge-Kutta method for thesame temporal accuracy. The compact third-order scheme in 3D can use a relative largeCFL number in the computations, and shows similar resolution as the fifth-order gas-kineticscheme with the non-compact WENO reconstruction. Overall, the time accurate evolutionmodel, the HWENO reconstruction, and the MSMD time marching technique make the finalscheme accurate, robust, and efficient for the compressible flow simulations with smooth anddiscontinuous solutions.Although the current scheme is constructed on structured mesh, it can be directly ex-tended to unstructured one. The compact third-order and fourth-order HGKS on three-dimensional unstructured mesh is on the development, which further enlarges the applicableregime of the high-order gas-kinetic schemes for flow computation with complex geome-try. The development of high-order compact scheme with implicit and other accelerationtechniques is on the investigation as well for the steady state solution.
Acknowledgment
The authors would like to thank Dr. Liang Pan for helpful discussion. The currentresearch is supported by National Numerical Windtunnel project, Hong Kong research grantcouncil 16206617, and National Science Foundation of China 11772281, 91852114.
ReferencesReferences [1] Dinshaw S Balsara, Sudip Garain, and Chi-Wang Shu. An efficient class of WENO schemes withadaptive order.
Journal of Computational Physics , 326:780–804, 2016.[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] Guiyu Cao, Hualin Liu, and Kun Xu. Physical modeling and numerical studies of three-dimensionalnon-equilibrium multi-temperature flows.
Physics of Fluids , 30(12):126104, 2018.[4] Guiyu Cao, Liang Pan, and Kun Xu. Three dimensional high-order gas-kinetic scheme for supersonicisotropic turbulence i: criterion for direct numerical simulation.
Computers & Fluids , 192:104273, 2019.
5] 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 Computa-tional Physics , 21(5):1231–1257, 2017.[6] James DeBonis. Solutions of the Taylor-Green vortex problem using high-resolution explicit finitedifference methods. In , page 382, 2013.[7] Zhifang Du and Jiequan Li. A Hermite WENO reconstruction for fourth order temporal accurateschemes based on the GRP solver for hyperbolic conservation laws.
Journal of Computational Physics ,355:385–396, 2018.[8] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor. Strong stability-preserving high-order time dis-cretization methods.
SIAM review , 43(1):89–112, 2001.[9] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes.
Journal of Computational Physics , 150(1):97–127, 1999.[10] Ht T Huynh. A flux reconstruction approach to high-order schemes including discontinuous Galerkinmethods.
AIAA paper , 4079:2007, 2007.[11] Xing Ji.
High-order non-compact and compact gas-kinetic schemes . PhD thesis, Hong Kong Univeristyof Science and Technology, 2019.[12] Xing Ji, Liang Pan, Wei Shyy, and Kun Xu. A compact fourth-order gas-kinetic scheme for the Eulerand Navier-Stokes equations.
Journal of Computational Physics , 372:446 – 472, 2018.[13] 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.[14] Xing Ji, Fengxiang Zhao, Wei Shyy, and Kun Xu. A family of high-order gas-kinetic schemes andits comparison with Riemann solver based high-order methods.
Journal of Computational Physics ,356:150–173, 2018.[15] 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.[16] PAN Jianhua, WANG Qian, Yusi Zhang, and REN Yuxin. 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.[17] Lilia Krivodonova and Marsha Berger. High-order accurate implementation of solid wall boundaryconditions in curved geometries.
Journal of computational physics , 211(2):492–512, 2006.[18] 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.[19] Wanai Li.
Efficient implementation of high-order accurate numerical methods on unstructured grids .Berlin, Heidelberg: Springer, 2014.[20] Hong Luo, Joseph D Baum, and Rainald L¨ohner. On the computation of steady-state compressible flowsusing a discontinuous galerkin method.
International Journal for Numerical Methods in Engineering ,73(5):597–623, 2008.[21] Liang Pan, Junxia Cheng, Shuanghu Wang, and Kun Xu. A two-stage fourth-order gas-kinetic schemefor compressible multicomponent flows.
Communications in Computational Physics , 22(4):1123–1149,2017.[22] Liang Pan and Kun Xu. A third-order compact gas-kinetic scheme on unstructured meshes for com-pressible Navier–Stokes solutions.
Journal of Computational Physics , 318:327–348, 2016.[23] Liang Pan and Kun Xu. Two-stage fourth-order gas-kinetic scheme for three-dimensional Euler andNavier-Stokes solutions.
International Journal of Computational Fluid Dynamics , 32:395–411, 2018.[24] 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.[25] Jianxian Qiu and Chi-Wang Shu. Hermite WENO schemes and their application as limiters forRunge–Kutta discontinuous Galerkin method: one-dimensional case.
Journal of Computational Physics , Handbook of Numerical Analysis , volume 17, pages 147–171. Elsevier, 2016.[27] Ravi Samtaney, Dale I Pullin, and Branko Kosovi´c. Direct numerical simulation of decaying compress-ible turbulence and shocklet statistics.
Physics of Fluids , 13(5):1415–1430, 2001.[28] Chi-Wang Shu. High order weighted essentially nonoscillatory schemes for convection dominated prob-lems.
SIAM review , 51(1):82–126, 2009.[29] Chi-Wang Shu. High order weno and dg methods for time-dependent convection-dominated pdes: Abrief survey of several recent developments.
Journal of Computational Physics , 316:598–613, 2016.[30] Y. Sun, Z. J. Wang, and Y. Liu. High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids.
Communications in Computational Physics ,2(2):310–333, 2007.[31] 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.[32] Qian Wang.
Compact High-Order Finite Volume Method on Unstructured Grids . PhD thesis, TsinghuaUniversity, 6 2017.[33] Qian Wang, Yu-Xin Ren, and Wanai Li. Compact high order finite volume method on unstructuredgrids II: Extension to two-dimensional Euler equations.
Journal of Computational Physics , 314:883–908,2016.[34] Zhijian Wang, Krzysztof Fidkowski, R´emi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary,Herman Deconinck, Ralf Hartmann, Koen Hillewaert, Hung T Huynh, et al. High-order CFD methods:current status and perspective.
International Journal for Numerical Methods in Fluids , 72(8):811–845,2013.[35] 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.[36] Bin Xie, Xi Deng, Ziyao Sun, and Feng Xiao. A hybrid pressure–density-based mach uniform algorithmfor 2d euler equations on unstructured grids by using multi-moment finite volume method.
Journal ofComputational Physics , 335:637–663, 2017.[37] Kun Xu. Gas-kinetic schemes for unsteady compressible flow simulations.
Lecture series-van KaremanInstitute for fluid dynamics , 3:C1–C202, 1998.[38] 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.[39] Kun Xu.
Direct modeling for computational fluid dynamics: construction and application of unifiedgas-kinetic schemes . World Scientific, 2014.[40] Kun Xu, Meiliang Mao, and Lei Tang. A multidimensional gas-kinetic BGK scheme for hypersonicviscous flow.
Journal of Computational Physics , 203(2):405–421, 2005.[41] Xiaoquan Yang, Jian Cheng, Hong Luo, and Qijun Zhao. A reconstructed direct discontinuous galerkinmethod for simulating the compressible laminar and turbulent flows on hybrid grids.
Computers &Fluids , 168:216–231, 2018.[42] 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.[43] Meilin Yu, Zhijian Wang, and Yen Liu. On the accuracy and efficiency of discontinuous Galerkin,spectral difference and correction procedure via reconstruction methods.
Journal of ComputationalPhysics , 259:70–95, 2014.[44] Fengxiang Zhao, Xing Ji, Wei Shyy, and Kun Xu. An acoustic and shock wave capturing compacthigh-order gas-kinetic scheme with spectral-like resolution. arXiv preprint arXiv:2001.01570 , 2019.[45] Fengxiang Zhao, Xing Ji, Wei Shyy, and Kun Xu. Compact higher-order gas-kinetic schemes withspectral-like resolution for compressible flow simulations.
Advances in Aerodynamics , 1(1):13, 2019.[46] Fengxiang Zhao, Liang Pan, and Shuanghu Wang. Weighted essentially non-oscillatory scheme onunstructured quadrilateral and triangular meshes for hyperbolic conservation laws.
Journal of Compu-tational Physics , in press, 2018.
47] Jun Zhu and Jianxian Qiu. Hermite WENO schemes and their application as limiters for Runge-Kuttadiscontinuous Galerkin method, III: unstructured meshes.
Journal of Scientific Computing , 39(2):293–321, 2009.[48] Jun Zhu and Jianxian Qiu. New finite volume weighted essentially nonoscillatory schemes on triangularmeshes.
SIAM Journal on Scientific Computing , 40(2):A903–A928, 2018.[49] 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.