High-Order Positivity-Preserving Well-Balanced Discontinuous Galerkin Methods for Euler Equations with Gravitation on Unstructured Meshes

In this paper, we propose a high-order accurate discontinuous Galerkin (DG) method for the compressible Euler equations under gravitational fields on unstructured meshes. The scheme preserves a general hydrostatic equilibrium state and provably guarantees the positivity of density and pressure at the same time. Comparing with the work on the well-balanced scheme for Euler equations with gravitation on rectangular meshes, the extension to triangular meshes is conceptually plausible but highly nontrivial. We first introduce a special way to recover the equilibrium state and then design a group of novel variables at the interface of two adjacent cells, which plays an important role in the well-balanced and positivity-preserving properties. One main challenge is that the well-balanced schemes may not have the weak positivity property. In order to achieve the well-balanced and positivity-preserving properties simultaneously while maintaining high-order accuracy, we carefully design DG spatial discretization with well-balanced numerical fluxes and suitable source term approximation. For the ideal gas, we prove that the resulting well-balanced scheme, coupled with strong stability preserving time discretizations, satisfies a weak positivity property. A simple existing limiter can be applied to enforce the positivity-preserving property, without losing high-order accuracy and conservation. Extensive one- and two-dimensional numerical examples demonstrate the desired properties of the proposed scheme, as well as its high resolution and robustness.     

A High Order Central DG Method of the Two-Layer Shallow Water Equations

In this paper, we focus on the numerical simulation of the two-layer shallow water equations over variable bottom topography. Although the existing numerical schemes for the single-layer shallow water equations can be extended to two-layer shallow water equations, it is not a trivial work due to the complexity of the equations. To achieve the well-balanced property of the numerical scheme easily, the two-layer shallow water equations are reformulated into a new form by introducing two auxiliary variables. Since the new equations are only conditionally hyperbolic and their eigenstructure cannot be easily obtained, we consider the utilization of the central discontinuous Galerkin method which is free of Riemann solvers. By choosing the values of the auxiliary variables suitably, we can prove that the scheme can exactly preserve the still-water solution, and thus it is a truly well-balanced scheme. To ensure the non-negativity of the water depth, a positivity-preserving limiter and a special approximation to the bottom topography are employed. The accuracy and validity of the numerical method will be illustrated through some numerical tests.     

A Compact High Order Space-Time Method for Conservation Laws

This paper presents a novel high-order space-time method for hyperbolic conservation laws. Two important concepts, the staggered space-time mesh of the space-time conservation element/solution element (CE/SE) method and the local discontinuous basis functions of the space-time discontinuous Galerkin (DG) finite element method, are the two key ingredients of the new scheme. The staggered space-time mesh is constructed using the cell-vertex structure of the underlying spatial mesh. The universal definitions of CEs and SEs are independent of the underlying spatial mesh and thus suitable for arbitrarily unstructured meshes. The solution within each physical time step is updated alternately at the cell level and the vertex level. For this solution updating strategy and the DG ingredient, the new scheme here is termed as the discontinuous Galerkin cell-vertex scheme (DG-CVS). The high order of accuracy is achieved by employing high-order Taylor polynomials as the basis functions inside each SE. The present DG-CVS exhibits many advantageous features such as Riemann-solver-free, high-order accuracy, point-implicitness, compactness, and ease of handling boundary conditions. Several numerical tests including the scalar advection equations and compressible Euler equations will demonstrate the performance of the new method.     

A Compact Third-Order Gas-Kinetic Scheme for Compressible Euler and Navier-Stokes Equations

