A novel stabilization method for high-order shock fitting with finite element methods

https://doi.org/10.1016/j.jcp.2020.110096Get rights and content

Highlights

  • We develop a moving-grid shock fitting method to accurately simulate flow problems with shock waves.

  • A novel stabilization term is introduced to the shock interface motion equation to ensure solution convergence.

  • High-order solution accuracy is demonstrated for quantities of interest in practical problems.

Abstract

A moving-grid, shock-fitting, finite element method has been implemented that can achieve high-order accuracy for flow simulations with shocks. In this approach, element edges in the computational mesh are fitted to the shock front and moved with the shock throughout the simulation. The Euler or Navier-Stokes equations are solved on the moving mesh in an arbitrary Lagrangian-Eulerian framework. The method is implemented in two-dimensions in the context of a streamwise upwind Petrov-Galerkin finite element discretization with unstructured triangular meshes and mesh adaptation. It is shown that the shock interface motion equation has a wave nature, and disturbances can propagate along the shock interface. A SUPG stabilization term is introduced to the interface motion equation that is critical for ensuring that interface disturbances do not lead to non-convergent solution behavior. The formal order of accuracy of the scheme is verified, and the performance of the proposed scheme is assessed for both inviscid and viscous problems. It was found that the present scheme predicts smooth and noise-free surface heating for hypersonic flow over a cylinder with purely irregular triangular elements.

Introduction

Finite element methods have the potential to significantly increase the reliability of flow predictions because they can attain high-order accuracy on unstructured meshes [1], [2], [3], [4], [5]. This allows accurate solutions to be obtained for complex geometries at a lower cost than second order accurate methods. One of the major challenges in applying these methods to hypersonic simulations is the accurate evaluation of shock waves.

To incorporate shocks into the solution, many researchers have implemented shock capturing techniques. In shock capturing, the shock discontinuity is spread across a number of mesh elements and a method must be found for maintaining the monotonicity of the solution in the shock region. It has been widely reported that, for strong shocks, many classical shock capturing schemes can provide unreliable results under some circumstances [6]. Some of the reported failings of shock capturing schemes include the carbuncle phenomenon [7] and odd-even decoupling along the shock wave [8]. Furthermore, in Ref. [9], a number of shock capturing flux functions were investigated for 1D and 2D hypersonic shock problems, and none were found to be universally stable across all dimensions and test cases. Additionally, as shown by Godunov [10], all linear, monotone schemes will necessarily be reduced to first order accuracy near the discontinuity. Among existing shock capturing schemes, high-order compact shock capturing techniques such as weighted essentially nonoscillatory (WENO) polynomial limiters are the preferred methods for high-order schemes as they preserve the order of accuracy in the smooth region, and are shown to be stable for strong and complicated shock structures [4].

Another approach for incorporating shocks into numerical simulations is the shock fitting technique. This approach aims to mitigate the problems associated with the shock capturing method by maintaining the shock wave as a true flow discontinuity. Mesh edges are aligned with the shock boundary throughout the simulation, and the Rankine-Hugoniot shock jump conditions are used to determine the shock behavior. The idea of shock fitting has been around for many years [11]. More recently, shock fitting techniques have been accurately applied to the solution of numerous shock problems [2], [12], [13], demonstrating the promise of utilizing these methods in practical applications.

There are several possible implementations of shock fitting including fixed grid (also called floating shock fitting) approaches [14], [15], [16], [17], [18], [19] and moving grid approaches [16], [20]. Fixed grid approaches allow the shock to move over a background mesh, and require local remeshing and/or solution interpolation in the shock region to maintain the shock as a flow discontinuity. Most fixed grid approaches have been limited to second-order accuracy, however recent work [16] has demonstrated higher-order accuracy by using one-sided finite difference schemes. In Ref. [16], moving grid shock fitting methods were found to be more robust than fixed grid approaches for problems with relatively simple shock geometries.

