An Immersed Discontinuous Galerkin Method for Compressible Navier-Stokes Equations on Unstructured Meshes
Hong Xiao Eky Febrianto Qiaoling Zhang Fehmi Cirak


H Xiao ET AL.Immersed discontinuous galerkin method on unstructured meshes






        1Department of Engineering, University of Cambridge, Cambridge, CB2 1PZ, UK 
 
        2School of Power and Energy, Northwestern Polytechnical University, Xi'an, 710072, China
    

f.cirak@eng.cam.ac.uk




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. 







An Immersed Discontinuous Galerkin Method for Compressible Navier-Stokes Equations on Unstructured Meshes
    Hong Xiao1,2, Eky Febrianto1, Qiaoling Zhang1, Fehmi Cirak1 
    
=========================================================================================================






§ 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 <cit.>.  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 <ref>. 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. <cit.>. 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.  <cit.>, and a few for the discontinuous Galerkin <cit.>, finite volume <cit.> and finite difference methods <cit.>.  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 <cit.>.  For integration the cut-elements  are triangulated into simplicial elements using the marching triangle and tetrahedra algorithms <cit.>. 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 <cit.>.  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 <cit.>. 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 <cit.>, 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 <cit.>. 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 <cit.> 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 <cit.>.

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.





§ MIXED MODAL DISCONTINUOUS GALERKIN METHOD FOR COMPRESSIBLE NAVIER-STOKES 






 §.§ Governing equations



The non-dimensional compressible Navier-Stokes equations, according to, e.g., Canuto et al. <cit.>, read


    ∂U⃗/∂t + ∇·F⃗_inv(U⃗)
    -∇·F⃗_vis(U⃗,∇U⃗)= 0⃗


with the conservation vector 


    U⃗  = [   ρ; ρu⃗; ρ E ]


and the inviscid and viscous flux vectors



    F⃗_inv(U⃗)     =
    [             ρu⃗; ρu⃗⊗u⃗+1/γ^2pI⃗;  (ρ E+1/γ^2p)u⃗ ]
    F⃗_vis(U⃗,∇U⃗)    = 1/[                   0;                  Π⃗; Π⃗u⃗-1/(γ-1) ^2Prq⃗ ]  ,



where ρ is the mass density, u⃗ the velocity vector, p the pressure, γ the specific heat ratio, I⃗ the identity tensor, E the total energy density, Π⃗ the deviatoric stress tensor and q⃗ the heat flux vector. 

The dimensionless Mach number , Reynolds number  and Prandtl number  are defined as


    = u_ref/(γ RT_ref)^1/2,    =ρ_ref u_ref L/η_ref,   and   = C_p_refη_ref/k_ref.


Here, the reference state variables are the velocity u_ref, temperature   T_ref, density ρ_ref, viscosity η_ref, heat capacity  at constant pressure C_p_ref and thermal conductivity k_ref. Furthermore, R is the universal gas constant and  L the characteristic length of the solid. For the viscosity parameter η and thermal conductivity k, a  power law model in non-dimensional form is assumed


    η=k=T^s   ,


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 <cit.>.  The constitutive equations for the Navier-Stokes equation are assumed to be





    Π⃗ = η[∇u⃗ + (∇u⃗)^] - 2/3η (∇·u⃗) I⃗, q⃗ =-k∇ T   .



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




 §.§ Review of the mixed discontinuous Galerkin method



A mixed formulation of the discontinuous Galerkin method was proposed by Bassi and Rebay <cit.> for the treatment of the viscous terms in the Navier-Stokes equation. 
To this end, the governing equation (<ref>), which depends on the second order derivatives of the velocity,  is reformulated as a first-order coupled system.  That is, the auxiliary variable S⃗ is introduced so that the coupled first order system  can be written as



    ∂U⃗/∂t+∇·F⃗_inv(U⃗)
    -∇·F⃗_vis(U⃗,S⃗) =0, 
       S⃗(U⃗) = T^s ∇U⃗.
 


The discretised field variables U⃗_h and S⃗_h are in each element ω_j represented with


    U⃗_h(x,t)=∑_i=1^K U⃗^i(t)φ_i(x),   S⃗_h(x,t)=∑_i=1^K S⃗^i(t)φ_i(x)   ,


