An efficient implementation of compact third-order implicit reconstruction solver with a simple WBAP limiter for compressible flows on unstructured meshes

Abstract

In this paper, we present the development of a third-order density-based solver within the OpenFOAM framework, tailored for handling compressible flows. The solver incorporates implicit variational reconstruction on three-dimensional unstructured meshes, as well as a novel technology coupling reconstruction and time integration for both steady and unsteady simulations. To address the challenge of achieving high-order accuracy for curved geometries, we introduce a new approach for curved wall boundary reconstruction, specifically designed for situations where high-order mesh information is not readily available in OpenFOAM. Furthermore, we propose a simple WBAP limiter capable of capturing shocks without necessitating the whole-domain successive limiting procedure. Numerical tests were conducted to assess the solver's performance. The results reveal that our established solver exhibits higher accuracy in smooth flow simulations, while maintaining an excellent balance between accuracy and robustness for problems involving strong shocks and other complex flow structures.

KEYWORDS:

• High-order scheme
• Curved wall boundary
• WBAP limiter
• Compressible flows
• Three-dimensional unstructured meshes

1. Introduction

Computational fluid dynamics (CFD) solvers, utilizing finite volume (FV) methods on unstructured meshes, play a crucial role in aerospace engineering research and aircraft design due to their capability to effectively handle complex geometries. Among the widely-used CFD libraries, OpenFOAM (Jasak, 1996) stands out as a prominent C++ based framework. Notably, several density-based solvers have been developed within OpenFOAM, including $rhoCentralFoam$ (Kurganov & Tadmor, 2000) in the standard package, $blastFoam$ (Heylmun et al., n.d.), and $RGDFoam$ (Qi et al., 2022). Despite the widespread use of the solvers developed in this framework (R. Li et al., 2022), it is crucial to acknowledge their reliance on second-order FV methods. This dependence gives rise to a notable limitation as second-order schemes are known to introduce considerable numerical dissipation, thereby compromising the accuracy of simulations when dealing with essential flow phenomena such as shocks, vortices, and turbulence.

To mitigate these limitations and improve simulation accuracy, there has been a surge of interest in the development of high-order FV methods. These advanced techniques aim to effectively capture intricate flow structures while preserving computational efficiency. Substantial further research and development in this domain are imperative to propel the capabilities of CFD solvers and enable more precise and dependable predictions in aerospace engineering and other fluid dynamics applications. Similar to finite difference (FD) methods (Dong et al., 2023; Zhou & Dong, 2021), FV methods are utilized to solve equations exclusively for the cell average of conservative variables or primitive variables. However, for attaining higher accuracy in simulations, both FD and FV methods necessitate the incorporation of larger stencils that encompass an increased number of neighbouring cells during the reconstruction process. One effective approach to attain this increased accuracy is through the implementation of weighted essentially non-oscillatory (WENO) schemes, as demonstrated in studies by (Ji et al., 2022, 2023; Tsoutsanis, 2019). Nonetheless, the extensive reconstruction stencil presents challenges in coding complexity and obstructs parallel computation, especially at inter-processor boundaries (Busto et al., 2020; Godunov, 1959; Tsoutsanis & Dumbser, 2021). Furthermore, the extensive reconstruction stencil cannot be directly applied to physical boundaries. Over the past two decades, research into high-order methods with compact stencils involving only face-neighbouring cells has emerged as a crucial area of interest, exemplified by the discontinuous Galerkin (DG) method (Proskurov et al., 2022).

Compared to the successful implementation of high-order time discretization (Mangani et al., 2015), achieving high-order spatial discretization within the OpenFOAM framework poses significant challenges. However, in recent years, notable advancements have been made through the development of high-order solvers, showcasing superior performance compared to traditional second-order solvers. Among these solvers, the first class (Vevek, 2020) utilizes the FD method, efficiently achieving high-order accuracy and demonstrating the ability to simulate problems with intricate flow structures. Nevertheless, these solvers have a limitation, as they are not directly applicable to unstructured meshes. An alternative class of solvers (Xu et al., 2017) relies on the DG method, which exhibits the capability to address flow problems on unstructured meshes. However, integrating the DG method into the OpenFOAM framework poses challenges, particularly in adapting pre-processing and post-processing interfaces. This article's primary objective is to develop a high-order density-based unstructured solver capable of simulating compressible flows.

Considering the interfaces of OpenFOAM, it is imperative to adopt reconstruction with a compact stencil for compressible solvers, as only direct face neighbouring information is accessible during communication across inter-processor boundaries. This limitation is the primary reason why high-order FV unstructured CFD solvers (Antoniadis et al., 2022) are seldom developed in OpenFOAM. In response to the challenges posed by the extensive stencil of high-order reconstruction, a series of compact implicit reconstruction procedures on unstructured meshes have been recently introduced (Wang et al., 2016a, 2016b, 2017). These implicit reconstruction procedures offer an efficient solution on the compact stencil. Furthermore, variational reconstruction procedures are specifically designed to address some of the key deficiencies associated with non-compact reconstruction procedures.

High-order accurate implementation of wall boundary in curved geometries is crucial, several approaches (Krivodonova & Berger, 2006; Lübon et al., 2006) have been presented for DG methods, which are particularly sensitive to the accuracy of curved wall boundaries. To achieve the desired solution accuracy, we employ the approach proposed by (Lübon et al., 2006) to reconstruct the curved boundary from the flat boundary for high-order FV methods. This step is necessary as high-order mesh information is currently unavailable in the OpenFOAM framework. Importantly, this approach does not require an exact description of the geometry surface.

A major concern with FV schemes that have accuracy higher than first order is that the polynomials of the solution cannot be directly used for convective flux evaluation due to the Gibbs phenomenon arising around the discontinuities. A common approach to address this concern is to introduce a limiter that reconstructs the polynomials of troubled cells. Traditionally, slope limiters like the Barth–Jespersen limiter (Barth & Jespersen, 1989) and Venkatakrishnan limiter (Venkatakrishnan, 1995) have been used in FV methods. However, these slope limiters are designed to reconstruct coefficients based on the bound of cell average around the troubled cell, which is not suitable for high-order schemes where the solution is non-monotonic in the control volume. Therefore, the limiting procedure needs to be carefully designed, and it has become a research focus to achieve the required robustness for high-order coefficients.

To ensure efficient computation, problem-independent troubled cell indicators (Fu & Shu, 2017; W. Li & Ren, 2012a; Qiu & Shu, 2005) are applied to cells near the discontinuities in high-order methods, while keeping the smooth domain with the original solution. However, many investigations have shown that achieving robustness near discontinuities and preserving the designed accuracy of cells in smooth extrema, which are often misidentified as troubled cells, is challenging with the limiting procedure. In this context, the weighted biased averaging procedure (WBAP) proposed in (W. Li et al., 2011; W. Li & Ren, 2012b) serves as an effective limiter for high-order FV methods and has successfully been applied to solve hyperbolic conservation laws. Similar to the WENO schemes, the WBAP limiter requires multiple candidate reconstructions, including primary reconstruction (PR) of the current cell and secondary reconstructions (SRs) of face-neighbouring cells proposed in (W. Li & Ren, 2012a). The successive limiting procedure, where the limit procedure is applied from high-order coefficients to low-order coefficients respectively, along with the WBAP limiter, effectively controls numerical oscillations.

For compressible flow simulations, the limiting procedure can significantly improve the robustness of high-order methods when coupled with characteristic-wise projection (W. Li & Ren, 2012b). This coupling is necessary because the solution within the control volume is non-monotonic, and direct application on conservation leads to numerical instability. However, the combination of projection and successive limiting procedures can result in high computational costs. To address this challenge, a new simple WBAP limiter is proposed, eliminating the need for the whole-domain successive limiting procedure while maintaining the same level of robustness as the original limiter. In this new limiter, the candidate reconstructions for a shared Jacobian matrix include the PR of the current cell and a SR of one of its face-neighbouring cells. Additionally, to prevent negative density and pressure solutions, a positivity-preserving limiter is applied after employing the new simple WBAP limiter. The limiting procedure introduced in this article successfully achieves a desirable balance between accuracy and robustness, making it a valuable tool for compressible flow simulations with high-order methods.

The remainder of this paper is organized as follows. Section 2.1 presents the fundamental idea of high-order FV methods on unstructured grids, while Section 2.2 introduces the optimized weighted variational reconstruction FV scheme. In Section 2.3, the implicit time integration is described. The implementation of curved wall boundary reconstruction is discussed in Section 3. Section 4 covers the limiting procedure, including a shock detector, a new simple WBAP limiter, and a positivity-preserving limiter. Section 5 presents the numerical results and discussions on the implementation of the CFD solver using the methods proposed in this article. Finally, the conclusion is given in Section 6.

2. High-order finite volume method

2.1. General framework

In this section, the framework for the high-order FV method is introduced below. The time-independent Euler equations can be written in conversation form (1) $\begin{array}{r}\frac{\mathrm{\partial }\mathit{U}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot {\mathit{F}}^c\left(\mathit{U}\right)=0,\end{array}$(1) where $\mathit{U}=[\rho ,\rho u,\rho v,\rho w,\rho E{]}^T$ and the Cartesian components ${\mathit{F}}^c=\left[{\mathit{F}}_1^c,{\mathit{F}}_2^c,{\mathit{F}}_3^c\right]$ are the components defined as $\begin{array}{rl}{\mathit{F}}_1^c& =\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uw\\ \rho u(E+p/\rho )\end{array}\right],{\mathit{F}}_2^c=\left[\begin{array}{c}\rho v\\ \rho vu\\ \rho v^2+p\\ \rho vw\\ \rho v(E+p/\rho )\end{array}\right],\\ {\mathit{F}}_3^c& =\left[\begin{array}{c}\rho w\\ \rho wu\\ \rho wv\\ \rho w^2+p\\ \rho w(E+p/\rho )\end{array}\right],\end{array}$where $\rho$ represents the density and $p$ is the pressure. The components $\left[u,v,w\right]$ correspond to the velocity in the $\left[x,y,z\right]$ directions within the Cartesian coordinate system.

