A characteristic-featured shock wave indicator for simulating high-speed inviscid flows on 3D unstructured meshes
Yiwei Feng1,2,
Tiegang Liu1,
Xiaofeng He1,
Bin Zhang1 &
Kun Wang ORCID: orcid.org/0000-0002-8293-78461,3
In this work, we extend the characteristic-featured shock wave indicator based on artificial neuron training to 3D high-speed flow simulation on unstructured meshes. The extension is achieved through dimension splitting. This leads to that the proposed indicator is capable of identifying regions of flow compression in any direction. With this capability, the indicator is further developed to combine with h-adaptivity of mesh refinement to improve resolution with less computational costs. The present indicator provides an attractive alternative for constructing high-resolution, high-efficiency shock-processing methods to simulate high-speed inviscid flows.
Recently, high-order, high-resolution schemes are widely used in computational fluid dynamics for their strength in capturing detailed features of wave structures. Specially, Discontinuous Galerkin (DG) methods [1–4] can treat wave propagation well due to its upwind numerical flux which solves a Riemann problem arising from the discontinuous representation of the solution at element interfaces. As compared to finite difference (FD) and finite volume (FV) methods, DG schemes are easier to achieve high order and easier to handle complex boundaries on unstructured meshes. However, DG methods have a higher probability of making solutions suffer from unphysical oscillations as compared with low-order methods. The oscillations are produced by over-fitting shock wave discontinuities using high-order degrees of freedom (DoFs), and usually cause robustness issues of simulating high-speed flows. In addition, DG is recognized as more expensive in terms of computational costs for that it requires more DoFs to achieve the same accuracy as compared to classical finite element (FE) methods.
To eliminate unphysical oscillations near shock/contact wave discontinuities and improve robustness of DG methods on simulations of high-speed flows, a number of shock-processing (correction) techniques are designed to correct high-moment DoFs between spatial discretization and temporal discretization, such as limiting techniques [3, 5], solution reconstruction [6–8], and artificial diffusion [9–11]. However, the use of correction techniques usually leads to inefficiencies of numerical schemes or even deterioration of accuracy. Therefore, it is necessary to employ a troubled-cell or shock wave indicator to guide accurate application of correction techniques. There are several practical troubled-cell indicators, such as the classical Total Variation Bounded (TVB) based indicator [3], KXRCF shock detector [12], Harten's troubled-cell indicator [13], the posteriori subcell-based indicator in [14], and a simple indicator in [15]. There are still very recent Artificial Neural Network (ANN) based troubled-cell indicators in [16–18] to improve the accuracy of discontinuity capturing. However, these indicators mentioned above mainly detect troubled-cells based on the posteriori errors of numerical oscillations, usually cannot theoretically avoid detecting noise such as wrongly marked extremums, and part of them need increased computational costs to evaluate their values.
In our recent work, we employ an Artificial Neuron (AN) with characteristic-compression-imbedded data to propose a characteristic-featured shock wave indicator. The resultant indicator is concise, compact and explicable, and we have proved in [19] that its output can measure the compression of characteristic curves, which evolve into shock wave discontinuities. The indicator has been extended to 2D unstructured meshes and its performance combined with various correction techniques was presented in [20]. In this work, we develop this indicator to 3D unstructured meshes via dimension by dimension. This extension technique makes the present indicator capable of detecting compressible waves in any direction. Such a feature enables the present indicator to work closely with h-adaptivity of mesh refinement to further improve resolution and efficiency of simulation of high-speed inviscid flows.
The rest of the paper is arranged as follows. In Section 2, we briefly introduce the governing equations, spatial Discontinuous Galerkin (DG) discretization, temporal discretization and correction techniques including artificial viscosity. The description of the characteristic-featured shock wave indicator, development on 3D unstructured meshes, and some corresponding properties are presented in Section 3. Mesh refinement strategy based on the present shock wave indicator is presented in Section 4. Numerical experiments are provided in Section 5 to show the performance of the present indicator combined with various correction techniques and mesh refinement on multi-dimensional unstructured meshes. We make a few conclusions in Section 6.
Numerical scheme framework
We consider a multi-dimensional conservation law
$$ \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \mathbf{F}(\mathbf{U})= 0, $$
where U∈Rm is the vector of conserved variable and F(U)=(Fd(U))∈Rm×Rd. m, d represent the dimension of conserved vector and spatial space, respectively.
In this work, we focus on the Euler equations, where the conservative variable vector U and the advection flux vector F(U) in 3D space are
$$ \mathbf{U}= \left(\begin{array}{c} \rho \\ \rho u_{1} \\ \rho u_{2} \\ \rho u_{3} \\ E \end{array}\right),\;\; \mathbf{F}_{1}= \left(\begin{array}{c} \rho u_{1} \\ \rho u_{1}^{2}+p \\ \rho u_{1}u_{2} \\ \rho u_{1}u_{3} \\ u_{1}(E+p) \end{array}\right), \;\; \mathbf{F}_{2}= \left(\begin{array}{c} \rho u_{2} \\ \rho u_{1}u_{2} \\ \rho u_{2}^{2}+p \\ \rho u_{2}u_{3} \\ u_{2}(E+p) \end{array}\right), \;\; \mathbf{F}_{3}= \left(\begin{array}{c} \rho u_{3} \\ \rho u_{1}u_{3} \\ \rho u_{2}u_{3} \\ \rho u_{3}^{2}+p \\ u_{3}(E+p) \end{array}\right), $$
here ρ, p, E denote the density, pressure, and total energy of the fluid, respectively. u=(u1,u2,u3) is the velocity vector of the flow in the coordinate direction. The pressure can be computed from the equation of state
$$ p=(\gamma-1)\left(E-\frac{1}{2}\rho u_{j}^{2}\right) $$
Here, γ=Cp/CV denotes the ratio of specific heats. We set γ=1.4 for all the test cases in this work.
General DG formulation
The Discontinuous Galerkin (DG) scheme is used to discretize the weak form of the governing Eq. (1) over the domain Ω
$$ \int_{\Omega}\frac{\partial \mathbf{U}}{\partial t}\;\phi\;d\Omega+ \int_{\Gamma}\mathbf{F}_{d}\cdot \mathbf{n}_{d} \; \phi\;d\Gamma - \int_{\Omega}\mathbf{F}_{d}\cdot \nabla \phi \; d\Omega =0, $$
where Γ(=∂Ω) denotes the boundary of Ω and nd is the unit outward normal vector to the interface boundary. We assume that the domain is subdivided into a collection of non-overlapping elements Ωe. We introduce the following broken space of vector-values polynomials with a degree of p
$$ \mathbf{V}_{h}^{p}=\{ v\in L^{2}(\Omega):v|_{\Omega_{e}}\in \mathbf{P}^{p}(\Omega_{e}) \}. $$
Provided the basis functions \(\{\phi _{i}^{e}(\mathbf {x})\}\) of \(\mathbf {V}_{h}^{p}\), the local numerical solution Uh on Ωe can be computed by the basis function
$$ \mathbf{U}_{h}(\mathbf{x},t^{n})=\sum_{i=1}^{p} \mathbf{U}_{i}^{n}\cdot \phi_{i}^{e}(\mathbf{x}). $$
The Taylor basis functions [21] are used in this work. For sake of simplicity, let us consider a quadratic polynomial solution Uh in 3D space which is expressed as
$$ \begin{aligned} \mathbf{U}_{h}=&\bar{\mathbf{U}}\;\phi_{0}(\mathbf{x})+\frac{\partial \mathbf{U}}{\partial x}|_{c}\Delta x\;\phi_{1}(\mathbf{x})+\frac{\partial \mathbf{U}}{\partial y}|_{c}\Delta y\;\phi_{2}(\mathbf{x})+\frac{\partial \mathbf{U}}{\partial z}|_{c}\Delta z\;\phi_{3}(\mathbf{x}) \\ &+\frac{\partial^{2} \mathbf{U}}{\partial x^{2}}|_{c}\Delta x^{2}\;\phi_{4}(\mathbf{x})+\frac{\partial^{2} \mathbf{U}}{\partial y^{2}}|_{c}\Delta y^{2}\;\phi_{5}(\mathbf{x})+\frac{\partial^{2} \mathbf{U}}{\partial z^{2}}|_{c}\Delta z^{2}\;\phi_{6}(\mathbf{x})\\ &+\frac{\partial^{2} \mathbf{U}}{\partial x \partial y}|_{c}\Delta x\Delta y\;\phi_{7}(\mathbf{x})+\frac{\partial^{2} \mathbf{U}}{\partial y \partial z}|_{c}\Delta y\Delta z\;\phi_{8}(\mathbf{x})+\frac{\partial^{2} \mathbf{U}}{\partial z \partial x}|_{c}\Delta x\Delta z\;\phi_{9}(\mathbf{x}), \end{aligned} $$
$$ \begin{aligned} &\phi_{0}(\mathbf{x})=1, \quad\phi_{1}(\mathbf{x})=\frac{x-x_{c}}{\Delta x}, \quad\phi_{2}(\mathbf{x})=\frac{y-y_{c}}{\Delta y},\quad \phi_{3}(\mathbf{x})=\frac{z-z_{c}}{\Delta z}\\ &\phi_{4}(\mathbf{x})=\frac{(x-x_{c})^{2}}{2\Delta x^{2}}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(x-x_{c})^{2}}{2\Delta x^{2}}d\Omega, \\ &\phi_{5}(\mathbf{x})=\frac{(y-y_{c})^{2}}{2\Delta y^{2}}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(y-y_{c})^{2}}{2\Delta y^{2}}d\Omega, \\ &\phi_{6}(\mathbf{x})=\frac{(z-z_{c})^{2}}{2\Delta z^{2}}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(z-z_{c})^{2}}{2\Delta z^{2}}d\Omega, \\ &\phi_{7}(\mathbf{x})=\frac{(x-x_{c})(y-y_{c})}{\Delta x\Delta y}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(x-x_{c})(y-y_{c})}{\Delta x\Delta y}d\Omega, \\ &\phi_{8}(\mathbf{x})=\frac{(y-y_{c})(z-z_{c})}{\Delta y\Delta z}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(y-y_{c})(z-z_{c})}{\Delta y\Delta z}d\Omega, \\ &\phi_{9}(\mathbf{x})=\frac{(x-x_{c})(z-z_{c})}{\Delta x\Delta z}-\frac{1}{|\Omega_{e}|}\int_{\Omega_{e}}\frac{(x-x_{c})(z-z_{c})}{\Delta x\Delta z}d\Omega, \\ \end{aligned} $$
where Δx=0.5(xmax−xmin), Δy=0.5(ymax−ymin), Δz=0.5(zmax−zmin), and (xc,yc,zc) is the centroid of the element Ωe.
Replacing test function ϕ in (4) by elemental basis function \(\phi _{i}^{e}(\mathbf {x})\), we can obtain the semi-discrete form by applying weak formulation on each element Ωe,
$$ \frac{d}{dt}\int_{\Omega_{e}}\mathbf{U}_{h}\phi_{i}^{e}\;d\Omega+\int_{\Gamma_{e}}\widehat{\mathbf{F}}_{d} (\mathbf{U}_{h}^{-},\mathbf{U}_{h}^{+};\mathbf{n}_{d})\phi_{i}^{e}d\Gamma - \int_{\Omega_{e}}\mathbf{F}_{d}(\mathbf{U}_{h})\cdot \nabla\phi_{i}^{e}\;d\Omega=\mathbf{0}, $$
where i is the index of polynomial basis on each element \(\Omega _{e}, \widehat {\mathbf {F}}_{d}\) is a consistent numerical flux.
Temporal discretization
The spatial discretization (9) leads to a system of ordinary differential equations,
$$ \mathbf{M}\frac{d\mathbf{U}}{dt}+\mathbf{R}(\mathbf{U})=\mathbf{0}, $$
where M denotes the mass matrix, U is the solution vector, and R(U) is the residual vector with the form of
$$ \begin{aligned} &\mathbf{M} = \left((\phi_{i}^{e},\phi_{j}^{e})_{\Omega_{e}}\right),\quad \mathbf{U}_{h} = \mathbf{U}\cdot (\phi_{1}^{e}, \phi_{2}^{e},.., \phi_{9}^{e})_{e\;loops} \\ &\mathbf{R}(\mathbf{U}) = \int_{\Gamma_{e}}\widehat{\mathbf{F}}_{d}\;\phi_{i}^{e}d\Gamma - \int_{\Omega_{e}}\mathbf{F}_{d}(\mathbf{U}_{h})\cdot \nabla\phi_{i}^{e}\;d\Omega \end{aligned} $$
A high-order and stable explicit temporal discretization SSP-RK3 [22] scheme is applied to treat unsteady problems. This work focuses more on steady problems, in order to speed up the convergence and improve the stability of simulation for steady problems, implicit temporal discretization is presented in this work. Using Euler implicit time-integration, the spatial discretized Eq. (9) can be linearized in the temporal direction and written as
$$ \begin{aligned} &\mathbf{M}\frac{\Delta\mathbf{U}^{n}}{\Delta t}+\mathbf{R}(\mathbf{U}^{n+1})=\mathbf{0}\\ \Longleftrightarrow &\mathbf{M}\frac{\Delta\mathbf{U}^{n}}{\Delta t}+\mathbf{R}(\mathbf{U}^{n})+(\frac{\partial \mathbf{R}}{\partial \mathbf{U}})^{n}\Delta \mathbf{U}^{n}=\mathbf{0}\\ \Longleftrightarrow &(\frac{\mathbf{M}}{\Delta t}\mathbf{I}+(\frac{\partial \mathbf{R}}{\partial \mathbf{U}})^{n})\Delta \mathbf{U}^{n}=-\mathbf{R}(\mathbf{U}^{n}), \end{aligned} $$
where \(\frac {\partial \mathbf {R}}{\partial \mathbf {U}}\) is the Jacobian matrix, Δt is the time increment, and ΔUn=Un+1−Un is the solution difference between time level n and n+1. This system of linear equations at each time step is solved by the GMRES method with a LU-SGS preconditioner [23].
Correction techniques
As mentioned in the introduction, there are a number of correction techniques to treat shock wave discontinuities. For example, we employ the Barth-Jespersen limiter [24] in an explicit temporal scheme to treat unsteady problems in [20], and employ the artificial viscosity with strong-residual term in an implicit temporal scheme to treat steady problems. Here, we describe artificial viscosity in details.
The artificial viscosity technique can treat shock waves through adding an adaptive artificial viscosity term to bring parabolic properties into hyperbolic conservation laws (1) to smear the shock wave discontinuities. Thus, the semi-DG scheme is modified into
$$ \begin{aligned} \frac{d}{dt}\int_{\Omega_{e}}\mathbf{U}_{h}\phi_{e}^{i}\;d\Omega &+\int_{\Gamma_{e}}\widehat{\mathbf{F}}_{k}(\mathbf{U}_{h}^{-},\;\mathbf{U}_{h}^{+};\mathbf{n}_{k})\phi_{e}^{i}d\Gamma - \int_{\Omega_{e}}\mathbf{F}_{k}(\mathbf{U}_{h})\cdot \nabla\phi_{e}^{i}\;d\Omega\\ & + \epsilon_{e}\;\int_{\Omega_{e}}\nabla\mathbf{U}_{h}\cdot \nabla\phi_{e}^{i}\;d\Omega=0. \end{aligned} $$
Parameter εe in this work is a Hartmann-type coefficient of strong residual
$$ \begin{aligned} \epsilon_{e}=&C_{\epsilon}h_{e}^{2-\beta}\frac{\int_{\Omega_{e}}\left|\nabla \cdot \mathbf{F}(\mathbf{u}_{h})\right|d\Omega}{|\Omega_{e}|}, \end{aligned} $$
where Cε and β are problem and mesh dependent parameters which are usually set to 0.0001−0.01,0.01−0.5, respectively.
The characteristic-featured shock wave indicator
In our recent work [19, 20], we proposed a characteristic-featured indicator based on training an artificial neuron (AN) to detect shock wave discontinuities. We firstly describe the indicator on 1D mesh, next develop the shock wave indicator on 3D unstructured meshes, then show some properties of the present indicator, and finally provide its application prospects as a shock wave indicator.
Indicator on 1D mesh
For 1D scalar conservation laws
$$\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial x}=0, $$
we consider a numerical solution uh(x,t) on a perturbed mesh, he denotes the mesh-size of element \(I_{e}, \bar {u}_{e}=\frac {1}{h_{e}}\int _{I_{e}} u_{h}(x,t)dx\) is the cell-avaerge of variable u inside Ie. The characteristic-featured shock wave indicator through training an artificial neuron (AN) [20] using the information of neighbours (Ie−1,Ie,Ie+1) is presented as follows,
$$ \widetilde{out}_{e}(\bar{\lambda}_{L},\bar{\lambda}_{R};h_{e-1},h_{e},h_{e+1})=\frac{1}{1+e^{-W(\bar{\lambda}_{L}-\bar{\lambda}_{R})+M_{1}(h_{e-1}+h_{e}+h_{e+1})+M_{2}}}. $$
Here, W=10,M1=4,M2=1, and \(\bar {\lambda }_{L}, \bar {\lambda }_{R}\) are the left and right side weighted averages of eigenvalue λ(u)=f′(u) which are convex combinations of cell-averages \(\bar {\lambda }_{e-1},\bar {\lambda }_{e}\) and \(\bar {\lambda }_{e+1}, \bar {\lambda }_{e}\), respectively (Refer to Fig. 1), in [20], \(\bar {\lambda }_{L}, \bar {\lambda }_{R}\) are set to integral-averages
$$ \left\{ \begin{array}{lr} \bar{\lambda}_{L} = \frac{h_{e-1}}{h_{e-1}+h_{e}}\cdot \bar{\lambda}_{e-1}+\frac{h_{e}}{h_{e-1}+h_{e}}\cdot \bar{\lambda}_{e} \\ \bar{\lambda}_{R} = \frac{h_{e+1}}{h_{e+1}+h_{e}}\cdot \bar{\lambda}_{e+1}+\frac{h_{e}}{h_{e+1}+h_{e}}\cdot \bar{\lambda}_{e}. \end{array} \right. $$
Geometry information for 1D scalar indicating
In practice, we found M1 in (15) is not sensitive to the convergence of AN. We retrain an AN with 3 arguments \(\bar {\lambda }_{L},\bar {\lambda }_{R},h_{M}\) (hM:=max{he−1,he,he+1}) and without hidden layer, and the AN convergents as well and maintains a similar working mechanism. In order to make the indicator more concise, we adopt a model with less input arguments and keep most coefficients in (15) unchanged. The expression of the modified shock wave indicator is present as follows,
$$ \widetilde{out}_{e}(\bar{\lambda}_{L},\bar{\lambda}_{R};h_{M})=\frac{1}{1+e^{-W(\bar{\lambda}_{L}- \bar{\lambda}_{R})+M_{1}h_{M}+M_{2}}}. $$
Here, W=10,M2=1 keep unchanged, \(M_{1}=12, \bar {\lambda }_{L},\bar {\lambda }_{R}\) are defined by (16) and hM=max{he−1,he,he+1} is the local maximum mesh-size. We drop the subscript "M" in the following text for convenience.
In [19], we have proved that the output of indicator (17) measures the admissible jump values of eigenvalues, the greater the value of \(\bar {\lambda }_{L}-\bar {\lambda }_{R}\) is, the stronger the admissible jump is (a big admissible jump value implies characteristic compressing and leads to the development of a shock wave), at the same time the output of the indicator (17) is closer to 1.
Extension to system
For 1D system of conservation laws (1), the indicator is extended by characteristic field splitting, that is to apply the scalar indicator (17) to each genuinely nonlinear characteristic field to detect the compressing or intersection of characteristic curves (refer to [25]). The expression is presented as follows,
$$ \widehat{out}_{e}={max}_{i}\left\{\widetilde{out}_{e}(\bar{\lambda}_{i,L},\bar{\lambda}_{i,R};h)\right\}. $$
Here, \(\bar {\lambda }_{i,L},\bar {\lambda }_{i,R}\) are the left and right side-weighted cell-averages of eigenvalue λi of i-th genuinely nonlinear characteristic field, respectively, and \(\widetilde {out}_{e}(\bar {\lambda }_{i,L},\bar {\lambda }_{i,R};h)\) is defined in (17). Especially, for the 1D Euler equations, the indicator (18) becomes
$$ \widehat{out}_{e}=max\left\{\widetilde{out}_{e}(\;(\bar{u}-\bar{a})_{L},(\bar{u}-\bar{a})_{R};h);\;\widetilde{out}_{e}(\;(\bar{u}+\bar{a})_{L},(\bar{u}+\bar{a})_{R};h)\right\}. $$
where \(a=\sqrt {\frac {\gamma p}{\rho }}\) represents the speed of sound. It should be noted that the contact discontinuity can be detected with the indicator (19), because there is always an admissible jump occurring to the 1- or 3- eigenvalues across the 2-contact discontinuity.
Development on 3D unstructured meshes
The indicator (17) is extended to 3D unstructured grids through dimension splitting, that is to construct side-weighted averages of eigenvalue \(\bar {\lambda }\) on each dimension, and then apply (17) in each spatial dimension to accomplish the indication. Note that compression of characteristic cones in any direction results in compression of the characteristic curves projected onto the x or y or z-direction. Therefore, multi-dimensional indicator through dimension splitting can capture shock discontinuities aligned to any directions, although it is only applied along x, y and z-directions.
Thus, careful evaluation of the weights of neighbours to construct \(\bar {\lambda }\) becomes crucial to avoid the influence of mesh shape and mesh size. Here, we present a relatively simple construction method shown in our recent work [20].
We set the 3D Euler equations on the following mesh (refer to Fig. 2) as an example to show the extension procedure. Firstly, we present the following notation to be used,
Geometry information for 3D Euler indicating on unstructured grids
Ω1,Ω2,Ω3 are the neighbours of Ω0,bk=(xk,yk,zk) is the centroid of Ωk,|Ωk| represents the volume of Ωk,Ak represents the volume ratio with the form of \(\frac {|\Omega _{k}|}{|\Omega _{0}|+|\Omega _{k}|}, \mathbf {n}_{k}=(n_{kx},n_{ky},n_{kz})\) represents the unit vector of centroid connection vector \(\frac {b_{0}b_{k}}{|b_{0}b_{k}|}, \lambda _{\Omega _{k}}\) represents cell-averages of eigenvalues in Ωk in the x,y,z direction, λ=u1±a,u2±a,u3±a.
Based on the methods shown in [20], in order to construct side-weighted averages of eigenvalue \(\bar {\lambda }\) and reduce the influence of mesh shape and mesh size, we need to solve two key issues, (a) one is how to select the left or right side element of Ω0 on each spatial dimension among all neighbours, such as right-side element of Ω0 in the z-direction (named as top-side element); (b) the other is how to select the "weight" for each side-element, such as how to evaluate the weights for each top-side element when constructing weighted average \(\bar {\lambda }_{T}\). The procedure of constructing side-weighted averages \(\bar {\lambda }\) on each spatial dimension is presented as follows.
Step-1: Selection of left, right side-element on each dimension. The direction of nk is used to determine the position of Ωk to Ω0. Ωk is recorded as the backward element of Ω0 if nkx<0, similarly,
$$\begin{aligned} &\Omega_{k}\; is \; left\; side \; element \qquad \quad\; if \;n_{ky} < 0 \\ &\Omega_{k}\; is \; bottom\; side \; element \qquad if \;n_{kz} < 0. \end{aligned} $$
For example, refer to Fig. 2, Ω1 is marked as a forward element in the x-direction, a left-side element in the y-direction, and a top-side element in the z-direction.
Step-2: Evaluation of weights of each side-element. As the position of Ωk to Ω0 on each dimension is determined, the next is to evaluate the weights that \(\lambda _{\Omega _{k}}\) contributes for constructing side-weighted averages \(\bar {\lambda }\) on each dimension. The weights are products of two parts, the direction weight and the size weight,
(a) direction weight. The direction weight represents the contribution of Ωk on each spatial dimension, for \(n_{kx}^{2}+n_{ky}^{2}+n_{kz}^{2}=1\), the greater \(n_{kx}^{2}\) is, the more contribution Ωk makes in the x-direction, while the less contribution Ωk makes in the y,z direction. Therefore, \(n_{kx}^{2}\) is selected as the direction weight for constructing backward-weighted average \(\bar {\lambda }_{b}\) or forward-weighted average \(\bar {\lambda }_{f}\), and \(n_{ky}^{2}, n_{kz}^{2}\) are selected as the direction weight for constructing \(\bar {\lambda }_{L}\) or \(\bar {\lambda }_{R}, \bar {\lambda }_{B}\) or \(\bar {\lambda }_{T}\), respectively.
(b) size weight. The mesh size needs to be considered as another factor for weights evaluating. The greater the volume |Ωk| is, the more contribution Ωk makes on each dimension. Therefore, we set the volume ratio \(A_{k}=\frac {|\Omega _{k}|}{|\Omega _{0}|+|\Omega _{k}|}\) as the size weight for constructing \(\bar {\lambda }\) on each dimension.
Based on Step-2 and refer to Fig. 2, we can construct the forward-weighted average \(\bar {\lambda }_{\Omega _{1}}^{f}\) contributed by Ω1, which is a convex combination of \(\bar {\lambda }_{\Omega _{1}}\) and \(\bar {\lambda }_{\Omega _{0}}\),
$$ \begin{aligned} &\bar{\lambda}^{f}_{\Omega_{1}}=w_{1}^{f} \bar{\lambda}_{\Omega_{1}}+w_{0}^{f} \bar{\lambda}_{\Omega_{0}}\quad (w_{0}^{f}+w_{1}^{f}=1),\\ &\tilde{w}_{1}^{f} = (n_{1x}^{2})\cdot (A_{1}), \qquad \tilde{w}_{0}^{f} = (1-n_{1x}^{2})\cdot (1-A_{1}),\\ &w_{1}^{f} = \tilde{w}_{1}^{f}/(\tilde{w}_{1}^{f}+\tilde{w}_{0}^{f}),\qquad w_{0}^{f}=1-w_{1}^{f}. \end{aligned} $$
where \(n_{1x}^{2}\) is the direction weight and A1 is the size weight for Ω1 in the x-direction, \(w_{1}^{f}\) is normalized weight combined by both direction and size weights.
Similarly, we have
$$ \begin{aligned} &\bar{\lambda}^{L}_{\Omega_{1}}=\frac{\tilde{w}_{0}^{L}}{\tilde{w}_{1}^{L}+\tilde{w}_{0}^{L}} \bar{\lambda}_{\Omega_{0}}+\frac{\tilde{w}_{1}^{L}}{\tilde{w}_{1}^{L}+\tilde{w}_{0}^{L}} \bar{\lambda}_{\Omega_{1}}\quad with \; \tilde{w}_{1}^{L} = n_{1y}^{2}\cdot A_{1},\; \tilde{w}_{0}^{T} = (1-n_{1y}^{2})\cdot (1-A_{1}) \\ &\bar{\lambda}^{T}_{\Omega_{1}}=\frac{\tilde{w}_{0}^{T}}{\tilde{w}_{1}^{T}+\tilde{w}_{0}^{T}} \bar{\lambda}_{\Omega_{0}}+\frac{\tilde{w}_{1}^{T}}{\tilde{w}_{1}^{T}+\tilde{w}_{0}^{T}} \bar{\lambda}_{\Omega_{1}}\quad with \; \tilde{w}_{1}^{T} = n_{1z}^{2}\cdot A_{1},\; \tilde{w}_{0}^{T} = (1-n_{1z}^{2})\cdot (1-A_{1}). \end{aligned} $$
At last, we assume all the same side-elements on each spatial dimension have an equivalent opportunity to contribute to the construction of side-weighted averages. For example, refer to Fig. 2, only Ω1 is determined as a left-side element in the y-direction, both Ω1 and Ω2 are determined as top-side elements in the z-direction, thus we have
$$ \bar{\lambda}_{L}=\bar{\lambda}^{L}_{\Omega_{1}},\quad\bar{\lambda}_{T}=\frac{1}{2}(\bar{\lambda}^{R}_{\Omega_{1}}+\bar{\lambda}^{R}_{\Omega_{2}}). $$
In addition, we define the characteristic mesh-size hx,hy,hz to be used as follows,
$$ h_{x}=max\left\{|x_{k}-x_{0}|\right\},\quad h_{y}=max\left\{|y_{k}-y_{0}|\right\},\quad h_{z}=max\left\{|z_{k}-z_{0}|\right\}. $$
According to the steps above and combined with formulae (20), (21) and (22), side-weighted averages of eigenvalue \(\bar {\lambda }_{b},\bar {\lambda }_{f}, \bar {\lambda }_{L},\bar {\lambda }_{R},\bar {\lambda }_{B},\bar {\lambda }_{T}\) can be constructed easily. Replacing h in (17) by hx in the x-direction, hy in the y-direction, hz in the z-direction, and \(\bar {\lambda }_{L},\;\bar {\lambda }_{R}\) in (17) by \(\bar {\lambda }_{b},\bar {\lambda }_{f}\) in the x-direction, \(\bar {\lambda }_{L},\bar {\lambda }_{R}\) in the y-direction, and \(\bar {\lambda }_{B},\bar {\lambda }_{T}\) in the z-direction, we can obtain a concise and generalized characteristic-based shock wave indicator on unstructured grids as follows,
$$ \widehat{out}_{e}^{3D}=max\left\{ \widehat{out}_{e}(\bar{\lambda}_{b},\;\bar{\lambda}_{f};\;h_{x}),\;\widehat{out}_{e}(\bar{\lambda}_{L},\;\bar{\lambda}_{R};\;h_{y}),\;\widehat{out}_{e}(\bar{\lambda}_{B},\;\bar{\lambda}_{T};\;h_{z}) \right\}. $$
Here, \(\widehat {out}_{e}\) is the indicator for system equations (18), \(\bar {\lambda }_{b},\bar {\lambda }_{f}, \bar {\lambda }_{L},\bar {\lambda }_{R},\bar {\lambda }_{B},\bar {\lambda }_{T}\) are defined by (20), (21), (22).
(24) is used as the shock wave indicator in this work. Based on statements in Section 3.1, we can conclude that (24) can detect shock wave discontinuities caused by compressing or intersecting of characteristic. Especially it can detect contact waves in the Euler equations.
Indicating procedure based on the present indicator
We present two approaches to detect troubled-cells having shock wave discontinuities based on the present indicator (24) on arbitrary grids. One is to flag cells by a critical value s, that is
$$ \begin{aligned} \left\{ \begin{array}{lr} \widehat{out}_{e}^{3D} \geq s \Longrightarrow \;\Omega_{e}\;is\;flagged\;as\;a\;troubled\;cell.\\ \widehat{out}_{e}^{3D} < s \Longrightarrow \;\Omega_{e}\;is\;not\;treated. \end{array} \right. \end{aligned} $$
The parameter s is set to 0.5, while s can be set smaller for identifying weaker admissible jump strength. This approach is suitable to be combined with TVB or WENO limiting techniques when using explicit temporal schemes to simulate unsteady flows.
The other is to flag cells with a fixed fraction θ. Only the top θ fraction of elements with output of (24) is flagged. This approach is suitable to be applied with artificial viscosity in implicit temporal schemes to treat steady flows, and suitable to be applied in adaptive mesh-refinement strategies as well. This use will be described in the next section.
Mesh adaptive method
Adaptive methods provide an attractive alternative for reducing computational costs of DG schemes. Appropriate adaptive methods are able to focus computational resources on needed regions to further capture the detailed features of wave structures (such as shock and contact waves, vortex) and avoid the additional increase in computational costs caused by global mesh refinement or global order enrichment.
The essence of adaptive methods is how to flag elements which need focusing on. As mentioned in the last section, the indicator (24) is able to measure the compressing of characteristics. This feature results in the indicator being naturally suitable as a special error estimator which measures the unidirectional variation of eigenvalues to capture specific wave structure such as shock, contact wave and vortex through detecting the compressing of characteristics.
There are two strategies for mesh refinement once the local error estimator is determined. One is adaptively refining the mesh with a fixed fraction (each time only the top or bottom fixed percent of mesh is refined or coarsened); the other is refining the mesh with a fixed number (each time only a fixed number of mesh is refined or coarsened). In this work we select the former strategy based on the following error estimator for mesh (quadrilateral grid) refinement,
$$ \eta_{e}:={max}_{d}\left\{\widehat{out}^{d}_{e}(\bar{\lambda}_{L,d},\bar{\lambda}_{R,d};h_{d})\right\}. $$
Here, the subscript d denotes the dimension of space, \(\widehat {out}^{d}_{e}\) is the 1D indicator (17), \(\bar {\lambda }_{L,d},\bar {\lambda }_{R,d}\) are side-weighted averages of eigenvalues in the d-direction defined by (20), (21), (22), and hd is the local characteristic mesh-size in the d-direction defined in (23).
The whole adaptive algorithm based on the above error estimator (26) can be briefly summarized in Algorithm 1 as follows:
The performance of the proposed adaptive method is presented and discussed in Section 5.
Numerical results
We now provide the performance of the present indicator (24) as employed to detect shock waves on multi-dimensional meshes. All the test cases presented in this section are governed by the Euler equations and discretized by the DG spatial discretization and the local Lax-Friedrichs type numerical flux.
2D test cases
2D Riemann problem
We first test the present indicator on a 2D unstructured mesh (S-\(\frac {1}{400}\)) to treat unsteady problems. The computational domain is [−0.5,0.5]×[−0.5,0.5], S-\(\frac {1}{400}\) mesh is obtained by a 400×400 structured quadrilateral mesh (mesh-size is \(\frac {1}{400}\)), and diving each quadrilateral element into two triangles. We apply the present indicator with a fixed value (s=0.5) and combined with the Barth-Jespersen limiter in an explicit temporal scheme to solve a 2D Riemann problem with the following initial conditions,
$$(\rho,\;u,\;v,\;p)_{0}(\mathbf{x})= \left\{ \begin{array}{lr} (1.5,\;0,\;0,\;1.5), \qquad \qquad \qquad \qquad if\;x>0,\;y>0, &\\ (0.5323,\;1.206,\; 0,\;0.3) \qquad \qquad \quad\; if\;x<0,\;y>0, &\\ (0.138,\; 1.206,\; 1.206,\; 0.029) \qquad \quad if\;x<0,\;y<0, &\\ (0.5323,\; 0,\; 1.206,\;0.3) \qquad \qquad \quad \; if\;x>0,\;y<0, &\\ \end{array} \right. $$
We simulate the results till T=0.3 under SSP-RKDG(p2) framework, and TVB-based indicators are presented for comparison. TVB-based indicators are sensitive to the parameter M, and the following TVB-3 indicator denotes the TVB-based indicator with M=1000 (refer to [20] for details of TVB-based indicators). The results obtained by different indicators are shown in Fig. 3. The results for this problem show a significant improvement with the present indicator, which leads to a better shock and contact wave resolution in the zoomed regions in Fig. 4. Furthermore, the results in Fig. 5 show the present indicator marks thinner shock and contact waves regions and fewer noises as observed in the central regions.
Solutions (density) of the Riemann problem case for the 2D Euler equations with various indicators with \(S-\frac {1}{400}\) grid at T=0.3 under RKDG(p2) framework
Zoomed solution of the Riemann problem case for 2D Euler equations with \(S-\frac {1}{400}\) grid at T=0.3 under RKDG(p2) framework
Troubled-cells flagged by various indicators in the Riemann problem case for 2D Euler equations with \(S-\frac {1}{400}\) grid at T=0.3 under RKDG(p2) framework
Readers can refer to [20] for more unsteady test cases.
2D steady flow
We now employ the present indicator on 2D unstructured meshes to treat steady problems. We apply the present indicator with a fixed fraction (θ=0.05) and combined with the artificial viscosity in an implicit temporal scheme to simulate a transonic inviscid flow (Ma=0.8,α=1.25∘) around the NACA0012 airfoil. A rather fine mesh is used as shown in Fig. 6, DG(p2) discretization is used in spatial direction, Cε and β in artificial viscosity (AV) term are set to 0.01 and 0.01, respectively. The convergent solution flow (density) is obtained by the present indicator (24), as presented in Fig. 7.(a), and the output values of the present indicator are presented in Fig. 7.(b). One can observe from these results that the present indicator can detect the shock regions with high resolution and low noise. Besides, to illustrate the present indicator can lead to a stably convergent status and maintain physics properties well, residual convergence history in Fig. 8 and pressure coefficient distributions on the surface of the airfoil as compared with experiment data from [26] in Fig. 9 are provided below as well,
A fine mesh for NACA0012 test case in the 2D Euler equations
Solution (density) and indicator in NACA0012 test case with Ma=0.8,α=1.25∘ in the 2D Euler equations on refined grids
Residual convergence history of the transonic flow over NACA0012 airfoil
Pressure coefficient distributions on the surface of the NACA0012 airfoil
In addition, indication with a fixed fraction performs efficiently as combined with artificial viscosity to treat high-speed steady flows. Readers can refer to [20] for efficiency comparison and more results.
This part is to extend the present indicator to 3D test cases and observe its performance when applied to treat steady flows.
Inviscid flow at Ma=0.8395, α=3.06∘ around the Onera M6 wings
We use a mesh file (Tetrahedral mesh with 621508 elements) as shown in Fig. 10, DG(p2) discretization is used in spatial direction, Cε and β in AV term are set to 0.1 and 0.01, respectively. The convergent solution flow (density) is obtained by the present indicator (24) with a fixed fraction of 0.05, as presented in Fig. 11.(a), and the output values of the present unstructured indicator are presented in Fig. 11.(b). Residual convergence history in Fig. 12 and chordwise pressure distributions at several spanwise stations on M6 wing as compared with experiment data from [27] in Fig. 13 are provided as follows.
Fig. 10
Mesh for M6 wing test case in the 3D Euler equations
Solution (density) and indicator in M6 wing test case with Ma=0.84,α=3.06∘ in the 3D Euler equations
Residual convergence history of the transonic flow over M6 wing
Chordwise pressure distributions at several spanwise stations on M6 wing
Inviscid flow at Ma=0.75, α=0∘ around the DLR-F4 body
We use a mesh file (Tetrahedral mesh with 550362 elements) as shown in Fig. 14, DG(p2) discretization is used in spatial direction, Cε and β in AV term are set to 0.01 and 0.2, respectively. The convergent solution flow (density) is obtained by the present indicator (24) with a fixed fraction of 0.05, as presented in Fig. 15.(a), and the output values of the present unstructured indicator are presented in Fig. 15.(b). Residual convergence history in Fig. 16 and chordwise pressure distributions at several spanwise stations on DLR-F4 body as compared with three groups of experiment data in Fig. 17 are provided as follows. Figure 17 shows the comparative results with experiment data. From results from CART3D and CFL3D shown in [28], we infer that the present methods by solving RANS equations or LES can obtain more perfect agreement with experimental results, which can be further developed in the future.
Mesh for DLR-F4 body test case in the 3D Euler equations
Solution (density) and indicator in DLR-F4 body test case
Residual convergence history of the transonic flow over DLR-F4 body
Chordwise pressure distributions at several spanwise stations on DLR-F4 body
From these two cases, one can observe the present indicator can not only detect the shock/compressing waves, but also identify the regions where characteristics compress much, such as the regions near the head and tail of the wing, which are regions that usually cause issues of convergence. As a result, employment of the present indicator in 3D can detect shock and compressing regions accurately, further leading to a stably convergent process and reduced computational costs in a steady problem.
Mesh adaptive refinement
The proposed h-adaptive method in Section 4 has been successfully applied to 2D test cases. Here, we present a simple case of transonic flow around NACA0012 airfoil to show its strength in terms of improving resolution with apparently reduced computational costs.
The initial mesh is rather coarse and is shown in Fig. 18.(a), and the initial solution (density) with low-resolution and the output of indicator are presented in Fig. 18.(b)(c). After 3 times of mesh refinement (θ=0.05), one can observe a significant improvement in terms of resolution in result (density), shown in Fig. 19.(b), and the reason can be explained from that the proposed error estimator (26) can capture shocks accurately (Fig. 19.(c)) and guide the computational resources (more refined mesh) focusing on the shock and compressing regions (Fig. 19.(c)).
Initial mesh and solution (density) of the transonic flow around the NACA0012 airfoil under DG(p2) framework
Final results of the transonic flow around the NACA0012 airfoil after 3 times of mesh refinement with θ=0.05
We described a characteristic-featured shock wave indicator for system of conservation laws and developed it on 3D unstructured meshes. The present indicator has been successfully applied as a local error estimator in h-adaptive (mesh refinement) method to further improve resolution with less computational costs.
The numerical results demonstrated that the present indicator has excellent and robust performance on shock wave detecting in different problems. It detects shock, contact waves with low noise and high efficiency, provides an attractive alternative as an error estimator to combine with h-adaptive method to design high-resolution and high-efficiency schemes.
All data generated or analysed during this study are included in this published article.
Cockburn B, Shu C-W (1998) The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J Num Anal 35(6):2440–2463.
Cockburn B, Lin S-Y, Shu C-W (1989) TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems. J Comput Phys 84(1):90–113.
Cockburn B, Shu C-W (1989) TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II: general framework. Math Comput 52(186):411–435.
MathSciNet MATH Google Scholar
Cockburn B, Hou S, Shu C-W (1990) The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV: the multidimensional case. Math Comput 54(190):545–581.
Van Leer B (1979) Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. J Comput Phys 32(1):101–136.
Harten A, Osher S (1997) Uniformly high-order accurate nonoscillatory schemes. I In: Upwind and High-Resolution Schemes, 187–217.. Springer, Berlin.
Shu C-W (1998) Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws In: Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, 325–432.. Springer, Berlin.
Qiu J, Shu C-W (2005) Runge–Kutta discontinuous Galerkin method using WENO limiters. SIAM J Sci Comput 26(3):907–929.
Jameson A, Schmidt W, Turkel E (1981) Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes In: 14th Fluid and Plasma Dynamics Conference, 1259. https://doi.org/10.2514/6.1981-1259.
Tadmor E (2003) Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numer 12:451–512.
Hesthaven JS, Warburton T (2007) Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer, Berlin.
Krivodonova L, Xin J, Remacle J-F, Chevaugeon N, Flaherty JE (2004) Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws. Appl Num Math 48(3-4):323–338.
Harten A (1989) ENO schemes with subcell resolution. J Comput Phys 83(1):148–184.
Dumbser M, Zanotti O, Loubère R, Diot S (2014) A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws. J Comput Phys 278:47–75.
Fu G, Shu C-W (2017) A new troubled-cell indicator for discontinuous Galerkin methods for hyperbolic conservation laws. J Comput Phys 347:305–327.
Wen X, Don WS, Gao Z, Hesthaven JS (2020) An edge detector based on artificial neural network with application to hybrid Compact-WENO finite difference scheme. J Sci Comput 83:1–21.
Ray D, Hesthaven JS (2018) An artificial neural network as a troubled-cell indicator. J Comput Phys 367:166–191.
Ray D, Hesthaven JS (2019) Detecting troubled-cells on two-dimensional unstructured grids using a neural network. J Comput Phys 397:108845.
Feng Y, Liu T, Wang K (2020) A characteristic-featured shock wave indicator for conservation laws based on training an artificial neuron. J Sci Comput 83(1):1–34.
Feng Y, Liu T (2021) A characteristic-featured shock wave indicator on unstructured grids based on training an artificial neuron. J Comput Phys:110446. https://doi.org/10.1016/j.jcp.2021.110446.
Luo H, Baum JD, Löhner R (2008) A discontinuous Galerkin method based on a taylor basis for the compressible flows on arbitrary grids. J Comput Phys 227(20):8875–8893.
Gottlieb S, Shu C-W, Tadmor E (2001) Strong stability-preserving high-order time discretization methods. SIAM Rev 43(1):89–112.
Luo H, Baum JD, Löhner R (1998) A fast, matrix-free implicit method for compressible flows on unstructured grids. J Comput Phys 146(2):664–690.
Barth T, Jespersen D (1989) The design and application of upwind schemes on unstructured meshes In: 27th Aerospace Sciences Meeting, 366. https://doi.org/10.2514/6.1989-366.
Toro EF (2013) Riemann solvers and numerical methods for fluid dynamics. Springer, Berlin. https://doi.org/10.1007/978-3-662-03490-3.
Yoshihara H, Norstrud H, Boerstoel J, Chiocchia G, Jones D (1985) Test cases for inviscid flow field methods. AGARD AR-211.
Schmitt V, Charpin F (1979) Pressure distributions on the Onera M6-wing at transonic mach numbers, experimental data base for computer program assessment. AGARD AR-138.
Rodriguez D, Sturdza P, Suzuki Y, Martins-Rivas H, Peronto A (2012) A rapid, robust, and accurate coupled boundary-layer method for cart3d In: 50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 302. https://doi.org/10.2514/6.2012-302.
The authors gratefully acknowledge the support of the National Numerical Wind Tunnel Project, the National Natural Science Foundation of China (No. 12001031), the Academic Excellence Foundation of BUAA for PhD Students, China Postdoctoral Science Foundation (No. 2020M680284). The authors are grateful to reviewers and editors for reviewing this work.
This work is supported by the National Numerical Wind Tunnel Project, the National Natural Science Foundation of China (No. 12001031), the Academic Excellence Foundation of BUAA for PhD Students, China Postdoctoral Science Foundation (No. 2020M680284).
School of Mathematical Sciences, Beihang University, Beijing, China
Yiwei Feng, Tiegang Liu, Xiaofeng He, Bin Zhang & Kun Wang
ShenYuan Honors College, Beihang University, Beijing, China
Yiwei Feng
School of Aeronautic Science and Engineering, Beihang University, Beijing, China
Kun Wang
Tiegang Liu
Xiaofeng He
Bin Zhang
The research output comes from a joint effort. All authors read and approved the final manuscript.
Correspondence to Kun Wang.
Feng, Y., Liu, T., He, X. et al. A characteristic-featured shock wave indicator for simulating high-speed inviscid flows on 3D unstructured meshes. Adv. Aerodyn. 3, 27 (2021). https://doi.org/10.1186/s42774-021-00084-9
High-order high-resolution method
Shock wave detector
H-adaptive method
Transonic supersonic flow