High-order Accurate Discontinuous Galerkin Research Articles

High-order accurate positivity-preserving and well-balanced discontinuous Galerkin schemes for ten-moment Gaussian closure equations with source terms

This paper proposes novel high-order accurate discontinuous Galerkin (DG) schemes for the one- and two-dimensional ten-moment Gaussian closure equations with source terms defined by a known potential function. Our DG schemes exhibit the desirable capability of being well-balanced (WB) for a known hydrostatic equilibrium state while simultaneously preserving positive density and positive-definite anisotropic pressure tensor. The well-balancedness is built on carefully modifying the solution states in the Harten–Lax–van Leer–contact (HLLC) flux, and appropriate reformulation and discretization of the source terms. Our novel modification technique overcomes the difficulties posed by the anisotropic effects, maintains the high-order accuracy, and ensures that the modified solution state remains within the physically admissible state set. We provide the rigorous positivity-preserving analyses of our WB DG schemes, based on several key properties of the admissible state set, the HLLC flux and the HLLC solver, as well as the geometric quasilinearization (GQL) approach in Wu and Shu (2023) [52], which was originally applied to analyze the admissible state set and the physical-constraints-preserving schemes for the relativistic magnetohydrodynamic equations in Wu and Tang (2017) [54], to address the difficulties arising from the nonlinear constraints on the pressure tensor. Moreover, the proposed WB DG schemes satisfy the weak positivity for the cell averages, implying the use of a simple scaling limiter to enforce the physical admissibility of the DG solution polynomials at certain points of interest. Extensive numerical experiments are conducted to validate the preservation of equilibrium states, accuracy in capturing small perturbations to such states, robustness in solving problems involving low density or low pressure, and high resolution for both smooth and discontinuous solutions.

Very High-Order Accurate Discontinuous Galerkin Method for Curved Boundaries with Polygonal Meshes

Preserving the optimal convergence order of discontinuous Galerkin (DG) discretisations in curved domains is a critical and well-known issue. The proposed approach relies on the reconstruction for off-site data (ROD) method developed originally within the finite volume framework. The main advantages are simplicity, since the PDE solver only considers polygonal domains, and versatility, since any type of boundary condition can be imposed. The developed DG–ROD method consists in splitting the boundary conditions treatment and the leading discrete equations from a classical DG formulation into two independent solvers coupled in a simple and efficient iterative procedure. A numerical benchmark is provided to assess the capability of the method with Dirichlet and Neumann boundary conditions prescribed on curved boundaries, demonstrating that the optimal convergence order is effectively achieved.

High-order accurate transient and free-vibration analysis of plates and shells

The limited availability of analytical solutions and the high cost associated with experimental testing motivate the use of computational tools to assess the dynamic behavior of load-bearing components, especially when a wide design space must be explored, as is often the case with composite structures. In this context, a novel high-order accurate discontinuous Galerkin formulation for transient and free-vibration analysis of multilayered plates and shells is presented and numerical validated. The starting point of the formulation is a generalized structural theory for multilayered shells with arbitrary curvature based on the expansion of the displacement covariant components throughout the shell thickness. The variational statement of three-dimensional elastodynamics allows deriving the strong form of the governing differential equations, which form the basis to obtain the corresponding discontinuous Galerkin weak statements. As the order of the through-the-thickness expansion and the order of the discontinuous Galerkin basis functions are free parameters, the proposed approach allows tuning the order of accuracy of the computed solution throughout both the shell thickness and the shell modeling domain. Numerical results are reported and discussed for several validation test cases in terms of h- and p-convergence analyses, demonstrating the high-order accuracy, robustness, and computational savings of the formulation.

BE-BDF2 Time Integration Scheme Equipped with Richardson Extrapolation for Unsteady Compressible Flows