The general FV formulation for Equation (1) can be obtained from the control volume balance equation as (2) $\begin{array}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }t}{\int }_{{\mathrm{\Omega }}_i}{\mathit{U}}_{i}d\mathrm{\Omega }+{\oint }_{\mathrm{\partial }{\mathrm{\Omega }}_i}\mathit{H}\cdot \mathit{n}d\mathrm{\Gamma }=0,\end{array}$(2) where ${\mathrm{\Omega }}_i$ is the current cell, $\mathrm{\partial }{\mathrm{\Omega }}_i$ is the cell boundary of ${\mathrm{\Omega }}_i$, and $\mathit{H}=[{\mathit{F}}_1^c,{\mathit{F}}_2^c,{\mathit{F}}_3^c{]}^T$ is the numerical flux on the cell face.

In numerical methods, we turn Equation (2) into the semi-discrete form as (3) $\begin{array}{r}{\overline{\mathrm{\Omega }}}_i\frac{\mathrm{\partial }{\overline{\mathit{U}}}_i}{\mathrm{\partial }t}=-\sum _{m=1}^M{\int }_{I_m}\mathit{H}\cdot \mathit{n}d\mathrm{\Gamma },\end{array}$(3) where ${\overline{\mathit{U}}}_i=\frac{1}{{\overline{\mathrm{\Omega }}}_i}{\int }_{{\mathrm{\Omega }}_i}{\mathit{U}}_{i}d\mathrm{\Omega }$ donates the average of solution on ${\mathrm{\Omega }}_i$ and ${\overline{\mathrm{\Omega }}}_i$ represents the cell volume. The right-hand side of Equation (3) can be computed according the form as (4) $\begin{array}{r}{\int }_{I_m}\mathit{H}\cdot \mathit{n}d\mathrm{\Gamma }=\left[\sum _{\phi }^{NG}{\omega }_{\phi }\mathit{H}\left(\mathit{U}\left({\mathit{x}}_{m,\phi }\right)\right)\right]\cdot {\mathit{n}}_m\mathrm{\Delta }{S}_m,\end{array}$(4) where $\mathit{H}\left(\mathit{U}\left({\mathit{x}}_{m,\phi }\right)\right)$ is evaluated on each Gaussian point of the shared interface by HLLC Riemann solver (Toro, 2009) in this paper., i.e. $\mathit{H}\left(\mathit{U}\left({\mathit{x}}_{m,\phi }\right)\right)\cdot {\mathit{n}}_m=\tilde{\mathit{H}}\left({\mathit{U}}^L\left({\mathit{x}}_{m,\phi }\right),{\mathit{U}}^R\left({\mathit{x}}_{m,\phi }\right),{\mathit{n}}_m\right),$where both sides of states ${\mathit{U}}^L\left({\mathit{x}}_{m,\phi }\right)$ and ${\mathit{U}}^R\left({\mathit{x}}_{m,\phi }\right)$ are evaluated respectively by the polynomials of the conservative variables on ${\mathrm{\Omega }}_i$ and ${\mathrm{\Omega }}_j$, which share the same interface $I_m$ with ${\mathrm{\Omega }}_i$. A semi-discrete form for time integration can be described as (5) $\begin{array}{r}\overline{\mathrm{\Omega }}\frac{\mathrm{\partial }\overline{\mathit{U}}}{\mathrm{\partial }t}=\mathit{R}\left(\mathit{U}\right),\end{array}$(5) where $\mathit{R}$ is computed from Equation (5). The time integration scheme to solve the ordinary differential equations (ODEs) to update the cell averages in this article will be described in Section 2.3.

2.2. Variational reconstruction

For third-order variational reconstruction described in this article, only the cells sharing the interfaces with the current cell are included in the compact stencil, for a typical hexahedral cell, denoted by $S_i=\left\{{\mathrm{\Omega }}_{j1},{\mathrm{\Omega }}_{j2},{\mathrm{\Omega }}_{j3},{\mathrm{\Omega }}_{j4},{\mathrm{\Omega }}_{j5},{\mathrm{\Omega }}_{j6}\right\}$ shown in Figure 1.

Figure 1. The compact stencils for interior cell (left) and boundary cell (right).

One of the components of conservative variables denoted by $u\left(\mathit{x}\right)$ is considered below. The conservation of the mean of the conservative variables requires (6) $\begin{array}{r}{\overline{u}}_i=\frac{1}{{\overline{\mathrm{\Omega }}}_i}{\int }_{{\mathrm{\Omega }}_{\text{i}}}u\left(\mathit{x}\right)d\mathrm{\Omega },\end{array}$(6) where ${\overline{u}}_i$ is the cell average of $u\left(\mathit{x}\right)$ and a degree $p$ polynomial (7) $\begin{array}{r}{u}_i\left(\mathit{x}\right)={\overline{u}}_i+\sum _{l=1}^{N\left(p\right)}{u}_i^l{\phi }_{l,i}\left(\mathit{x}\right),\end{array}$(7) is evaluated on ${\mathrm{\Omega }}_{\text{i}}$, the degree of the coefficients of degree $p$ polynomial is $N\left(p\right)=\frac{\left(p+1\right)\left(p+2\right)\left(p+3\right)}{6}-1$ in three-dimensional cases. For the presented third-order FV method in three-dimension, $N=9$. The non-dimensional Taylor basis functions are taken as $\begin{array}{rl}{\phi }_1\left(\mathit{x}\right)& =\frac{x-x_c}{\mathrm{\Delta }x},{\phi }_2\left(\mathit{x}\right)=\frac{y-y_c}{\mathrm{\Delta }y},{\phi }_3\left(\mathit{x}\right)=\frac{z-z_c}{\mathrm{\Delta }z}\\ {\phi }_4\left(\mathit{x}\right)& =\frac{{\phi }_1^2\left(\mathit{x}\right)-{\overline{\phi }}_{1,1}}{2},{\phi }_5\left(\mathit{x}\right)=\frac{{\phi }_2^2\left(\mathit{x}\right)-{\overline{\phi }}_{2,2}}{2},\\ {\phi }_6\left(\mathit{x}\right)& =\frac{{\phi }_3^2\left(\mathit{x}\right)-{\overline{\phi }}_{3,3}}{2},{\phi }_7\left(\mathit{x}\right)={\phi }_1\left(\mathit{x}\right){\phi }_2\left(\mathit{x}\right)-{\overline{\phi }}_{1,2},\\ {\phi }_8\left(\mathit{x}\right)& ={\phi }_1\left(\mathit{x}\right){\phi }_3\left(\mathit{x}\right)-{\overline{\phi }}_{1,3},\\ {\phi }_9\left(\mathit{x}\right)& ={\phi }_2\left(\mathit{x}\right){\phi }_3\left(\mathit{x}\right)-{\overline{\phi }}_{2,3},\\ \mathrm{\Delta }x& =\left(x_{max}-x_{min}\right)/2,\mathrm{\Delta }y=\left(y_{max}-y_{min}\right)/2,\\ \mathrm{\Delta }z& =\left(z_{max}-z_{min}\right)/2,\\ {\overline{\phi }}_{m,n}& =\frac{1}{{\overline{\mathrm{\Omega }}}_i}{\int }_{{\mathrm{\Omega }}_{\text{i}}}{\phi }_m\left(\mathit{x}\right){\phi }_n\left(\mathit{x}\right)d\mathrm{\Omega }.\end{array}$

Let the centroid of cell be ${\mathit{x}}_c=\left(x_c,y_c,z_c\right)$. The length scales (Luo et al., 2008) in the functions are defined to make to basis non-dimensional and lower the condition number for the iteration solver for reconstruction. $x_{max}$, $y_{max}$, $z_{max}$, $x_{min}$, $y_{min}$ and $z_{min}$ are computed from coordinates of the cell ${\mathrm{\Omega }}_i$ in $x$-, $y$- and $z$-direction. This formulation is simple and keeps the cell average unchanged. The following degree $p$ polynomial of the solution on ${\mathrm{\Omega }}_{\text{i}}$ can be described as (8) $\begin{array}{r}{u}_i\left(\mathit{x}\right)={\overline{u}}_i+\sum _{l=1}^9{u}_i^l{\phi }_{l,i}\left(\mathit{x}\right).\end{array}$(8)