where φ_i(x) are the basis functions, and U⃗^i and S⃗^i are the corresponding coefficients. The number of basis functions per element K depends on the polynomial degree of approximation, see Section <ref>. In this work the same basis functions φ_i are used in each local fluid element ω_j, even though it is straightforward to relax this restriction. To obtain the discretised weak form of the governing equations (<ref>), the discretised coupled system is element-wise multiplied with the test functions φ_i and then integrated by parts over each element ω_j yielding



    d/d t∫_ω_jU⃗_h φ_i  v
    -∫_ω_j (  F⃗_inv -  F⃗_vis) ·∇φ_i   v
    +∫_∂ω_jφ_i ( F⃗_inv  -  F⃗_vis)  ·n⃗ s = 0⃗
       ∫_ω_jS⃗_h φ_i  v
    +∫_ω_j T^s U⃗_h ⊗∇φ_i   v
    -∫_∂ω_j T^s φ_i U⃗_h  ⊗n⃗ s=0⃗  ,
  


where ∂ω_j denotes the boundaries of the element ω_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 F⃗_inv·n⃗ with the Lax-Friedrichs flux 


    h⃗_inv(U⃗_h^-,U⃗_h^+ ; n⃗)=
    1/2[F⃗_inv(U⃗_h^-)+F⃗_inv(U⃗_h^+) ] ·n⃗  - C/2(U⃗_h^+-U⃗_h^-)


with  


    C = max(| u⃗^- ·n⃗|+c^-/,| u⃗^+ ·n⃗|+c^+/)   ,


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  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 <cit.>. As a flux limiter, in case of linear basis functions the technique presented in the work of Cockburn et al. <cit.> is applied direction-by-direction and in case of quadratic basis functions a coupled slope limiter similar to the one in Le et al. <cit.> is used.

We approximate the viscous flux F⃗_vis·n⃗ with the central flux 


    h⃗_vis(U⃗_h^-,S⃗_h^-,U⃗_h^+,S⃗_h^+; n⃗) = 
    1/2[F⃗_vis(U⃗_h^-,S⃗_h^-)+
    F⃗_vis(U⃗_h^+,S⃗_h^+)]·n⃗,


and the term T^sU⃗_h  ⊗n⃗ in the auxiliary equation (<ref>) with


    h⃗_aux(U⃗_h^-,U⃗_h^+; n⃗) = 
    1/2[T^s-U⃗_h^-+
    T^s+U⃗_h^+] ⊗n⃗  .


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

 

    LU/ t   =R_ U(U), 
    LS   = R_ S(U)   ,
 


where L is the mass matrix and U and S are arrays containing the coefficients of the basis functions.  The corresponding residual vectors are denoted with R_ U and R_ S. We solve this system using a third-order five-step Runge-Kutta time integration <cit.>.  The explicit time integration scheme allows us to decouple the solution of (<ref>) in two steps. Firstly, equation (<ref>) for the auxiliary unknown S⃗_h  is solved to compute the derivatives ∇U⃗_h of the conservative variables. Subsequently, equation (<ref>) is solved with the already computed  ∇U⃗_h.  As will be discussed in Section <ref>, we choose orthogonal basis functions so that the mass matrix L is readily invertible.  

The time step size Δ t is chosen according to Cockburn and Shu <cit.> with 


    Δ t=min{Δ t_j },   Δ t_j=1/2m+1 h_j  /|u|+c /+η/(ρ h_j )  ,


where the index j refers to the element ω_j with the characteristic element size h_j, m is the polynomial degree of the basis functions and < 1 is the Courant-Friedrichs-Lewy (CFL) number. 




 §.§ Boundary and initial conditions 



The boundary conditions for the compressible Navier-Stokes equations (<ref>) depend critically on the Knudsen number , which is the ratio of molecular mean free path λ  to the characteristic length of the solid L and can be expressed with 


    = λ/L = √(γπ/2)/  .


Especially, for rarefied (large λ) 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  <cit.>. 

