An Immersed Discontinuous Galerkin Method for Compressible Navier-Stokes Equations on Unstructured Meshes

TL;DR

This paper presents an immersed high-order discontinuous Galerkin method for solving compressible Navier-Stokes equations on unstructured meshes, enabling accurate simulations on complex geometries without boundary-fitted meshes.

Contribution

The paper introduces novel techniques for handling cut-elements and boundary representation in an immersed DG framework, improving stability and accuracy for complex flow simulations.

Findings

Validated on subsonic and hypersonic flows

Achieved high-order accuracy with curved boundary representation

Successfully simulated flows over complex geometries

Abstract

We introduce an immersed high-order discontinuous Galerkin method for solving the compressible Navier-Stokes equations on non-boundary-fitted meshes. The flow equations are discretised with a mixed discontinuous Galerkin formulation and are advanced in time with an explicit time marching scheme. The discretisation meshes may contain simplicial (triangular or tetrahedral) elements of different sizes and need not be structured. On the discretisation mesh the fluid domain boundary is represented with an implicit signed distance function. The cut-elements partially covered by the solid domain are integrated after tessellation with the marching triangle or tetrahedra algorithms. Two alternative techniques are introduced to overcome the excessive stable time step restrictions imposed by cut-elements. In the first approach the cut-basis functions are replaced with the extrapolated basis functions from the nearest largest element. In the second approach the cut-basis functions are simply scaled proportionally to the fraction of the cut-element covered by the solid. To achieve high-order accuracy additional nodes are introduced on the element faces abutting the solid boundary. Subsequently, the faces are curved by projecting the introduced nodes to the boundary. The proposed approach is verified and validated with several two- and three-dimensional subsonic and hypersonic low Reynolds number flow applications, including the flow over a cylinder, a space capsule and an aerospace vehicle.

Tables

Table 1. Table 1: Comparison of wake separation length L w / d subscript 𝐿 𝑤 𝑑 L_{w}/d .

$\mathrm{𝑅𝑒}$	20	40
Dennis and Chang [55]	0.94	2.35
Fornberg [56]	0.91	2.24
Ye et al. [53]	0.92	2.27
Emblemsvag et al. [54]	0.93	2.28
$P^2$ - scaled basis	0.920	2.280
$P^2$ - extended basis	0.929	2.284

Table 2. Table 2: Comparison of the computed drag coefficients C d subscript 𝐶 𝑑 C_{d} with the experiments.

	$\mathrm{𝑀𝑎}$=5.92, $\mathrm{𝐾𝑛}$=0.016	$\mathrm{𝑀𝑎}$=5.92, $\mathrm{𝐾𝑛}$=0.08
Exp. by Maslach et al. [57]	1.55	1.63
$P^2$ - scaled basis (error)	1.561 (0.71%)	1.647 (1.04%)
$P^2$ - extended basis (error)	1.563 (0.84%)	1.649 (1.17%)

Equations

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\text{inv}}(\mathbf{U})-\nabla \cdot {\mathbf{F}}_{\text{vis}}(\mathbf{U},\nabla \mathbf{U})=0$$

$$\mathbf{U}=\left(\begin{array}{c}\rho \\ \rho \mathbf{u}\\ \rho E\end{array}\right)$$

$${\mathbf{F}}_{\text{inv}}(\mathbf{U})$$

$${\mathbf{F}}_{\text{vis}}(\mathbf{U},\nabla \mathbf{U})$$

$$\mathrm{Ma}=\frac{u_{\text{ref}}}{(\gamma R{T}_{\text{ref}}{)}^{1/2}},\mathrm{Re}=\frac{{\rho }_{\text{ref}}{u}_{\text{ref}}L}{{\eta }_{\text{ref}}},\text{and}\mathrm{Pr}=\frac{C_{p_{\text{ref}}}{\eta }_{\text{ref}}}{k_{\text{ref}}}.$$

$$\eta =k=T^s,$$

$$\mathbf{\Pi }=\eta [\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\mathsf{T}}]-\frac{2}{3}\eta (\nabla \cdot \mathbf{u})\mathbf{I},\mathbf{q}=-k\nabla T.$$

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\text{inv}}(\mathbf{U})-\nabla \cdot {\mathbf{F}}_{\text{vis}}(\mathbf{U},\mathbf{S})=0,$$

$$\mathbf{S}(\mathbf{U})=T^s\nabla \mathbf{U}.$$

$${\mathbf{U}}_h(\mathbf{x},t)=\sum _{i=1}^K{\mathbf{U}}^i(t){\phi }_i(\mathbf{x}),{\mathbf{S}}_h(\mathbf{x},t)=\sum _{i=1}^K{\mathbf{S}}^i(t){\phi }_i(\mathbf{x}),$$

$$\frac{d}{dt}{\int }_{{\omega }_j}{\mathbf{U}}_h{\phi }_{i}dv-{\int }_{{\omega }_j}({\mathbf{F}}_{\text{inv}}-{\mathbf{F}}_{\text{vis}})\cdot \nabla {\phi }_{i}dv+{\int }_{\partial {\omega }_j}{\phi }_i({\mathbf{F}}_{\text{inv}}-{\mathbf{F}}_{\text{vis}})\cdot \mathbf{n}ds=0$$

$${\int }_{{\omega }_j}{\mathbf{S}}_h{\phi }_{i}dv+{\int }_{{\omega }_j}{T}^s{\mathbf{U}}_h\otimes \nabla {\phi }_{i}dv-{\int }_{\partial {\omega }_j}{T}^s{\phi }_i{\mathbf{U}}_h\otimes \mathbf{n}ds=0,$$

$$\bm{h}_{\text{inv}}(\bm{U}_{h}^{-},\bm{U}_{h}^{+};\bm{n})=\frac{1}{2}\Big{[}\bm{F}_{\text{inv}}(\bm{U}_{h}^{-})+\bm{F}_{\text{inv}}(\bm{U}_{h}^{+})\Big{]}\cdot\bm{n}-\frac{C}{2}(\bm{U}_{h}^{+}-\bm{U}_{h}^{-})$$