The PR of the target cell and several SRs are necessary in limiting procedure. For example, a SR can be transformed from its PR of one of face-neighbouring cells of stencil $S_i$ and is denoted by (9) $\begin{array}{r}{u}_{j\to i}\left(\mathit{x}\right)={\overline{u}}_i+\sum _{l=1}^{N\left(p\right)}{u}_{j\to i}^l{\phi }_{l,i}\left(\mathit{x}\right),j\in S_i.\end{array}$(9) The specific formulation of $u_{j\to i}^l$ for the solution ($p=2$) are computed according to $\begin{array}{rl}{u}_{j\to i}^1& =\frac{\mathrm{\Delta }{x}_i}{\mathrm{\Delta }{x}_j}\left(u_j^1+\frac{\mathrm{\Delta }{x}_{ij}}{\mathrm{\Delta }{x}_j}{u}_j^4+\frac{\mathrm{\Delta }{y}_{ij}}{\mathrm{\Delta }{y}_j}{u}_j^7+\frac{\mathrm{\Delta }{z}_{ij}}{\mathrm{\Delta }{z}_j}{u}_j^8\right),\\ u_{j\to i}^2& =\frac{\mathrm{\Delta }{y}_i}{\mathrm{\Delta }{y}_j}\left(u_j^2+\frac{\mathrm{\Delta }{x}_{ij}}{\mathrm{\Delta }{x}_j}{u}_j^7+\frac{\mathrm{\Delta }{y}_{ij}}{\mathrm{\Delta }{y}_j}{u}_j^5+\frac{\mathrm{\Delta }{z}_{ij}}{\mathrm{\Delta }{z}_j}{u}_j^9\right),\\ u_{j\to i}^3& =\frac{\mathrm{\Delta }{z}_i}{\mathrm{\Delta }{z}_j}\left(u_j^3+\frac{\mathrm{\Delta }{x}_{ij}}{\mathrm{\Delta }{x}_j}{u}_j^8+\frac{\mathrm{\Delta }{y}_{ij}}{\mathrm{\Delta }{y}_j}{u}_j^9+\frac{\mathrm{\Delta }{z}_{ij}}{\mathrm{\Delta }{z}_j}{u}_j^6\right),\\ u_{j\to i}^4& =\frac{\mathrm{\Delta }{x}_i^2}{\mathrm{\Delta }{x}_j^2}{u}_j^4,u_{j\to i}^5=\frac{\mathrm{\Delta }{y}_i^2}{\mathrm{\Delta }{y}_j^2}{u}_j^5,u_{j\to i}^6=\frac{\mathrm{\Delta }{z}_i^2}{\mathrm{\Delta }{z}_j^2}{u}_j^6,\\ u_{j\to i}^7& =\frac{\mathrm{\Delta }{x}_i\mathrm{\Delta }{y}_i}{\mathrm{\Delta }{x}_j\mathrm{\Delta }{y}_j}{u}_j^7,u_{j\to i}^8=\frac{\mathrm{\Delta }{x}_i\mathrm{\Delta }{z}_i}{\mathrm{\Delta }{x}_j\mathrm{\Delta }{z}_j}{u}_j^8,\\ u_{j\to i}^9& =\frac{\mathrm{\Delta }{y}_i\mathrm{\Delta }{z}_i}{\mathrm{\Delta }{y}_j\mathrm{\Delta }{z}_j}{u}_j^9,\end{array}$in which $\mathrm{\Delta }{x}_{ij}=x_i-x_j$, $\mathrm{\Delta }{y}_{ij}=y_i-y_j$ and $\mathrm{\Delta }{z}_{ij}=z_i-z_j$ computed from the two centroids of ${\mathrm{\Omega }}_i$ and ${\mathrm{\Omega }}_j$.

In high-order FV schemes, we opt to reconstruct the conservative variables instead of using primitive variable, as done in conventional second-order FV schemes. This choice is driven by the fact that reconstructing the primitive variables can only achieve second-order accuracy, even when employing a high-order scheme. Moreover, the uncertainty in the characteristic direction across the entire domain renders reconstruction in characteristic space unsuitable.

We can evaluate the coefficients of polynomials of solution by a weighted variational reconstruction. The cost functions of the approach are flexible, and the cost function in this article is defined by (10) $\begin{array}{rl}{\text{I}}_{LR}& =\frac{1}{d_{LR}}{\int }_{{\mathrm{\Gamma }}_{LR}}\\ & \times \sum _{p=0}^2{\left(w_p\frac{{\mathrm{\partial }}^p{u}_L}{\mathrm{\partial }{\tilde{x}}^p}{\left(d_{LR}\right)}^p-w_p\frac{{\mathrm{\partial }}^p{u}_R}{\mathrm{\partial }{\tilde{x}}^p}{\left(d_{LR}\right)}^p\right)}^{2}d\mathrm{\Gamma },\end{array}$(10) interface ${\mathrm{\Gamma }}_{LR}$ is the shared face between the cells on both sides, $d_{LR}$ donates the distance between the centroids and $\tilde{x}$ is the spatial derivatives along the interface normal vector. In the presented work, the weights $w_p$ are introduced to control the relative importance of each term in the weighted IJI function. Based on the findings in Section 5.4, we have selected optimized weights $\left[w_0,w_1,w_2\right]=\left[1.0,0.44,0.18\right]$ to achieve a more favourable balance between accuracy and efficiency in this study.

The cost function is to minimize the sum of the weighted IJI of all of the interfaces and boundaries, i.e. and we can obtain the following linear equations (11) $\begin{array}{r}{\mathit{A}}_{\text{i}}{\mathit{u}}_i-\sum _{j\in S_i,j\ne i}{\mathit{B}}_i^j{\mathit{u}}_j={\mathit{b}}_i,i=1,2,\dots ,N,\end{array}$(11) where the specific forms of the matrices are (12) $\begin{array}{rl}{\mathit{A}}_{\text{i}}& ={\left[\sum _{j\in S_i,j\ne i}{\int }_{{\mathrm{\Gamma }}_{\text{ij}}}\sum _{p=0}^2{w}_p^2\frac{{\mathrm{\partial }}^p{\phi }_{l,i}}{\mathrm{\partial }{\tilde{x}}^p}\frac{{\mathrm{\partial }}^p{\phi }_{m,i}}{\mathrm{\partial }{\tilde{x}}^p}{\left(d_{ij}\right)}^{2p-1}d\mathrm{\Gamma }\right]}_{9\times 9},\\ {\mathit{B}}_i^j& ={\left[{\int }_{{\mathrm{\Gamma }}_{\text{ij}}}\sum _{p=0}^2{w}_p^2\frac{{\mathrm{\partial }}^p{\phi }_{l,i}}{\mathrm{\partial }{\tilde{x}}^p}\frac{{\mathrm{\partial }}^p{\phi }_{n,j}}{\mathrm{\partial }{\tilde{x}}^p}{\left(d_{ij}\right)}^{2p-1}d\mathrm{\Gamma }\right]}_{9\times 9},\\ {\mathit{b}}_i& ={\left[\sum _{j\in S_i,j\ne i}\frac{1}{d_{ij}}{\int }_{{\mathrm{\Gamma }}_{\text{ij}}}{w}_0^2{\phi }_{l,i}\left({\overline{u}}_j-{\overline{u}}_i\right)d\mathrm{\Gamma }\right]}_9.\end{array}$(12)

We can calculate the matrices (12) at the beginning of the computation, once we obtain the mesh information, and these matrices remain constant if the meshes are fixed. Consequently, the coefficients of the polynomial are obtained as the solution to the system of linear equations (13) $\begin{array}{rl}\left(\mathit{D}-\mathit{L}-\mathit{U}\right)\mathit{u}& =\mathit{b},\\ \mathit{D}& =\left\{{\mathit{A}}_{\text{i}}\right\},\mathit{L}=\left\{{\mathit{B}}_{\text{i}}^{\text{j}},j<i\right\},\mathit{U}=\left\{{\mathit{B}}_{\text{i}}^{\text{j}},j>i\right\},\\ \mathit{u}& =\left\{{\mathit{u}}_{\text{i}}\right\},\mathit{b}=\left\{{\mathit{b}}_i\right\}.\end{array}$(13)

Obviously, because the weighted IJI function only includes both sides of the interface, the stencil $S_i$ for the reconstruction is compact.

For the boundary cell ${\mathrm{\Omega }}_i$ in Figure 2, this type of situation often poses significant challenges when coding for the reconstruction of classical high-order schemes. In contrast, this reconstruction procedure allows for easier handling of boundary treatment without sacrificing accuracy, unlike other high-order schemes. The sum of cost function of $N_{LR}$ internal faces and $N_{BC}$ boundary faces of the computational domain can be described as $\text{I}=\sum _{LR=1}^{N_{LR}}{\text{I}}_{LR}+\sum _{BC=1}^{N_{BC}}{\text{I}}_{BC}$and the boundary cost function need to treated specifically.

Figure 2. A curved triangle boundary face (left) and vertex normal vectors around sharp geometry (right).

For the slip or symmetry boundary conditions where spatial derivatives are not available, the cost function of a solution $q$ is used as the form ${\text{I}}_{BC}=\frac{1}{d_{BC}}{\int }_{{\mathrm{\Gamma }}_{BC}}{w}_0^2(q_L-q_{BC}{)}^{2}d\mathrm{\Gamma }.$The solutions of conservative variables on boundary faces are computed according to $\begin{array}{rl}{\rho }_{BC}& ={\rho }_L,\rho E_{BC}=\rho E_L,\\ [\rho u,\rho v,\rho w{]}_{BC}\cdot \mathit{n}& =-[\rho u,\rho v,\rho w{]}_L\cdot \mathit{n},\\ [\rho u,\rho v,\rho w{]}_{BC}\cdot \mathit{t}& =-[\rho u,\rho v,\rho w{]}_L\cdot \mathit{t}\end{array}$where $\mathit{n}$ is the unit outward normal vector of face ${\mathrm{\Gamma }}_{BC}$, $\mathit{t}$ is the unit tangential vector of face ${\mathrm{\Gamma }}_{BC}$.

For inter-processor or periodic boundary conditions, we need to update the $patchNeighbourField$ after the reconstruction iteration step described below and the limiting procedure. In the context of OpenFOAM, for example, the second-order coefficients of the solution in Equation (8) are stored in the container $volSymmTensorField$, which allows for efficient communication of information across inter-processor boundaries through the function $correctBoundaryConditions$. This function is responsible for updating the values of ghost cells, ensuring accurate parallel computing treatment during the simulation. Other information about the high-order solution can also be treated in a similar manner. The iteration step is similar to internal faces, with the neighbouring processor cell number $j$ set as$j<i$, where $i$ is the number of the current internal cell.

