I Introduction
Increasing engineering demands for accuracy and efficiency have catalyzed the development of highorder methods in the field of computational fluid dynamics. Complex geometries are often involved in engineering problems, necessitating the use of unstructured meshes, which are more flexible and more easily generated than structured meshes. Several flow problems require low dissipation and dispersion errors, such as large eddy simulation and direct numerical simulation of turbulence and problems in aeroacoustics. Loworder methods are usually too dissipative to resolve smallscale structures. In contrast, highorder methods can potentially achieve increased accuracy with low computational cost. However, existing highorder methods are generally less robust and more complicated than loworder methods. Therefore, it is necessary to develop new highorder methods for unstructured meshes to improve robustness, simplicity, and efficiency.
This study was performed in the finite volume (FV) framework. Compared to finite difference (FD) methods, finite volume methods, and finite element (FE) methods are more suitable for unstructured meshes. FV methods are usually simpler and more robust than FE methods, especially for highspeed flows containing discontinuities (e.g., shocks). Many highorder FV methods were developed on unstructured meshes, such as the kexact method T.J Barth and P.O. Frederickson , the essentially nonoscillatory (ENO) method Abgrall (1992), and the weighted ENO (WENO) method Friedrich (1998)
. Because only cell averages are available in common FV methods, solution polynomials with higher degrees of freedom (DoFs) are required to represent the flow structures, which may reduce the compactness, efficiency, and robustness. For example, the noncompactness resulting from a large stencil for reconstruction leads to problems with numerical dissipation, highorder numerical approximations of boundary conditions, reduced parallel efficiency, and caches missing.
To achieve compactness for highorder FV methods, the compact leastsquare reconstruction and variational reconstruction Wang et al. (2016, 2017) were developed, where a large linear equation system should be solved and implicit timestepping methods are necessary. By subdividing each cell into a set of subcells, the spectral volume (SV) achieves compactness as wellZ.J. Wang et al. (2004), and subcellaveraged solutions are updated under the FV framework. Compared to traditional FV methods, the SV method can achieve higher resolution with the same mesh size owing to the use of subcells, which helps to capture flow structures with small scales more accurately. However, the SV method requires the number of subcells to be equal to that of the unknowns in the solution polynomial. The complex subdivision of the main cells thus becomes a bottleneck, and no stable subdivision has yet been found for quadrilateral or hexahedral meshes. As an extension of the SV method, the subcell finite volume (SCFV) method Pan and Ren (2017); Pan et al. (2017) overcomes the difficulty of subdivision by combining the idea of the method, avoiding the restriction of the number of subcells in SV. The subcells from faceneighboring cells are involved in the stencil for reconstruction. As a result, the subdivision is much easier and more flexible. Theoretically, the SCFV method can achieve arbitrarily highorder accuracy with a compact stencil. In summary, the SCFV method has the following three advantages over existing highorder methods:
(i) Compact reconstruction is more easily implemented compared to traditional FV methods. Generally, the arbitrary distribution and number of unstructured meshes make it difficult to choose and search neighboring cells to form the stencil. The SCFV method combines the features of both unstructured and structured meshes. First, a set of coarser unstructured meshes are generated, and are called main cells. Each main cell is further subdivided into a set of subcells that are stored in a structured way. In fact, these subcells are the real control volumes (CVs), and the subcell averages are updated in the FV framework. In this way, the complexity of choosing and searching neighbors is reduced greatly. As only faceneighboring main cells are considered, there are usually enough DOFs, i.e., subcells, for a compact reconstruction.
(ii) The reconstruction is efficient compared to traditional FV methods. Assume that each main cell is subdivided into subcells. Similar to finite element (FE) methods, the subcells (CVs) are treated as inner DOFs, and a common solution polynomial is reconstructed for them. With the same number of CVs, reconstruction is conducted only times compared to traditional FV methods for which the reconstruction needs to be conducted for each CV separately. Moreover, the memory space for the corresponding coefficient matrix is also of that in traditional FV. Therefore, the reconstruction in the SCFV method is more efficient in terms of both computational cost and memory storage.
(iii) The subcell resolution can be well preserved compared to FE methods. There are many FE or hybrid methods that also achieve compactness easily, such as the DG method Cockburn and Shu (1998), the CPR method Huynh ; Z.J. Wang and Gao (2009), and the method Dumbser et al. (2008). However, the solution distribution needs to remain continuous inside each main cell, which results in the difficulty with respect to capturing shocks. It is difficult to fully maintain the advantage of inner DOFs, i.e., the subcell resolution. In contrast, for the SCFV method, discontinuities can even exist inside each main cell. The resolution and robustness for shock capturing can thus be enhanced effectively. Besides, the SCFV method can take larger CFL numbers, and the formulation is much simpler compared to FE methods.
The SCFV method was originally developed on quadrilateral meshes to solve the Euler equations Pan and Ren (2017); Pan et al. (2017). Then, it was further developed on triangular meshes and extended to the NS equations by combining the singlestage thirdorder gaskinetic flux solver (GKS) Zhang and Q.B. Li (2021). Although other solvers, such as the GRP solver, BenArtzi et al. (2006); BenArtzi and Li (2007) can also be adopted, this study still uses the GKS solver but with secondorder accuracy Xu (2001); Q.B. Li and Fu (2006) for the flux construction because the inviscid and viscous flux are coupled and computed simultaneously L.Q. Zhang et al. (2019); Zhan et al. (2021); L.M. Yang et al. (2021). Besides, a series of unified gaskinetic scheme (UGKS) has been developed for the entire flow regime from the continuum flow to the highly rarefied flowXu and J.C. Huang (2010); Zhang et al. (2019); Chen et al. (2020); Zhong et al. (2021). During the past decade, highorder GKS has been developed systematically. By taking a secondorder Taylor expansion, a thirdorder multidimensional GKS was developed within a single stage Q.B. Li et al. (2010). The timedependent gas distribution function can even be used to evolve the solutions at cell interfaces, which provides additional DoFs to help construct a compact thirdorder FVGKS Pan and Xu (2016); Ji et al. (2020). Other extensions of singlestage thirdorder GKS can also be found for the DG method Ren et al. (2015) and the CPR method Zhang et al. (2018). As mentioned above, a singlestage thirdorder SCFVGKS has been proposed recently Zhang and Q.B. Li (2021). However, when fourthorder accuracy is considered, the singlestage flux function is quite complicated, and the computational cost increases significantly Liu and Tang (2014); the robustness may also be weakened. Fortunately, the twostage fourthorder timestepping method Li and Du (2016) provides an efficient way for LaxWendrofftype solvers to achieve fourthorder time accuracy. Only the secondorder gaskinetic flux solver is needed, which is much simpler and more robust Zhao et al. (2019); Yang et al. (2021); F. Zhao et al. (2021). The computational cost of the flux evaluation can be reduced effectively compared to the singlestage highorder GKSPan et al. (2016).
Hence, this study developed a twostage temporalspatial fourthorder SCFVGKS to solve compressible flows over triangular meshes, with the aim of enhancing the compactness, efficiency, and robustness. Compared to the singlestage SCFVGKS, both accuracy and efficiency are improved significantly. By adopting the twostage temporal discretization, the flux evaluation under the SCFV framework can be simplified greatly, making it easier to implement. The robustness is also enhanced because of the use of the secondorder gaskinetic flux. These improvements make the SCFV method more promising for engineering applications.
This paper is organized as follows. Section 2 presents the compact fourthorder reconstruction based on the SCFV method and the twostage gaskinetic flux evolution. Section 3 presents numerical tests with several benchmark cases to demonstrate the high accuracy, efficiency, and robustness of the current method for simulating compressible flows. The last section presents the conclusions.
Ii Gas kinetic solverbased subcell finite volume method
In this section, we propose a GKSbased SCFV method for compressible fluid flows. The primary steps are the compact subcell technique and the twostage fourthorder temporal advancing based on the secondorder GKS.
ii.1 Subcell finite volume method
First, we briefly review the SCFV method. Consider the twodimensional (2D) conservation law
(1) 
where are the conservative variables, where are the macroscopic velocities and
is the flux vector. The computational domain is divided into
nonoverlapping triangular cells , which are also referred to as main cells. As shown in Figure 1, each main cell is uniformly subdivided into four similar subcells . In the SCFV method, subcellaveraged solutions are stored and updated under the FV framework. More clearly, a subcell can also be called a CV when compared to that for traditional FV methods. The primary task is to reconstruct solution polynomials over each main cell. Then, the numerical flux at subcell interfaces can be determined and used to update subcell averages. The reconstruction is implemented for each component of conserved variables . To achieve fourthorder space accuracy, over each main cell is approximated by a solution polynomial of degree , where unknown coefficients need to be determined. As shown in Figure 1, by involving faceneighboring cells, there are 16 subcells available to sufficiently determine the unknowns.In comparison, Figure 2 also shows the stencil used in a GKSbased FV method and the classical kexact FV Zhao et al. . Note that for the GKSbased FV, cellaveraged slopes are updated and participate in the reconstruction, which helps to avoid the large stencil for the kexact FV. Nonetheless, with only face neighbors, it is still not enough for a fourthorder reconstruction, which results in the enlargement of the stencil. In contrast, the stencil of SCFV is most compact.
For the convenience of further description, all subcells involved in the stencil are denoted as . Besides, is specifically used to indicate subcells from the target main cell . The solution polynomial on is expressed as
(2) 
where , indicates the centroid of . The zeromean basis functions are adopted so that the averaged solution can be conserved automatically on . Thus, there are unknowns to be determined for each . The definition of is
(3) 
where and satisfy . is the volume of . Integrating over gives
(4) 
where ,
(5) 
and is the averaged solution on the subcell . Note that because has been used in Eq.(2) and it has the relation with as
(6) 
the target main cell can only provide three more DOFs for reconstruction besides . Here, the central subcell is excluded from the stencil. Therefore, 15 equations are left in Eq.(4), which forms an overdetermined system. The unknowns are solved by the weighted leastsquare technique
(7) 
where the weight , indicates the distance between the centroid of and . The resulting solution polynomial is shared by the four subcells . However, the averaged solutions of the four subcells are generally not conserved directly by the above weighted leastsquare reconstruction, i.e.,
(8) 
Studies showed that Pan and Ren (2017); Pan et al. (2017), it is necessary to conserve the subcell averages on the four subcells, otherwise, the scheme may not be stable. A simple correction is adopted by shifting on each subcell , i.e., where
(9) 
Then, is used for the flux evaluation. The correction introduces jumps artificially across interfaces between subcells, as shown in Figure 3. The computational cost for the flux evaluation increases at interfaces between subcells.
Some remarks are made below. Because the secondorder gas distribution function is much simpler than the thirdorder one, and the computational cost for numerical flux construction has been reduced significantly, we simply take the weighted leastsquare (WLS) technique to determine the unknowns, rather than the constrained leastsquare (CLS) technique adopted in the singlestage SCFVGKS Zhang and Q.B. Li (2021). Both the techniques have advantages and disadvantages,. With the CLS technique, the subcell averages can be conserved directly and the distribution inside each main cell can remain continuous so that the computational cost of flux is reduced. Thus, it introduces more benefits for the thirdorder gaskinetic flux. However, the technique itself is more complicated and timeconsuming, and it may reduce the robustness of SCFV. In contrast, the WLS technique is simpler and more efficient at the expense of increasing the computational cost of flux because of the correction. However, the secondorder gaskinetic flux is cheaper, and by introducing jumps, the robustness of SCFV is enhanced, and the resolution for discontinuities can be improved, especially for shock capturing.
For smooth flows, the solution polynomial obtained by the above reconstruction can be used directly for flux evaluation. However, for shock capturing, it may lead to numerical oscillations near shock waves. It is necessary to apply an effective limiting procedure or limiter to the solution polynomial. To reduce the computational cost, a shock detector Li and Ren (2012) is adopted to mark troubled cells near shock waves. For other unmarked cells, no limiter is needed and the high accuracy in smooth flow regions can be maintained. For troubled cells, a limiting procedure based on hierarchical reconstruction (or HR limiter) is implemented to obtain a new solution polynomial, which is able to suppress oscillations effectivelyXu et al. (2009). The limiter can keep the designed order of accuracy, and only faceneighboring cells are needed, which remains the compactness of the current scheme. For the details of the limiting procedure, we refer to Ref.Zhang and Q.B. Li, 2021.
ii.2 Twostage gaskinetic flux evolution
By integrating Eq.(1) over each subcell , the semidiscrete form of the SCFV framework can be obtained as follows.
(10) 
where is the area of , is the interfaces surrounding , and is the outward unit normal vector. To maintain fourthorder space accuracy, the line integral on the righthand side is discretized using the Gaussian quadrature rule with two points
(11) 
where is the length of subcell interfaces, and are the quadrature weights of the two Gaussian points. These Gaussian points are also referred to as flux points for clarity. For the temporal discretization of Eq.(10), by fully exploiting the timeevolving flux function provided by the secondorder gaskinetic flux solver, fourthorder time accuracy can be achieved more efficiently with the help of the twostage temporal discretization Li and Du (2016), which can be expressed with the current notations as
(12) 
where is the flux at the intermediate stage in the time interval . The aforementioned twostage temporal discretization has been proved to achieve fourthorder time accuracy for the hyperbolic conservation law Li and Du (2016). Assuming that the computational mesh does not change with time, is a linear function of , giving . Thus, the problem remaining is to determine the flux and its time derivative at each flux point.
Based on the gaskinetic theory, both the macroscopic conservative variables and the numerical flux can be obtained from the gas distribution function with the following relations
(13) 
where is a function of physical space , time , particle velocity , and the internal DoFs ,
is the vector of moments, and
is the element of the phase space. For convenience, the summation convention is adopted in this subsection, such as . The governing equation of the gas distribution function is the 2D BGK equation P.L. Bhatnagar et al. (1954)(14) 
where is the collision time related to the viscosity and pressure . The local equilibrium state corresponds to the macroscopic variables, and can be expressed as
(15) 
where is the total number of . is the gas constant. The conservation law Eq.(1) can be recovered by taking moments of Eq.(14), where the collision term vanishes automatically owing to the conservation of mass, moments, and total energy during collisions, i.e., the compatibility condition. In particular, the NaiverStokes equations can be recovered through the firstorder ChapmanEnskog expansion Xu (2002, 2015)
(16) 
By eliminating the firstorder term in Eq.(16), the Euler equations can be recovered as well. To compute the numerical flux through Eq.(13), a timedependent gas distribution function is constructed at each flux point based on the analytical solution of Eq.(14)
(17) 
where is the initial distribution function at the beginning of each time step (), and is the local equilibrium state. For convenience, the subcell interface is assumed to be perpendicular to the axis, and the flux point is assumed to be the origin in a local coordinate system. Using a firstorder Taylor expansion of and around the flux point and combined with Eq.(16), the timedependent gas distribution function can be obtained as follows.
(18) 
where the coefficients are related to the Taylor expansion of the corresponding Maxwellian functions, i.e., the derivatives of , and the coefficients and correspond to and , respectively. is the Heaviside function. These coefficients are determined by the reconstructed conservative variables and the slopes, as well as the compatibility condition. For details about the evaluation of these coefficients, we refer to Ref.Zhang et al., 2018. The gaskinetic flux solver is intrinsically multidimensional by involving both normal and tangential spatial derivatives in the construction of the gas distribution function. Besides, compared to the thirdorder gaskinetic flux solver, the construction of the above gas distribution function Eq.(18) is much simpler, and the computational cost of flux can be reduced significantly.
Now, the flux and its time derivative are determined as below. With Eq.(18), a timedependent flux function can be constructed according to Eq.(13), denoted as , which is a nonlinear function of . A simple fitting method is adopted to obtain the approximated flux, and its firstorder time derivative Pan et al. (2016), where is approximated by a linear function within the time interval . Denote the time integration of within as
(19) 
Then, an equation set can be obtained
(20) 
where the lefthand side is the time integration of within and , while the righthand side is obtained according to Eq.(19). By solving the equation set Eq.(20), we have
(21) 
Similarly, the approximated flux and its time derivative can be obtained by simply replacing the superscript in Eq.(21) with . Finally, the twostage time stepping method for updating subcell averages can be summarized as
(22) 
where the flux at each flux point is respectively computed according to
(23) 
Although there are two stages in the current scheme, the flux evaluation can still be more efficient than the previous singlestage SCFVGKS because the secondorder gaskinetic flux is much cheaper. Besides, the robustnes of SCFVGKS can also be enhanced.
Iii Numerical tests
In this section, several benchmark cases are tested to validate the performance of the current scheme. For inviscid flows, the collision time is computed by
(24) 
where is set as for included in the exponential term in Eq.(18) to provide the required numerical dissipation; otherwise, is set as to better approximate the inviscid assumption. The second term in Eq.(24) represents the artificial viscosity, where and are the pressure at the left and right sides of the subcell interface. The coefficient is set to 10.0. For viscous flows, the first term should be replaced by the molecular viscosity , where is the dynamic viscous coefficient, and is the pressure at the subcell interface.
The CFL number is set to for all test cases when the time step is computed with reference to the mesh size of main cells, denoted as . Note that it is more reasonable to consider the size of subcells (or CVs) when compared to traditional FV methods. Thus, the corresponding CFL number used in the current scheme is actually with reference to the size of subcells, i.e., . According to our tests, the upper limit of the CFL number is empirically 0.35 with the size of main cells (or 0.7 with the size of subcells) for a stable time marching, which is much higher than the DG and CPR methods, and which is comparable with traditional highorder FV methods. The ratio of specific heat is set to . In the accuracy test, the error and error are computed by
(25) 
where and indicate the averaged numerical solution and the analytical solution on subcell , respectively, and is the number of main cells. The twodimensional (2D) contours of flow fields are plotted based on subcell averaged solutions.
iii.1 Isentropic vortex propagation
The isentropic vortex propagation Shu (1998) was simulated to verify the accuracy and efficiency of the current scheme. An isentropic vortex was added to the mean flow by introducing perturbations in the velocity and temperature , but without perturbations in the entropy , i.e.,
(26) 
where the vortex strength , . The exact solution is that the vortex propagates with constant velocity . The computational domain is , and the periodic boundary condition is adopted for all boundaries. A sample mesh is shown in Figure 4. The mesh is refined by splitting each cell into four similar finer cells. The computational time is . Note that for such a smooth flow problem, no troubled cells are marked by the shock detector, and the HR limiter is not activated. Nevertheless, we also test this case by applying the limiter on all cells artificially to validate its accuracy. The density errors and convergence orders are presented in Table 1 and 2. The expected order of accuracy is achieved by the current scheme whether or not the limiter is applied. With the limiter, the errors only slightly increase, which means that the high accuracy can still be maintained in smooth regions even if there are cells that are incorrectly marked as troubled cells for shock capturing problems. Thus, the dependence on the shock detector can be reduced effectively.
h  error  Order  error  Order 

