Increasing engineering demands for accuracy and efficiency have catalyzed the development of high-order 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. Low-order methods are usually too dissipative to resolve small-scale structures. In contrast, high-order methods can potentially achieve increased accuracy with low computational cost. However, existing high-order methods are generally less robust and more complicated than low-order methods. Therefore, it is necessary to develop new high-order 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 high-speed flows containing discontinuities (e.g., shocks). Many high-order FV methods were developed on unstructured meshes, such as the k-exact method T.J Barth and P.O. Frederickson , the essentially non-oscillatory (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 non-compactness resulting from a large stencil for reconstruction leads to problems with numerical dissipation, high-order numerical approximations of boundary conditions, reduced parallel efficiency, and caches missing.
To achieve compactness for high-order FV methods, the compact least-square reconstruction and variational reconstruction Wang et al. (2016, 2017) were developed, where a large linear equation system should be solved and implicit time-stepping 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 subcell-averaged 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 face-neighboring 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 high-order accuracy with a compact stencil. In summary, the SCFV method has the following three advantages over existing high-order 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 face-neighboring 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 single-stage third-order gas-kinetic flux solver (GKS) Zhang and Q.B. Li (2021). Although other solvers, such as the GRP solver, Ben-Artzi et al. (2006); Ben-Artzi and Li (2007) can also be adopted, this study still uses the GKS solver but with second-order 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 gas-kinetic 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, high-order GKS has been developed systematically. By taking a second-order Taylor expansion, a third-order multi-dimensional GKS was developed within a single stage Q.B. Li et al. (2010). The time-dependent gas distribution function can even be used to evolve the solutions at cell interfaces, which provides additional DoFs to help construct a compact third-order FV-GKS Pan and Xu (2016); Ji et al. (2020). Other extensions of single-stage third-order 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 single-stage third-order SCFV-GKS has been proposed recently Zhang and Q.B. Li (2021). However, when fourth-order accuracy is considered, the single-stage flux function is quite complicated, and the computational cost increases significantly Liu and Tang (2014); the robustness may also be weakened. Fortunately, the two-stage fourth-order time-stepping method Li and Du (2016) provides an efficient way for Lax-Wendroff-type solvers to achieve fourth-order time accuracy. Only the second-order gas-kinetic 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 single-stage high-order GKSPan et al. (2016).
Hence, this study developed a two-stage temporal-spatial fourth-order SCFV-GKS to solve compressible flows over triangular meshes, with the aim of enhancing the compactness, efficiency, and robustness. Compared to the single-stage SCFV-GKS, both accuracy and efficiency are improved significantly. By adopting the two-stage 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 second-order gas-kinetic flux. These improvements make the SCFV method more promising for engineering applications.
This paper is organized as follows. Section 2 presents the compact fourth-order reconstruction based on the SCFV method and the two-stage gas-kinetic 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 solver-based subcell finite volume method
In this section, we propose a GKS-based SCFV method for compressible fluid flows. The primary steps are the compact subcell technique and the two-stage fourth-order temporal advancing based on the second-order GKS.
ii.1 Subcell finite volume method
First, we briefly review the SCFV method. Consider the two-dimensional (2D) conservation law
where are the conservative variables, where are the macroscopic velocities and
is the flux vector. The computational domain is divided intonon-overlapping 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, subcell-averaged 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 fourth-order 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 face-neighboring cells, there are 16 subcells available to sufficiently determine the unknowns.
In comparison, Figure 2 also shows the stencil used in a GKS-based FV method and the classical k-exact FV Zhao et al. . Note that for the GKS-based FV, cell-averaged slopes are updated and participate in the reconstruction, which helps to avoid the large stencil for the k-exact FV. Nonetheless, with only face neighbors, it is still not enough for a fourth-order 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
where , indicates the centroid of . The zero-mean 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
where and satisfy . is the volume of . Integrating over gives
and is the averaged solution on the subcell . Note that because has been used in Eq.(2) and it has the relation with as
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 over-determined system. The unknowns are solved by the weighted least-square technique
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 least-square reconstruction, i.e.,
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
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 second-order gas distribution function is much simpler than the third-order one, and the computational cost for numerical flux construction has been reduced significantly, we simply take the weighted least-square (WLS) technique to determine the unknowns, rather than the constrained least-square (CLS) technique adopted in the single-stage SCFV-GKS 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 third-order gas-kinetic flux. However, the technique itself is more complicated and time-consuming, 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 second-order gas-kinetic 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 face-neighboring 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 Two-stage gas-kinetic flux evolution
By integrating Eq.(1) over each subcell , the semi-discrete form of the SCFV framework can be obtained as follows.
where is the area of , is the interfaces surrounding , and is the outward unit normal vector. To maintain fourth-order space accuracy, the line integral on the right-hand side is discretized using the Gaussian quadrature rule with two points
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 time-evolving flux function provided by the second-order gas-kinetic flux solver, fourth-order time accuracy can be achieved more efficiently with the help of the two-stage temporal discretization Li and Du (2016), which can be expressed with the current notations as
where is the flux at the intermediate stage in the time interval . The aforementioned two-stage temporal discretization has been proved to achieve fourth-order 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 gas-kinetic theory, both the macroscopic conservative variables and the numerical flux can be obtained from the gas distribution function with the following relations
where is a function of physical space , time , particle velocity , and the internal DoFs ,
is the vector of moments, andis 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)
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
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 Naiver-Stokes equations can be recovered through the first-order Chapman-Enskog expansion Xu (2002, 2015)
By eliminating the first-order term in Eq.(16), the Euler equations can be recovered as well. To compute the numerical flux through Eq.(13), a time-dependent gas distribution function is constructed at each flux point based on the analytical solution of Eq.(14)
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 first-order Taylor expansion of and around the flux point and combined with Eq.(16), the time-dependent gas distribution function can be obtained as follows.
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 gas-kinetic 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 third-order gas-kinetic 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 time-dependent flux function can be constructed according to Eq.(13), denoted as , which is a non-linear function of . A simple fitting method is adopted to obtain the approximated flux, and its first-order time derivative Pan et al. (2016), where is approximated by a linear function within the time interval . Denote the time integration of within as
Then, an equation set can be obtained
Similarly, the approximated flux and its time derivative can be obtained by simply replacing the superscript in Eq.(21) with . Finally, the two-stage time stepping method for updating subcell averages can be summarized as
where the flux at each flux point is respectively computed according to
Although there are two stages in the current scheme, the flux evaluation can still be more efficient than the previous single-stage SCFV-GKS because the second-order gas-kinetic flux is much cheaper. Besides, the robustnes of SCFV-GKS 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
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 high-order FV methods. The ratio of specific heat is set to . In the accuracy test, the error and error are computed by
where and indicate the averaged numerical solution and the analytical solution on subcell , respectively, and is the number of main cells. The two-dimensional (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.,
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.
To estimate the efficiency of the current two-stage SCFV-GKS, Figure5 shows the relations between the error and the CPU time. For comparison, the result obtained by the single-stage SCFV-GKS 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 two-stage fourth-order SCFV-GKS is much more efficient than the two kinds of third-order SCFV-GKS. The superiority becomes more distinct as the error decreases. Besides, the single-stage third-order SCFV-GKS is more efficient than the two-stage third-order SCFV-GKS 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 single-stage SCFV-GKS. Nevertheless, in the current study, we focus on developing the temporal-spatial fourth-order SCFV-GKS, which is much more efficient. Owing to the two-stage temporal discretization, a simpler and more efficient way of achieving fourth-order 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 Tirarev-Toro problem
To confirm the ability of the current scheme to capture high-frequency waves, the Tirarev-Toro problem Titarev and Toro (2014) was tested. The initial condition is
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 high-frequency 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 6th-order and 8th-order compact GKS Zhao et al. (2019). Even the result with is comparable with that obtained by the WENO-7JS 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
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 one-dimensional high-order 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 right-moving shock wave with the Mach number , initially located at ,
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 post-shock condition is applied to the left boundary and the bottom boundary fromto . 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 FV-GKS, a recently proposed fourth-order FV-GKS 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 FV-GKS 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 FV-GKS with . It can be observed that the vortex in the slip line captured by SCFV-GKS has a larger scale than that captured by FV-GKS. 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 FV-GKS, which means that the number of DOFs is twice that for FV-GKS, which can help to improve the resolution. However, in FV-GKS, 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 fourth-order reconstruction, while the stencil in FV-GKS is beyond face neighbors. Therefore, the current scheme is more compact than FV-GKS.
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 non-slip and adiabatic conditions are applied for other boundaries. The initial condition is
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 fourth-order 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.
A compact two-stage fourth-order gas-kinetic 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 single-stage third-order SCFV-GKS, both the accuracy and efficiency are improved significantly by combining the fourth-order compact reconstruction with the second-order accurate flux evolution. With the two-stage fourth-order temporal discretization, we only need to construct the second-order 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 SCFV-GKS is enhanced. More importantly, the two-stage temporal discretization provides a more efficient approach to achieve fourth-order time accuracy. Compared to the fourth-stage Runge–Kutta method, one half of the stages can be saved. With the gas-kinetic 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 two-stage fourth-order SCFV-GKS 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 Sino-German 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.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
- On essentially non-oscillatory 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 three-dimensional unified gas-kinetic wave-particle solver for flow computation in all regimes. Phys. Fluids 32 (9), pp. 096108. Cited by: §I.
- The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J. Comput. Phys. 141, pp. 199–224. Cited by: §I.
- High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations. J. Comput. Phys. 193, pp. 563–594. Cited by: §III.5.
- A unified framework for the construction of one-step 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 high-order compact gas-kinetic scheme on unstructured mesh. Phys. Fluids 33, pp. 083613. Cited by: §I.
Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids. J. Comput. Phys. 144, pp. 194–212. Cited by: §I.
-  A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. AIAA Paper 2007-4079. Cited by: §I.
- A HWENO reconstruction based high-order compact gas-kinetic scheme on unstructured mesh. J. Comput. Phys. 410, pp. 109367. Cited by: §I.
- High-order gas kinetic flux solver for simulation of two dimensional incompressible flows. Phys. Fluids 33, pp. 017107. Cited by: §I.
- An improved discrete gas-kinetic scheme for two-dimensional viscous incompressible and compressible flows. Phys. Fluids 31, pp. 066103. Cited by: §I.
- A two-stage fourth order time-accurate discretization for Lax-Wendroff type flow solvers, I: hyperbolic conservation laws. SIAM J. Sci. Comput. 38 (5), pp. A3046–A3069. Cited by: §I, §II.2.
- High-order k-exact 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 high-order accurate gas-kinetic scheme for one-and two-dimensional 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 one-component systems. Phys. Rev. 94, pp. 511–525. Cited by: §II.2.
- High order sub-cell finite volume schemes for solving hyperbolic conservation laws II: extension to two-dimensional systems on unstructured grids. J. Comput. Phys. 338, pp. 165–198. Cited by: §I, §I, §II.1.
- High order sub-cell finite volume schemes for solving hyperbolic conservation laws I: basic formulation and one-dimensional analysis. Sci. China, Phys. Mech. Astron. 60 (8), pp. 084711. Cited by: §I, §I, §II.1.
- An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Euler and Navier-Stokes equations. J. Comput. Phys. 326, pp. 197–221. Cited by: §I, §II.2.
- A third-order compact gas-kinetic scheme on unstructured meshes for compressible Navier-Stokes solutions. J. Comput. Phys. 318, pp. 327–348. Cited by: §I.
- On the multidimensional gas-kinetic BGK scheme. J. Comput. Phys. 220, pp. 532–548. Cited by: §I.
- A high-order gas-kinetic Navier-Stokes flow solver. J. Comput. Phys. 229, pp. 6715–6731. Cited by: §I, §III.3.
- A multi-dimensional high-order discontinuous Galerkin method based on gas kinetic theory for viscous flow computations. J. Comput. Phys. 292, pp. 176–193. Cited by: §I.
- Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Lecture Notes in Mathematics. Springer. Cited by: §III.1.
-  Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA Paper 1990-0013. Cited by: §I.
- Finite volume WENO schemes for three-dimensional 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 two-dimensional 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 two-dimensional fluid flow with strong shocks. J. Comput. Phys. 54, pp. 115–173. Cited by: §III.3, §III.4.
- A unified gas-kinetic scheme for continuum and rarefied flows. J. Comput. Phys. 229, pp. 7747–7764. Cited by: §I.
- A gas-kinetic BGK scheme for the Navier-Stokes 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 Gas-kinetic Schemes. World Scientific. Cited by: §II.2.
- Hierarchical reconstruction for discontinuous Galerkin methods on unstructured grids with a WENO-type linear reconstruction and partial neighboring cells. J. Comput. Phys. 228, pp. 2194–2212. Cited by: §II.1.
- High-order gas-kinetic scheme on three-dimensional 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 two-dimensional systems. J. Comput. Phys. 194 (2), pp. 716–741. Cited by: §I.
- Discrete gas-kinetic scheme-based arbitrary Lagrangian–Eulerian method for moving boundary problems. Phys. Fluids 33, pp. 067101. Cited by: §I.
- A third-order gas-kinetic CPR method for the Euler and Navier-Stokes equations on triangular meshes. J. Comput. Phys. 363, pp. 329–353. Cited by: §I, §II.2.
- A third-order subcell finite volume gas-kinetic scheme for the Euler and Navier-Stokes 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.
-  Compact high-order gas-kinetic scheme on unstructured mesh for acoustic and shock wave computations. arXiv preprint arXiv:2010.05717v2 (2020). Cited by: §II.1, §III.4.
- Compact higher-order gas-kinetic schemes with spectral-like 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.
- Grid-converged solution and analysis of the unsteady viscous flow in a two-dimensional shock tube. Phys. Fluids 30, pp. 016102. Cited by: §III.5.