$$C=\max\Big{(}\big{|}\bm{u}^{-}\cdot\bm{n}\big{|}+\frac{c^{-}}{\operatorname{\mathit{M\kern-1.49994pta}}},\big{|}\bm{u}^{+}\cdot\bm{n}\big{|}+\frac{c^{+}}{\operatorname{\mathit{M\kern-1.49994pta}}}\Big{)}\,,$$

$$\bm{h}_{\text{vis}}(\bm{U}_{h}^{-},\bm{S}_{h}^{-},\bm{U}_{h}^{+},\bm{S}_{h}^{+};\bm{n})=\frac{1}{2}\Big{[}\bm{F}_{\text{vis}}(\bm{U}_{h}^{-},\bm{S}_{h}^{-})+\bm{F}_{\text{vis}}(\bm{U}_{h}^{+},\bm{S}_{h}^{+})\Big{]}\cdot\bm{n},$$

$$\bm{h}_{\text{aux}}(\bm{U}_{h}^{-},\bm{U}_{h}^{+};\bm{n})=\frac{1}{2}\Big{[}T^{s-}\bm{U}_{h}^{-}+T^{s+}\bm{U}_{h}^{+}\Big{]}\otimes\bm{n}\,.$$

$$\mathsf{L}\frac{d\mathsf{U}}{dt}$$

$$\mathsf{L}\mathsf{S}$$

$$\Delta t=\min \{\Delta t_j\},\Delta t_j=\frac{1}{2m+1}\frac{h_j\mathrm{CFL}}{|u|+c/\mathrm{Ma}+\eta /(\rho h_j\mathrm{Re})},$$

$$\mathrm{Kn}=\frac{\lambda }{L}=\sqrt{\frac{\gamma \pi }{2}}\frac{\mathrm{Ma}}{\mathrm{Re}}.$$