0.5  7.05E05  1.20E04  
0.25  3.95E06  4.16  7.39E06  4.02 
0.125  2.69E07  3.87  5.23E07  3.82 
0.0625  1.87E08  3.84  3.64E08  3.85 
h  error  Order  error  Order 

0.5  1.27E04  2.80E04  
0.25  6.84E06  4.22  1.88E05  3.89 
0.125  3.37E07  4.34  9.66E07  4.29 
0.0625  2.47E08  3.77  5.07E08  4.25 
To estimate the efficiency of the current twostage SCFVGKS, Figure
5 shows the relations between the error and the CPU time. For comparison, the result obtained by the singlestage SCFVGKS is also presented, in which the CLS reconstruction is adopted. Because the computational cost introduced by the limiter has been reduced as much as possible by using the shock detector, for convenience, only the results without the limiter are considered here. It can be seen that the twostage fourthorder SCFVGKS is much more efficient than the two kinds of thirdorder SCFVGKS. The superiority becomes more distinct as the error decreases. Besides, the singlestage thirdorder SCFVGKS is more efficient than the twostage thirdorder SCFVGKS due to the use of the CLS reconstruction, rather than the WLS reconstruction adopted in this study. With the CLS reconstruction, the computational cost of flux can also be reduced effectively for the singlestage SCFVGKS. Nevertheless, in the current study, we focus on developing the temporalspatial fourthorder SCFVGKS, which is much more efficient. Owing to the twostage temporal discretization, a simpler and more efficient way of achieving fourthorder accuracy is constructed for the SCFV method. More detailed investigations about the influence of different reconstruction techniques are still needed in the future.iii.2 TirarevToro problem
To confirm the ability of the current scheme to capture highfrequency waves, the TirarevToro problem Titarev and Toro (2014) was tested. The initial condition is
(27) 
The computational domain is with the mesh obtained using the simple triangulation of a rectangular mesh. Figure 6 presents the density distribution at along the horizontal centerline with the mesh size and . The reference data is computed by the current scheme with . The highfrequency waves are captured by the current scheme accurately, especially with , which the density profile matches very well with the reference data. With the same mesh size , the result is slightly better than that obtained by the 6thorder and 8thorder compact GKS Zhao et al. (2019). Even the result with is comparable with that obtained by the WENO7JS Zhao et al. (2019) with . The results demonstrate the high resolution of the current scheme.
iii.3 Blast wave problem
Next, the blast wave problem Woodward and Colella (1984) was tested to validate the robustness of numerical schemes for capturing extremely strong shock waves. The initial condition is
(28) 
The computational domain is with the mesh obtained by a simple triangulation of a rectangular mesh. The periodic boundary condition is applied to upper and lower boundaries, and the reflecting boundary condition is applied to the left and right boundaries. The reference data are computed by the onedimensional highorder GKS (HGKS) Q.B. Li et al. (2010) with the mesh size . Figure 7 shows the density distribution at along the horizontal centerline with the mesh size and . The interaction of strong shock waves is captured by the current scheme with no oscillations. The results with agree well with the reference data.
iii.4 Double Mach reflection
The double Mach reflection is a 2D benchmark case that is used to validate the robustness and resolution of numerical schemes for shock capturing Woodward and Colella (1984). A rightmoving shock wave with the Mach number , initially located at ,
(29) 
which impinges on a wedge and leads to the double Mach reflection. Figure 8 shows the computational domain along with a sample mesh with , where
indicates the size of main cells uniformly distributed in the region near the wedge, and coarser cells are set elsewhere. The reflecting boundary condition is applied to the wedge and upper boundary. The exact postshock condition is applied to the left boundary and the bottom boundary from
to . The density contours are shown in Figure 9 with the mesh size . The shock waves are captured sharply and the slip line is resolved with high resolution. Besides, Figure 10 shows the density contours near the Mach stem, where the solid black lines indicate main cells, and the dashed grey lines indicate subcells. It can be observed that the thickness of the shock wave has nearly the same size as the main cell. Because discontinuities can exist inside main cells, the subcell resolution is fully maintained for shock capturing. In contrast, for DG and CPR, the distributions inside main cells need to remain continuous, and the thickness of the shock waves is usually larger than the size of main cells. It is difficult to achieve the high resolution for shock waves as in the SCFV method.In addition, for comparison with traditional FVGKS, a recently proposed fourthorder FVGKS based on WENO reconstruction is considered Zhao et al. . As shown in Figure 11, with the same mesh size , the result obtained by the current scheme is much more accurate than that obtained by FVGKS owing to the subcell resolution of the current scheme. With a mesh size of , the size of CVs, i.e., subcells, can be considered to be . Therefore, the result is also compared with that obtained by FVGKS with . It can be observed that the vortex in the slip line captured by SCFVGKS has a larger scale than that captured by FVGKS. Note that only the averaged solutions are updated on each subcell (CV) in the current scheme, while both the averaged solutions and its slopes are updated on each CV in the FVGKS, which means that the number of DOFs is twice that for FVGKS, which can help to improve the resolution. However, in FVGKS, the reconstruction needs to be implemented for each CV, while in the current scheme, a common reconstruction is implemented for four CVs (or subcells), which is more efficient. Besides, only face neighbors are involved in the current scheme for the fourthorder reconstruction, while the stencil in FVGKS is beyond face neighbors. Therefore, the current scheme is more compact than FVGKS.
iii.5 Viscous shock tube flow
The viscous shock tube problem is a benchmark case for supersonic viscous flows Daru and Tenaud (2004), where the flow is bounded by a unit square, and complex unsteady interactions occur between the shock wave and the boundary layer, which requires numerical schemes with both a strong robustness and high resolution. The Reynolds number is set as with a constant dynamic viscosity , and the Prandtl number is . The computational domain is set as . The symmetrical condition is applied for the upper boundary while the nonslip and adiabatic conditions are applied for other boundaries. The initial condition is
(30) 
The density contours at are presented in Figure 12 with the mesh size and . The complex flow structures, including the lambda shock and the vortex configurations, are well captured by the current scheme. However, the result with is still too dissipative, especially for capturing the primary vortex. By reducing the mesh size to , the resolution is improved significantly. For additional verification, Figure 13 shows the density profile along the bottom wall. Table 3 presents the estimation of the height of the primary vortex. The results are compared with the reference data Zhou et al. (2018) provided by a HGKS with . As can be seen, the result with clearly deviates from the reference data in terms of both the density profile and the height of the primary vortex. In contrast, the result with is much more accurate and agrees well with the reference data. Moreover, this flow case is also computed by the compact fourthorder CLSFV and VFV methods Wang et al. (2017) with . The current result with is comparable with that computed by the CLSFV and VFV methods with . The predicted height of the primary vortex is even closer to the reference data. Considering the advantage of the subcell resolution, the current scheme can achieve accurate solutions with a much coarser mesh compared to traditional FV methods. The results here demonstrate the good performance of the current scheme in compressible viscous flows.
Scheme  h  Height 