In this paper, a compact third-order gas-kinetic scheme is proposed for the compressible Euler and Navier-Stokes equations. The main reason for the feasibility to develop such a high-order scheme with compact stencil, which involves only neighboring cells, is due to the use of a high-order gas evolution model. Besides the evaluation of the time-dependent flux function across a cell interface, the high-order gas evolution model also provides an accurate time-dependent solution of the flow variables at a cell interface. Therefore, the current scheme not only updates the cell averaged conservative flow variables inside each control volume, but also tracks the flow variables at the cell interface at the next time level. As a result, with both cell averaged and cell interface values, the high-order reconstruction in the current scheme can be done compactly. Different from using a weak formulation for high-order accuracy in the Discontinuous Galerkin method, the current scheme is based on the strong solution, where the flow evolution starting from a piecewise discontinuous high-order initial data is precisely followed. The cell interface time-dependent flow variables can be used for the initial data reconstruction at the beginning of next time step. Even with compact stencil, the current scheme has third-order accuracy in the smooth flow regions, and has favorable shock capturing property in the discontinuous regions. It can be faithfully used from the incompressible limit to the hypersonic flow computations, and many test cases are used to validate the current scheme. In comparison with many other high-order schemes, the current method avoids the use of Gaussian points for the flux evaluation along the cell interface and the multi-stage Runge-Kutta time stepping technique. Due to its multidimensional property of including both derivatives of flow variables in the normal and tangential directions of a cell interface, the viscous flow solution, especially those with vortex structure, can be accurately captured. With the same stencil of a second order scheme, numerical tests demonstrate that the current scheme is as robust as well-developed second-order shock capturing schemes, but provides more accurate numerical solutions than the second order counterparts.     

A Slope Constrained 4th Order Multi-Moment Finite Volume Method with WENO Limiter

This paper presents a new and better suited formulation to implement the limiting projection to high-order schemes that make use of high-order local reconstructions for hyperbolic conservation laws. The scheme, so-called MCV-WENO4 (multi-moment Constrained finite Volume with WENO limiter of 4th order) method, is an extension of the MCV method of Ii & Xiao (2009) by adding the 1st order derivative (gradient or slope) at the cell center as an additional constraint for the cell-wise local reconstruction. The gradient is computed from a limiting projection using the WENO (weighted essentially non-oscillatory) reconstruction that is built from the nodal values at 5 solution points within 3 neighboring cells. Different from other existing methods where only the cell-average value is used in the WENO reconstruction, the present method takes account of the solution structure within each mesh cell, and thus minimizes the stencil for reconstruction. The resulting scheme has 4th-order accuracy and is of significant advantage in algorithmic simplicity and computational efficiency. Numerical results of one and two dimensional benchmark tests for scalar and Euler conservation laws are shown to verify the accuracy and oscillation-less property of the scheme.     

A High-Order Cell-Centered Discontinuous Galerkin Multi-Material Arbitrary Lagrangian-Eulerian Method

In this paper, a high-order cell-centered discontinuous Galerkin (DG) multi-material arbitrary Lagrangian-Eulerian (MMALE) method is developed for compressible fluid dynamics. The MMALE method utilizes moment-of-fluid (MOF) interface reconstruction technology to simulate multi-materials of immiscible fluids. It is an explicit time-marching Lagrangian plus remap type. In the Lagrangian phase, an updated high-order discontinuous Galerkin Lagrangian method is applied for the discretization of hydrodynamic equations, and Tipton's pressure relaxation closure model is used in the mixed cells. A robust moment-of-fluid interface reconstruction algorithm is used to provide the information of the material interfaces for remapping. In the rezoning phase, Knupp's algorithm is used for mesh smoothing. For the remapping phase, a high-order accurate remapping method of the cell-intersection-based type is proposed. It can be divided into four stages: polynomial reconstruction, polygon intersection, integration, and detection of problematic cells and limiting. Polygon intersection is based on the "clipping and projecting" algorithm, and detection of problematic cells depends on a troubled cell marker, and a posteriori multi-dimensional optimal order detection (MOOD) limiting strategy is used for limiting. Numerical tests are given to demonstrate the robustness and accuracy of our method.     

Comparison of Fifth-Order WENO Scheme and Finite Volume WENO-Gas-Kinetic Scheme for Inviscid and Viscous Flow Simulation

