A high order kinetic flow solver based on flux reconstruction framework
aa r X i v : . [ phy s i c s . c o m p - ph ] J u l A high order kinetic flow solver based on flux reconstructionframework
Ji Li, ∗ Chengwen Zhong, † and Sha Liu ‡ National Key Laboratory of Science and Technology on Aerodynamic Design and Research,Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China. (Dated: August 7, 2020)
Abstract
The goal of this paper is to develop a high order numerical method based on Kinetic Inviscid Flux(KIF) method and Flux Reconstruction (FR) framework. The KIF aims to find a balance betweenthe excellent merits of Gas-Kinetic Scheme (GKS) and the lower computational costs. The ideaof KIF can be viewed as an inviscid-viscous splitting version of the gas-kinetic scheme, and Shuand Ohwada have made the fundamental contribution. The combination of Totally ThermalizedTransport (TTT) scheme and Kinetic Flux Vector Splitting (KFVS) method are achieved in KIF.Using a coefficient which is related to time step δt and averaged collision time τ , KIF can adjust theweights of TTT and KFVS flux in the simulation adaptively. By doing the inviscid-viscous splitting,KIF is very suitable and easy to integrate into the existing framework. The well understood FRframework is used widely for the advantages of robustness, economical costs and compactness. Thecombination of KIF and FR is originated by three motivations. The first purpose is to develop ahigh order method based on the gas kinetic theory. The second reason is to keep the advantagesof GKS. The last aim is that the designed method should be more efficient. In present work, weuse the KIF method to replace the Riemann flux solver applied in the interfaces of elements. Thecommon solution at the interface is computed according to the gas kinetic theory, which makes thecombination of KIF and FR scheme more reasonable and available. The accuracy and performanceof present method are validated by several numerical cases. The Taylor-Green vortex problem hasbeen used to verify its potential to simulate turbulent flows. ∗ Ji Li: [email protected] † Corresponding author: Chengwen Zhong, [email protected] ‡ Sha Liu: [email protected] . INTRODUCTION The Gas-Kinetic Scheme (GKS) developed by Xu [1, 2] is based on the idea of Bhatnagaret al. [3]. In recent years, GKS is on the way to become the preferred numerical method inthe fluid dynamics. Compared with traditional Navier-Stokes numerical method, GKS is ofhigh spatial and temporal accuracy. The advantages of GKS have been recognized in thesimulation of turbulent flows [4–8], shock-boundary interaction, hypersonic flows [2, 9] andnon-equilibrium simulations [10, 11]. A series studies based on the GKS have been advanced,such as immersed boundary method [12–14], implicit temporal marching [15], and dual-timestrategy [16] for unsteady flows.In the field of Computational Fluid Dynamics (CFD), a numerical algorithm can beclassified as a low or high order method according to the numerical accuracy approached.For the features of robustness and economical costs, the low order methods are still popularin the aeronautical industry. Compared with the low order method, the high order methodis more accurate, which has the potential of providing more details of the flow fields [17, 18].However, using the high order method in the real industry is still a challenge, which attractsthe interests of many researchers from all over the world. The development of high orderGKS can be traced back to the study of Q. Li [19]. J. Luo [20], G. Zhou [21], L. Pan [22–24],X. Ji [25, 26], F. Zhao [27] have made the great contribution to the high order gas-kineticscheme.Although high order gas-kinetic scheme has been studied quite well, the development ofhigh order gas-kinetic scheme is never stopped. In present work, we focus on the combinationof Kinetic Inviscid Flux (KIF) and Flux Reconstruction method (FR). The combination ofKIF and FR is originated by three motivations. The first purpose is to develop a high ordermethod based on the gas kinetic theory. The second reason is to keep the advantages ofGKS. The last aim is that the designed method should be more efficient. In order to find abalance between the advantages of gas-kinetic scheme and lower computational costs, KIF isproposed by S. Liu[28]. The KIF scheme is a combination of Totally Thermalized Transport(TTT) scheme and Kinetic Flux Vector Splitting method (KFVS). The TTT scheme, whichdoes not introduce extra artificial viscosity in smooth flow area, can approach the boundarylayer accurately. TTT has the property similar to central schemes, so it also can not capturethe discontinuity properly. The KFVS is a shock capturing scheme with good robustness.2he combination is a good idea, which means that we use TTT scheme where the flow issmooth and use KFVS method where discontinuity exists. The kernel of KIF method is toadjust the weights of TTT and KFVS in the simulation automatically.To develop a high order kinetic flux solver, it is critical to adopt the advantages fromthe traditional high order method based on the Navier-Stokes equation. Over the last20 years, the high order numerical method is one of the research hotspots in the field ofCFD. A great many of researchers have devoted their attentions to such a challenge, and alarge numbers of high order numerical methods have been developed under the frameworksof finite volume method, finite difference method and finite element method et al. Someof the schemes are of notable features and have been used widely. For example, k-exactmethod [29], Essentially Non-Oscillatory (ENO) method [30–33], Weighted ENO (WENO)method [34–36], Discontinuous Galerkin (DG) method [37], radial basis function method [38],Compact Least-Squares (CLS) reconstruction method [39, 40] and variational reconstructionmethod [41]. An excellent review of the high order methods can be referred to the work ofZ. Wang [42].Flux reconstruction method, first proposed by H. T. Huynh [43, 44], is aimed to be morepopular in both of the research and real industry fields. The designed features of robustness,economical costs and compactness make it well understood and available. A particular FRscheme depends on three factors [45], namely the distribution of solution points, the Riemannflux solver applied at the interfaces, and the choice of the correction functions G and H .It has been proved that the flux reconstruction method can recover the simplified DG andstaggered grid scheme with specific factors, and the conservation also has been proved inRef. [43]. Based on the study of Jameson [46], a class of energy stable flux reconstructionmethod was proposed by Vincent, Castonguay and Jameson (VCJH) [45]. And then, theVCJH scheme was used for triangular elements [47]. Up to now, the VCJH correct functionhas played an important role in the FR framework.In present work, the combination of KIF and FR is achieved by (a) replacing the Riemannsolver applied on the interface of elements with KIF, (b) using the gas kinetic theory tocompute the common solution on the interface, and (c) implementing the inviscid-viscoussplitting strategy in the simulation. The present paper is organized as follows. In Sec. II,KIF method and the flux reconstruction framework are introduced. Several numerical testsare set up in the Sec. III, and the numerical accuracy of present method is validated. The3ast section of paper is a short conclusion. II. NUMERICAL METHOD
The FR method, which takes advantages of DG and staggered grid scheme [48, 49], isfirst developed by H. T. Huynh [43, 44]. It focuses on the features of robustness, economicalcosts and compactness. Benefiting from the merits of well understood and available, theflux reconstruction method has attracted a great many of attentions. In this section, thehigh order kinetic flux solver based on flux reconstruction framework will be introduced.
A. Governing equation
For one-dimensional problem, the Bhatnagar-Gross-Krook (BGK) model [3] in x-directionis f t + uf x = g − fτ , (1)where u is the particle velocity, f represents the gas distribution function, g denotes theequilibrium state approached by f , and τ is related to the averaged collision time. Theequilibrium state is known as a Maxwellian distribution reads g = ρ ( λπ ) K +12 e − (( u − U ) + ξ ) , (2)where ρ is the density, U is the macroscopic velocity. λ , which reads λ = m/ (2 kT ), is relatedto the temperature T of gas, m represents the molecular mass, and k denotes the Boltzmannconstant. The total number of degrees of freedom K in ξ equals to (5 − γ ) / ( γ −
1) + 2, γ is the ratio of specific heat, and ξ = P Ki =1 ξ i .According to the kinetic theory of gases, both the distribution function f and the equilib-rium state g are functions of particle velocities u , space x and time t . Taking the moments ofdistribution function f , the macroscopic conservative variable w can be obtained as follow w = ρρUE = Z ψ f d Ξ , ψ = (cid:18) , u, (cid:0) u + ξ (cid:1)(cid:19) T , (3)where d Ξ = du (cid:18) K Q i =1 dξ i (cid:19) . In order to obtain the spatial discretization in the flux recon-struction framework, we take moments of ψ in Eq. (1) and integrate it with d Ξ in phase4pace, Z ( f t + uf x ) ψ d Ξ = − Z f − gτ ψ d Ξ . (4)Using the compatibility condition Z g − fτ ψ d Ξ = 0 . (5)We can get the following formula w t + G x = 0 , G = Z uf ψd Ξ , (6)where G is the flux corresponding to conservative variables w along the x-direction. Then,we can solve the Eq. (6) within the flux reconstruction framework, and the flux G can becomputed using KIF. B. Kinetic Inviscid Flux
The Riemann flux solver employed at the interfaces is one of the three critical factorsof flux reconstruction method. In the present work, we implement the kinetic inviscidflux to determine the common flux at the interface of elements. The motivation of ourwork is to reach a compromise between good performance of gas-kinetic scheme and lowercomputational costs. The KIF is a kind combination of TTT scheme and KFVS method.TTT scheme has been discussed by Xu in Ref. [50]. The first step of TTT scheme is toget the Maxwell distribution on both sides of the surface, g l = [ ρ ( λπ ) K +12 e − (( u − U ) + ξ ) ] l , g r = [ ρ ( λπ ) K +12 e − (( u − U ) + ξ ) ] r . (7)The subscripts l and r represent the state based on the macroscopic variables from left andright sides of interface respectively. The second step is to make particles crossing the interfacecollide sufficiently. The new Maxwell distribution g is assumed to have the following form g = [ ρ ( λπ ) K +12 e − (( u − U ) + ξ ) ] , w = Z g ψd Ξ = Z ((1 − H ( x )) g l + H ( x ) g r ) ψd Ξ , (8)where H ( x ) is the Heaviside function H ( x ) = , x < , , x ≥ . (9)5inally, the flux of TTT scheme reads F T T T = Z ug ψd Ξ . (10)The TTT scheme has the similar property to central scheme, which does not introduceextra artificial viscosity in smooth flow area and can capture the boundary layer accurately.However, it can not deal with shock wave, because of lacking artificial viscosity.Equilibrium Flux Method (EFM) [51] and KFVS [52] are similar, and we just call it KFVShere. KFVS is another kinetic scheme and is a shock capturing method. After getting theMaxwell distribution beside the interface using Eq. (7), KFVS gives flux by calculatingparticles across the interface as F KF V S = Z u> ug l ψd Ξ + Z u< ug r ψd Ξ . (11)KFVS has the properties as good robustness and positivity-preserving [53], which make thescheme suitable for capturing the discontinuity. However, it introduces enormous artificialviscosity, and the essential problem, which is analyzed in Ref. [50], is that the equationsolved at interface is the collisionless Boltzmann equation.Combination is a good idea, which means that we use TTT in the smooth flow area andKFVS in the flow field where shock wave exists. The idea of KIF can be viewed as aninviscid-viscous splitting version of the gas-kinetic scheme, and Shu [54] and Ohwada [55]have made the fundamental contribution. The most significant difference between KIF [28]and the works of Shu and Ohwada is the weight of TTT and KFVS. Ref. [28] adopted thephilosophy of direct modeling [56], and constructed two kinds of KIF method (namely KIF1method from GKS strategy [1] and KIF2 from DUGKS strategy [57]). In present work, weadopt the KIF1 in the simulations. The KIF1 can be expressed as F = { τδt (1 − e − δt/τ ) } F KF V S + { − τδt (1 − e − δt/τ ) } F T T T , (12) δt = rτ = τ max Ω h | p l − p r || p l + p r | , max( M a l , M a r ) i . (13) δt is the observation time scale and is measured in mean collision time (or the relaxationtime τ ) in discontinuities. Since KIF is an inviscid-viscous splitting version of GKS, the6iscous fluxes can be computed using the traditional central viscosity scheme. The detailsof KIF1 can be referred to Ref. [28].In another point of view, KIF is a kind of balance between kinetic scheme and traditionalmacroscopic numerical method. In recent years, kinetic schemes have a significant develop-ment [57–61], which mainly aim at nonequilibrium flow. With a view at equilibrium state,KIF replaces the complicated nonequilibrium part by traditional central viscosity scheme.One motivation is to find a balance between advantages and efficiency, while another is tobe suitable and easy to integrate into the existing framework. C. Spatial discretization in flux reconstruction framework
In this section, the FR framework used to solve Eq. (6) is introduced. For the one-dimensional problem, the computational domain Ω can be divided into N subdomainsΩ = { Ω i | i = 0 , , · · · , N − } , Ω i = [ x i , x i +1 ] , x < x < · · · < x N . (14)Within the element Ω i , the solution points are set as x i,k ( k = 0 , , , · · · , P ). It is obviousthat the number of solution points within a standard elements is P + 1, and P is relatedwith the accuracy order of the numerical method. The set x i,k can be chosen as Gauss,Radau, Lobatto or equidistant points. It has been proved in Refs. [43, 44] by H. T. Huynhthat Fourier stability and accuracy analysis of FR framework are independent of the typeof solution points.Since dealing with every elements Ω i is very tedious, all the element Ω i should be mappedinto the same standard element Ω s = { ξ | ξ ∈ [ − , } to simplify the implementation of thealgorithm. The mapping function θ ( ξ ) can be expressed as x = θ i ( ξ ) = (cid:18) − ξ (cid:19) x i + (cid:18) ξ (cid:19) x i +1 . (15)In order to be consistent with the existed literature of FR framework, the denotation offlux G in Eq. (6) is replaced by f . With the mapping expressed as Eq. (15), the evolutionof macroscopic variables w within each Ω i can be transformed as the Eq. (16) within thestandard element. ˆ w t + 1 J n ˆ f ξ = 0 , (16)7 IG. 1. The discontinuous solution polynomials and interface common solutions within elementsΩ i − , Ω i and Ω i +1 . where ˆ w = w ( θ i ( ξ ) , t ) in Ω i , (17)ˆ f = f ( θ i ( ξ ) , t ) in Ω i , (18) J n = ∂x∂ξ in Ω i . (19)The FR framework for solving the Eq. (16) within the standard element Ω s consists ofseven subsequent steps. In the first step, the solution polynomial ˆ w δi ( x ) can be obtainedthrough the macroscopic variable ˆ w i,k at the solution points ξ k , ˆ w δi ( x ) = k = P X k =0 ˆ w i,k φ k ( ξ k ) , (20)where the symbol δ denotes the solution polynomial is always discontinuous at the elementinterface. φ k ( ξ k ) is the 1D Lagrange polynomial equal to 1 at the kth solution point and 0at the others, φ k ( ξ k ) = P Y l =0 ,l = k ξ − ξ l ξ k − ξ l . (21)Fig. 1 shows the solution polynomials at element Ω i and the neighbors in the physicalspace. Take the interface x i +1 / as an example, w δi ( x x +1 / ) and w δi +1 ( x x +1 / ) represent themacroscopic variables at the interface from Ω i (left) and Ω i +1 (right) respectively. Generallyspeaking, w δi ( x x +1 / ) and w δi +1 ( x x +1 / ) are not equal. Since ξ belongs to the interval [ − , s , w δi ( x x +1 / ) and w δi +1 ( x x +1 / ) equal to ˆ w δ i (1) and ˆ w δ i +1 ( − ˆ w δ i (1) does not equal to ˆ w δi +1 ( −
1) in most of cases, which isthe “Jump” or “Discontinuous” at the boundaries of element.8n the second step, it is turn to determine the common solution ˆ w CI at the boundaries ofstandard element, i.e. ξ = ±
1. The common solution ˆ w CI at interface is used to make thesolution within the standard element to feel the effect of the boundaries, so the superscript C also has the meanings “Corrected” and “Continuous”. In present scheme, the commonsolution ˆ w CI is computed using the following expression ˆ w CIi +1 / = Z u> g l ψ d Ξ + Z u< g r ψ d Ξ , (22)where g l and g r , which are corresponding to ˆ w i (1) and ˆ w i +1 ( − ˆ w δi (1) and ˆ w δi +1 ( −
1) can be obtainedeasily using Eq. (20). The demonstration of ˆ w CI is shown in the Fig. 1.It must be noticed that if the jumps at boundaries of element are ignored, the solutionwithin the element can not feel the effect of the boundaries and the evolution of scheme mustget an erroneous result. Thus, the third step is to construct the corrected (or continuous)solution polynomial in the standard element. As shown in Fig. 1, the common solutionsat the two end-points of the element of Ω i are ˆ w CIi,L and ˆ w CIi,R respectively. The correctedsolution polynomial within the element is named as ˆ w Ci ( ξ ), and must has features as ˆ w Ci ( −
1) = ˆ w CIi − / = ˆ w CIi,L , ˆ w Ci (1) = ˆ w CIi +1 / = ˆ w CIi,R . (23)The corrected solution polynomial ˆ w Ci ( ξ ) is assumed to have the following form ˆ w Ci ( ξ ) = ˆ w δi ( ξ ) + ( ˆ w Ci ( − − ˆ w δi ( − G L ( ξ ) + ( ˆ w Ci (1) − ˆ w δi (1)) G R ( ξ ) . (24) G L ( ξ ) and G R ( ξ ) are the correct functions related with the left and the right end-points ofthe element, and G L ( ξ ) and G R ( ξ ) should satisfy the following conditions G L ( −
1) = 1 , G L (1) = 0 , G R ( −
1) = 0 , G R (1) = 1 . (25)The correct function is one of the most important factor of the FR Framework. More detailsand introductions of correct function can be referred to the works of H. T. Huynh[43, 44].Now that the corrected solution polynomial is computed, the corrected solution derivativespolynomial ˆ w Ci,ξ ( ξ ) can be obtained directly using the Eq. (24). The correction procedurecan be seen in Fig. 2 briefly.The fourth step is to compute the flux ˆ f i,k at solution points. The flux at the solutionspoint can be evaluated by the macroscopic variables ˆ w i,k and the corrected derivatives9 a) (b) FIG. 2. The demonstration of solution correction: (a) Corrected solution beside the interface and(b) The solution correction within element Ω i . ˆ w Ci,ξ ( ξ k ). The flow is considered as continuous within the element, so the flux solver forsmooth flow is used at the solution points. After the fluxes at solution points have beencomputed, the flux polynomial within the element can be obtained, ˆ f δi ( x ) = k = P X k =0 ˆ f i,k φ k ( ξ k ) . (26) ˆ f δi means that the flux polynomial is always discontinuous at the boundaries of elements.The fifth step focuses on the common flux at the two-end points of element. The fluxeswithin the whole computational domain are assumed as a piecewise function, which has theform as Eq. (26) in each individual element. The fluxes across the boundaries of elementsare always discontinuous. To make the fluxes within element feel the effect of boundaries,it is very important to compute the continuous (or common) fluxes at the interface to getthe accurate results.The common fluxes at boundaries of elements is another critical factor of the FR frame-work. The discontinuity is considered in the computation of common fluxes, and differentRiemann solvers are used according to the numerical method. In our present work, the KIFmethod is applied to solve the fluxes across the interface.The sixth step is to correct the fluxes using the common fluxes at the interfaces. Theprocedure of flux correction is very similar to solution correction, and the Fig. 3 exhibits theprocedure primitively. The corrected flux polynomial reads ˆ f Ci ( ξ ) = ˆ f δi ( ξ ) + ( ˆ f Ci ( − − ˆ f δi ( − H L ( ξ ) + ( ˆ f Ci (1) − ˆ f δi (1)) H R ( ξ ) , (27)where H is the correction function, which is similar to the G used in the correction procedure10 a) (b) FIG. 3. The demonstration of flux correction: (a) Corrected flux beside the interface and (b) Theflux correction within element Ω i . of solution polynomial, and H L ( ξ ) and H R ( ξ ) should satisfy the following conditions H L ( −
1) = 1 , H L (1) = 0 , H R ( −
1) = 0 , H R (1) = 1 . (28)The final step is to compute the divergence of the corrected fluxes at the solution points.Since the corrected flux polynomial is expressed as Eq. (27), the divergence of correctedflux can be obtained directly. By now, all the preconditions of using Eq. (16) to update themacroscopic variables at the solution points are completed for 1D advection problem.The correction function is critical for the flux reconstruction method. G and H have agreat effect on the accuracy and stability. The VCJH scheme developed by Vincent[45] isused in our work. The VCJH scheme has been proved as an energy stable scheme, and itcan be recovered to a particular existing scheme, such as nodal DG and Spectral Difference(SD) methods, with a specific parameter.In terms of time integration, an explicit adaptive Runge-Kutta 45 (RK45) method is usedin present work. D. Extension to multidimensional problem
Extension to quadrilateral and hexahedral elements are straight forward [43, 44]. Fortriangles [47, 62] and tetrahedra [63], they are more complicated in algorithm and imple-mentation. But, the procedure is analogous to FR in one dimension.11 . Shock capturing method
The robust shock capturing is a main difficulty for high-order FE-type CFD method. Inthe vicinity of discontinuities, the smooth indicator [64–66] is used to detect the discontinuity.Once the shock has been sensed, the shock capturing method is applied on the elements. Inthe present work, we follow the idea of Jameson [66]. We use two parameters s and κ todecide whether the shock capturing method should be applied. The determination of values s and κ can be referred to works [67, 68]. In the present paper, we set s + κ around thevalue 0 .
01. We find this setting can keep the scheme robust and accurate.
III. NUMERICAL TEST CASES
In this section, numerical tests are set up for the validation of present method. The accu-racy order, shock capturing method, viscous flow problem and various boundary conditionsare all validated in the section. Finally, the potential of present scheme to simulate theturbulent flow is verified in the Taylor-Green Vortex problem.Our algorithmic code is deployed on the HiFiLES open-source platform [69], thanks fortheir great works. It should also be noticed that the polynomial order p = 3 is used inthis section. Several one-dimensional problems are simulated using multidimensional codein present paper. The upper and bottom bounds of the computational domain are treatedas periodic boundaries in these cases. A. Accuracy tests
In this case, the advection of density perturbation problem [19] is presented to validatethe accuracy of our method on Cartesian grid. The initial condition is given as ρ ( x ) = 1 + 0 . sin ( πx ) , u ( x ) = 1 , v ( x ) = 0 , p ( x ) = 1 , (29)and the analytic solution at the time t can be expressed as ρ ( x, t ) = 1 + 0 . sin ( π ( x − t )) , u ( x, t ) = 1 , v ( x, t ) = 0 , p ( x, t ) = 1 . (30)12he case is a one-dimensional problem, and we simulate it using a two-dimensional solveron the Cartesian grid. The computational domain is { ( x, y ) | x ∈ [0 , , y ∈ [0 , h ] } , (31)where the length scale h equals to 1 /N . N is the number of elements along the x-direction.Table I gives the L normal of density distribution. The numerical results are obtained at t = 2, and the time step is set as ∆ t = 0 . /N . It can be concluded from the Table I thatthe accuracy order is reached quite well. TABLE I. L normal of density for the advection of density problem N L order10 2 . E − −
20 1 . E −
06 4 . . E −
08 4 . . E −
09 4 . . E −
10 3 . . E −
11 3 . . E −
13 3 . The second case is the isentropic vortex problem, which is a two-dimensional problemalways used to validate the accuracy of high order method. The computational domain is a[ − , × [ − ,
5] square. The periodic boundary is applied on the four bounds of the square.The diagonal uniform form flow, ( ρ, u, v, p ) = (1 , , , δu, δv ) = ǫ π e (0 . − r )) ( − ¯ y, ¯ x ) , (32) δT = − ( γ − ǫ γπ e (1 − r ) , δS = 0 . (33)where (¯ x, ¯ y ) = ( x, y ) . r = (¯ x + ¯ y ) . Since the vortex is moved along the diagonal line with the time marching in the case, thenumerical results is obtained at t = 10. The vortex is just back to the origin position at13he moment. The triangular grid, which is similar to the grid shown in Fig. 11, is used inthe simulation. The time step is set as ∆ t = 0 . /N . N is the number of elements on thebound. Table II gives the L normal of density. The designed accuracy order can be clearlyseen in the table. TABLE II. L N L order10 4 . E − −
20 2 . E −
04 4 . . E −
05 4 . . E −
07 3 . . E −
08 3 . B. One dimensional Riemann problem
The first one-dimensional Riemann problem is the Sod problem [70], which is always usedto validate the ability of numerical schemes to capture the discontinuity. The computationaldomain is ( x, y ) ∈ [0 , × [0 , h ], and the Cartesian grid with different length scale h is usedin the approach. The initial condition reads( ρ, u, v, p ) = (1 , , , , < x < . , (0 . , , , . , . ≤ x ≤ . (34)Fig. 4 shows the density, velocity, pressure, and temperature distributions at t = 0 .
2. Thenumerical results have a good accordance with the exact solution, and it can be obviouslyseen that the accuracy is improved with the h criterion grid refinement. For the shockcapturing method used in present work, the value of s + κ has a great effect to the numericalaccuracy. It is known to all that the larger value of s + κ is set, the less accuracy lost isobtained. Fig. 5 plots the simulation with different values of s + κ . It is evident that thehigher accuracy can be obtained with the larger value of s + κ .The second one-dimensional Riemann problem is the Lax problem[71]. Compared withSod problem, Lax problem has a much stronger discontinuity. The computational domainis ( x, y ) ∈ [0 , × [0 , h ], and the Cartesian grid with different length scale h is used in the14 .0 0.2 0.4 0.6 0.8 1.0 x ρ ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (a) x u ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (b) x p ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (c) x T ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (d)
FIG. 4. The h criterion grid refinement tests for Sod problem: (a) density, (b) u-velocity, (c)pressure, and (d) temperature distributions at t = 0 . approach. The value of s + κ is set as 0 .
01 in the computation.The initial condition is expressed as( ρ, u, v, p ) = (0 . , . , . , . , < x < . , (0 . , . , . , . , . ≤ x ≤ . (35)Fig. 6 shows the density, velocity, pressure and temperature distributions at t = 0 . .0 0.2 0.4 0.6 0.8 1.0 x ρ Exacts +κ=0.01s +κ=0 (a) x u Exacts0+κ=0.01s0+κ=0 (b) x p Exacts +κ=0.01s +κ=0 (c) x T Exacts +κ=0.01s +κ=0 (d) FIG. 5. The results of Sod problem with different values of s + κ : (a) density, (b) u-velocity, (c)pressure, and (d) temperature distributions at t = 0 .
2. The length scale h is set as 1 / C. Shu-Osher problem
The problem of Shu-Osher [71] describes the interaction of a sinusoidal density-wave witha Mach 3 normal shock. The purpose of this case is to validate the behavior of our methodon the shock-wave interaction problem. The shock capturing method is also examined inthe case, and s + κ is set as 0 .
01 in the simulation. The computational domain used in thesimulation is taken as [0 , × [0 , h ]. The Cartesian grid is used in the simulations, and thelength scale of the grid is h = 1 / ρ, u, v, p ) = (3 . , . , . , . , ≤ x < , (1 + 0 . sin (5( x − , . , . , . , ≤ x ≤ . (36)16 .0 0.2 0.4 0.6 0.8 1.0 x ρ ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (a) x u ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (b) x p ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (c) x T ExactPresent, h=1/64Present, h=1/32Present, h=1/16 (d)
FIG. 6. The h criterion grid refinement tests for Lax problem: (a) density, (b) u-velocity, (c)pressure, and (d) temperature distributions at t = 0 . Fig. 7 shows the density distribution at t = 1 .
8, and the zoom in view near the highfrequency wave is also exhibited. Because the exact solution of this problem can not becomputed directly, the solution of fourth order WENO method with 10000 grid points inone dimension is taken as the exact result. It can be seen in Fig. 7 that the performanceis improved with increasing mesh resolution. The discontinuity is captured well, and theability of present scheme to capture the frequency wave is also be verified in the case.
D. Shock vortex interaction problem
The shock vortex interaction problem [35] is always used to validate the performanceof high order method. Compared to lower order method, the high order scheme have theadvantages on resolving the vortex and interaction.In the simulation, a stationary normal shock and a small perturbation are initialed in the17 x ρ ExactKIF 4th: h=1/20KIF 4th: h=1/40
FIG. 7. Shu-Osher problem: the distribution of density-wave at t = 1 . flow field. A Mach 1.1 normal shock wave is located at the position x = 0 .
5. The left sidestate (
M a = 1 .
1) of the shock wave is given as follows( ρ, u, v, p ) = (
M a , √ γ, . , . , T = p/ρ, S = ln ( p/ρ γ ) . (37)A small and weak vortex is superposed to the left side of the normal shock. The center ofthe vortex is ( x c , y c ) = (0 . , . δu, δv ) = κηe µ (1 − η ) ( sinθ, − cosθ ) , (38) δT = − ( γ − κ µγ e µ (1 − η ) ( sinθ, − cosθ ) , δS = 0 , (39)18here κ = 0 . ,µ = 0 . ,η = r/r c ,r c = 0 . ,r = p ( x − x c ) + ( y − y c ) . (40)The computational domain and boundary conditions are exhibited in the Fig. 8. Inthe simulation, the Cartesian grid is used and the grid size h is 1 / s + κ is set as 0 .
01. The initial status at t = 0 and the contourplots at t = 0 . t = 0 . t = 0 . t = 0, and then the vortex moves from left to right across the shock. The profileof vortex varies with the movement, and the interaction of shock and vortex can be seenobviously in the Fig. 9 and Fig. 10. The plots show that our present method can capturethe shock and vortex interaction with enough resolution, and the vortex is recovered well.It also shows clearly in Fig. 10 that the shock bifurcations reaches to the top boundary, andthe reflection is evident. left side right sideshock out (cid:1) owin (cid:0) ow symmetry FIG. 8. Shock vortex interaction problem: the computational domain and the boundary conditions.In the simulation, Cartesian grid is used, and the element length scale of mesh h equals to 1 / E. Lid-driven cavity flow
The lid-driven cavity flow [72] is one of the benchmarks for validating the performanceof the viscous flow solver, and the aim of this case is also to examine performance of presentmethod on viscous solid wall. An incompressible flow is initialed in the computational19 a) (b)
FIG. 9. Shock vortex interaction problem: (a) the initial status and (b) the contour of pressure at t = 0 . (a) (b) FIG. 10. Shock vortex interaction problem: (a) the contours of pressure at t = 0 . t = 0 . domain, and the Mach number of the lid is set as M a = 0 .
1. The Reynolds number are Re = 400 , , , × [0 , h = 1 /
16. The u-velocityprofiles along the vertical center-line and the v-velocity profiles along the horizontal center-line are all compared with the reference data in Figs. 12-14. The figures exhibit that thepresent results match quite well with the data from Ghia [72]. The numerical data extractedfrom Ref. [22] at Re = 1000 , ×
65 Cartesian grid are also shown in the figures.It is obvious that the present method can reach the same accuracy with fewer mesh nodes.
F. Blasius incompressible laminar flat plate
The incompressible boundary layer flow over a flat plate is simulated. The Mach numberis
M a ∞ = 0 .
15, and the Reynolds number based on the length of plate is Re ∞ = 1 × . The20 IG. 11. The triangular grid used in the simulation of lid-driven cavity flow (816 triangularelements). u / U ∞ GhiaKIF, Cartesian, h=1/32KIF, Cartesian, h=1/16KIF, Triangular, h=1/16 (a) v / U ∞ GhiaKIF, Cartesian, h=1/32KIF, Cartesian, h=1/16KIF, Triangular, h=1/16 (b)
FIG. 12. The velocity profile of the lid driven cavity flow at Re = 400: (a) u-velocity profiles at x = 0 . y = 0 . subscript ∞ indicates the state of free stream flow. The length of the flat plate is L = 100,and the leading edge of the flat plate is located at x = 0. The computational domain is[ − , × [0 , .0 0.2 0.4 0.6 0.8 1.0y−0.4−0.20.00.20.40.60.81.0 u / U ∞ GhiaPan, Cartesian, h=1/65KIF, Cartesian, h=1/32KIF, Cartesian, h=1/16KIF, Triangular, h=1/16 (a) v / U ∞ GhiaPan, Cartesian, h=1/65KIF, Cartesian, h=1/32KIF, Cartesian, h=1/16KIF, Triangular, h=1/16 (b)
FIG. 13. The velocity profile of the lid driven cavity flow at Re = 1000: (a) u-velocity profiles at x = 0 . y = 0 . u / U ∞ GhiaPan, Cartesian, h=1/65KIF, Carteisian, h=1/32KIF, Triangular, h=1/16 (a) v / U ∞ GhiaPan, Cartesian, h=1/65KIF, Carteisian, h=1/32KIF, Triangular, h=1/16 (b)
FIG. 14. The velocity profile of the lid driven cavity flow at Re = 3200: (a) u-velocity profiles at x = 0 . y = 0 . coarser grid is shown in Fig. 15, and the details of the two grids are shown in Table III.“Grid 2” has fewer elements than “Grid 1” in the boundary layer.The velocity profiles versus η = y q U ∞ νx are shown in Figs. 16-18. It can be seen clearlyin the pictures that the velocity profiles have a good accordance with the Blasius solutionat every positions. It is obvious that the present method can approach boundary layer withvery few elements. The skin friction coefficient is shown in Fig. 19. The numerical resultsshows an excellent performance compared with Blasius solution.22 a) (b) FIG. 15. The “Grid 2” used in the simulation of laminar flow over a flat plate.TABLE III. The details of the grids used in the laminar boundary layer case.Grid size Number of elementsin boundary layerGrid 1 4876 25Grid 2 2336 6 u/U ∞ η BlasiusGrid 1Grid 2 (a) v √ ∞ η BlasiusGrid 1Grid 2 (b)
FIG. 16. Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocity profiles at x/L = 0 . G. Taylor-Green vortex at Re = 1 , The Taylor-Green Vortex (TGV) is a simple test case for the resolution of the small scalesof a turbulent flow by a numerical method. The compressible TGV at Re = 1600 was one ofthe benchmark problems in the 1st and 2nd International Workshops on High-Order CFDMethods. The reference solution used in current paper was obtained by Debonis [73] using23 .0 0.2 0.4 0.6 0.8 1.0 u/U ∞ η BlasiusGrid 1Grid 2 (a) u profile v √ ∞ η BlasiusGrid 1Grid 2 (b) v profile FIG. 17. Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocity profiles at x/L = 0 . u/U ∞ η BlasiusGrid 1Grid 2 (a) u profile v √ ∞ η BlasiusGrid 1Grid 2 (b) v profile FIG. 18. Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocity profiles at x/L = 0 . x/L C f BlasiusGrid 1Grid 2
FIG. 19. Blasius incompressible laminar flat plate: the skin friction coefficient distribution alongthe flat plate.
24 high-order dispersion-relation-preserving (DRP) scheme on a mesh of 512 elements. Thecomputational domain is a cubic box of dimensions [0 , π ] , and the periodic boundary isapplied on the faces of the cube. In the case, both of the 32 and 64 grids are used.The initial condition is set as u ( t ) = u sin ( x/L ) cos ( y/L ) cos ( z/L ) ,v ( t ) = − u cos ( x/L ) sin ( y/L ) cos ( z/L ) ,w ( t ) = 0 ,p ( t ) = p + ρ V
16 [ cos (2 x/L ) + cos (2 y/L )] [ cos (2 z/L ) + 2] , (41)where ρ = 1, p = 100, u = 1, and L = 1. The Mach number is set to 0 .
08, and the initialtemperature is 300 K . The volume-averaged kinetic energy and the dissipation rate of thekinetic energy are computed. The volume-averaged kinetic energy is read as E k = < k > = 1 ρ Ω Z Ω ρ u i u i d Ω , (42)and the dissipation rate of the kinetic energy is give by ǫ ( E k ) = − dE k dt . (43)The numerical results of averaged kinetic energy and the dissipation rate of kinetic energyare compared with reference data from Debonis [73]. The results have a good accordancewith the reference data. The iso-surfaces of Q criterions colored by velocity magnitude attime 3, 5, 7, and 9 are shown in Fig. 21. The evolution of flow structure from large scalevortices to small vortices can be clearly seen in the figure. IV. CONCLUSION
In present paper, a high order numerical scheme is proposed based on the flux reconstruc-tion framework and Kinetic Inviscid Flux. KIF, which aims to find a balance between theadvantages of gas-kinetic scheme and lower computational costs, is a combination of TTTscheme and KFVS method. FR framework is well understood and available. It has theproperties of robustness and compactness, which are of great importance for the high ordermethod. The accuracy order have been verified using the advection of density perturba-tion problem and isentropic vortex problem. The results show that the accuracy of present25 .0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 t(s) < k > / < k > Debonis, DNS, 512KIF, DNS, 32KIF, DNS, 64 (a) t(s) ε Debonis, DNS, 512KIF, DNS, 32KIF, DNS, 64 (b)
FIG. 20. The kinetic energy and the dissipation rate: (a) E k and (b) ǫ ( E k ). The integer numberdenotes the grid size. For example, the number 32 indicates that the grid size is 32 . method reaches to the designed order. KIF also can be viewed as an inviscid-viscous splittingversion of the gas-kinetic scheme, and the excellent performance of proposed method can beseen in the simulations of lid-driven cavity flow and Blasius incompressible laminar bound-ary layer. The good numerical results have shown the success of inviscid-viscous strategy ofgas-kinetic scheme. The Taylor-Green Vortex problem has been used to verify the potentialof present method to simulate turbulent flow, and the excellent results are obtained. Theshock capturing method used in current paper is still an open question which is needed infurther study. The ability to capture strong discontinuity should be improved further. ACKNOWLEDGMENTS
The present work is supported by National Natural Science Foundation of China (GrantNo. 11702223, No. 11902266 and No. 11902264), National Numerical Wind Tunnel Projectof China, and 111 Project of China (Grant No. B17037), as well as the ATCFD Project(2015-F-016). 26 a) (b)(c) (d)
FIG. 21. Iso-surfaces of Q ( Q = 0 .
5) criterion colored by velocity magnitude: (a) t = 3, and (b) t = 5, (c) t = 7, and (d) t = 9. EFERENCES [1] K. Xu, Journal of Computational Physics , 289 (2001).[2] K. Xu, M. Mao, and L. Tang, Journal of Computational Physics , 405 (2005).[3] P. L. Bhatnagar, E. P. Gross, and M. Krook, Physical Review , 511 (1954).[4] Q. Li, S. Tan, S. Fu, and K. Xu, in (2010) pp. 12–17.[5] S. Xiong, C. Zhong, C. Zhuo, K. Li, X. Chen, and J. Cao, International Journal for NumericalMethods in Fluids , 1833 (2011).[6] D. Pan, C. Zhong, J. Li, and C. Zhuo, International Journal for Numerical Methods in Fluids , 748 (2016).[7] J. Li, C. Zhong, D. Pan, and C. Zhuo, Computers & Mathematics with Applications , 1227 (2019).[8] G. Cao, H. Su, J. Xu, and K. Xu, Aerospace Science and Technology , 958 (2019).[9] Q. Li, S. Fu, and K. Xu, AIAA Journal , 2170 (2005).[10] H. Liu, K. Xu, T. Zhu, and W. Ye, Computers & Fluids , 115 (2012).[11] T. Zhu and W. Ye, Numerical Heat Transfer, Part B: Fundamentals , 203 (2010).[12] R. Yuan, C. Zhong, and H. Zhang, Journal of Computational Physics , 184 (2015).[13] R. Yuan and C. Zhong, Applied Mathematical Modelling , 417 (2018).[14] H. Dong and L. M. Yang, Advances in Applied Mathematics and Mechanics , 1177 (2019).[15] W. Li, M. Kaneda, and K. Suga, Computers & Fluids , 100 (2014).[16] J. Li, C. Zhong, Y. Wang, and C. Zhuo, Physical Review E , 053307 (2017).[17] K. Fujii, Progress in Aerospace Sciences , 455 (2005).[18] M. R. Visbal and D. V. Gaitonde, Journal of Computational Physics , 155 (2002).[19] Q. Li, K. Xu, and S. Fu, Journal of Computational Physics , 6715 (2010).[20] J. Luo and K. Xu, Science China Technological Sciences , 2370 (2013).[21] G. Zhou, K. Xu, and F. Liu, Journal of Computational Physics , 146 (2017).[22] L. Pan, K. Xu, Q. Li, and J. Li, Journal of Computational Physics , 197 (2016).[23] L. Pan and K. Xu, Computers & Fluids , 250 (2015).[24] L. Pan and K. Xu, Communications in Computational Physics , 985 (2015).
25] X. Ji, L. Pan, W. Shyy, and K. Xu, Journal of Computational Physics , 446 (2018).[26] X. Ji, F. Zhao, W. Shyy, and K. Xu, Journal of Computational Physics , 109367 (2020).[27] F. Zhao, X. Ji, W. Shyy, and K. Xu, Advances in Aerodynamics , 13 (2019).[28] S. Liu, J.-Z. Cao, and C. Zhong, arXiv: Computational Physics (2020).[29] T. J. Barth and P. O. Frederickson, AIAA paper , 0013 (1990).[30] R. Abgrall, Journal of Computational Physics , 45 (1994).[31] L. J. Durlofsky, B. Engquist, and S. Osher, Journal of Computational Physics , 64 (1992).[32] C. F. Ollivier-Gooch, Journal of Computational Physics , 6 (1997).[33] T. Sonar, Computer Methods in Applied Mechanics and Engineering , 157 (1997).[34] X.-D. Liu, S. Osher, and T. Chan, Journal of Computational Physics , 200 (1994).[35] C.-W. Shu, in Advanced Numerical Approximation of Nonlinear Hyperbolic Equations (Springer, 1998) pp. 325–432.[36] C. Hu and C.-W. Shu, Journal of Computational Physics , 97 (1999).[37] B. Cockburn and C.-W. Shu, Journal of Scientific Computing , 173 (2001).[38] Y. Liu, W. Zhang, Y. Jiang, and Z. Ye, Computers & Mathematics with Applications , 1096 (2016).[39] Q. Wang, Y.-X. Ren, and W. Li, Journal of Computational Physics , 863 (2016).[40] Q. Wang, Y.-X. Ren, and W. Li, Journal of Computational Physics , 883 (2016).[41] Q. Wang, Y.-X. Ren, J. Pan, and W. Li, Journal of Computational Physics , 1 (2017).[42] Z. Wang, Progress in Aerospace Sciences , 1 (2007).[43] H. T. Huynh, AIAA paper , 4079 (2007).[44] H. T. Huynh, AIAA paper , 403 (2009).[45] P. E. Vincent, P. Castonguay, and A. Jameson, Journal of Scientific Computing , 50 (2011).[46] A. Jameson, Journal of Scientific Computing , 348 (2010).[47] P. Castonguay, P. E. Vincent, and A. Jameson, Journal of Scientific Computing , 224(2012).[48] D. A. Kopriva and J. H. Kolias, Journal of Computational Physics , 244 (1996).[49] D. A. Kopriva, Journal of Computational Physics , 475 (1996).[50] K. Xu, Numerical Hydrodynamics from Gas-Kinetic Theory , Ph.D. thesis, COLUMBIA UNI-VERSITY. (1993).[51] D. I. Pullin, Journal of Computational Physics , 231 (1980).[52] J. C. Mandal and S. M. Deshpande, Computers & Fluids , 447 (1994).
53] T. Tao and K. Xu, Zeitschrift fr angewandte Mathematik und Physik , 258 (1999).[54] Y. Sun, C. Shu, L. M. Yang, and C. J. Teo, Advances in Applied Mathematics and Mechanics , 703 (2016).[55] T. Ohwada, Y. Shibata, T. Kato, and T. Nakamura, Journal of Computational Physics ,131 (2018).[56] K. Xu, Direct modeling for computational fluid dynamics: construction and application ofunified gas-kinetic schemes (World Scientific, 2015).[57] Z. Guo, K. Xu, and R. Wang, Physical Review E , 033305 (2013).[58] K. Qu, C. Shu, and Y. T. Chew, Physical Review E , 036706 (2007).[59] K. Xu and J.-C. Huang, Journal of Computational Physics , 7747 (2010).[60] L. Wu, J. Zhang, J. M. Reese, and Y. Zhang, Journal of Computational Physics , 602(2015).[61] C. Liu, Y. Zhu, and K. Xu, Journal of Computational Physics , 108977 (2020).[62] D. Williams, P. Castonguay, P. Vincent, and A. Jameson,Journal of Computational Physics , 53 (2013).[63] D. M. Williams and A. Jameson, Journal of Scientific Computing , 721 (2014).[64] P.-O. Persson and J. Peraire, AIAA paper , 112 (2006).[65] P.-O. Persson, AIAA paper , 3061 (2013).[66] A. Sheshadri and A. Jameson, AIAA paper , 2688 (2014).[67] M. L. Yu, F. X. Giraldo, M. Peng, and Z. J. Wang, Monthly Weather Review , 4823(2015).[68] R. Vandenhoeck and A. Lani, Computer Physics Communications , 1 (2019).[69] M. R. Lpez, A. Sheshadri, J. R. Bull, T. D. Economon, J. Romero, J. E. Watkins, D. M.Williams, F. Palacios, A. Jameson, and D. E. Manosalvas, AIAA paper , 3168 (2014).[70] G. A. Sod, Journal of Computational Physics , 1 (1978).[71] C.-W. Shu and S. Osher, Journal of Computational Physics , 439 (1988).[72] U. Ghia, K. Ghia, and C. Shin, Journal of Computational Physics , 387 (1982).[73] J. R. DeBonis, AIAA paper , 0382 (2013). IST OF FIGURES i − , Ω i and Ω i +1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 The demonstration of solution correction: (a) Corrected solution beside theinterface and (b) The solution correction within element Ω i . . . . . . . . . . . . . . . . 103 The demonstration of flux correction: (a) Corrected flux beside the interfaceand (b) The flux correction within element Ω i . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 The h criterion grid refinement tests for Sod problem: (a) density, (b) u-velocity, (c) pressure, and (d) temperature distributions at t = 0 .
2. . . . . . . . . . 155 The results of Sod problem with different values of s + κ : (a) density, (b)u-velocity, (c) pressure, and (d) temperature distributions at t = 0 .
2. Thelength scale h is set as 1 /
64. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 The h criterion grid refinement tests for Lax problem: (a) density, (b) u-velocity, (c) pressure, and (d) temperature distributions at t = 0 .
14. . . . . . . . . 177 Shu-Osher problem: the distribution of density-wave at t = 1 .
8. . . . . . . . . . . . . 188 Shock vortex interaction problem: the computational domain and the bound-ary conditions. In the simulation, Cartesian grid is used, and the elementlength scale of mesh h equals to 1 / t = 0 .
3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2010 Shock vortex interaction problem: (a) the contours of pressure at t = 0 . t = 0 .
8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2011 The triangular grid used in the simulation of lid-driven cavity flow (816 tri-angular elements). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2112 The velocity profile of the lid driven cavity flow at Re = 400: (a) u-velocityprofiles at x = 0 . y = 0 .
5. . . . . . . . . . . . . . . . . . . . 2113 The velocity profile of the lid driven cavity flow at Re = 1000: (a) u-velocityprofiles at x = 0 . y = 0 .
5. . . . . . . . . . . . . . . . . . . . 2214 The velocity profile of the lid driven cavity flow at Re = 3200: (a) u-velocityprofiles at x = 0 . y = 0 .
5. . . . . . . . . . . . . . . . . . . . 2215 The “Grid 2” used in the simulation of laminar flow over a flat plate. . . . . . . . 23316 Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocityprofiles at x/L = 0 .
1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2317 Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocityprofiles at x/L = 0 .
3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2418 Blasius incompressible laminar flat plate: (a) u-velocity and (b) v-velocityprofiles at x/L = 0 .
5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2419 Blasius incompressible laminar flat plate: the skin friction coefficient distri-bution along the flat plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2420 The kinetic energy and the dissipation rate: (a) E k and (b) ǫ ( E k ). The integernumber denotes the grid size. For example, the number 32 indicates that thegrid size is 32 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2621 Iso-surfaces of Q ( Q = 0 .