The fraction of the boundary surface covered by gas molecules is given by the ratio α, with  and can be expressed in dimensional form as 


    α=
    β p1+β p   for  monatomic
    √(β p)1+√(β p)   for  diatomic,


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


    β=√(π/32)π/c_i^2T_r/T_wexp(D_e/R T_w)1/p_r ,


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 α as follows  

 

    u⃗    =αu⃗_w+(1-α)u⃗_g 
    
    T    =α T_w+(1-α) T_g,
 


where u⃗_w and u⃗_g denote  wall velocity and far-field gas velocity respectively; and, T_g is far-field gas temperature. As can be inferred from (<ref>) and (<ref>), the ratio α approaches  α = 1 with the decrease in Knudsen number  so that the boundary condition according (<ref>) becomes no-slip.

In the mixed discontinuous Galerkin method with explicit time marching  it is straightforward to apply the wall boundary conditions <cit.>. The inviscid and auxiliary fluxes in (<ref>) and (<ref>) are computed with the  slip velocity (<ref>), the temperature (<ref>) and other variables taken from inside the fluid, that is, 



    h⃗_inv(U⃗^-,U⃗^+; n⃗)|_w   =
    F⃗_inv(U⃗ |_w)·n⃗  ,
    
    h⃗_aux(U⃗^-,U⃗^+;n⃗)|_w   =
    T^s |_wU⃗ |_w⊗n⃗  .



The viscous flux in  (<ref>) is evaluated using the value of the unknown auxiliary variable S⃗,


    h⃗_vis(U⃗^-,S⃗^-,U⃗^+,S⃗^+; n⃗)|_w=
    F⃗_vis(U⃗_w,S⃗_w)·n⃗.


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 (<ref>) are initialised as follows


    U⃗^1 |_t=0 = U⃗_∞,     U⃗^i|_t=0 = 0   for  i = 2,...,K   .


Note that for the considered modal basis functions we have φ_1(x⃗)=1. The components of auxiliary unknown are set to zero, 


    S⃗^i|_t=0 = 0    for  i=1,2,...,K.





 §.§ Modal basis functions 



We discretise the two field variables U⃗ and S⃗, c.f. (<ref>), using the modal Dubiner basis functions <cit.>, which are derived from the Jacobi polynomial


    P_n^α,β(x)=(-1)^n/2^nn!(1-x)^-α(1+x)^-β^n/ x^n[(1-x)^n+α(1+x)^n+β]


with {α, β}∈ℕ_0, not to be confused with the same symbols used in Section <ref>. 

In addition, the following three functions are defined 



    Θ_i^a(x)    =P_i^0,0(x),
    Θ_ij^b(x)    =(1-x/2)^iP_j^(2i+1),0(x),
    Θ_ijk^c(x)    =(1-x/2)^i+jP_k^(2i+2j+1),0(x).





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


    ℛ^2 = { (η_1,η_2) ∈ℝ^2; -1 ≤η_1 ≤ 1; -1 ≤η_2 ≤ 1 }


and then map it to the reference triangle 


    𝒯^2 = { (ξ_1,ξ_2) ∈ℝ^2; ξ_1 ≥ -1; ξ_2 ≥ -1; ξ_1 + ξ_2 ≤ 0 }  .


To this end, the linear mapping 


    η_1 =  1 + 2 ξ_1 + ξ_2 / 1 - ξ_2   and  η_2 = ξ_2


is used, which collapses one of the edges of the square. The orthogonal basis on 𝒯^2 is given by  


    φ_(i,j)( ξ_1, ξ_2)=Θ_i^a(η_1)·Θ_ij^b(η_2)   .




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


    ℛ^3 = { (η_1, η_2, η_3) ∈ℝ^3; 0 ≤η_1 ≤ 1; 0 ≤η_2 ≤ 1 ; 0 ≤η_3 ≤ 1 }


and then mapped into a reference tetrahedron 


    𝒯^3 = { (ξ_1,ξ_2,ξ_3) ∈ℝ^3; ξ_1 ≥ 0; ξ_2 ≥ 0; ξ_3 ≥ 0; ξ_1 + ξ_2 + ξ_3 ≤ 1 }  .