The development of high-order schemes has been mostly concentrated on the limiters and high-order reconstruction techniques. In this paper, the effect of the flux functions on the performance of high-order schemes will be studied. Based on the same WENO reconstruction, two schemes with different flux functions, i.e., the fifth-order WENO method and the WENO-Gas-Kinetic scheme (WENO-GKS), will be compared. The fifth-order finite difference WENO-SW scheme is a characteristic variable reconstruction based method which uses the Steger-Warming flux splitting for inviscid terms, the sixth-order central difference for viscous terms, and three stages Runge-Kutta time stepping for the time integration. On the other hand, the finite volume WENO-GKS is a conservative variable reconstruction based method with the same WENO reconstruction. But it evaluates a time dependent gas distribution function along a cell interface, and updates the flow variables inside each control volume by integrating the flux function along the boundary of the control volume in both space and time. In order to validate the robustness and accuracy of the schemes, both methods are tested under a wide range of flow conditions: vortex propagation, Mach 3 step problem, and the cavity flow at Reynolds number 3200. Our study shows that both WENO-SW and WENO-GKS yield quantitatively similar results and agree with each other very well provided a sufficient grid resolution is used. With the reduction of mesh points, the WENO-GKS behaves to have less numerical dissipation and present more accurate solutions than those from the WENO-SW in all test cases. For the Navier-Stokes equations, since the WENO-GKS couples inviscid and viscous terms in a single flux evaluation, and the WENO-SW uses an operator splitting technique, it appears that the WENO-SW is more sensitive to the WENO reconstruction and boundary treatment. In terms of efficiency, the finite volume WENO-GKS is about 4 times slower than the finite difference WENO-SW in two dimensional simulations. The current study clearly shows that besides high-order reconstruction, an accurate gas evolution model or flux function in a high-order scheme is also important in the capturing of physical solutions. In a physical flow, the transport, stress deformation, heat conduction, and viscous heating are all coupled in a single gas evolution process. Therefore, it is preferred to develop such a scheme with multi-dimensionality, and unified treatment of inviscid and dissipative terms. A high-order scheme does prefer a high-order gas evolution model. Even with the rapid advances of high-order reconstruction techniques, the first-order dynamics of the Riemann solution becomes the bottleneck for the further development of high-order schemes. In order to avoid the weakness of the low order flux function, the development of high-order schemes relies heavily on the weak solution of the original governing equations for the update of additional degree of freedom, such as the non-conservative gradients of flow variables, which cannot be physically valid in discontinuous regions.     

Three Discontinuous Galerkin Methods for One- and Two-Dimensional Nonlinear Dirac Equations with a Scalar Self-Interaction

This paper develops three high-order accurate discontinuous Galerkin (DG) methods for the one-dimensional (1D) and two-dimensional (2D) nonlinear Dirac (NLD) equations with a general scalar self-interaction. They are the Runge-Kutta DG (RKDG) method and the DG methods with the one-stage fourth-order Lax-Wendroff type time discretization (LWDG) and the two-stage fourth-order accurate time discretization (TSDG). The RKDG method uses the spatial DG approximation to discretize the NLD equations and then utilize the explicit multistage high-order Runge-Kutta time discretization for the first-order time derivatives, while the LWDG and TSDG methods, on the contrary, first give the one-stage fourth-order Lax-Wendroff type and the two-stage fourth-order time discretizations of the NLD equations, respectively, and then discretize the first- and higher-order spatial derivatives by using the spatial DG approximation. The $L^2$ stability of the 2D semi-discrete DG approximation is proved in the RKDG methods for a general triangulation, and the computational complexities of three 1D DG methods are estimated. Numerical experiments are conducted to validate the accuracy and the conservation properties of the proposed methods. The interactions of the solitary waves, the standing and travelling waves are investigated numerically and the 2D breathing pattern is observed.     

A Colocalized Scheme for Three-Temperature Grey Diffusion Radiation Hydrodynamics