In this work we investigate the effectiveness of the Backward Euler-Backward Differentiation Formula (BE-BDF2) in solving unsteady compressible inviscid and viscous flows. Furthermore, to improve its accuracy and its order of convergence, we have equipped this time integration method with the Richardson Extrapolation (RE) technique. The BE-BDF2 scheme is a second-order accurate, A-stable, L-stable and self-starting scheme. It has two stages: the first one is the simple Backward Euler (BE) and the second one is a second-order Backward Differentiation Formula (BDF2) that uses an intermediate and a past solution. The RE is a very simple and powerful technique that can be used to increase the order of accuracy of any approximation process by eliminating the lowest order error term(s) from its asymptotic error expansion. The spatial approximation of the governing Navier–Stokes equations is performed with a high-order accurate discontinuous Galerkin (dG) method. The presented numerical results for canonical test cases, i.e., the isentropic convecting vortex and the unsteady vortex shedding behind a circular cylinder, aim to assess the performance of the BE-BDF2 scheme, in its standard version and equipped with RE, by comparing it with the ones obtained by using more classical methods, like the BDF2, the second-order accurate Crank–Nicolson (CN2) and the explicit third-order accurate Strong Stability Preserving Runge–Kutta scheme (SSP-RK3).

On the entropy conserving/stable implicit DG discretization of the Euler equations in entropy variables

The aim of this paper is to investigate the behavior of a high-order accurate Discontinuous Galerkin entropy conserving/stable scheme in space for unsteady compressible inviscid flows. In order to ensure an entropy conserving/stable scheme in space, several entropy conserving/stable numerical fluxes are considered. For the time discretization, high-order accurate linearly implicit Rosenbrock-type Runge–Kutta schemes are used. These schemes cannot be provably entropy conserving/stable, but, for an enough small time step size, the error due to the time integration has negligible contribution, and entropy conserving/stable properties can be fulfilled. The properties of the fully discrete system of equations are assessed in a series of numerical experiments of growing complexity, i.e.: (i) the isentropic vortex convection problem, (ii) the double shear layer, (iii) the Sod shock tube and (iv) the Taylor–Green vortex. For each one of these test-case the accuracy and the related order of convergence of the time integration schemes and of the numerical fluxes employed will be investigated, as well as the comparison between the accuracy provided by the entropy and the primitive set of variables and their stability properties. Furthermore, same relevant issues related to the use of entropy conserving fluxes are addressed, like the necessity to use higher accurate quadrature rules, the spatial sub-optimal order of convergence, and the behavior of the fully discrete system of equations when long-time simulations are performed.

High order weight‐adjusted discontinuous Galerkin methods for wave propagation on moving curved meshes

AbstractThis article presents high order accurate discontinuous Galerkin (DG) methods for wave problems on moving curved meshes with general choices of basis and quadrature. The proposed method adopts an arbitrary Lagrangian–Eulerian formulation to map the wave equation from a time‐dependent moving physical domain onto a fixed reference domain. For moving curved meshes, weighted mass matrices must be assembled and inverted at each time step when using explicit time‐stepping methods. We avoid this step by utilizing an easily invertible weight‐adjusted approximation. The resulting semi‐discrete weight‐adjusted DG scheme is provably energy stable up to a term that (for a fixed time interval) converges to zero with the same rate as the optimal error estimate. Numerical experiments using both polynomial and B‐spline bases verify the high order accuracy and energy stability of proposed methods.

An artificial Equation of State based Riemann solver for a discontinuous Galerkin discretization of the incompressible Navier–Stokes equations

In this work we present a novel exact Riemann solver for an artificial Equation of State (EoS) based modification of the incompressible Euler equations. Differently from the well known artificial compressibility method, this new approach overcomes the lack of the pressure-velocity coupling by using a suitably designed artificial EoS. The modified set of equations fits into the framework of first order hyperbolic conservation laws. Accordingly, an exact Riemann solver is derived. This new solver has the advantage to avoid a wave pattern violation issue which can affect the exact Riemann problem solution obtained by the standard artificial compressibility approach. The new artificial EoS based Riemann solver can be used as a tool for the definition of the advective Godunov fluxes in a Finite Volume or a discontinuous Galerkin discretization of the incompressible Navier–Stokes (INS) equations. We assess and analyse the new exact Riemann solver on 1D Riemann problems. Its capability and effectiveness are then shown in the context of a high-order accurate discontinuous Galerkin discretization of the INS equations. In particular, we verify the convergence properties, both in space and time, considering two test cases in 2D, namely, the Kovasznay test case and a damped travelling waves test case. Finally, we perform an implicit large eddy simulation of the incompressible turbulent flow over periodic hills where the classic forcing term formulation is modified to deal with variable time steps. Solution comparisons with respect to numerical and experimental results from the literature are given.