The constitutive relations of the variational reconstruction are derived by minimizing the total $I$ with respect to the coefficients of the reconstruction polynomial. This leads to a system of linear Equations (13), which is then solved using the symmetric Gauss-Seidel (SGS) method as (14) $\begin{array}{rl}{\mathit{u}}_i^{\left(s+1\right)}& =({\mathit{A}}_i{)}^{-1}\left[\sum _{j\in L\left(i\right)}{\mathit{B}}_i^j{\mathit{u}}_j^{\left(s+1\right)}+\sum _{j\in U\left(i\right)}{\mathit{B}}_i^j{\mathit{u}}_j^{\left(s\right)}-{\mathit{b}}_i\right],\\ i& =1,2,\dots ,N,\end{array}$(14) where $s$ denotes the current number of iteration step. $L\left(i\right)=\{j|j\in S_i,j<i\}$ and $U\left(i\right)=\left\{j|j\in S_i,j〉i\right\}$ represent the lower and upper cells of compact stencil $S_i$. Solving the large system of linear Equations (13) in each time step is computationally very expensive. However, the coupled iterative procedure of reconstruction and implicit time integration can significantly improve efficiency, and the SGS iterative method is utilized once in every inner iteration of implicit time integration.

To achieve third-order accuracy approximation, a sufficient number of Gauss integration points are required to evaluate values related to quadratic polynomials in Equation (8). In the context of OpenFOAM, the integration of the numerical flux for each interface in Equation (4) must meet the second-order approximating demand, necessitating three Gauss integration points for triangle faces and four for quadrangle faces. When computing the constant terms for quadratic Taylor polynomials, tetrahedron cells require four Gauss integration points, prism cells require six, and hexahedron cells require eight. Additionally, the matrices of the linear system in Equation (12) have a maximum polynomial degree of $p=4$, which mandates the use of the quadrature rule for polynomials with a degree of $p>=4$. The coordinates of integration points can be determined by converting roots of Legendre polynomials from the reference framework to the global framework, based on the vertices of faces and volumes. Similarly, the weights in the reference framework must also be converted into the global framework. The computation of integration points and weights can be performed at the beginning and saved for future use.

2.3. Time integration

The variational reconstruction Equation (13) needs to be solved at each step for explicit time integration, and this can result in high computational costs for three-dimensional problems. However, in implicit time integration, there is no need to prioritize accuracy before achieving full convergence of the solution, allowing for synchronous convergence. By coupling implicit reconstruction with time stepping, the high-order scheme can achieve computational efficiency. Specifically, the reconstruction iteration step is employed once in every single inner step of implicit time integration for unsteady simulations or in a single time step for steady-state simulations. This approach efficiently balances accuracy and computational costs, making it well-suited for practical simulations.

For steady-state simulation, the GMRES with a matrix-free LU-SGS preconditioner method (Luo et al., 1998) is used to solve the equation as (15) $\begin{array}{r}\left(\frac{\mathrm{\Omega }}{\mathrm{\Delta }t}\mathit{I}-\frac{\mathrm{\partial }\mathit{R}}{\mathrm{\partial }\mathit{U}}\right)\mathrm{\Delta }\mathit{U}=\mathit{R},\end{array}$(15) which is linearized from Equation (5). Instead of computing and storing the full matrix, the matrix-vector products can be approximated as $\frac{\mathrm{\partial }\mathit{R}}{\mathrm{\partial }\mathit{U}}\mathrm{\Delta }\mathit{U}\approx \frac{\mathit{R}\left(\mathit{U}+h\mathrm{\Delta }\mathit{U}\right)-\mathit{R}\left(\mathit{U}\right)}{h}$in which the evaluation of $\mathit{R}\left(\mathit{U}+h\mathrm{\Delta }\mathit{U}\right)$ is performed using the first-order upwind FV method which is efficient and friendly for parallel computing. $h$ is a small scalar (Nielsen et al., 1995) as $h=\frac{\sqrt{ϵ}}{\mathrm{\Delta }{\mathit{U}}_2},ϵ=1\times {10}^{-16}.$Hence, for $m$ inner iterations in the GMRES method, $m+1$ evaluations of the right-hand-side R are required. In this article, a relative tolerance of $0.1$ is set, and the number of search directions is $8$, with a maximum inner iteration limit of $40$.

To achieve high-order accuracy, the two-stage, third-order accurate singly-diagonally implicit Runge–Kutta (SDIRK) method (Ferracina & Spijker, 2008) is employed. This method provides third-order numerical accuracy, which matches the order of accuracy of the spatial discretization used in this article. To solve Equation (5), the time integration can be expressed as (16) $\begin{array}{r}\overline{\mathrm{\Omega }}\left({\overline{\mathit{U}}}^{n+1}-{\overline{\mathit{U}}}^n\right)=\mathrm{\Delta }t\left[{\text{β}}_1\mathit{R}\left({\mathit{U}}^1\right)+{\text{β}}_2\mathit{R}\left({\mathit{U}}^2\right)\right],\end{array}$(16) where (17) $\begin{array}{r}\overline{\mathrm{\Omega }}\left({\overline{\mathit{U}}}^{\alpha }-{\overline{\mathit{U}}}^n\right)=\mathrm{\Delta }t\sum _{\text{β}=1}^{\alpha }{\mathit{\alpha }}_{\alpha \text{β}}\mathit{R}\left({\mathit{U}}^{\text{β}}\right).\end{array}$(17)

The specific values of ${\text{β}}_{\alpha }$ and ${\mathit{\alpha }}_{\alpha \text{β}}$ are given according to (KvarnO et al., n.d.). A pseudo-time variable $\tau$ is introduced in the dual stepping inner iteration of the SDIRK method, and the inner iteration step is (18) $\begin{array}{r}\left(\frac{\overline{\mathrm{\Omega }}}{\mathrm{\Delta }\tau }+\frac{\overline{\mathrm{\Omega }}}{\mathrm{\Delta }t}-{\mathit{\alpha }}_{\alpha \alpha }\frac{\mathrm{\partial }\mathit{R}}{\mathrm{\partial }\overline{\mathit{U}}}\right)\mathrm{\Delta }\overline{\mathit{U}}={\mathit{R}}_t\end{array}$(18) in which (19) $\begin{array}{rl}\mathrm{\Delta }\overline{\mathit{U}}& ={\overline{\mathit{U}}}^{\left(s+1\right)}-{\overline{\mathit{U}}}^{\left(s\right)},\\ {\mathit{R}}_t& =\frac{\overline{\mathrm{\Omega }}}{\mathrm{\Delta }t}\left({\overline{\mathit{U}}}^n-{\overline{\mathit{U}}}^{\left(s\right)}\right)+{\alpha }_{\alpha \alpha }\mathit{R}\left({\mathit{U}}^{\left(s\right)}\right)+\sum _{\text{β}=1}^{\alpha -1}{\alpha }_{\alpha \text{β}}\mathit{R}\left({\mathit{U}}^{\text{β}}\right).\end{array}$(19)

The reconstruction iteration step (14) is solved once before an inner iteration step using the matrix-free LU-SGS approach (Yoon & Jameson, 1988) for Equation (18). The implicit scheme allows for the use of larger time steps compared to the maximum allowable time step imposed by an explicit stability requirement. However, using excessively large time steps can result in numerical errors (Luo et al., 2001), particularly around discontinuities. Conversely, using smaller time steps can lead to a reduced number of inner iteration steps and higher accuracy. In this article, we carefully select time steps such that the maximum Courant–Friedrichs-Lewy (CFL) numbers are in the range of $2$ to $5$. By doing so, we observe a significant reduction in the residual of the inner iteration, typically by three orders of magnitude within five steps. This approach strikes a well-balanced compromise between accuracy and computational efficiency.

3. Curved wall boundary reconstruction

High-order FV methods can achieve higher accuracy by increasing the order of reconstruction polynomials. However, the flat mesh boundaries of curved geometry may limit the accuracy of the solution, especially for solid wall boundaries. While some mesh generator software can provide high-order mesh information, the framework of OpenFOAM currently only supports information about linear meshes. In this article, we simplify the strategy presented in (Lübon et al., 2006), which was originally designed for DG methods with tetrahedral elements, and extend it to unstructured meshes with triangle and quadrangle faces.

The key strategy involves curved wall boundary reconstruction, and we obtain the information of the curved wall boundary by solving linear equations. The difference between flat and curved boundaries is given by $h=\overrightarrow{P_{flat}{P}_{curved}}\cdot {\mathit{n}}_f$. Notably, $h=0$ for each vertex of the current boundary face, and the normals at these vertices are computed by weighted averaging the normal vectors of neighbouring boundary faces, as shown in Figure 2. To ensure numerical stability, the weighted average of normal vectors should exclude the normal vector that satisfies ${\mathit{n}}_s\cdot {\mathit{n}}_f\le 0$ where ${\mathit{n}}_s$ and ${\mathit{n}}_f$ are the normal vectors of a neighbouring face and the current face, respectively.

Spatial transformation from physical space to the $\xi -\eta -h$ system involves the equation of the curved boundary of the current face, which is represented as $F\left(\xi ,\eta ,h\right)=0$. For each vertex $p$ of the current face, the equation satisfies: $\begin{array}{rl}F(\xi ,\eta ,0{)}_p& =0,\\ [F_x,F_y,F_z{]}_p& =\frac{{\mathit{n}}_p}{{\mathit{n}}_p\cdot {\mathit{n}}_f}.\end{array}$

According to the number of vertices, we can compute the unknown coefficients of the curved face equation. The equation for a triangle face is represented as Equation (20), and the equation for a quadrangle face is given by Equation (21). (20) $\begin{array}{rl}F\left(\xi ,\eta ,h\right)& =h-a_0-a_1\xi -a_2\eta -a_3{\xi }^2-a_4\xi \eta \\ & -a_5{\eta }^2-a_6{\xi }^3-a_7{\xi }^2\eta -a_8{\xi }^3\end{array}$(20) (21) $\begin{array}{rl}F\left(\xi ,\eta ,h\right)& =h-a_0-a_1\xi -a_2\eta -a_3{\xi }^2-a_4\xi \eta \\ & -a_5{\eta }^2-a_6{\xi }^3-a_7{\xi }^2\eta -a_8\xi {\eta }^2\\ & -a_9{\eta }^3-a_{10}{\xi }^4-a_{11}{\eta }^4\end{array}$(21)