SCFVGKS  1/60  0.142 
SCFVGKS  1/120  0.165 
Reference  1/1500  0.166 
Iv Conclusions
A compact twostage fourthorder gaskinetic SCFV method is developed to solve compressible flows on triangular meshes. The difficulty of compactness faced by traditional FV methods is overcome with the SCFV method by subdividing each cell into a set of subcells in order to increase the number of DoFs for the representation of solution polynomials. The reconstruction is also more efficient than traditional FV because a set of subcells (CVs) share a common reconstruction. Compared to the singlestage thirdorder SCFVGKS, both the accuracy and efficiency are improved significantly by combining the fourthorder compact reconstruction with the secondorder accurate flux evolution. With the twostage fourthorder temporal discretization, we only need to construct the secondorder gas distribution function, where the flux evaluation is simplified significantly and the computational cost of flux is reduced to a great extent. Besides, the robustness of SCFVGKS is enhanced. More importantly, the twostage temporal discretization provides a more efficient approach to achieve fourthorder time accuracy. Compared to the fourthstage Runge–Kutta method, one half of the stages can be saved. With the gaskinetic flux, there is no need to compute the viscous flux for solving viscous flows. Several benchmark cases are tested to verify the performance of the current scheme in compressible flows. The high accuracy and efficiency are validated. For shock capturing, the current scheme suppresses oscillations near shock waves effectively; meanwhile, high accuracy can be remained in smooth flow regions. Because discontinuities can exist inside each cell, the subcell resolution is fully maintained, and thus shock waves can be resolved sharply. Hence, the current twostage fourthorder SCFVGKS is very promising for the practical simulation of compressible flows. The extension of the proposed method to other types of 2D and 3D unstructured meshes is under consideration in the near future.
Acknowledgements.
Qibing Li was supported by the National Natural Science Foundation of China (Nos. 11672158, 91852109) and the National Key Basic Research and Development Program (No. 2014CB744100). Peng Song was supported by the National Natural Science Foundation of China (No. 12031001) and the CAEP foundation (No. CX20200026). Jiequan Li was supported by the National Natural Science Foundation of China (Nos. 11771054, 12072042, 91852207), the SinoGerman Research Group Project (No. GZ1465), and the National Key Project (GJXM92579). The authors would like to thank the technical support of PARATERA and the “Explorer 100” cluster system of Tsinghua National Laboratory for Information Science and Technology.Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
references
References
 On essentially nonoscillatory schemes on unstructured meshes: analysis and implementation. J. Comput. Phys. 114, pp. 45–58. Cited by: §I.
 A direct Eulerian GRP scheme for compressible fluid flows. J. Comput. Phys. 218, pp. 19–34. Cited by: §I.
 Hyperbolic balance laws: Riemann invariants and the generalized Riemann problem. Numer. Math. 106, pp. 369–425. Cited by: §I.
 A threedimensional unified gaskinetic waveparticle solver for flow computation in all regimes. Phys. Fluids 32 (9), pp. 096108. Cited by: §I.
 The RungeKutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J. Comput. Phys. 141, pp. 199–224. Cited by: §I.
 High order onestep monotonicitypreserving schemes for unsteady compressible flow calculations. J. Comput. Phys. 193, pp. 563–594. Cited by: §III.5.
 A unified framework for the construction of onestep finite volume and discontinuous Galerkin schemes on unstructured meshes. J. Comput. Phys. 227, pp. 8209–8253. Cited by: §I.
 The study of shallow water flow with bottom topography by highorder compact gaskinetic scheme on unstructured mesh. Phys. Fluids 33, pp. 083613. Cited by: §I.