Provably physical-constraint-preserving discontinuous Galerkin methods for multidimensional relativistic MHD equations

We propose and analyze a class of robust, uniformly high-order accurate discontinuous Galerkin (DG) schemes for multidimensional relativistic magnetohydrodynamics (RMHD) on general meshes. A distinct feature of the schemes is their physical-constraint-preserving (PCP) property, i.e., they are proven to preserve the subluminal constraint on the fluid velocity and the positivity of density, pressure, and internal energy. This is the first time that provably PCP high-order schemes are achieved for multidimensional RMHD. Developing PCP high-order schemes for RMHD is highly desirable but remains a challenging task, especially in the multidimensional cases, due to the inherent strong nonlinearity in the constraints and the effect of the magnetic divergence-free condition. Inspired by some crucial observations at the PDE level, we construct the provably PCP schemes by using the locally divergence-free DG schemes of the recently proposed symmetrizable RMHD equations as the base schemes, a limiting technique to enforce the PCP property of the DG solutions, and the strong-stability-preserving methods for time discretization. We rigorously prove the PCP property by using a novel “quasi-linearization” approach to handle the highly nonlinear physical constraints, technical splitting to offset the influence of divergence error, and sophisticated estimates to analyze the beneficial effect of the additional source term in the symmetrizable RMHD system. Several two-dimensional numerical examples are provided to further confirm the PCP property and to demonstrate the accuracy, effectiveness and robustness of the proposed PCP schemes.

Uniformly High-Order Structure-Preserving Discontinuous Galerkin Methods for Euler Equations with Gravitation: Positivity and Well-Balancedness

This paper presents a class of novel high-order accurate discontinuous Galerkin (DG) schemes for the compressible Euler equations under gravitational fields. A notable feature of these schemes is that they are well-balanced for a general hydrostatic equilibrium state, and at the same time, provably preserve the positivity of density and pressure. In order to achieve the well-balanced and positivity-preserving properties simultaneously, a novel DG spatial discretization is carefully designed with suitable source term reformulation and a properly modified Harten-Lax-van Leer contact (HLLC) flux. Based on some technical decompositions as well as several key properties of the admissible states and HLLC flux, rigorous positivity-preserving analyses are carried out. It is proven that the resulting well-balanced DG schemes, coupled with strong stability preserving time discretizations, satisfy a weak positivity property, which implies that one can apply a simple existing limiter to effectively enforce the positivity-preserving property, without losing high-order accuracy and conservation. The proposed methods and analyses are applicable to the Euler system with general equation of state. Extensive one- and two-dimensional numerical tests demonstrate the desired properties of these schemes, including the exact preservation of the equilibrium state, the ability to capture small perturbation of such state, the robustness for solving problems involving low density and/or low pressure, and good resolution for smooth and discontinuous solutions.

Discontinuous Galerkin solution of the RANS and [formula omitted] equations for natural and bypass transition

Transition modelling represents a key ingredient for improving the performance predictions of many industrial applications. Among transition models, local formulations seem to guarantee better robustness, accuracy and easiness of implementation in modern CFD solvers. These models have been proposed in the finite volume context to predict the laminar-turbulent transition, but only few attempts have been made in the high-order framework. In this paper a new phenomenological transition model based on the concept of laminar kinetic energy has been proposed and implemented in a high-order accurate Discontinuous Galerkin code, named MIGALE. The transition model is validated and assessed by computing the transitional flow around flat plates with zero/adverse pressure gradients and through different turbine nozzles (T106A and LS89) for different values of Reynolds number and turbulence intensity. The computed results have been compared with experimental data and reference numerical solutions.