The details of the curved wall boundary reconstruction are shown in Algorithm 1. After the reconstruction procedure, the new point $P_{curved}$ and normal vectors $\mathit{n}=\left[n_x,n_y,n_z\right]$ of the curved faces in physical space can be computed by $\begin{array}{rl}& \overrightarrow{P_{flat}{P}_{curved}}=h\cdot {\mathit{n}}_f,\\ & \{\begin{array}{rl}& n_x{x}_{\xi }+n_y{y}_{\xi }+n_z{z}_{\xi }=F_{\xi }\\ & n_x{x}_{\eta }+n_y{y}_{\eta }+n_z{z}_{\eta }=F_{\eta }\\ & n_x{n}_{fx}+n_y{n}_{fy}+n_z{n}_{fz}=1\end{array}.\end{array}$

Table

Algorithm 1. Curved wall boundary reconstruction.
For $f=1,2,\dots ,nBoundaryFaces$ do
For $p=1,2,\dots ,nFaceVertices$ do
For $s=1,2,\dots ,nVertexFaces$ do
If $\mathit{S}{\mathit{f}}_s\cdot \mathit{S}{\mathit{f}}_f>0$ then
$\mathit{S}{\mathit{f}}_p+=\mathit{S}{\mathit{f}}_s,\left|\mathit{S}{\mathit{f}}_p\right|+=\left|\mathit{S}{\mathit{f}}_s\right|$
End if
${\mathit{n}}_p=\mathit{S}{\mathit{f}}_p/\left|\mathit{S}{\mathit{f}}_p\right|$
End for
Add $F\left({\xi }_p,{\eta }_p,0\right)=0,[F_x,F_y,F_z{]}_p=\frac{{\mathit{n}}_p}{{\mathit{n}}_p\cdot {\mathit{n}}_f}$ to linear equations
End for
Compute coefficients of $F\left(\xi ,\eta ,h\right)=0$ by solving linear equations
End for

4. Limiting procedure

4.1. The problem-independent shock detector

Applying the limiter function to all the cells of the computational domain leads to high computational cost and low accuracy. Therefore, we detect troubled cells that need to be limited using the indicator $I{S}_i$ defined by (22) $\begin{array}{r}I{S}_i=\frac{\sum _{j\in S_i,j\ne i}|u_i\left({\mathit{x}}_i\right)-u_{j\to i}\left({\mathit{x}}_i\right)|}{N_j{h}_i^{\left(k+1\right)/2}\text{max}\left\{{\overline{u}}_i,{\overline{u}}_j\right\}},j\in S_i,\end{array}$(22) in which $N_j=4$ for tetrahedral cells and $N_j=6$ for hexahedral cells. ${\mathit{x}}_i$ is the centroid of ${\mathrm{\Omega }}_i$. The solution of the degree $p$ polynomial satisfies (23) $\begin{array}{r}{u}_i\left({\mathit{x}}_i\right)-u_{j\to i}\left({\mathit{x}}_i\right)=\{\begin{array}{ll}O\left(h_i^{p+1}\right),& \text{in smooth region}\\ O\left(1\right),& \text{near discontinuity}.\end{array}\end{array}$(23)

Therefore, $I{S}_i\to 0$ as either $h_i\to 0$ or $p\to 0$ in smooth regions, whereas $I{S}_i\to \mathrm{\infty }$ near discontinuities, so the shock detector is also an error-based smoothness indicator. (24) $\begin{array}{r}(I{S}_i<{\overline{S}}_{dis})=\{\begin{array}{ll}true,& \text{smooth region}\\ false,& \text{shock region}.\end{array}\end{array}$(24) In this article, we use ${\overline{S}}_{dis}=1$ for one-dimensional problems and ${\overline{S}}_{dis}=2$ for multi-dimensional problems. Furthermore, we employ the solutions of density as the variables to be computed with the shock detector. According to the zero-mean basis, we can compute the left side of Equation (16) using (25) $\begin{array}{rl}& u_i\left({\mathit{x}}_i\right)-u_{j\to i}\left({\mathit{x}}_i\right)\\ & =\left(u_i^4-u_{j\to i}^4\right){\overline{\phi }}_{1,1}\mathrm{\setminus }2+\left(u_i^5-u_{j\to i}^5\right){\overline{\phi }}_{2,2}\mathrm{\setminus }2\\ & +\left(u_i^6-u_{j\to i}^6\right){\overline{\phi }}_{3,3}\mathrm{\setminus }2+\left(u_i^7-u_{j\to i}^7\right){\overline{\phi }}_{1,2}\\ & +\left(u_i^8-u_{j\to i}^8\right){\overline{\phi }}_{1,3}+\left(u_i^9-u_{j\to i}^9\right){\overline{\phi }}_{2,3},\end{array}$(25) and there is no need to evaluate the special values of ${\phi }_{l,i}\left({\mathit{x}}_i\right)$ with the quadrature-free approach, making the method computationally efficient.

4.2. The simple WBAP limiter

In this article, a simple WBAP limiter is presented to achieve high computational efficiency, derived from the original WBAP limiter (W. Li & Ren, 2012b). The limiter's fundamental function is described as (26) $\begin{array}{r}L\left(a_0,a_1,\dots ,a_{N_j}\right)=a_0\cdot W\left(1,\frac{a_1}{a_0},\dots ,\frac{a_{N_j}}{a_0}\right),\end{array}$(26) where $a_0,a_1,\dots ,a_{N_j}$ represent the candidates, with $a_0$ being the coefficient that needs to be limited. The candidates consist of the PR and SRs coefficients of the current troubled cell. The function $W$ is defined by (27) $\begin{array}{rl}& W\left(1,{\theta }_1,\dots ,{\theta }_{N_j}\right)\\ & =\{\begin{array}{cc}{\displaystyle \frac{n_p+\sum _{m=1}^{N_j}1/{\theta }_m^3}{n_p+\sum _{m=1}^{N_j}1/{\theta }_m^4},}& \text{if }{\theta }_1,\dots ,{\theta }_{N_j}>0\\ 0,& \text{else}\text{.}\end{array}\end{array}$(27)

The limited reconstruction polynomials for Equation (1) are computed in the form of (28) $\begin{array}{r}{\tilde{U}}_i\left(\mathit{x}\right)={\overline{U}}_i+\sum _{l=1}^{N\left(p\right)}{\tilde{U}}_i^l{\phi }_{l,i}\left(\mathit{x}\right),\end{array}$(28) where the coefficients of conservative variables ${\tilde{U}}_i^l=\left[u_{i,\rho }^l,u_{i,\rho u}^l,u_{i,\rho v}^l,u_{i,\rho w}^l,u_{i,\rho E}^l\right]$ are computed according to (29) $\begin{array}{r}{\tilde{U}}_i^l=L\left({\hat{U}}_{ij1}^l,,{\hat{U}}_{ij1}^l,\dots ,{\hat{U}}_{ij{N}_j}^l\right).\end{array}$(29)

To illustrate the limiting procedure, troubled cells have been detected, as shown in Figure 3.

Figure 3. Troubled cells marked by red in limiting procedure.

For the original WBAP, Figure 4 and Algorithm 2 shows that the successive limiting procedure is performed with a characteristic-wise projection procedure for each troubled cell. It involves transforming the second-order scaled derivatives to first-order scaled derivatives, which are computed using unlimited first-order scaled derivatives and limited second-order scaled derivatives. However, this characteristic-wise projection on PR and SRs coefficients is computationally very expensive. In the limiting procedure, the projection function corresponding to $V_i^l=R^{-1}\left({\overline{U}}_i,{\mathit{n}}_{ij}\right)\cdot U_i^l$ is required for $p\times N_{trouble}\times \left[N_j\times \left(N_j+1\right)\right]$ times for $N_{trouble}$ troubled cells. Additionally, the evaluations of $R\left({\overline{U}}_i,{\mathit{n}}_{ij}\right)$ and $R^{-1}\left({\overline{U}}_i,{\mathit{n}}_{ij}\right)$ are required $p\times N_{trouble}\times N_j$ times for $N_{trouble}$ troubled cells. This computational cost can become substantial as the number of troubled cells increases.

Figure 4. Original WBAP successive limiting procedure.

Table

Algorithm 2. Original WBAP successive limiting procedure.
For $k=p,p-1,\dots ,1$ do
For $i=1,2,\dots ,\text{nCells}$ do
If $I{S}_i>{\overline{S}}_{dis}$ then
For $j=j_1,j_2,\dots ,j_{N_j}$ do
${\hat{U}}_{i,j}^k=R\left({\overline{U}}_i,{\mathit{n}}_{ij}\right)\cdot L\left(R^{-1}\left({\overline{U}}_i,{\mathit{n}}_{ij}\right)\cdot \left(U_i^k,U_{j_1\to i}^k,U_{j_2\to i}^k,\dots ,U_{j_{N_j}\to i}^k\right)\right)$
End for
${\tilde{U}}_i^k=L\left({\hat{U}}_{i,j_1}^k,{\hat{U}}_{i,j_2}^k,\dots ,{\hat{U}}_{i,j_{N_j}}^k\right)$
End if
End for
For $i=1,2,\dots ,\text{nCells}$ do
If $I{S}_i>{\overline{S}}_{dis}$ then
$U_i^k={\tilde{U}}_i^k$
End if
End for
End for

Obviously, the successive limiting procedure leads to higher computational cost as the degree $p$ increases. The key to designing the new simple limiter is to achieve similar robustness to the original WBAP limiter while eliminating the whole-domain succussive limiting procedure. To accomplish this, the successive limiting procedure for a direction of characteristic-wise projection can be computed within an edge-based mini stencil. The procedure is performed between the left and right cells of the shared interface, instead of involving the whole computational domain. The details of the new simple WBAP limiting procedure are shown in Figure 5 and Algorithm 3.

Figure 5. Simple WBAP successive limiting procedure.

Table

Algorithm 3. Simple WBAP successive limiting procedure.
For $f=1,2,\dots ,nFaces$ do
If $\text{max}\left(I{S}_i,I{S}_j\right)>{\overline{S}}_{dis}$ then
$i$ is the left cell of $f$, and $j$ is the right cell of $f$
$\left(V_i,V_j\right)=R^{-1}\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)\cdot \left(U_i,U_j\right)$
For $k=p,p-1,\dots ,1$ do
${\hat{V}}_i^k=L\left(V_i^k,V_{j\to i}^k\right),{\hat{V}}_j^k=L\left(V_j^k,V_{i\to j}^k\right)$
$V_i^k={\hat{V}}_i^k,V_j^k={\hat{V}}_j^k$
End for
${\hat{U}}_{i,j}=R\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)\cdot V_i,{\hat{U}}_{j,i}=R\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)\cdot V_j$
End if
End for
For $i=1,2,\dots ,\text{nCells}$ do
If $I{S}_i>{\overline{S}}_{dis}$ then
${\tilde{U}}_i=L\left({\hat{U}}_{i,j1},{\hat{U}}_{i,j2},\dots ,{\hat{U}}_{i,j{N}_j}\right)$
End if
End for