The orthogonal basis functions in the reference tetrahedron are given by 


    φ_(i,j,k)(ξ_1, ξ_2, ξ_3) = Θ_i^a(η_1) ·Θ_ij^b(η_2) ·Θ_ijk^c (η_3)


with the mappings


    η_1 = -1 + 2ξ_1 + ξ_2 + ξ_3/1 - ξ_2 - ξ_3,   η_2=-1 + 2ξ_2 + ξ_3/1-ξ_3  and  η_3=-1 + 2ξ_3    .
















§ IMMERSED DISCONTINUOUS GALERKIN METHOD ON NON-BOUNDARY-FITTED MESHES








 §.§ Implicit geometry description



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



    ϕ (x⃗) = 
        
         -   ( x⃗, Γ)      if x⃗∈Ω
    -0    if x⃗∈Γ
    -(x⃗,Γ)    otherwise ,


where (x⃗,Γ) = min_y⃗∈Γ |x⃗ - y⃗|  is  the shortest distance between the point x⃗ and the fluid-solid interface Γ. By definition, the level set function ϕ(x⃗) is negative inside the solid domain Ω and positive inside the fluid domain. The unit normal to the boundary is given by the normalised gradient of the level set function


    n⃗(x⃗) = ∇ϕ (x⃗)/∇ϕ (x⃗)   ,


where x⃗∈Γ is a point on the boundary.


In problems where an implicit function ϕ(x⃗) is not available a sufficiently accurate signed distance function needs to be computed. This can be accomplished either with numerical signed distance computations <cit.> or analytic implicitisation techniques <cit.>.
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 ϕ(x⃗), see Figure <ref>. 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. 









 §.§ Cut-element  integration 



For elements entirely inside the fluid domain the element and boundary integrals appearing in the weak form of the governing equations (<ref>) are evaluated with Gaussian quadrature, see  Section <ref>. 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 <cit.> 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 <ref>. 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 <ref>. 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 <ref> 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.












 §.§ 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 L in the system of  equations (<ref>) are proportional to the physically active support size of the corresponding basis functions. Consequently, the condition number of  the matrix 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 Δ 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 (<ref>) 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 (<ref>). To decide which cut-elements are susceptible to stability issues, the following coverage ratio is computed


    f_j=Vol(ω_j^F)/Vol(ω_j)  ,