High order ADER-DG schemes for the simulation of linear seismic waves induced by nonlinear dispersive free-surface water waves

In this paper, we propose a unified and high order accurate fully-discrete one-step ADER Discontinuous Galerkin method for the simulation of linear seismic waves in the sea bottom that are generated by the propagation of free surface water waves. In particular, a hyperbolic reformulation of the Serre-Green-Naghdi model for nonlinear dispersive free surface flows is coupled with a first order velocity-stress formulation for linear elastic wave propagation in the sea bottom. To this end, Cartesian non-conforming meshes are defined in the solid and fluid domains and the coupling is achieved by an appropriate time-dependent pressure boundary condition in the three-dimensional domain for the elastic wave propagation, where the pressure is a combination of hydrostatic and non-hydrostatic pressure in the water column above the sea bottom. The use of a first order hyperbolic reformulation of the nonlinear dispersive free surface flow model leads to a straightforward coupling with the linear seismic wave equations, which are also written in first order hyperbolic form. It furthermore allows the use of explicit time integrators with a rather generous CFL-type time step restriction associated with the dispersive water waves, compared to numerical schemes applied to classical dispersive models that contain higher order derivatives and typically require implicit solvers. Since the two systems that describe the seismic waves and the free surface water waves are written in the same form of a first order hyperbolic system they can also be efficiently solved in a unique numerical framework. In this paper we choose the family of arbitrary high order accurate discontinuous Galerkin finite element schemes, which have already shown to be suitable for the numerical simulation of wave propagation problems. The developed methodology is carefully assessed by first considering several benchmarks for each system separately, i.e. in the framework of linear elasticity and non-hydrostatic free surface flows, showing a good agreement with exact and numerical reference solutions. Finally, also coupled test cases are addressed. Throughout this paper we assume the elastic deformations in the solid to be sufficiently small so that their influence on the free surface water waves can be neglected.

Space-time adaptive ADER discontinuous Galerkin schemes for nonlinear hyperelasticity with material failure

We are concerned with the numerical solution of a unified first order hyperbolic formulation of continuum mechanics that goes back to the work of Godunov, Peshkov and Romenski [65,68,97] (GPR model) and which is an extension of nonlinear hyperelasticity that is able to describe simultaneously nonlinear elasto-plastic solids at large strain, as well as viscous and ideal fluids. The proposed governing PDE system also contains the effect of heat conduction and can be shown to be symmetric and thermodynamically compatible, as it obeys the first and second law of thermodynamics. In this paper we extend the GPR model to the simulation of nonlinear dynamic rupture processes, which can be achieved by adding an additional scalar to the governing PDE system. This extra parameter describes the material damage and is governed by an advection-reaction equation, where the stiff and highly nonlinear reaction mechanisms depend on the ratio of the local equivalent stress to the yield stress of the material. The stiff reaction mechanisms are integrated in time via an efficient exponential time integrator. Due to the multiple spatial and temporal scales involved in the problem of crack generation and propagation, the model is solved on space–time adaptive Cartesian meshes using high order accurate discontinuous Galerkin finite element schemes endowed with an a posteriori subcell finite volume limiter. A key feature of our new model is the use of a twofold diffuse interface approach that allows the cracks to form anywhere and at any time, independently of the chosen computational grid, which is simply adaptive Cartesian (AMR). This is substantially different from many fracture modeling approaches that need to resolve discontinuities explicitly, such as for example dynamic shear rupture models used in computational seismology, where the geometry of the rupture fault needs to prescribed a priori. We furthermore make use of a scalar volume fraction function α that indicates whether a given spatial point is inside the solid (α=1) or outside (α=0), thus allowing the description of solids of arbitrarily complex shape. We show extensive numerical comparisons with experimental results for stress-strain diagrams of different real materials and for the generation and propagation of fracture in rocks and pyrex glass at low and high velocities. Overall, a very good agreement between numerical simulations and experiments is obtained. The proposed model is also naturally able to describe material fatigue.