This procedure only limits $V_i^l$ and $V_j^l$ by the SR of each other, and both of the two cells share the same $R^{-1}\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)$, where ${\overline{U}}_{ij}=\left({\overline{U}}_i+{\overline{U}}_j\right)/2$, instead of limiting by all the direct face neighbours of ${\mathrm{\Omega }}_i$. Additionally, the successive limiting procedure is only used in computing $L\left(V_i^l,V_{j\to i}^l\right)$ to avoid the complicated limiting steps in the original WBAP limiter. The need of projection function in the new simple WBAP limiting procedure is $N_{trouble}\times N_j$, which is approximately $1/\left(p\times (N_j+1\right))$ of the original computation requirements. The need for the evaluations of $R\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)$ and $R^{-1}\left({\overline{U}}_{ij},{\mathit{n}}_{ij}\right)$ is required for $N_{trouble}\times N_j/2$ times. Additionally, the requirement for parallel communication is $1/p$ of the original procedure. The parameter $n_p=2$ is used to emphasize the PR of ${\mathrm{\Omega }}_i$ for the edge-based mini stencil, and $n_p=1$ is taken for Equation (29) as we treat every direction equally. The numerical result in the following section shows similar robustness performance to the original WBAP limiter.

4.3. The positivity-preserving limiter

The computational cost is high to detect negative pressure or density for each cell on each Gauss point of each interface and volume. To efficiently detect cells with negative solutions of density and pressure while maintaining high-order accuracy, a new positivity-preserving limiter is proposed to perform a positivity correction procedure on the troubled cells indicated by Equation (22). Using the constant attribute of basis function polynomial, we compute (30) $\begin{array}{r}\overline{u}-\sum _{l=1}^9|u^l|{\tilde{\phi }}_l<0\end{array}$(30) in which $\begin{array}{rl}{\tilde{\phi }}_4& =\frac{1}{2}\left({\tilde{\phi }}_1^2-{\overline{\phi }}_{1,1}\right),{\tilde{\phi }}_{5,i}=\frac{1}{2}\left({\tilde{\phi }}_2^2-{\overline{\phi }}_{2,2}\right),\\ {\tilde{\phi }}_6& =\frac{1}{2}\left({\tilde{\phi }}_5^2-{\overline{\phi }}_{3,3}\right),\\ {\tilde{\phi }}_7& ={\tilde{\phi }}_1{\tilde{\phi }}_2+|{\overline{\phi }}_{1,2}|,{\tilde{\phi }}_8={\tilde{\phi }}_1{\tilde{\phi }}_3+|{\overline{\phi }}_{1,3}|,\\ {\tilde{\phi }}_9& ={\tilde{\phi }}_2{\tilde{\phi }}_3+|\overline{{\overline{\phi }}_{2,3}}|,\\ {\tilde{\phi }}_1& =max\left(x_{max}-x_c,x_c-x_{min}\right)/\mathrm{\Delta }x,\\ {\tilde{\phi }}_2& =max\left(y_{max}-y_c,y_c-y_{min}\right)/\mathrm{\Delta }y,\\ {\tilde{\phi }}_2& =max\left(y_{max}-y_c,y_c-y_{min}\right)/\mathrm{\Delta }y.\end{array}$