A positivity-preserving, conservative and entropic numerical scheme is presented for the three-temperature grey diffusion radiation hydrodynamics model. More precisely, the dissipation matrices of the colocalized semi-Lagrangian scheme are defined in order to enforce the entropy production on each species (electron or ion) proportionally to its mass as prescribed in [34]. A reformulation of the model is then considered to enable the derivation of a robust convex combination based scheme. This yields the positivity-preserving property at each sub-iteration of the algorithm while the total energy conservation is reached at convergence. Numerous pure hydrodynamics and radiation hydrodynamics test cases are carried out to assess the accuracy of the method. The question of the stability of the scheme is also addressed. It is observed that the present numerical method is particularly robust.     

High-Order Gas-Kinetic Scheme in Curvilinear Coordinates for the Euler and Navier-Stokes Solutions

The high-order gas-kinetic scheme (HGKS) has achieved success in simulating compressible flows with Cartesian meshes. To study the flow problems in general geometries, such as the flow over a wing-body, the development of HGKS in general curvilinear coordinates becomes necessary. In this paper, a two-stage fourth-order gas-kinetic scheme is developed for the Euler and Navier-Stokes solutions in the curvilinear coordinates from one-dimensional to three-dimensional computations. Based on the coordinate transformation, the kinetic equation is transformed first to the computational space, and the flux function in the gas-kinetic scheme is obtained there and is transformed back to the physical domain for the update of flow variables inside each control volume. To achieve the expected order of accuracy, the dimension-by-dimension reconstruction based on the WENO scheme is adopted in the computational domain, where the reconstructed variables are the cell averaged Jacobian and the Jacobian-weighted conservative variables. In the two-stage fourth-order gas-kinetic scheme, the point values as well as the spatial derivatives of conservative variables at Gaussian quadrature points have to be used in the evaluation of the time dependent flux function. The point-wise conservative variables are obtained by ratio of the above reconstructed data, and the spatial derivatives are reconstructed through orthogonalization in physical space and chain rule. A variety of numerical examples from the accuracy tests to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme for both inviscid and viscous flows. The precise satisfaction of the geometrical conservation law in non-orthogonal mesh is also demonstrated through the numerical example.     

Unconditional Positivity-Preserving and Energy Stable Schemes for a Reduced Poisson-Nernst-Planck System

The Poisson-Nernst-Planck (PNP) system is a widely accepted model for simulation of ionic channels. In this paper, we design, analyze, and numerically validate a second order unconditional positivity-preserving scheme for solving a reduced PNP system, which can well approximate the three dimensional ion channel problem. Positivity of numerical solutions is proven to hold true independent of the size of time steps and the choice of the Poisson solver. The scheme is easy to implement without resorting to any iteration method. Several numerical examples further confirm the positivity-preserving property, and demonstrate the accuracy, efficiency, and robustness of the proposed scheme, as well as the fast approach to steady states.     

A Decoupled and Positivity-Preserving DDFVS Scheme for Diffusion Problems on Polyhedral Meshes

We propose a decoupled and positivity-preserving discrete duality finite volume (DDFV) scheme for anisotropic diffusion problems on polyhedral meshes with star-shaped cells and planar faces. Under the generalized DDFV framework, two sets of finite volume (FV) equations are respectively constructed on the dual and primary meshes, where the ones on the dual mesh are derived from the ingenious combination of a geometric relationship with the construction of the cell matrix. The resulting system on the dual mesh is symmetric and positive definite, while the one on the primary mesh possesses an M-matrix structure. To guarantee the positivity of the two categories of unknowns, a cutoff technique is introduced. As for the local conservation, it is conditionally maintained on the dual mesh while strictly preserved on the primary mesh. More interesting is that the FV equations on the dual mesh can be solved independently, so that the two sets of FV equations are decoupled. As a result, no nonlinear iteration is required for linear problems and a general nonlinear solver could be used for nonlinear problems. In addition, we analyze the well-posedness of numerical solutions for linear problems. The properties of the presented scheme are examined by numerical experiments. The efficiency of the Newton method is also demonstrated by comparison with those of the fixed-point iteration method and its Anderson acceleration.     