Interactions of two bubbles along a gaseous interface undergoing the Richtmyer–Meshkov instability in two dimensions

The vorticity-driven interactions involving bubbles and spikes along a shock-accelerated interface significantly affect the dynamics of the flow and likely contribute to an eventual turbulent transition. In this study, the interaction of two adjacent bubbles of different sizes undergoing the Richtmyer–Meshkov process is investigated numerically in two dimensions. In particular, we compare two-bubble configurations with different initial amplitude ratios and wavelength ratios. Simulations are performed using an in-house, high-order accurate Discontinuous Galerkin code solving the Euler equations for a shock-accelerated interface initially perturbed with a periodic array of alternating sinusoids of different sizes. Two types of departures from existing bubble-merging models are observed, caused by a vorticity imbalance on the sides of the bubbles leading to spikes that either diverge or converge beneath the larger bubble, thus resulting in the generation of vortex dipoles escaping the confines of the mixing region or in the reacceleration of the larger bubbles. A vorticity-based criterion for predicting the generation of vortex dipoles and the reacceleration of larger bubbles from the initial conditions is developed.

Entropy stable numerical approximations for the isothermal and polytropic Euler equations

In this work we analyze the entropic properties of the Euler equations when the system is closed with the assumption of a polytropic gas. In this case, the pressure solely depends upon the density of the fluid and the energy equation is not necessary anymore as the mass conservation and momentum conservation then form a closed system. Further, the total energy acts as a convex mathematical entropy function for the polytropic Euler equations. The polytropic equation of state gives the pressure as a scaled power law of the density in terms of the adiabatic index gamma . As such, there are important limiting cases contained within the polytropic model like the isothermal Euler equations (gamma {=}1) and the shallow water equations (gamma {=}2). We first mimic the continuous entropy analysis on the discrete level in a finite volume context to get special numerical flux functions. Next, these numerical fluxes are incorporated into a particular discontinuous Galerkin (DG) spectral element framework where derivatives are approximated with summation-by-parts operators. This guarantees a high-order accurate DG numerical approximation to the polytropic Euler equations that is also consistent to its auxiliary total energy behavior. Numerical examples are provided to verify the theoretical derivations, i.e., the entropic properties of the high order DG scheme.

Well-Balanced Discontinuous Galerkin Method for Shallow Water Equations with Constant Subtraction Techniques on Unstructured Meshes

The classical Saint–Venant shallow water equations on complex geometries have wide applications in many areas including coastal engineering and atmospheric modeling. The main numerical challenge in simulating Saint–Venant equations is to maintain the high order of accuracy and well-balanced property simultaneously. In this paper, we propose a high-order accurate and well-balanced discontinuous Galerkin (DG) method on two dimensional (2D) unstructured meshes for the Saint–Venant shallow water equations. The technique used to maintain well-balanced property is called constant subtraction and proposed in Yang et al. (J Sci Comput 63:678–698, 2015). Hierarchical reconstruction limiter with a remainder correction technique is introduced to control numerical oscillations. Numerical examples with smooth and discontinuous solutions are provided to demonstrate the performance of our proposed DG methods.

Positivity Preserving Limiters for Time-Implicit Higher Order Accurate Discontinuous Galerkin Discretizations

Currently, nearly all positivity preserving discontinuous Galerkin (DG) discretizations of partial differential equations are coupled with explicit time integration methods. Unfortunately, for many problems this can result in severe time-step restrictions. The techniques used to develop explicit positivity preserving DG discretizations cannot, however, easily be combined with implicit time integration methods. In this paper, we therefore present a new approach. Using Lagrange multipliers, the conditions imposed by the positivity preserving limiters are directly coupled to a DG discretization combined with a diagonally implicit Runge--Kutta time integration method. The positivity preserving DG discretization is then reformulated as a Karush--Kuhn--Tucker (KKT) problem, which is frequently encountered in constrained optimization. Since the limiter is only active in areas where positivity must be enforced, it does not affect the higher order DG discretization elsewhere. The resulting nonsmooth nonlinear algebraic equations have, however, a different structure compared to most constrained optimization problems. We therefore develop an efficient active set semismooth Newton method that is suitable for the KKT formulation of time-implicit positivity preserving DG discretizations. Convergence of this semismooth Newton method is proven using a specially designed quasi-directional derivative of the time-implicit positivity preserving DG discretization. The time-implicit positivity preserving DG discretization is demonstrated for several nonlinear scalar conservation laws, which include the advection, Burgers, Allen--Cahn, Barenblatt, and Buckley--Leverett equations.