$$\alpha =\{\begin{array}{ll}\frac{\beta p}{1+\beta p}& \text{for monatomic}\\ \frac{\sqrt{\beta p}}{1+\sqrt{\beta p}}& \text{for diatomic},\end{array}$$

$$\beta=\sqrt{\frac{\pi}{32}}\frac{\pi}{c_{i}^{2}}\frac{T_{r}}{T_{w}}\exp\bigg{(}\frac{D_{e}}{RT_{w}}\bigg{)}\frac{1}{p_{r}\operatorname{\mathit{K\kern-1.99997ptn}}},$$

$$\mathbf{u}$$

$$T$$

$$\displaystyle\bm{h}_{\text{inv}}(\bm{U}^{-},\bm{U}^{+};\bm{n})\Big{|}_{\text{w}}$$

$$\displaystyle\bm{h}_{\text{aux}}(\bm{U}^{-},\bm{U}^{+};\bm{n})\Big{|}_{\text{w}}$$

$$\bm{h}_{\text{vis}}(\bm{U}^{-},\bm{S}^{-},\bm{U}^{+},\bm{S}^{+};\bm{n})\Big{|}_{\text{w}}=\bm{F}_{\text{vis}}(\bm{U}_{\text{w}},\bm{S}_{\text{w}})\cdot\bm{n}.$$

$$\bm{U}^{1}\Big{|}_{t=0}=\bm{U}_{\infty},\qquad\bm{U}^{i}\Big{|}_{t=0}=0\quad\text{for }i=2,...,K\,.$$

$$\displaystyle\bm{S}^{i}\Big{|}_{t=0}=0\quad\text{for }i=1,2,...,K.$$

$$P_{n}^{\alpha,\beta}(x)=\frac{(-1)^{n}}{2^{n}n!}(1-x)^{-\alpha}(1+x)^{-\beta}\frac{\operatorname{d}\!^{n}}{\operatorname{d}\!x^{n}}\big{[}(1-x)^{n+\alpha}(1+x)^{n+\beta}\big{]}$$

$${\Theta }_i^a(x)$$

$${\Theta }_{ij}^b(x)$$

$${\Theta }_{ijk}^c(x)$$

$${\mathcal{R}}^2=\{({\eta }_1,{\eta }_2)\in {\mathbb{R}}^2;-1\le {\eta }_1\le 1;-1\le {\eta }_2\le 1\}$$

Full text

1 Introduction

The discontinuous Galerkin method has a number of appealing properties, including accuracy and robustness, when applied to fluid dynamics problems. In industrial applications with complex geometries a key limitation of computational fluid dynamics using the discontinuous Galerkin method is the need for costly mesh generation. This is common with the finite element and finite volume methods which also require boundary-fitted meshes. As also identified in the recent NASA CFD Vision 2030 Study, the generation of boundary fitted meshes is a significant impediment to the implementation of autonomous and reliable computational simulation workflows [1]. Mesh generation, however, becomes straightforward when the element boundaries need not conform to the domain boundary as in immersed methods, also referred to as embedded or fictitious domain methods, see Figure 1. Meshes can still be unstructured and have elements of different sizes in order, for instance, to adapt to critical flow features and boundary layers. In cut-elements intersected by the domain boundary alternative techniques are required to enforce the boundary conditions and evaluate the element integrals. The development of these techniques is greatly aided by the variational structure of the discontinuous Galerkin method.

Immersed methods, which do not require a boundary-fitted mesh, have a long history in computational fluid dynamics, especially in the context of low order finite difference and finite volume methods. These techniques are often restricted to structured Cartesian grids and can represent the boundaries usually only in a diffuse sense, see e.g. [2, 3, 4, 5]. Critically, their accuracy cannot be systematically improved. Hence immersed methods and high-order accurate schemes appeared until recently to be intrinsically incompatible. Contrary to this perception, lately a number of high-order immersed methods have been introduced mostly for the finite element method, see e.g.  [6, 7, 8, 9, 10, 11, 12], and a few for the discontinuous Galerkin [13, 14, 15, 16], finite volume [17, 18, 19] and finite difference methods [20, 21]. Evidently, the techniques for the finite element and discontinuous Galerkin methods must share a number of features due to their common origins. Indeed, they all have a method for representing the domain boundaries on the computational mesh, a method for evaluating integrals on elements partially covered by the fluid and a method for enforcing boundary conditions.

The immersed discontinuous Galerkin method that we have developed builds on earlier work on immersed discontinuous Galerkin and finite element methods. On the computational mesh the domain is described with an implicit signed distance function, similar to that, for instance, in [22, 7, 13, 16]. For integration the cut-elements are triangulated into simplicial elements using the marching triangle and tetrahedra algorithms [23, 8, 24]. After the triangulation, element integrals are evaluated with standard Gauss integration. There are alternative techniques available which do not require a cut-element triangulation, but have other limitations [16, 25, 26]. To achieve high-order accuracy the domain boundaries must be approximated with curved cut-element boundaries of the same polynomial order as the discontinuous Galerkin method basis functions. To this end, we use the curved boundary treatment recently introduced by Zhang [27]. The element-specific domain integrals are evaluated on the simplicial elements with straight boundaries, but the boundary integrals are evaluated on curved boundaries. In our mixed discontinuous Galerkin implementation, according to Bassi and Rebay [28], the imposition of the Dirichlet and Neumann boundary conditions on the cut-elements is straightforward. The fluxes over the curved cut-element boundaries are treated in the same way as in the original discontinuous Galerkin method with boundary fitted meshes.

In most immersed methods numerical accuracy and stability are compromised when only an exceedingly small portion of a cut-element is covered by the fluid. When the discretised equations are solved with an explicit time marching scheme, as in our work, this leads to severe time step size restrictions. In earlier work on immersed discontinuous Galerkin methods element agglomeration strategies have often been used to deal with this problem. Agglomeration techniques are, however, difficult to implement on unstructured meshes and especially in three dimensions. Borrowing from immersed finite elements, we consider in our work two alternative techniques to stabilise the small cut-elements. In the first approach the cut-basis functions are replaced with the extrapolated basis functions from the nearest largest element. This is a consequent extension of the extrapolation technique introduced earlier in connection with B-spline basis functions [29, 7, 8]. In the second approach the basis functions of the cut-element are scaled proportionally to the fraction of the cut-element covered by the solid. This scaling is motivated by a stabilisation technique proposed in the work of Mößner and Reif [30] for B-spline basis functions. As our numerical computations confirm, both approaches provide an easy to implement and effective means of dealing with cut-element stability issues. Alternatively, as recently proposed, the adaptation of the ghost penalty cut-element stabilisation technique known from immersed, or cut, finite elements to the discontinuous Galerkin setting is also feasible [31].

The paper is organised as follows. In Section 2 the mixed modal discontinuous Galerkin method for compressible Navier-Stokes equation is reviewed. Subsequently, in Section 3, the developed high-order immersed discontinuous Galerkin method is introduced. In Section 4 the application of the immersed discontinuous Galerkin method to problems in subsonic and hypersonic flow regimes is demonstrated. The stability of the two proposed cut-element stabilisation methods is numerically evaluated and the results obtained are compared with those in literature and experiments. Finally, conclusions and perspectives are drawn in Section 5.

2 Mixed modal discontinuous Galerkin method for compressible Navier-Stokes

2.1 Governing equations

The non-dimensional compressible Navier-Stokes equations, according to, e.g., Canuto et al. [32], read

[TABLE]

with the conservation vector

[TABLE]

and the inviscid and viscous flux vectors

[TABLE]

where $\rho$ is the mass density, $\bm{u}$ the velocity vector, $p$ the pressure, $\gamma$ the specific heat ratio, $\bm{I}$ the identity tensor, $E$ the total energy density, $\bm{\Pi}$ the deviatoric stress tensor and $\bm{q}$ the heat flux vector.

The dimensionless Mach number $\operatorname{\mathit{M\kern-1.49994pta}}$, Reynolds number $\operatorname{\mathit{R\kern-0.39993pte}}$ and Prandtl number $\operatorname{\mathit{P\kern-0.29999ptr}}$ are defined as

[TABLE]

Here, the reference state variables are the velocity $u_{\text{ref}}$, temperature $T_{\text{ref}}$, density $\rho_{\text{ref}}$, viscosity $\eta_{\text{ref}}$, heat capacity at constant pressure $C_{p_{\text{ref}}}$ and thermal conductivity $k_{\text{ref}}$. Furthermore, $R$ is the universal gas constant and $L$ the characteristic length of the solid. For the viscosity parameter $\eta$ and thermal conductivity $k$, a power law model in non-dimensional form is assumed

[TABLE]

where $s$ is a gas constant associated with the inverse power law of gas molecules. In our examples we take $s=0.78$ for diatomic gas molecules [33]. The constitutive equations for the Navier-Stokes equation are assumed to be

[TABLE]

Finally, the ideal gas equation of state $p=\rho T$ is used in the non-dimensional form.

2.2 Review of the mixed discontinuous Galerkin method

A mixed formulation of the discontinuous Galerkin method was proposed by Bassi and Rebay [28] for the treatment of the viscous terms in the Navier-Stokes equation. To this end, the governing equation (1), which depends on the second order derivatives of the velocity, is reformulated as a first-order coupled system. That is, the auxiliary variable $\bm{S}$ is introduced so that the coupled first order system can be written as

[TABLE]

The discretised field variables $\bm{U}_{h}$ and $\bm{S}_{h}$ are in each element $\omega_{j}$ represented with

[TABLE]

where $\varphi_{i}(\bm{x})$ are the basis functions, and $\bm{U}^{i}$ and $\bm{S}^{i}$ are the corresponding coefficients. The number of basis functions per element $K$ depends on the polynomial degree of approximation, see Section 2.4. In this work the same basis functions $\varphi_{i}$ are used in each local fluid element $\omega_{j}$, even though it is straightforward to relax this restriction. To obtain the discretised weak form of the governing equations (6), the discretised coupled system is element-wise multiplied with the test functions $\varphi_{i}$ and then integrated by parts over each element $\omega_{j}$ yielding

[TABLE]

where $\partial\omega_{j}$ denotes the boundaries of the element $\omega_{j}$. Due to the discontinuous approximation at element boundaries, the flux functions appearing in the boundary integrals of each element are replaced by numerical flux functions. We approximate the inviscid term $\bm{F}_{\text{inv}}\cdot\bm{n}$ with the Lax-Friedrichs flux

[TABLE]

with

[TABLE]

where $c=T^{1/2}$ is the non-dimensional speed of sound. The superscripts $-$ and $+$ denote the variables belonging to the elements on the left and right of an edge. Note that the appearance of $\operatorname{\mathit{M\kern-1.49994pta}}$ in the definition of $C$ is due to the used dimensionless formulation. For other more advanced flux approximation schemes see, e.g., the work of Toro [34]. As a flux limiter, in case of linear basis functions the technique presented in the work of Cockburn et al. [35] is applied direction-by-direction and in case of quadratic basis functions a coupled slope limiter similar to the one in Le et al. [36] is used.

We approximate the viscous flux $\bm{F}_{\text{vis}}\cdot\bm{n}$ with the central flux

[TABLE]

and the term $T^{s}\bm{U}_{h}\otimes\bm{n}$ in the auxiliary equation (8b) with

[TABLE]

All the integrals in the preceding discretised equations are evaluated with Gaussian quadrature. This leads to a semi-discretised system of equations, which can be compactly written as

[TABLE]

where $\bm{\mathsf{L}}$ is the mass matrix and $\bm{\mathsf{U}}$ and $\bm{\mathsf{S}}$ are arrays containing the coefficients of the basis functions. The corresponding residual vectors are denoted with $\bm{\mathsf{R}}_{\bm{\mathsf{U}}}$ and $\bm{\mathsf{R}}_{\bm{\mathsf{S}}}$. We solve this system using a third-order five-step Runge-Kutta time integration [37]. The explicit time integration scheme allows us to decouple the solution of (13) in two steps. Firstly, equation (13b) for the auxiliary unknown $\bm{S}_{h}$ is solved to compute the derivatives $\nabla\bm{U}_{h}$ of the conservative variables. Subsequently, equation (13a) is solved with the already computed $\nabla\bm{U}_{h}$. As will be discussed in Section 2.4, we choose orthogonal basis functions so that the mass matrix $\bm{\mathsf{L}}$ is readily invertible.

The time step size $\Delta t$ is chosen according to Cockburn and Shu [37] with

[TABLE]

where the index $j$ refers to the element $\omega_{j}$ with the characteristic element size $h_{j}$, $m$ is the polynomial degree of the basis functions and $\operatorname{\mathit{C\kern-1.99997ptF\kern-1.99997ptL}}<1$ is the Courant-Friedrichs-Lewy (CFL) number.

2.3 Boundary and initial conditions

The boundary conditions for the compressible Navier-Stokes equations (1) depend critically on the Knudsen number $\operatorname{\mathit{K\kern-1.99997ptn}}$, which is the ratio of molecular mean free path $\lambda$ to the characteristic length of the solid $L$ and can be expressed with

[TABLE]

Especially, for rarefied (large $\lambda$) or microscale gas flows (small $L$) the typical non-slip boundary conditions are not suitable. To cover all Knudsen number regimes, we consider the Langmuir boundary conditions that are derived from the Langmuir adsorption isotherm at the gas and boundary surface interface [38, 39, 40].

The fraction of the boundary surface covered by gas molecules is given by the ratio $\alpha$, with $0\leq\alpha\leq 1$, and can be expressed in dimensional form as

[TABLE]

where $p$ is the pressure at the boundary and $\beta$ is an equilibrium constant that describes the interaction between the boundary and gas molecules. The equilibrium constant may be derived as

[TABLE]

where $c_{i}$ is the exponent of the inverse power law for the gas particle interaction potential, $p_{r}$ the reference pressure, $T_{r}$ the reference temperature, $T_{w}$ the wall temperature, $D_{e}$ the heat of adsorption and $R$ the universal gas constant. The two constitutive constants in this model for air and nitrogen are $c_{i}=1.1908$ and $D_{e}=5255\;J/mol$.

The boundary conditions for the velocity and temperature are given in dependence of the ratio $\alpha$ as follows

[TABLE]

where $\bm{u}_{w}$ and $\bm{u}_{g}$ denote wall velocity and far-field gas velocity respectively; and, $T_{g}$ is far-field gas temperature. As can be inferred from (16) and (17), the ratio $\alpha$ approaches $\alpha=1$ with the decrease in Knudsen number $\operatorname{\mathit{K\kern-1.99997ptn}}$ so that the boundary condition according (18) becomes no-slip.

In the mixed discontinuous Galerkin method with explicit time marching it is straightforward to apply the wall boundary conditions [28]. The inviscid and auxiliary fluxes in (9) and (12) are computed with the slip velocity (18a), the temperature (18b) and other variables taken from inside the fluid, that is,

[TABLE]

The viscous flux in (11) is evaluated using the value of the unknown auxiliary variable $\bm{S}$,

[TABLE]

The initial condition for the conservation variables throughout the domain is specified using values at the far-field. The coefficients of the basis functions in the interpolation equation (7) are initialised as follows

[TABLE]

Note that for the considered modal basis functions we have $\varphi_{1}(\bm{x})=1$. The components of auxiliary unknown are set to zero,

[TABLE]

2.4 Modal basis functions

We discretise the two field variables $\bm{U}$ and $\bm{S}$, c.f. (6), using the modal Dubiner basis functions [41], which are derived from the Jacobi polynomial

[TABLE]

with $\{\alpha,\beta\}\in\mathbb{N}_{0}$, not to be confused with the same symbols used in Section 2.3. In addition, the following three functions are defined

[TABLE]

In two-dimensions, as illustrated in Figure 2, we firstly construct the basis in a reference square

[TABLE]

and then map it to the reference triangle

[TABLE]

To this end, the linear mapping

[TABLE]

is used, which collapses one of the edges of the square. The orthogonal basis on $\mathcal{T}^{2}$ is given by

[TABLE]

Similarly, in three-dimensions as shown in Figure 3 the basis functions are first constructed in a reference hexahedron

[TABLE]

and then mapped into a reference tetrahedron

[TABLE]

The orthogonal basis functions in the reference tetrahedron are given by

[TABLE]

with the mappings

[TABLE]

3 Immersed discontinuous Galerkin method on non-boundary-fitted meshes

3.1 Implicit geometry description

To begin with, we assume that the embedded solid domain $\Omega$ is described with a scalar-valued implicit (or, level set) function

[TABLE]

where $\operatorname{dist}(\bm{x},\Gamma)=\min_{\bm{y}\in\Gamma}|\bm{x}-\bm{y}|$ is the shortest distance between the point $\bm{x}$ and the fluid-solid interface $\Gamma$. By definition, the level set function $\phi(\bm{x})$ is negative inside the solid domain $\Omega$ and positive inside the fluid domain. The unit normal to the boundary is given by the normalised gradient of the level set function

[TABLE]

where $\bm{x}\in\Gamma$ is a point on the boundary.

In problems where an implicit function $\phi(\bm{x})$ is not available a sufficiently accurate signed distance function needs to be computed. This can be accomplished either with numerical signed distance computations [42, 43, 8] or analytic implicitisation techniques [44]. The use of an implicit instead of a parametric representation greatly simplifies all geometric operations. For instance, the classification of elements as solid, fluid or cut involves only the evaluation of the signed distance function $\phi(\bm{x})$, see Figure 4. In our present implementation elements are classified according to the signed distance function value at their nodes. In the discontinuous Galerkin computations only the fluid and cut-elements are considered.

3.2 Cut-element integration

For elements entirely inside the fluid domain the element and boundary integrals appearing in the weak form of the governing equations (8) are evaluated with Gaussian quadrature, see Section 4. The cut-elements, which are only partly covered by the fluid domain, are first triangulated before evaluating the domain and boundary integrals. Note that this triangulation is only for integration. The field variables are still approximated by the modal basis functions defined over the entire cut-element.

In our implementation the cut-elements are triangulated according to Rüberg and Cirak [8] using marching triangle in 2D and marching tetrahedra algorithms in 3D. In 2D there are $2^{3}$ combinations for the sign of the level set function at the three nodes of the triangle. Considering symmetry these reduce to the two canonical cases shown in Figure 5. In 3D, there are $2^{4}$ combinations for the sign of the level set function, which reduce to the three canonical cases shown in Figure 6. In evaluating the element and flux integrals, we use in each sub-triangle or sub-tetrahedron the same number of Gauss integration points as indicated in Section 4 for the elements fully covered by the fluid. Hence, the total number of integration points for cut-elements is usually higher than for other elements. This does however not lead to a significant time overhead because only very few elements in a mesh are cut-elements. Moreover, while the elements are classified as fluid, solid or cut according to the signed distance function value at their nodes, the intersection of the boundary with the level set function is found by bisectioning.

3.3 Cut-element stabilisation

When an exceedingly small portion of a cut-element is covered by the fluid domain, numerical accuracy and stability are usually compromised. The diagonal entries of the mass matrix $\bm{\mathsf{L}}$ in the system of equations (13) are proportional to the physically active support size of the corresponding basis functions. Consequently, the condition number of the matrix $\bm{\mathsf{L}}$ can become arbitrarily large if there are cut-elements with only a negligible part covered by the fluid domain. This results in significantly smaller stable time step sizes $\Delta t$ when an explicit time marching scheme is used. Although no estimates exist for such cut-elements, it can be argued that the characteristic element size $h_{j}$ in the estimate (14) has to be chosen in accordance with the physically active support size of each basis function, i.e. the area in 2D or the volume in 3D of the triangulated part of the cut-element.

We propose two alternative stabilisation techniques to overcome the excessive stable time restriction imposed by cut-elements. Both stabilisation techniques can also be interpreted as preconditioning of the system of discretised discontinuous Galerkin equations (13). To decide which cut-elements are susceptible to stability issues, the following coverage ratio is computed

[TABLE]

where $\text{Vol}(\omega_{j}^{F})$ is the size, area in 2D or volume in 3D, of the triangulated part of the cut-element and the $\text{Vol}(\omega_{j})$ is the size of the entire cut-element. Although, in principle, all basis functions belonging to cut-elements with $f_{j}<1$ may lead to stability issues, we process in our implementation for numerical efficiency reasons only cut-elements with $f_{j}<0.3$.

In the first proposed stabilisation approach, we replace the basis functions in the critical cut-element with the extrapolated basis functions of the neighbouring fluid element of largest size, see Figure 7. The neighbourhood of a cut-element element $\omega_{j}$ is defined as the set of fluid elements $\{\omega_{i}\}$ which share a common node with $\omega_{j}$. This approach is motivated by the extended B-splines proposed in [29]. It ensures that the polynomial reproduction property of the basis functions is automatically retained so that optimal convergence rates can be achieved irrespective of the extension process. In effect, this technique enlarges the size of the element from which the basis functions are extrapolated akin to classical element agglomeration strategies. The basis functions with a very small physically active support size and their coefficients do not appear in the system of discretised discontinuous Galerkin equations (13). In a cut-element $\omega_{j}$ with the largest neighbouring fluid element $\omega_{j}^{\ast}$, see Figure 7, the interpolation equations are modified as follows

[TABLE]

In terms of implementation, if the coverage ratio of a cut-element $\omega_{j}$ is below $f_{j}<0.3$, the discretised conservation vector $\bm{U}_{h}$ and auxiliary vector $\bm{S}_{h}$ on $\omega_{j}$ are extrapolated from the basis functions and coefficients of the element $\omega_{j}^{*}$. According to (35), this is accomplished by simply evaluating the interpolation equations for the element $\omega_{j}^{\ast}$ at $\bm{x}\in\omega_{j}$. Similarly, the test functions in the weak form (8) are replaced with the basis functions of the element $\omega_{j}^{*}$. In computing the fluxes (9), (11) and (12) the extrapolated field variables are used. The explicit time integration of the stabilised system of equations  (13) proceeds as usual.

In the second proposed stabilisation approach, as introduced in the work of Mößner and Reif [30] for immersed B-spline finite elements, the cut-basis functions are rescaled in proportion to the ratio $1/f_{j}$. This leads in a cut-element $\omega_{j}$ to the following modification of the interpolation equations

[TABLE]

The implementation of this stabilisation approach is slightly easier than the extrapolation approach. In the field discretisation (7) and the weak form (8) the basis functions are replaced with the scaled basis functions $\varphi_{i}(\bm{x})/f_{j}$. This leads to a scaling of the mass matrix $\bm{\mathsf{L}}$ and the coefficient arrays $\bm{\mathsf{U}}$ and $\bm{\mathsf{S}}$ in the semi-discretised system of equations (13), but not of the right-handside vectors $\bm{\mathsf{R}}_{\bm{\mathsf{U}}}$ and $\bm{\mathsf{R}}_{\bm{\mathsf{S}}}$. Again, the explicit time integration of the stabilised system of equations (13) proceeds as usual.

3.4 High-order boundary reconstruction

The cut-element integration using (straight-sided) simplices leads to large errors when the boundary is highly curved. Even when high-order basis functions are used, the domain approximation errors will dominate resulting in a low-order discretisation scheme [45]. In practical terms, an additional complication concerns the interpretation of the boundary conditions which are usually prescribed on the curved boundary.

It is straightforward to achieve high-order accuracy by suitably curving the boundaries of the simplicial integration elements introduced in Section 3.2 and evaluating all element integrals on the curved simplices. As shown in the work of Zhang [27], high-order accuracy can also be achieved by evaluating the domain integrals on straight-sided simplicial elements and evaluating only the boundary, i.e. flux, integrals along curved boundaries, see Figure 8. That is, only the integration of the boundary integrals needs special treatment in case of curved domain boundaries. According to [27], for nonconvex domains it is assumed that the solution can be smoothly extended outside the domain and that the distance between the straight and the curved boundary is sufficiently small for the sake of stability. Aside from these restrictions, to obtain an optimally convergent discretisation scheme the curved domain boundary has to be approximated with the same polynomial order like the basis functions used for field interpolation, i.e. (7). To this end, we modify the boundaries of the integration simplices abutting the domain boundary as follows. Firstly, additional nodes are introduced on the straight-sided simplices. Subsequently, all boundary nodes are projected to the prescribed domain boundary with the help of the gradient of the signed distance function (33). The gradient of the signed distance function is well defined close to a smooth domain boundary. When evaluating the boundary integrals the nonlinear mapping between the straight and curved element boundaries has to be taken into account. To establish this, the curved boundary is assumed to be parametrised with Lagrange polynomials. Throughout the mapping the parametric coordinates remain the same while their physical coordinates change.

As an example, in Figure 8 the process of approximating the domain boundary geometry with quadratic polynomials ($P^{2}$) is illustrated. First additional mid-nodes are introduced on the edges of the integration simplices. Subsequently, the introduced mid-nodes are pushed along the gradient of the signed distance function $\nabla\phi$ until the boundary is reached. The intersection of the line given by the gradient with the boundary $\phi(\bm{x})=0$ is found by bisectioning.

4 Results and Discussion

The theoretical convergence properties and accuracy of the discontinuous Galerkin method are well understood for hyperbolic and elliptic systems. However, the verification and validation of the discontinuous Galerkin solutions for complex engineering flow problems remain challenging. We present several examples of increasing complexity which are all discretised with unstructured triangular or tetrahedral meshes. As basis functions linear ($P^{1}$) or quadratic ($P^{2}$) Dubiner basis functions are used. The first set of examples verifies the convergence rate of the method and demonstrates the efficacy of the two implemented cut-element stabilisation techniques. In the later more complex examples, our results are both verified and validated mostly by comparison with results from the literature. The last example concerning the hypersonic flow around a vehicle is validated by our experimental results. The examples included cover the subsonic and hypersonic flow regimes. For stationary flows the discontinuous Galerkin equations are advanced in time with explicit time marching until a stationary state is reached. We choose $\operatorname{\mathit{C\kern-1.99997ptF\kern-1.99997ptL}}=0.8$ for subsonic flows and $\operatorname{\mathit{C\kern-1.99997ptF\kern-1.99997ptL}}=0.1$ for hypersonic flows. The stationary solution is assumed to be reached, depending on the example, either when the increment of the density or drag coefficient becomes smaller than $10^{-7}$.

The surface and volume integrals in the weak form (8) are evaluated numerically. In the present work, the Gauss-Legendre quadrature rule is employed within the elements and element interfaces. In two dimensions, we use for the element interfaces 3 Gauss points for linear basis functions ($P^{1}$), and 5 Gauss points for quadratic basis functions ($P^{2}$). For the element integrals we use 9 and 25 Gauss points for $P^{1}$ and $P^{2}$, respectively. In three dimensions, boundary integrals are evaluated with 9 and 25 Gauss points for $P^{1}$ and $P^{2}$, and volume integrals are computed with 11 and 21 Gauss points for $P^{1}$ and $P^{2}$, respectively.

4.1 Flow around a sphere and convergence

We verify the convergence order of the proposed approach with an inviscid flow around a sphere at $Ma=0.2$. The sphere has a radius 1 and is placed at the centre of a cuboidal domain with an edge length of 10. The sphere is defined with the signed distance function

[TABLE]

The cuboidal domain is meshed in turn with $64^{3}$, $96^{3}$, $128^{3}$, $160^{3}$ and $192^{3}$ hexahedrons each of which is later split into six tetrahedra.

The flow is isentropic since it is inviscid and no shocks are present. The numerically computed entropy differs from the free-stream entropy and converges with increasing mesh refinement. The error in the computed entropy is defined as

[TABLE]

where $\Delta s=p/\rho^{\gamma}-s_{\infty}$ and the free-stream entropy is $s_{\infty}=p_{\infty}/\rho^{\gamma}_{\infty}=1$. The convergence of the entropy error with increasing mesh refinement is shown in Figure 9. As is to be expected, the entropy error converges with increasing mesh refinement approximately at a rate of three for quadratic ($P^{2}$) basis functions and a rate of two for linear ($P^{1}$) basis functions. The cut-element stabilisation by extension gives slightly larger errors than by scaling of the basis functions. This may be due to the slight increase of the characteristic element size in the extension approach.

4.2 Shock trains and cut-element stability

In this example the axisymmetric, two-dimensional compressible Navier-Stokes equations are solved to study the shock train formation inside a confined duct, see Figures 10 and 11. The conditions for the flow are $\operatorname{\mathit{M\kern-1.49994pta}}=5.0$ and $\operatorname{\mathit{R\kern-0.39993pte}}=800$. The Reynolds number is defined with respect to the radius of the solid cone placed inside the duct. The domain is discretised with a structured mesh consisting of 8000 triangular elements of uniform size and quadratic ($P^{2}$) Dubiner basis functions. To study the influence of the proposed cut-element stabilisation techniques the discretisation mesh has been slightly displaced to obtain successively smaller minimum coverage ratios $f_{j}\in\{1/2,\,1/8,\,1/32,\,1/128\}$.

Stationary shock trains consisting of a series of successive shock and rarefaction waves are frequently used for verification and validation of high-order schemes  [46]. Although such problems have been extensively studied, see e.g. [47, 48], they provide valuable information about the capability of the proposed immersed discontinuous Galerkin method in dealing with intricate gas dynamics problems. The structure of the shock train can be complex even for simple inflow and wall conditions and the flow remains hypersonic throughout the duct. For instance, in a symmetric duct with uniform, steady inflow conditions, asymmetric flow structures may develop and self-excited oscillations may appear [48, 49]. These features are not present in our simplified axisymmetric, two-dimensional set-up. The shock trains are usually classified as either normal or oblique shock trains depending on the presence or absence of bifurcated normal shock waves [47]. Lower upstream Mach numbers typically lead to normal shock trains and higher Mach numbers, as in our computations, lead to oblique shock waves.

The pressure isocontours of the present computations are shown in Figures 10 and 11 for the scaled and extended basis stabilisation techniques, respectively. As the gas flow impinges on the cone, first a planar oblique shock is generated. After this shock wave impinges on the duct wall, a reflected shock wave is generated. This process continues and a series of reflected shock waves are created. The results for four sets of meshes and two cut-element stabilisation techniques show that the size, structure and locations of the incident and reflected shock waves are visually identical.

4.3 Flow around a two-dimensional cylinder

There is extensive literature on the study of sub- and supersonic flows around a cylinder, which makes it an ideal application to verify and validate the developed immersed discontinuous Galerkin method. In all the presented computations, a cylinder with diameter $d=1$ is placed at the centre of a fluid domain with radius $30d$. For meshing, $100$ nodes are placed along the circumference of the outer fluid domain boundary and across the diameter, on average, $300$ nodes are present. In the computations the boundary-fitted and non-boundary-fitted meshes depicted in Figures 12a and 12b are used. Each of the meshes consists of about $60000$ triangular elements. As basis functions the quadratic ($P^{2}$) Dubiner basis functions are used. Along the circumference of the outer fluid domain boundary far-field boundary conditions are prescribed.

4.3.1 Subsonic flow

For very low Reynolds numbers, with $\operatorname{\mathit{R\kern-0.39993pte}}\ll 1$ (based on cylinder’s diameter), the flow is symmetrical upstream and downstream of the cylinder. As Reynolds number increases, the upstream-downstream symmetry breaks. When the Reynolds number exceeds $\operatorname{\mathit{R\kern-0.39993pte}}>4$ two attached eddies emerge behind the cylinder, see e.g. [50, 51, 52]. These eddies increase in size with increasing Reynolds numbers before the flow becomes nonstationary at about $\operatorname{\mathit{R\kern-0.39993pte}}\gtrsim 46$ and periodic vortex shedding begins [53].

In the first set of computations, we use the non-boundary-fitted mesh depicted in Figure 12b to verify and validate the accuracy of the developed techniques. This mesh contains $62$ cut-elements and nine of them require stabilisation. In Figure 13 the computed average drag coefficients $C_{d}$ for Reynolds numbers up to $\operatorname{\mathit{R\kern-0.39993pte}}=500$ are plotted. The obtained values are in close agreement with the computational results reported by Emblemsvag et al. [54] and the experimental results reported by Wieselsberger [50] and Tritton [51].

In addition to the average drag coefficient, the wake separation length (for $\operatorname{\mathit{R\kern-0.39993pte}}\lesssim 40$) is also a feature of the flow that is often compared. The wake separation length is defined as the distance from the rear of the cylinder to the tail end of the wake. Table 1 shows that our results are in very good agreement with computations by others. Moreover, the introduced two cut-element stabilisation techniques lead to almost the same results.

In the second set of computations, we compare the vorticity fields obtained with the boundary-fitted and non-boundary fitted meshes. The vorticity field is highly sensitive to the accurate resolution of the velocity gradients in the vicinity of the cylinder boundary. In Figure 14 the vorticity isocontours for the boundary-fitted and non-boundary fitted meshes and three different Reynolds numbers $\operatorname{\mathit{R\kern-0.39993pte}}=20$, $\operatorname{\mathit{R\kern-0.39993pte}}=40$ and $\operatorname{\mathit{R\kern-0.39993pte}}=80$ are plotted. For the two lower Reynolds numbers $\operatorname{\mathit{R\kern-0.39993pte}}=20$ and $\operatorname{\mathit{R\kern-0.39993pte}}=40$ the flow is stationary. As evident from the vorticity isocontours, the boundary-fitted and non-boundary fitted meshes yield very similar results. For $\operatorname{\mathit{R\kern-0.39993pte}}=20$ the three isocontour plots are essentially indistinguishable. For $\operatorname{\mathit{R\kern-0.39993pte}}=40$ there is some difference in the back of the cylinder, which appears to be due to the slightly larger recirculation velocity in case of the immersed method. For the instationary $\operatorname{\mathit{R\kern-0.39993pte}}=80$ the essential features of the vorticity field are the same while there are some minor differences in the wake flow. In making these comparisons it has to be noted that the element sizes and connectivity of the boundary-fitted and non-boundary-fitted meshes are, especially close to the cylinder, different.

4.3.2 Hypersonic flow

Next, the hypersonic flow around a two-dimensional cylinder is considered to validate the proposed immersed discontinuous Galerkin method. Shock detection and capturing in high-order discontinuous Galerkin methods is challenging and remains an active area of research. For finite difference and finite volume methods many ingenious shock capturing methods are known. Although the extension of these methods to the one-dimensional discontinuous Galerkin method works well, their extension to two- and three-dimensional unstructured meshes is not always satisfactory.

We consider the two experiments reported by Maslach et al. [57]. Between the two experiments only the Knudsen number $\operatorname{\mathit{K\kern-1.99997ptn}}=0.016$ and $\operatorname{\mathit{K\kern-1.99997ptn}}=0.08$ is different. The other flow parameters are $\operatorname{\mathit{M\kern-1.49994pta}}=5.92$, $\operatorname{\mathit{P\kern-0.29999ptr}}=0.707$, $T_{\infty}=36.6K$ and $T_{w}=293K$. The dimensionless parameters correspond to a physical problem with cylinder radius $R=1.55\,mm$ and $R=0.31\,mm$ with air as the working gas.

In Figure 15 the computed pressure isocontours for the two experiments are shown. Both computations have been in turn computed with the scaled and extended basis stabilisation techniques. As evident from Figure 15 the two alternative cut-element stabilisation techniques yield indistinguishable results in terms of bow shock standoff distance and pressure contours in the proximity of the cylinder. In Table 2 the computed drag coefficients are reported which are in excellent agreement with the experiments, with the scaled basis stabilisation giving slightly better results.

4.4 Hypersonic gas flow around the Apollo 6 command module

As a more complex engineering problem, the three-dimensional gas flow around the Apollo 6 command module [58] is investigated. The assumed freestream flow and wall conditions are $\operatorname{\mathit{M\kern-1.49994pta}}=5.0$, $\operatorname{\mathit{R\kern-0.39993pte}}=1500$, $T_{0}=288\,K$ and $T_{w}=288\,K$. The angle of attack is $25^{\circ}$. The axisymmetric command module geometry and the unstructured tetrahedral discretisation mesh are shown in Figure 16. The command module is placed in the centre of a spherical fluid domain with a diameter of $30\,d$, where $d$ is the maximum diameter of the module. The mesh consists of 386600 elements in total and across the diameter of the fluid domain, on average, 720 elements are present. As basis functions the quadratic ($P^{2}$) Dubiner basis functions are used and the total number of degrees of freedom is approximately 19.3 million. The surface mesh generated during cut-element integration consists of 54576 triangular elements. The computed pressure isocontours are shown in Figure 17 with the clearly captured shockwave. Unfortunately, we do not consider the reacting flow so that it is not possible to compare our computational results with experimental data.

4.5 Hypersonic gas flow around a complex vehicle

As the final example, the complex three-dimensional gas flow around a hypersonic aerospace vehicle is considered, see Figure 18. The freestream flow and wall conditions are $\operatorname{\mathit{M\kern-1.49994pta}}=4.98,\operatorname{\mathit{R\kern-0.39993pte}}=1477,\operatorname{\mathit{K\kern-1.99997ptn}}=0.005,T_{\infty}=184.6K,T_{w}=1100K$. The hypersonic vehicle is placed in the centre of a spherical fluid domain with a diameter of $30\,L$, where $L$ is the maximum length of the vehicle. In the computations, five meshes of decreasing element size are used, namely with 2105672, 5201671, 8519752 and 11007986 elements. In all computations the quadratic (P2) Dubiner basis functions are used. The finest mesh results in about 550 million degrees of freedom. For this example the boundaries have only been resolved linearly because our current high-order boundary approximation does not allow for sharp features, such as corners and edges. The computed pressure distribution on the vehicle surface and the cross-section of the flow are shown in Figure 19 for the coarsest mesh with 2105672 elements. The pressure coefficient $C_{p}$ distribution along the symmetry axis on the bottom surface of the vehicle is plotted in Figure 20. In the same figure, the experimentally obtained pressure distribution is also shown. The experiment and the computations exhibit an excellent agreement.

5 Conclusions

We introduced a robust immersed discontinuous Galerkin method for solving the compressible Navier-Stokes equations on unstructured non-boundary fitted meshes. Despite their many advantages, discontinuous Galerkin methods suffer from their reliance on boundary-fitted meshes. This requirement has become a major bottleneck towards the realisation of automated design and simulation workflows in industrial applications with complex geometries and moving boundaries. As discussed throughout the paper, in discontinuous Galerkin methods there is no intrinsic need for boundary-fitted meshes. It is straightforward to apply the discontinuous Galerkin method on non-boundary fitted meshes by suitably adapting the evaluation of the element integrals and boundary fluxes in cut-elements. The stability or conditioning issues arising from cut-elements partially covered by the fluid can be effectively sidestepped with the introduced two cut-element stabilisation techniques. As our numerical studies confirm the proposed cut-element stabilisation techniques also work well in combination with the employed limiter for the fluxes. There are no spurious oscillations present close to the boundaries and there is no need to reduce the time step size in cut-elements. In addition, the comparison of the computed flow parameters with experiments and previous computations shows good agreement. Amongst the two investigated cut-element stabilisation techniques, especially the scaling of the basis functions is straightforward to implement.

Our present high-order boundary reconstruction approach relies on implicitly defined smooth geometries. Most industrial CAD geometries consist however of trimmed NURBS (Non-Uniform Rational B-Spline) patches or other emerging geometry representations, like subdivision surfaces [59, 60]. In addition, industrial geometries usually have many sharp features in form of corners and edges and contain many small geometry features, including chamfers, fillets and holes, that cannot be realistically resolved with a discretisation mesh. There is already extensive amount of work in computer graphics on the conversion between different geometry representations, the treatment of sharp features and multiresolution geometry representations, see e.g. [61, 62]. The application of these techniques in immersed discretisation methods aiming an achievable trade-off between accuracy and robustness is presently an active area of research [24, 63].

Acknowledgement

This study is partly supported by the Ministry of Education of the People’s Republic of China as part of the China 111 project under grant number of B17037.