A novel stabilization method for high-order shock fitting with finite element methods
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)
- et al.
High-order accurate discontinuous finite element solution of the 2D Euler equations
J. Comput. Phys.
(1997) - et al.
High-order adaptive arbitrary-Lagrangian–Eulerian (ALE) simulations of solidification
Comput. Fluids
(2018) - et al.
Efficient high-order discontinuous Galerkin schemes with first-order hyperbolic advection-diffusion system approach
J. Comput. Phys.
(2016) - et al.
An optimization-based approach for high-order accurate discretization of conservation laws with discontinuous solutions
J. Comput. Phys.
(2018) - et al.
A shock-fitting technique for 2D unstructured grids
Comput. Fluids
(2009) - et al.
On high-order shock-fitting and front-tracking schemes for numerical simulation of shock–disturbance interactions
J. Comput. Phys.
(2010) - et al.
Shock interaction computations on unstructured, two-dimensional grids using a shock-fitting technique
J. Comput. Phys.
(2011) - et al.
An unstructured, three-dimensional, shock-fitting solver for hypersonic flows
Comput. Fluids
(2013) - et al.
A shock-fitting technique for cell-centered finite volume methods on unstructured dynamic meshes
J. Comput. Phys.
(2017) - et al.
Computation of three dimentional flows about aircraft configurations
Comput. Fluids
(1973)
A triangular spectral element method; application to the incompressible Navier-Stokes equations
Comput. Methods Appl. Mech. Eng.
Continuous finite elements on tetrahedra with local mass matrix inversion to solve the preconditioned compressible Navier–Stokes equations
Comput. Fluids
Preconditioning for dual-time-stepping simulations of the shallow water equations including Coriolis effects
J. Comput. Phys.
A two-fluid spectral element method
Comput. Methods Appl. Mech. Eng.
Implicit method for the computation of unsteady flows on unstructured grids
J. Comput. Phys.
A four-stage index 2 diagonally implicit Runge-Kutta method
Appl. Numer. Math.
Tetrahedral hp finite elements: algorithms and flow simulations
J. Comput. Phys.
A moving discontinuous Galerkin finite element method for flows with interfaces
Int. J. Numer. Methods Fluids
Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: triangular elements
J. Comput. Phys.
A contribution to the great Riemann solver debate
Int. J. Numer. Methods Fluids
Cited by (12)
High-order implicit shock tracking boundary conditions for flows with parametrized shocks
2023, Journal of Computational PhysicsA robust, high-order implicit shock tracking method for simulation of complex, high-speed flows
2022, Journal of Computational PhysicsCitation 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 PhysicsCitation 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 EngineeringCitation 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 CommunicationsCitation 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.