A stable and high-order accurate discontinuous Galerkin based splitting method for the incompressible Navier–Stokes equations

In this paper we consider discontinuous Galerkin (DG) methods for the incompressible Navier–Stokes equations in the framework of projection methods. In particular we employ symmetric interior penalty DG methods within the second-order rotational incremental pressure correction scheme. The major focus of the paper is threefold: i) We propose a modified upwind scheme based on the Vijayasundaram numerical flux that has favourable properties in the context of DG. ii) We present a novel postprocessing technique in the Helmholtz projection step based on H(div) reconstruction of the pressure correction that is computed locally, is a projection in the discrete setting and ensures that the projected velocity satisfies the discrete continuity equation exactly. As a consequence it also provides local mass conservation of the projected velocity. iii) Numerical results demonstrate the properties of the scheme for different polynomial degrees applied to two-dimensional problems with known solution as well as large-scale three-dimensional problems. In particular we address second-order convergence in time of the splitting scheme as well as its long-time stability.

Artificial compressibility Godunov fluxes for variable density incompressible flows

In this work we present and compare three Riemann solvers for the artificial compressibility perturbation of the 1D variable density incompressible Euler equations. The goal is to devise an artificial compressibility flux formulation to be used in Finite Volume or discontinuous Galerkin discretizations of the variable density incompressible Navier–Stokes equations. Starting from the constant density case, two Riemann solvers taking into account density jumps at fluid interfaces are first proposed. By enforcing the divergence free constraint in the continuity equation, these approximate Riemann solvers deal with density as a purely advected quantity. Secondly, by retaining the conservative form of the continuity equation, the exact Riemann solver is derived. The variable density solution is fully coupled with velocity and pressure unknowns. The Riemann solvers are compared and analysed in terms of robustness on harsh 1D Riemann problems. The extension to multidimensional problems is described. The effectiveness of the exact Riemann solver is demonstrated in the context of an high-order accurate discontinuous Galerkin discretization of variable density incompressible flow problems. We numerically validate the implementation considering the Kovasznay test case and the Rayleigh–Taylor instability problem.

A Front Tracking Method Based on Runge-Kutta Discontinuous Galerkin Methods

In this paper, a high-resolution front tracking method was presented for interface tracking simulation with Runge-Kutta discontinuous Galerkin methods. An interface treating method of the discontinuous methods is presented. This method don’t construct the ghost fluid and the flow information on both sides next to the interface is used to solve the interfacial status. The limiter adopted the combination of the shock detection and monotonicity-preserving limiter and level set method is used for tracking the interface. Result shown that the front tracking of the high-order accurate Runge-Kutta discontinuous Galerkin method exhibits very good agreement with exact solution in the interface condition that contain strong shock.

High-order Discontinuous Galerkin Method for Solving Elliptic Interface Problems

In this study, we develop a high-order accurate discontinuous Galerkin scheme using curvilinear quadrilateral elements for solving elliptic interface problems. To maintain accuracy for curvilinear quadrilateral elements with $Q^k$-polynomial basis functions, we select Legendre-Gauss-Lobatto quadrature rules with $k+2$ integration points, including end points, on each edge. Numerical experiments show quadrature rules are appropriate and the numerical solution converges with the order $k+1$. Moreover, we implement the Uzawa method to rewrite the original linear system into smaller linear systems. The resulting method is three times faster than the original system. We also find that high-order methods are more efficient than low-order methods. Based on this result, we conclude that under the condition of limited computational resources, the best approach for achieving optimal accuracy is to solve a problem using the coarsest mesh and local spaces with the highest degree of polynomials.