We apply the positivity-preserving limiter on the solutions of density and total energy, and compute the smaller of the two factors according to $\theta =\{\begin{array}{cc}\overline{u}/\mathrm{\Delta }u,& \overline{u}-\mathrm{\Delta }u<0\\ 1,& \text{else}\text{.}\end{array}$

The coefficients of all the conservative variables need to be reduced by the factor above. This positivity-preserving limiter can maintain positive solutions for pressure and density without significantly increasing computational cost.

5. Numerical results and discussion

In this section, we present the results of various problems to demonstrate the performance of the methods proposed in this article. We compare our solver with two typical density-based solvers implemented in OpenFOAM. The first one is $blastFoam$ (Synthetik Applied Technologies, LLC, 2020), which is a classical second-order unstructured solver. The second improved second-order solver is $ROUNDA$ (Deng, 2023), which is based on the newly-developed low dissipative, structure-preserving ROUND (Reconstruction Operators on Unified Normalized-variable Diagram) scheme in a three-cell compact stencil and has shown significant improvements in resolution compared to conventional second-order schemes (Cheng et al., 2023; Deng et al., 2023). The open-source code for both solvers is available on the website and can be used within a similar framework as the $VR3$ solver presented in this study. The schemes employed by these three solvers are summarized in Table 1. The first two solvers use the explicit TVD Runge–Kutta method RK3SSP (Gottlieb et al., 2001) for time integration, while the $VR3$ solver utilizes the SDIRK3 method as discussed previously. To ensure accuracy and mitigate the impact of the time step on the results, we set $CFL=0.4$ for the explicit time schemes and $CFL=2$ for the implicit time scheme of $VR3$ in unsteady flow simulations.

Table 1. Information of solvers for numerical tests.

Solver	$blastFoam$	$ROUNDA$	$VR3$
Time scheme	RK3SSP	RK3SSP	SDIRK3 GMRES + LUSGS
Grad scheme	Gauss linear	Gauss linear	VR quadratic
Interpolation scheme	Linear vanAlbada	Non-linear low dissipative ROUND	Quadratic polynomial with WBAP limiters
Mesh	Unstructured	Uniform structured	Unstructured

5.1. Accuracy test for the scalar transport equation

In the accuracy test, we consider an initial smooth distribution, as given by $u\left(x\right)=\text{sin}\left(\pi x\right),x\in \left[-1,1\right],$and the computation is carried out for one period (at $t=2.0$). The results of the accuracy test are presented in Figure 6, clearly demonstrating that the $VR3$ scheme outperforms the other two solvers as well as the typical third-order upwind scheme. The spectral properties shown in Figure 7 indicate that the third-order scheme without limiters can achieve better resolution in relatively smooth flow region. Especially, the spectral property of ROUNDA is close to the third-order upwind scheme.

Figure 6. The accuracy test for the schemes of the blastFoam, ROUNDA and VR3.

Figure 7. The spectral properties of the schemes of the blastFoam, ROUNDA and VR3.

5.2. Subsonic flow past a sphere

This test case involves a subsonic flow past a sphere at a Mach number of $M{a}_{\mathrm{\infty }}=0.38$, and it is designed to evaluate the order of error convergence rate on hexahedral meshes. Four meshes with dimensions $8\times 8\times 24$, $16\times 16\times 48$, $24\times 24\times 72$ and $32\times 32\times 96$ (as shown in Figure 8) are used for comparison. The sphere has a radius of $r_{\text{sphere}}=0.5$, and the computational domain is bounded by $r_{\text{farField}}=20$. The focus of the comparison is on the $32\times 32\times 96$ mesh. As illustrated by the Mach contours in Figure 9, the $VR3$ demonstrates superior performance compared to $blastFoam$ and $ROUNDA$. Notably, the curved $VR3$ variant exhibits the best results. One of the reasons for the relatively poor performance of $ROUNDA$ is its limitation to uniform meshes for the property of second-order non-linear schemes. Consequently, to further evaluate the solvers on unstructured meshes, the following tests are conducted using $blastFoam$ and $VR3$.

Figure 8. Four successively refined hexahedral meshes for subsonic flow past a sphere.

Figure 9. Subsonic flow past a sphere. Computed 31 equally spaced Mach number contours from 0.019 to 0.589 by the three solvers on 32 × 32 × 96 mesh.

To assess the accuracy of the schemes, we present the $L_1$ and $L_2$ entropy errors, along with their rates of convergence, and these comparisons are graphically depicted in Figure 10. The findings indicate that the methods employed in $VR3$ yield more accurate solutions compared to the classical second-order FV method in $blastFoam$. Furthermore, the incorporation of the curved wall boundary reconstruction approach significantly enhances the achievement of the desired design accuracy for the third-order schemes. The entropy error defined by ${ϵ}_{\text{entropy}}=\frac{p/p_{\mathrm{\infty }}}{{\left(\rho /{\rho }_{\mathrm{\infty }}\right)}^{\gamma }}-1.$

Figure 10. Order of error convergence for subsonic flow past a sphere.

Figure 11. Convergence history comparison for subsonic flow past a sphere with mesh 32 × 32 × 96.

Figure 12. Surface mesh used for transonic flows past a ONERA M6 wing.

Figure 13. Transonic flows past a ONERA M6 wing. Computed pressure coefficient contours from −1.1 to 1.4.

Figure 14. Comparison of experiment and computed surface pressure coefficient of the three methods at different semi-span locations.

Figure 15. Comparison of entropy production of the three methods at different semi-span locations.

We conductea comparison of the time integration methods using a $32\times 32\times 96$ mesh. For the GMRES-LUSGS method, the local CFL number was set to $40$, while $blastFoam$ and $ROUNDA$ employed the explicit time scheme with a maximum CFL number of $1$. The convergence history results are depicted in Figure 11 Notably, the implicit time integration proves to be highly efficient for the third-order solver in steady flow simulations, despite the computational expense involved in reconstruction and flux evaluation.

5.3. Transonic flows past a ONERA M6 wing

This test cast is transonic flows past a ONERA M6 wing which has a leading-edge sweep angle of ${30}^{\circ }$, an aspect of ${3.8}^{\circ }$ and a taper ratio of $0.562$. . The Mach number of the flow is $0.8398$ and $\text{AOA}={3.06}^{\circ }$. The mesh shown in Figure 12 consists of 341,797 tetrahedral cells, and the solutions are listed in Figure 13. Figure 14 shows the comparison of experimental and computed surface pressure coefficient at $20\mathrm{\%}$, $44\mathrm{\%}$ and $90\mathrm{\%}$ semi-span locations. The effectiveness of the simple WBAP limiter used with third-order metho is evident in its ability to suppress spurious oscillations and accurately capture shocks on the surface and the third-order methods produce sharper shock profiles compared to the second-order method in $blastFoam$, especially at the $90\mathrm{\%}$ semi-span location. Figure 15 illustrates the entropy production at three different semi-span locations, revealing that the third-order method with curved boundary exhibits lower entropy production, particularly in the smooth domain before encountering shocks.

5.4. Isentropic vortex problem (Hu & Shu, 1999)

This is a two-dimensional smooth flow problem without shocks, and we use the scheme presented in this article to demonstrate its advantages. The computational domain is $[0,10]\times [0,10]$. The mean flow condition is $\left(\rho ,u,v,p\right)=\left(1,1,1,1\right)$, and we introduce the following perturbations in the form of an isentropic vortex described as (31) $\begin{array}{rl}\left(\delta u,\delta v\right)& =\frac{\chi }{2\pi }{e}^{0.5\left(1-r^2\right)}\left(-\overline{y},\overline{x}\right),\\ \delta T& =-\frac{\left(\gamma -1\right){ϵ}^2}{8\gamma {\pi }^2}{e}^{1-r^2},\delta \left(S=p/{\rho }^{\gamma }\right)=0,\end{array}$(31) where $\left(\overline{x},\overline{y}\right)=\left(x-5,y-5\right),r^2={\overline{x}}^2+{\overline{y}}^2$, and the vortex strength $\chi =5$. All the boundary conditions are periodic. The four types of meshes used, including regular prismatic mesh, irregular prismatic mesh, regular hexahedral mesh and irregular hexahedral mesh, are shown in Figure 16. The tests are performed until $t=2.0$ to demonstrate that the third-order variational reconstruction FV method in $VR3$ achieves the intended accuracy for the four types of meshes as listed in Table 2.

To optimize the weights of weighted IJI function in variational reconstruction, we designed a numerical test to compare different weights and find the optimal ones. Figure 17 illustrates the influence of weights on average inner iteration steps on a regular hexahedral mesh at $h=1/4$, and the physical time step set to $1/10$. All the weights can achieve intended accuracy, and we set $w_0=w_1$ and $w_1=w_2$ in the left result and right result, respectively. Based on the results, we choose $\left[w_0,w_1,w_2\right]=\left[1.0,0.44,0.18\right]$ to achieve high computational efficiency for numerical tests in this article. Table 3 presents a comparison of average inner iteration steps between the original weights $\left[w_0,w_1,w_2\right]=\left[1.0,1.0,0.5\right]$ as presented in (Wang et al., 2017) and the optimized weights in the four types of meshes used in the accuracy tests above. The optimized weights achieve the intended accuracy in all tests with lower computational cost than the original weights.

Figure 16. Four types of meshes for isentropic vortex problem.

Figure 17. Influence of weights on average inner iteration steps.

Table 2. Errors and convergence rates for the isentropic vortex problem on four types of meshes.

Mesh type	Size	${\mathit{L}}_1$error	Order	${\mathit{L}}_2$error	Order	${\mathit{L}}_{\mathrm{\infty }}$error	Order
Regular prismatic mesh	$1$	$3.33\text{E}-3$		$7.36\text{E}-3$		$5.91\text{E}-2$
	$1/2$	$2.17\text{E}-4$	$3.94$	$5.24\text{E}-4$	$3.81$	$5.59\text{E}-3$	$3.40$
	1/4	2.61E-5	3.05	$6.36\text{E}-5$	$3.05$	$5.08\text{E}-4$	$3.46$
	$1/8$	2.09E-6	$3.65$	$5.65\text{E}-6$	$3.49$	$6.26\text{E}-5$	$3.02$
Irregular prismatic mesh	$1$	$1.72\text{E}-3$		$3.73\text{E}-3$		$2.49\text{E}-2$
	$1/2$	$2.17\text{E}-4$	$2.99$	$5.24\text{E}-4$	$2.83$	$5.59\text{E}-3$	$2.16$
	$1/4$	1.34E-5	$4.01$	$3.24\text{E}-5$	$4.01$	$2.56\text{E}-4$	$4.45$
	$1/8$	1.35E-6	$3.31$	$4.23\text{E}-6$	$2.94$	$6.42\text{E}-5$	$2.00$
Regular hexahedral mesh	$1$	$9.35\text{E}-3$		$2.22\text{E}-2$		$1.33\text{E}-1$
	$1/2$	$1.39\text{E}-3$	$2.75$	$3.39\text{E}-3$	$2.71$	$2.27\text{E}-2$	$2.55$
	$1/4$	1.55E-4	$3.16$	$3.76\text{E}-4$	$3.17$	$3.53\text{E}-3$	$2.68$
	$1/8$	8.11E-6	$4.26$	$1.98\text{E}-5$	$4.25$	$1.44\text{E}-4$	$4.62$
Irregular hexahedral mesh	$1$			$1.89\text{E}-2$		$1.40\text{E}-1$
	$1/2$	$1.25\text{E}-3$	$2.72$	$2.98\text{E}-3$	$2.67$	$2.32\text{E}-2$	$2.60$
	$1/4$	$1.32\text{E}-4$	$3.25$	$3.19\text{E}-4$	$3.22$	$2.70\text{E}-3$	$3.10$
	$1/8$	$7.02\text{E}-6$	$4.23$	$1.70\text{E}-5$	$4.23$	$1.52\text{E}-4$	$4.15$

Table 3. Comparison of average inner iteration steps of the original weights and the optimized weights in the four types of meshes.

	Original weights	Optimized weights
Regular prismatic mesh	16.925	9.05
Irregular prismatic mesh	20.9	13.025
Regular hexahedral mesh	24.175	16.075
Irregular hexahedr mesh	19.225	21.65

To investigate the numerical dissipation of schemes and the two types of WBAP limiters, Figure 18 shows the comparisons of density solutions along centre x axis ($y=5$) computed by different schemes on a regular hexahedral mesh at $h=1/8$ after 5 and 10 time periods, respectively. It shows that ROUND schemes can better preserve the vortex strength than conventional second-order schemes due to the improved numerical dissipation property. The solver with an unlimited third-order scheme reduces the errors of numerical dissipation, and the advantage of high-order schemes becomes more evident with increasing time of performance. We also set ${\overline{S}}_{dis}=0$ to apply the WBAP limiters on the whole computation domain, and the simple WBAP limiter ($n_p=2$) produces less numerical dissipation, compared to the original WBAP limiter ($n_p=10$).

Figure 18. Isentropic vortex evolution after 5 time periods (left) and 10 time periods (right).

The comparison of wall-clock time listed in Table 4 shows that the new simple WBAP limiter is more efficient than the original WBAP limiter. The numerical tests in this article were performed on an Ubuntu operating system with an Intel Core i7-10700 2.9Ghz CPU and 16 GB memory.

Table 4. Comparison of wall-clock time for different limiters for the smooth isentropic vortex problem on a regular hexahedral mesh at h = 1/8 after 10 time periods.

Limiters	Unlimited	Simple WBAP	Original WBAP
wall-clock time (s)	888	1963	$5189$

5.5. Shock tube problem

This is a class of one-dimensional benchmark tests, and we can use the comparisons of numerical and exact solutions to assess the performance of the solvers in flow simulation. The initial conditions are $\left(\rho ,u,p\right)\left(\mathit{x}\right)=\{\begin{array}{cc}\left({\rho }_L,u_L,p_L\right)& \text{if }x\le 0.5\\ \left({\rho }_R,u_R,p_R\right)& \text{if }x>\mathrm{0.5.}\end{array}$The cases include Sod shock tube problem (Sod, 1978) and Lax problem (Lax, 1954) whose initial conditions are show in Table 5. The mesh size is $h=1/200$.

Table 5. Two shock tube cases.

Problem	${\rho }_L$	$u_L$	$p_L$	${\rho }_R$	$u_R$	$p_R$	$t_{end}$
Sod	1.0	0.75	1.0	0.125	0.0	0.1	0.2
Lax	0.445	0.698	3.528	0.5	0.0	0.571	0.1

As noted by (Jiang & Shu, 1996), these two cases serve as challenging tests to assess the characteristic-based limiters for high-order schemes. The tests were conducted until $t_{end}$ using $blastFoam$, $ROUNDA$, and $VR3$ with two types of limiters. The results distribution in Figure 19 shows that the simple WBAP limiter can capture high-resolution flow structures similar to the scheme of $ROUNDA$ solver, without requiring the computationally expensive successive limiting procedure. On the other hand, $blastFoam$ exhibits lower resolution due to its higher numerical dissipation.

Figure 19. Density distribution of Sod case at t = 0.2 (left) and Lax case at t = 0.1 (right).

5.6. Shu–Osher problem

The test is the one-dimensional Shu and Osher (1989) problem which is a simple model for shock-turbulence interactions. The flow structures caused by the interaction of left moving shock and the right smooth sine wave is complicated and they are demanding for numerical schemes to capture the high-frequency waves. Both of ends are reflective boundary. The initial conditions are $\begin{array}{rl}& \left(\rho ,u,p\right)\left(\mathit{x}\right)\\ & =\{\begin{array}{ll}\left(3.857143,2.629369,10.33333\right)& \text{if }0\le x<1\\ \left(1+0.2\text{sin}\left(5x\right),0,1\right)& \text{if }x\ge 1.\end{array}\end{array}$The density distribution results at $t=1.8$ are depicted in Figure 20 with $400$ cells, and the ‘exact’ solution is obtained using the one-dimensional fifth-order WENO scheme (Jiang & Shu, 1996) with $5000$ cells. The comparisons clearly demonstrate the significant advantages of high-order schemes over the second-order scheme of $blastFoam$ in simulations involving high-frequency waves and excellent shock capturing capability. The $VR3$ solver with WBAP limiters exhibits smaller numerical dissipation compared to the $ROUNDA$ solver in this particular test. This achievement demonstrates the effectiveness of the simple WBAP limiter in maintaining high accuracy and resolving flow features efficiently.

Figure 20. The Shu-Osher problem. Density distribution at t = 1.8.

5.7. Two interacting blast waves problem

This test case is to simulate interactions of two blast waves problem (Lax & Liu, 1998), and uniform cells with $h=1/400$ are used. The dramatic interactions are hardly to capture high-resolution flow structures without numerical oscillation, and the solutions of density and pressure are near zero near discontinuities, so positivity-preserving limiter is necessary for the achievement of robustness. Especially, the situation become more severe for high-order schemes because of their non-monotonic polynomials. The initial conditions are $\left(\rho ,u,p\right)\left(\mathit{x}\right)=\{\begin{array}{ll}\left(1,0,1000\right)& \text{for }0\le x<0.1\\ \left(1,0,0.01\right)& \text{for }0.1\le x<0.9\\ \left(1,0,100\right)& \text{otherwise}\text{.}\end{array}$The density distribution results at $t=0.038$ are depicted in Figure 21, and the ‘exact’ solution is obtained using the one-dimensional fifth-order WENO scheme (Jiang & Shu, 1996) with $5000$ cells. Through the comparisons of results computed by the solvers introduced in this article, we observe that the smearing of contact discontinuities appears more pronounced in the second-order scheme of $blastFoam$, while the third-order scheme achieves better performance. Notably, the VR3 solvers exhibit more smearing due to higher numerical dissipation caused by the WBAP limiter compared to the ROUNDA scheme around strong shocks.

Figure 21. Two blast wave problem. Density distribution at t = 0.038.

Figure 22. Mach 3 wind tunnel with a step. Thirty equally spaced density contour lines from ρ = 0.3 to 6.15.

5.8. A Mach 3 wind tunnel with a step (Lax & Liu, 1998)

The problem involves strong shock interactions induced by a step located at the corner $\left(x,y\right)=\left(0.6,0.2\right)$, within the computational domain $\left[0,3\right]\times \left[0,1\right]$. The initial condition and inflow are specified as $\left(\rho ,u,v,p\right)=\left(1.4,3,0,1\right)$, and the wind tunnel's solid wall is treated as a slip boundary. The simulation is conducted using a uniform mesh with $h=1/320$.

Both $ROUNDA$ and $VR3$ with the simple WBAP limiter are capable of capturing the wavy contact discontinuity originating from the triple point with schemes beyond second order. It is well known that the corner of the step can lead to an erroneous entropy layer at the downstream bottom wall. In this scenario, the $VR3$ solver with the simple WBAP limiter outperformed the other solvers in handling the spurious Mach stem at the bottom wall. The distribution of troubled cells identified is shown in Figure 23.

Figure 23. Forward step problem. Troubled cells are marked black.

5.9. Double Mach reflection problem (Lax & Liu, 1998)

The problem consists of interactions caused by a right-moving shock wave and a reflecting wall, and the computational domain $\left[0,3\right]\times \left[0,1\right]$. The initial conditions are $\begin{array}{rl}& \left(\rho ,u,v,p\right)\left(\mathit{x}\right)\\ & =\{\begin{array}{ll}\left(1.4,0,0,1.0\right)& \text{if }x>1/6+\sqrt{3}/3y,\\ \left(1+0.2\text{sin}\left(5x\right),0,1\right)& \text{if }x\le 1/6+\sqrt{3}/3y.\end{array}\end{array}$In Figure 24, the density contours at $t=0.2$ computed show that the $ROUNDA$ and $VR3$ solver with the simple WBAP limiter can capture more intricate details of the complicated interactions around the right shock's intersection point. The simple WBAP limiter can introduce larger numerical dissipation than the $ROUNDA$ scheme around strong shock waves, which helps to obtain non-oscillatory results due to its characteristic. However, this level of numerical dissipation may not be sufficient for the cells around vortex structures. The interaction between the rolled-up vortex structure and the right shock differs due to their numerical dissipative properties. The troubled cells at the last time step are shown in Figure 25. The troubled cells on the upper boundary are identifiedp} by the indicator function because the information of the moving shock is present in the boundary condition.

Figure 24. Double Mach reflection. Thirty equally spaced density contour lines from ρ = 1.75 to 21.5.

Figure 25. Double Mach reflection problem. Troubled cells are marked black.

5.10. 3d explosion problem (Zanotti et al., 2015)

The 3D explosion problem describes the evolution of a strong spherical shock. We use an eighth sphere with $r=1$ as the computation domain. The initial conditions about the spherical shock are described as $\left(\rho ,u,v,w,p\right)=\{\begin{array}{ll}\left(1,0,0,0,1\right)& \text{if }r\le 0.5\\ \left(0.125,0,0,0,0.1\right)& \text{if }r>0.5,\end{array}$where $r=\sqrt{x^2+y^2+z^2}$ is the radius. The exact solutions can be computed by solver of simple one-dimensional model with geometric source terms. The exact solutions are computed using a solver for a simple one-dimensional model with geometric source terms.

To test the resolution properties on three-dimensional unstructured meshes, we solve the problem using $blastFoam$ and the $VR3$ solvers with two types of WBAP limiters. The mesh consists of uniform tetrahedral cells with a characteristic size of $1/160$., which corresponds to $9.2$ million cells. The results in Figure 26 display the density profiles and the troubled cells marked in the radial direction at $t=0.2$. For comparison, we plot the numerical solutions in Figure 27 obtained with the original WBAP limiter and the new simple WBAP limiter. The results show that the solutions of the third-order solvers exhibit high resolution compared to the conventional second-order solver.

Figure 26. Numerical results for 3D explosion problem in the radial direction at time t = 0.2. Density profile (left) and troubled cells marked (right).

Figure 27. The results distribution along the radius direction.

The parallel scaling efficiency of the introduced $VR3$ solver is presented in the Figure 28. The compute times are measured as the average time over $20$ inner iteration steps. The VR3 solver with the simple WBAP limiter demonstrates good parallel scaling efficiency when the number of inter-processor faces accounts for a small percentage. This scaling study, which uses the same mesh as the 3D explosion problem but varies the number of cores, is conducted on a computing cluster composed of core nodes with two Intel Xeon Gold 6240 processors each. The maximum number of cores for each node is $32$ in this test. For the computing task corresponding to $128$ cores, the following information is listed: (1) average number of processor faces = $72387$, with the fluctuation of cell numbers under 2%; (2) max number of processor patches = $18$ ($59.0643\mathrm{\%}$ above the average of $10.6875$); (3) max number of faces between processors = $6705$ ($29.0641\mathrm{\%}$ above the average of $5195.09$). In this typical problem, the max number of troubled cells is within $4\mathrm{\%}$, which has a minimal effect on the parallel scaling efficiency.

Figure 28. Parallel scaling efficiency of the VR3 solver with a simple WBAP limiter.

6. Conclusion

In this article, we present a compact third-order density-based solver developed within the framework of OpenFOAM. The solver utilizes implicit variational reconstruction on three-dimensional unstructured meshes and employs technology coupling reconstruction and time integration for achieving high-order accuracy. Our numerical results demonstrate that the third-order solver with curved wall boundary reconstruction successfully attains the desired accuracy in smoother flow simulations without requiring high-order mesh information. Furthermore, the new simple WBAP limiter effectively suppresses spurious oscillations without necessitating a whole-domain successive limiting procedure, while still maintaining similar robustness to the original WBAP limiter in flow simulations with shocks.

Comparing the third-order solver to two typical second-order solvers reveals the clear advantage of high-order schemes, as our solver provides higher resolution results. Additionally, the third-order solver exhibits outstanding properties that make it a superior tool for problems involving complex flow structures. The compactness of reconstruction and limiting procedures also makes the solver highly amenable to parallelization, enhancing its computational efficiency. It is worth noting that all the tests performed in this study are inviscid flow problems, and future work will focus on extending the solver to tackle turbulent compressible flows.

Acknowledgements

The authors thank Weicheng Pei for his careful instruction of writing and coding. Author Contributions: Conceptualization, M.Y. and S.L.; methodology, M.Y.; software, M.Y.; validation, M.Y.; formal analysis, M.Y.; investigation, M.Y.; resources, M.Y.; data curation, M.Y.; writing – original draft preparation, M.Y.; writing – review and editing, M.Y. and S.L.; visualization, M.Y.; supervision, S.L.; All authors have read and agreed to the published version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was funded by the National High-tech R&D Program (863 Program) [grant number 2012AA112201].