where Vol(ω_j^F) is the size, area in 2D or volume in 3D, of the triangulated part of the cut-element and the Vol(ω_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 <ref>. The neighbourhood of a cut-element element ω_j is defined as the set of fluid elements {ω_i} which share a common node with ω_j. This approach is motivated by the extended B-splines proposed in <cit.>.  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 (<ref>).   In a cut-element ω_j with the largest neighbouring fluid element ω_j^∗, see Figure <ref>,  the interpolation equations are modified as follows 


    U⃗_h(x,t)|_ω_j = ∑_i=1^K   [ U^i(t) φ_i(x)  ]__ω_j^∗,   S⃗_h(x,t)|_ω_j = ∑_i=1^K  [ S^i(t)  φ_i(x)  ]__ω_j^∗.



In terms of implementation, if the coverage ratio of a cut-element ω_j is below f_j < 0.3, the discretised conservation vector U⃗_h and auxiliary vector S⃗_h on ω_j are extrapolated from the basis functions and coefficients of the element ω_j^*. According to (<ref>), this is accomplished by simply evaluating the  interpolation equations for the element ω_j^∗ at x⃗∈ω_j. Similarly, the test functions in the weak form (<ref>) are replaced with the basis functions of the  element ω_j^*. In computing the fluxes (<ref>), (<ref>) and (<ref>) the extrapolated field variables are used.  The explicit time integration of the stabilised system of equations  (<ref>) proceeds as usual. 







In the second proposed stabilisation approach, as introduced in the work of Mößner and Reif <cit.> 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 ω_j  to the following modification of the interpolation equations  


    U⃗_h(x,t)=1/f_j∑_i=1^K U^i(t)φ_i(x),   S⃗_h(x,t)=1/f_j∑_i=1^K S^i(t)φ_i(x).


The implementation of this stabilisation approach is slightly easier than the extrapolation approach. In the field discretisation (<ref>) and the weak form (<ref>) the basis functions are replaced with the scaled basis functions φ_i(x⃗) / f_j. This leads to a scaling of the mass matrix L and the coefficient arrays U and S in the semi-discretised  system of equations (<ref>), but not of the right-handside vectors R_ U and R_ S. 

Again, the explicit time integration of the stabilised system of equations (<ref>) proceeds as usual.




 §.§ 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 <cit.>. 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 <ref> and evaluating all element integrals on the curved simplices.  As shown in the work of Zhang <cit.>, 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 <ref>.  That is, only the integration of the boundary integrals needs special treatment in case of curved domain boundaries. According to <cit.>, 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. (<ref>). 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 (<ref>).  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 <ref> 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 ∇ϕ  until the boundary is reached. The intersection of the line given by the gradient with the boundary ϕ(x⃗)=0 is found by bisectioning.




§ 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 =0.8 for subsonic flows and  =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  (<ref>) 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.   



 §.§ 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 


    ϕ(x)=x^2+y^2+z^2-1


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


    e_s=√(∫_Ω(Δ s)^2   v)  ,

where Δ s=p/ρ^γ-s_∞ and the free-stream entropy is s_∞=p_∞/ρ^γ_∞=1.

The convergence of the entropy error with increasing mesh refinement is shown in Figure <ref>. 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. 







 §.§ 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 <ref> and <ref>. The conditions for the flow are =5.0 and =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 





Stationary shock trains consisting of a series of successive shock and rarefaction waves  are frequently used for verification and validation of high-order schemes  <cit.>. Although such problems have been extensively studied, see e.g. <cit.>, 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 <cit.>.  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 <cit.>. 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 <ref>  and  <ref> 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.




 §.§ 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 30 d. 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 <ref> and <ref> 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. 
 
 





  §.§.§ Subsonic flow



For very low Reynolds numbers, with ≪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 > 4 two attached eddies emerge behind the cylinder, see e.g. <cit.>. These eddies increase in size with increasing Reynolds numbers before the flow becomes nonstationary at about ≳ 46 and periodic vortex shedding begins <cit.>. 

In the first set of computations, we use the non-boundary-fitted mesh depicted in Figure <ref> to verify and validate the accuracy of the developed techniques. This mesh contains 62 cut-elements and nine of them require stabilisation. In Figure <ref> the computed average drag coefficients C_d for Reynolds numbers up to = 500 are plotted. The obtained values are in close agreement with the computational results reported by Emblemsvag et al. <cit.> and the experimental results reported by  Wieselsberger <cit.> and Tritton <cit.>. 



In addition to the average drag coefficient, the wake separation length (for ≲ 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 <ref> 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 <ref> the vorticity isocontours for the boundary-fitted and non-boundary fitted meshes and three different Reynolds numbers =20, =40 and =80 are plotted. 

For the two lower Reynolds numbers =20 and =40 the flow is stationary. As evident from the vorticity isocontours, the boundary-fitted and non-boundary fitted meshes yield very similar results. For =20 the three isocontour plots are essentially indistinguishable. For =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 =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. 






  §.§.§ 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. <cit.>. Between the two experiments only the Knudsen number =0.016 and =0.08  is different. The other flow parameters are =5.92,  = 0.707, T_∞=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 <ref> 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 <ref> 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 <ref> the computed drag coefficients are reported which are in excellent agreement with the experiments, with the scaled basis stabilisation giving slightly better results.  










 §.§ 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 <cit.> is investigated. The assumed freestream flow and wall conditions are =5.0,  = 1500, T_0=288  K and T_w= 288   K.  The angle of attack is 25^∘. The axisymmetric command module geometry and the unstructured tetrahedral discretisation mesh are shown in Figure <ref>. 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 <ref> 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. 








 §.§ 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 <ref>.  The freestream flow and wall conditions are = 4.98, =1477,  =0.005, T_∞=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 <ref> 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 <ref>. In the same figure, the experimentally obtained pressure distribution is also shown. The experiment and the computations exhibit an excellent  agreement. 











§ 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 <cit.>. 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. <cit.>. 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 <cit.>. 

 


§ 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.


wileyj