In the moving grid approach, the entire computational mesh moves with the fitted shock so that a single mesh boundary remains aligned with the shock interface throughout the simulation. This approach has been used previously for shocks to obtain high order accurate results [16], [21] by using either spectral or high-order finite difference approaches. However, the high-order solutions were typically limited to fairly simple shock and physical geometries. More complex geometries have been evaluated using low order methods, for example a simulation of an aircraft with interacting shocks is presented in Ref. [22] using a moving grid approach. More recently, complex shock interactions including triple junctions, reflected shocks, and contact discontinuities have been simulated using a cell-centered moving-mesh approach on moving shock-fitted grids [20].

In both moving grid and fixed grid shock fitting approaches, the Rankine-Hugoniot shock jump conditions are used to determine the shock interface motion and the corresponding mesh motion. A different approach to shock fitting was developed in [2], [12] which follows an optimization based approach to maintaining mesh alignment with the shock discontinuity. Although they follow slightly different procedures, these methods generally enforce interface alignment by incorporating an optimization step into the finite element solution process. The objective function used during optimization is developed such that optimal solutions are attained for an interface-aligned mesh; thus, the optimization process ensures that all finite-element solutions are obtained on discontinuity-aligned meshes.

The use of unstructured meshes, which has become more common in the CFD community since the 1990s [23], presents two major assets to the development of modern shock fitting techniques. Unstructured meshes are able to describe complex geometries including complex shock structures, and they are well-suited for applications which require mesh deformation and adaptation. As demonstrated in recent shock fitting advances [14], [20], these advantages greatly improve the ability to track shock waves with complex geometries that undergo significant motion and deformation. The development of robust, high-order finite element methods for unstructured meshes inspires our effort to develop a shock fitting method which can achieve high-order accuracy for complex shock problems.

The goal of the present work is to develop a robust numerical technique for obtaining high-order accurate results for supersonic and hypersonic flow problems with shock waves. To that end, we combine a continuous Galerkin finite element method with an arbitrary Lagrangian-Eulerian mesh movement simulation which tracks the position of the shock following the moving grid approach. A novel stabilization term is developed for the equation governing the shock interface motion to ensure that convergent solution behavior is maintained. Unstructured mesh adaptation and deformation techniques are used to maintain mesh quality for problems with significant shock interface motion and deformation. The method is implemented in two spatial dimensions and high-order results are demonstrated for a supersonic nozzle problem and a bow shock problem.

Section snippets

Mathematical formulation

This work focuses on the accurate numerical simulation of high-speed, compressible flows in two spatial dimensions. The compressible Navier-Stokes or Euler equations govern the flow behavior, and the Rankine-Hugoniot conditions govern the behavior of the shock wave. In the shock fitting approach, the computational domain is divided into two regions separated by the shock discontinuity. The mesh edges which are aligned with the shock interface move with the shock to remain aligned throughout the

Numerical formulation

Both the upstream and downstream domain regions are discretized using unstructured meshes made up of triangular elements. The two meshes coincide along the shock interface, and the mesh interface moves to track the position of the shock wave. The interior grid points move to reduce mesh deformation while the outer grid boundaries remain fixed in space, and mesh adaptation is used to ensure that mesh quality is maintained for problems with significant shock motion. The flow solution is evaluated

Results

In this section the proposed method is implemented to solve multiple test cases. The tests are used to evaluate the accuracy of the method as well as its applicability to realistic shock problems. For each problem, the initial shock location and topology are assumed, and a mesh boundary is placed at this initial location. The shock evaluation requires the initial shock location to be chosen such that the upstream flow is supersonic. An unsteady simulation allows the shock boundary to move and

Conclusions

An arbitrary Lagrangian-Eulerian moving grid shock fitting finite element method has been developed for solving hypersonic flow problems. The method is implemented for a two dimensional continuous Galerkin finite element code with mesh motion and adaptation. It is shown that high-wave-number numerical instabilities can propagate along the shock interface due to the wave-like nature of the equation governing the normal motion of the shock interface. A novel SUPG stabilization term is introduced

CRediT authorship contribution statement

Luke M. D'Aquila: Investigation, Software, Validation, Writing – original draft. Brian T. Helenbrook: Conceptualization, Methodology, Supervision, Writing – review & editing. Alireza Mazaheri: Conceptualization, Methodology, Supervision, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Funding: This research was funded by the NASA Fellowship Activity 2018 cohort [grant number 80NSSC18K1698].