Weighted essentially nonoscillatory schemes for the interpolation of mean values on unstructured grids
. J. Comput. Phys. 144, pp. 194–212. Cited by: §I.  [10] A flux reconstruction approach to highorder schemes including discontinuous Galerkin methods. AIAA Paper 20074079. Cited by: §I.
 A HWENO reconstruction based highorder compact gaskinetic scheme on unstructured mesh. J. Comput. Phys. 410, pp. 109367. Cited by: §I.
 Highorder gas kinetic flux solver for simulation of two dimensional incompressible flows. Phys. Fluids 33, pp. 017107. Cited by: §I.
 An improved discrete gaskinetic scheme for twodimensional viscous incompressible and compressible flows. Phys. Fluids 31, pp. 066103. Cited by: §I.
 A twostage fourth order timeaccurate discretization for LaxWendroff type flow solvers, I: hyperbolic conservation laws. SIAM J. Sci. Comput. 38 (5), pp. A3046–A3069. Cited by: §I, §II.2.
 Highorder kexact WENO finite volume schemes for solving gas dynamic Euler equations on unstructured grids. Int. J. Numer. Methods Fluids 70 (6), pp. 742–763. Cited by: §II.1.
 A highorder accurate gaskinetic scheme for oneand twodimensional flow simulation. Commun. Comput. Phys. 15 (4), pp. 911–943. Cited by: §I.
 A model for collision processes in gases I: small amplitude processes in charged and neutral onecomponent systems. Phys. Rev. 94, pp. 511–525. Cited by: §II.2.
 High order subcell finite volume schemes for solving hyperbolic conservation laws II: extension to twodimensional systems on unstructured grids. J. Comput. Phys. 338, pp. 165–198. Cited by: §I, §I, §II.1.
 High order subcell finite volume schemes for solving hyperbolic conservation laws I: basic formulation and onedimensional analysis. Sci. China, Phys. Mech. Astron. 60 (8), pp. 084711. Cited by: §I, §I, §II.1.
 An efficient and accurate twostage fourthorder gaskinetic scheme for the Euler and NavierStokes equations. J. Comput. Phys. 326, pp. 197–221. Cited by: §I, §II.2.
 A thirdorder compact gaskinetic scheme on unstructured meshes for compressible NavierStokes solutions. J. Comput. Phys. 318, pp. 327–348. Cited by: §I.
 On the multidimensional gaskinetic BGK scheme. J. Comput. Phys. 220, pp. 532–548. Cited by: §I.
 A highorder gaskinetic NavierStokes flow solver. J. Comput. Phys. 229, pp. 6715–6731. Cited by: §I, §III.3.
 A multidimensional highorder discontinuous Galerkin method based on gas kinetic theory for viscous flow computations. J. Comput. Phys. 292, pp. 176–193. Cited by: §I.
 Essentially nonoscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws, Lecture Notes in Mathematics. Springer. Cited by: §III.1.
 [26] Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA Paper 19900013. Cited by: §I.
 Finite volume WENO schemes for threedimensional conservation laws. J. Comput. Phys. 201, pp. 238–260. Cited by: §III.2.
 Compact high order finite volume method on unstructured grids II: extension to twodimensional Euler equations. J. Comput. Phys. 314, pp. 883–908. Cited by: §I.
 Compact high order finite volume method on unstructured grids III: variational reconstruction. J. Comput. Phys. 337, pp. 1–26. Cited by: §I, §III.5.
 The numerical simulation of twodimensional fluid flow with strong shocks. J. Comput. Phys. 54, pp. 115–173. Cited by: §III.3, §III.4.
 A unified gaskinetic scheme for continuum and rarefied flows. J. Comput. Phys. 229, pp. 7747–7764. Cited by: §I.
 A gaskinetic BGK scheme for the NavierStokes equations and its connection with artificial dissipation and Godunov method. J. Comput. Phys. 171, pp. 289–335. Cited by: §I.
 Regularization of the Chapman–Enskog expansion and its description of shock structure. Phys. Fluids 14, pp. L17. Cited by: §II.2.
 Direct Modeling for Computational Fluid Dynamics: Construction and Application of Unified Gaskinetic Schemes. World Scientific. Cited by: §II.2.
 Hierarchical reconstruction for discontinuous Galerkin methods on unstructured grids with a WENOtype linear reconstruction and partial neighboring cells. J. Comput. Phys. 228, pp. 2194–2212. Cited by: §II.1.
 Highorder gaskinetic scheme on threedimensional unstructured meshes for compressible flows. Phys. Fluids 33, pp. 096102. Cited by: §I.
 A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. J. Comput. Phys. 228, pp. 8161–8186. Cited by: §I.
 Spectral (finite) volume method for conservation laws on unstructured grids IV: extension to twodimensional systems. J. Comput. Phys. 194 (2), pp. 716–741. Cited by: §I.
 Discrete gaskinetic schemebased arbitrary Lagrangian–Eulerian method for moving boundary problems. Phys. Fluids 33, pp. 067101. Cited by: §I.
 A thirdorder gaskinetic CPR method for the Euler and NavierStokes equations on triangular meshes. J. Comput. Phys. 363, pp. 329–353. Cited by: §I, §II.2.
 A thirdorder subcell finite volume gaskinetic scheme for the Euler and NavierStokes equations on triangular meshes. J. Comput. Phys. 436, pp. 110245. Cited by: §I, §II.1, §II.1.
 Discrete unified gas kinetic scheme for flows of binary gas mixture based on the McCormack model. Phys. Fluids 31, pp. 017101. Cited by: §I.
 [43] Compact highorder gaskinetic scheme on unstructured mesh for acoustic and shock wave computations. arXiv preprint arXiv:2010.05717v2 (2020). Cited by: §II.1, §III.4.
 Compact higherorder gaskinetic schemes with spectrallike resolution for compressible flow simulations. Adv. Aerodyn. 1, pp. 13. Cited by: §I, §III.2.
 A simplified discrete unified gas–kinetic scheme for compressible flow. Phys. Fluids 33, pp. 036103. Cited by: §I.
 Gridconverged solution and analysis of the unsteady viscous flow in a twodimensional shock tube. Phys. Fluids 30, pp. 016102. Cited by: §III.5.
Comments
There are no comments yet.