An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes

We propose an all regime Lagrange-Projection like numerical scheme for the gas dynamics equations. By all regime, we mean that the numerical scheme is able to compute accurate approximate solutions with an under-resolved discretization with respect to the Mach number M, i.e. such that the ratio between the Mach number M and the mesh size or the time step is small with respect to 1. The key idea is to decouple acoustic and transport phenomenon and then alter the numerical flux in the acoustic approximation to obtain a uniform truncation error in term of M. This modified scheme is conservative and endowed with good stability properties with respect to the positivity of the density and the internal energy. A discrete entropy inequality under a condition on the modification is obtained thanks to a reinterpretation of the modified scheme in the Harten Lax and van Leer formalism. A natural extension to multi-dimensional problems discretized over unstructured mesh is proposed. Then a simple and efficient semi-implicit scheme is also proposed. The resulting scheme is stable under a CFL condition driven by the (slow) material waves and not by the (fast) acoustic waves and so verifies the all regime property. Numerical evidences are proposed and show the ability of the scheme to deal with tests where the flow regime may vary from low to high Mach values.     

A Positivity-Preserving Second-Order BDF Scheme for the Cahn-Hilliard Equation with Variable Interfacial Parameters

We present and analyze a new second-order finite difference scheme for the Macromolecular Microsphere Composite hydrogel, Time-Dependent Ginzburg-Landau (MMC-TDGL) equation, a Cahn-Hilliard equation with Flory-Huggins-deGennes energy potential. This numerical scheme with unconditional energy stability is based on the Backward Differentiation Formula (BDF) method in time derivation combining with Douglas-Dupont regularization term. In addition, we present a pointwise bound of the numerical solution for the proposed scheme in the theoretical level. For the convergent analysis, we treat three nonlinear logarithmic terms as a whole and deal with all logarithmic terms directly by using the property that the nonlinear error inner product is always non-negative. Moreover, we present the detailed convergent analysis in $ℓ^∞$(0,$T$;$H_h^{-1}$)∩$ℓ^2$(0,$T$;$H_h^1$) norm. At last, we use the local Newton approximation and multigrid method to solve the nonlinear numerical scheme, and various numerical results are presented, including the numerical convergence test, positivity-preserving property test, spinodal decomposition, energy dissipation and mass conservation properties.     

An Efficient Positivity-Preserving Finite Volume Scheme for the Nonequilibrium Three-Temperature Radiation Diffusion Equations on Polygonal Meshes

This paper develops an efficient positivity-preserving finite volume scheme for the two-dimensional nonequilibrium three-temperature radiation diffusion equations on general polygonal meshes. The scheme is formed as a predictor-corrector algorithm. The corrector phase obtains the cell-centered solutions on the primary mesh, while the predictor phase determines the cell-vertex solutions on the dual mesh independently. Moreover, the flux on the primary edge is approximated with a fixed stencil and the nonnegative cell-vertex solutions are not reconstructed. Theoretically, our scheme does not require any nonlinear iteration for the linear problems, and can call the fast nonlinear solver (e.g. Newton method) for the nonlinear problems. The positivity, existence and uniqueness of the cell-centered solutions obtained on the corrector phase are analyzed, and the scheme on quasi-uniform meshes is proved to be $L^2$- and $H^1$-stable under some assumptions. Numerical experiments demonstrate the accuracy, efficiency and positivity of the scheme on various distorted meshes.     

Mortar DG Method with Staggered Hybridization for Rayleigh Waves Simulation

The simulation of Rayleigh waves is important in a variety of geophysical applications. The computational challenge is the fact that very fine mesh is necessary as the waves are concentrated at the free surface and decay exponentially away from the free surface. To overcome this challenge and to develop a robust high order scheme for the simulation of Rayleigh waves, we develop a mortar discontinuous Galerkin method with staggered hybridization. The use of the mortar technique allows one to use fine mesh in only a local region near the free surface, and use coarse mesh in most of the domain. This approach reduces the computational cost significantly. The staggered hybridization allows the preservation of the strong symmetry of the stress tensor without complicated construction of basis functions. In particular, the basis functions are piecewise polynomial without any continuity requirement, and the coupling of the basis functions is performed by using carefully chosen hybridized variables. The resulting scheme is explicit in time, and only local saddle point system is solved for each time step. We will present several benchmark problems to demonstrate the performance of the proposed method.     