References (40)

  • S. Sherwin et al.

    A triangular spectral element method; application to the incompressible Navier-Stokes equations

    Comput. Methods Appl. Mech. Eng.

    (1995)
  • M.J. Brazell et al.

    Continuous finite elements on tetrahedra with local mass matrix inversion to solve the preconditioned compressible Navier–Stokes equations

    Comput. Fluids

    (2013)
  • B.T. Helenbrook et al.

    Preconditioning for dual-time-stepping simulations of the shallow water equations including Coriolis effects

    J. Comput. Phys.

    (2008)
  • B.T. Helenbrook

    A two-fluid spectral element method

    Comput. Methods Appl. Mech. Eng.

    (2001)
  • V. Venkatakrishnan et al.

    Implicit method for the computation of unsteady flows on unstructured grids

    J. Comput. Phys.

    (1996)
  • R. Williams et al.

    A four-stage index 2 diagonally implicit Runge-Kutta method

    Appl. Numer. Math.

    (2002)
  • S.J. Sherwin et al.

    Tetrahedral hp finite elements: algorithms and flow simulations

    J. Comput. Phys.

    (1996)
  • A. Corrigan et al.

    A moving discontinuous Galerkin finite element method for flows with interfaces

    Int. J. Numer. Methods Fluids

    (2018)
  • A. Mazaheri et al.

    Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: triangular elements

    J. Comput. Phys.

    (2019)
  • J.J. Quirk

    A contribution to the great Riemann solver debate

    Int. J. Numer. Methods Fluids

    (1994)
  • Cited by (12)

    • A robust, high-order implicit shock tracking method for simulation of complex, high-speed flows

      2022, Journal of Computational Physics
      Citation Excerpt :

      These methods are not easily applicable to discontinuities whose topologies are not known a priori. While interest in shock tracking/fitting has seen somewhat of a resurgence in recent years [8,4,15,12], shock tracking is not widely used for practical applications. A new class of high-order numerical methods has recently emerged, implicit shock tracking, that includes the High-Order Implicit Shock Tracking (HOIST) method [45,48,47] and the Moving Discontinous Galerkin Method with Interface Condition Enforcement (MDG-ICE) [10,24,25].

    • Implicit shock tracking for unsteady flows by the method of lines

      2022, Journal of Computational Physics
      Citation Excerpt :

      These methods are not easily applicable to discontinuities whose topologies not known a priori. While interest in shock tracking/fitting methods has seen somewhat of a resurgence in recent years [7,5,18,12], shock tracking is largely not used in practical CFD today. To our knowledge, the only other approach to implicit shock tracking is the Moving Discontinuous Galerkin Method with Interface Condition Enforcement (MDG-ICE), proposed in [9,11], where the authors enforce a DG discretization with unconventional numerical fluxes and the Rankine-Hugoniot interface conditions in a minimum-residual sense.

    • Extrapolated DIscontinuity Tracking for complex 2D shock interactions

      2022, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      Most of the aforementioned developments have been included in an open-source platform [13] with the goal of fostering even more the interest in such methods. Thanks to these new developments the interest in shock-fitting/tracking methods has seen a resurgence with the proposal of several implicit and explicit algorithms [14–22]. A limitation of this technique is that it heavily relies on the flexibility of triangular and tetrahedral grids to locally produce a fitted unstructured grid around the discontinuities.

    • UnDiFi-2D: An unstructured discontinuity fitting code for 2D grids

      2022, Computer Physics Communications
      Citation Excerpt :

      The use of shape-functions that are continuous across the element interfaces, which is the case with SUPG, or discontinuous, such as in DG, has implications on how discontinuities are fitted. Similarly to the algorithm described in this paper, in the SUPG-FEM of [41] the discontinuities are internal boundaries of zero thickness: by doing so, a finite jump in the dependent variables can take place while crossing the discontinuity. On the other hand, numerical methods that employ a data representation which is discontinuous across the cell interfaces, which is the case with DG-FEM, but also with cell-centered FV methods, allow to fit discontinuities as a collection of edges of the mesh, without introducing internal boundaries.

    View all citing articles on Scopus
    View full text