A High-Order Discontinuous Galerkin Solver for Helically Symmetric Flows

We present a high-order discontinuous Galerkin (DG) scheme to solve the system of helically symmetric Navier-Stokes equations which are discussed in [28]. In particular, we discretize the helically reduced Navier-Stokes equations emerging from a reduction of the independent variables such that the remaining variables are: $t$, $r$, $ξ$ with $ξ=az+bϕ$, where $r$, $ϕ$, $z$ are common cylindrical coordinates and $t$ the time. Beside this, all three velocity components are kept non-zero. A new non-singular coordinate $η$ is introduced which ensures that a mapping of helical solutions into the three-dimensional space is well defined. Using that, periodicity conditions for the helical frame as well as uniqueness conditions at the centerline axis $r=0$ are derived. In the sector near the axis of the computational domain a change of the polynomial basis is implemented such that all physical quantities are uniquely defined at the centerline.For the temporal integration, we present a semi-explicit scheme of third order where the full spatial operator is split into a Stokes operator which is discretized implicitly and an operator for the nonlinear terms which is treated explicitly. Computations are conducted for a cylindrical shell, excluding the centerline axis, and for the full cylindrical domain, where the centerline is included. In all cases we obtain the convergence rates of order $\mathcal{O}(h^{k+1})$ that are expected from DG theory.In addition to the first DG discretization of the system of helically invariant Navier-Stokes equations, the treatment of the central axis, the resulting reduction of the DG space, and the simultaneous use of a semi-explicit time stepper are of particular novelty.     

A Computational Model for Simulation of Shallow Water Waves by Elastic Deformations in the Topography

We propose a coupled model to simulate shallow water waves induced by elastic deformations in the bed topography. The governing equations consist of the depth-averaged shallow water equations including friction terms for the water free-surface and the well-known second-order elastostatics formulation for the bed deformation. The perturbation on the free-surface is assumed to be caused by a sudden change in the bottom beds. At the interface between the water flow and the bed topography, transfer conditions are implemented. Here, the hydrostatic pressure and friction forces are considered for the elastostatic equations whereas bathymetric forces are accounted for in the shallow water equations. The focus in the present study is on the development of a simple and accurate representation of the interaction between water waves and bed deformations in order to simulate practical shallow water flows without relying on complex partial differential equations with free boundary conditions. The effects of location and magnitude of the deformation on the flow fields and free-surface waves are investigated in details. Numerical simulations are carried out for several test examples on shallow water waves induced by sudden changes in the bed. The proposed computational model has been found to be feasible and satisfactory.     

The Error-Minimization-Based Strategy for Moving Mesh Methods

The typical elements in a numerical simulation of fluid flow using moving meshes are a time integration scheme, a rezone method in which a new mesh is defined, and a remapping (conservative interpolation) in which a solution is transferred to the new mesh. The objective of the rezone method is to move the computational mesh to improve the robustness, accuracy and eventually efficiency of the simulation. In this paper, we consider the one-dimensional viscous Burgers' equation and describe a new rezone strategy which minimizes the L₂ norm of error and maintains mesh smoothness. The efficiency of the proposed method is demonstrated with numerical examples.     

A Moving Mesh Method for the Euler Flow Calculations Using a Directional Monitor Function

This paper is concerned with the adaptive grid method for computations of the Euler equations in fluid dynamics. The new feature of the present moving mesh algorithm is the use of a dimensional-splitting type monitor function, which is to increase grid concentration in regions containing shock waves and contact discontinuities or their interactions. Several two–dimensional flow problems are computed to demonstrate the effectiveness of the present adaptive grid algorithm.