An analysis and comparison of conservative upwind difference schemes for ideal and non-ideal gases

In a recent paper [P. Glaister, Conservative upwind difference schemes for the Euler equations, Comput. Math. Appl. 45 (2003) 1673–1682] a number of numerical schemes were presented for the Euler equations governing compressible flows of an ideal gas, the principal one of which is based on a conservative linearisation approach. This scheme was subsequently extended to encompass compressible flows of real gases where the equation of state allows for non-ideal gases [P. Glaister, Conservative upwind difference schemes for compressible flows of a real gas, Comput. Math. Appl. 48 (2004) 469–480]. These schemes use different parameter vectors in their construction and, consequently, the scheme in [P. Glaister, Conservative upwind difference schemes for compressible flows of a real gas, Comput. Math. Appl. 48 (2004) 469–480] when applied to the special case of an ideal gas is not identical to the principal ideal gas scheme in [P. Glaister, Conservative upwind difference schemes for the Euler equations, Comput. Math. Appl. 45 (2003) 1673–1682]. In this paper it is shown how these schemes are related, followed by a numerical comparison when each is applied to two standard test problems.     

Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method II: Two dimensional case

A class of fifth-order weighted essentially non-oscillatory (WENO) schemes based on Hermite polynomials, termed HWENO (Hermite WENO) schemes, for solving one dimensional non-linear hyperbolic conservation law systems, was developed and applied as limiters for the Runge-Kutta discontinuous Galerkin (RKDG) methods in [J. Comput. Phys. 193 (2003) 115]. In this paper, we extend the method to solve two dimensional non-linear hyperbolic conservation law systems. The emphasis is again on the application of such HWENO finite volume methodology as limiters for RKDG methods to maintain compactness of RKDG methods. Numerical experiments for two dimensional Burgers’ equation and Euler equations of compressible gas dynamics are presented to show the effectiveness of these methods.     

Computation of incompressible bubble dynamics with a stabilized finite element level set method

This paper presents a stabilized finite element method for the three dimensional computation of incompressible bubble dynamics using a level set method. The interface between the two phases is resolved using the level set approach developed by Sethian [Level Set Methods and Fast Marching Methods, Cambridge University Press, 1999], Sussman et al. [J. Comput. Phys. 114 (1994) 146], and Sussman et al. [J. Comput. Phys. 148 (1999) 81–124]. In this approach the interface is represented as a zero level set of a smooth function. The streamline-upwind/Petrov–Galerkin method was used to discretize the governing flow and level set equations. The continuum surface force (CSF) model proposed by Brackbill et al. [J. Comput. Phys. 100 (1992) 335–354] was applied in order to account for surface tension effects. To restrict the interface from moving while re-distancing, an improved re-distancing scheme proposed in the finite difference context [J. Comput. Phys. 148 (1999) 81–124] is adapted for finite element discretization. This enables us to accurately compute the flows with large density and viscosity differences, as well as surface tension. The capability of the resultant algorithm is demonstrated with two and three dimensional numerical examples of a single bubble rising through a quiescent liquid, and two bubble coalescence.     

Hybrid Finite Difference Weighted Essentially Non-oscillatory Schemes for the Compressible Ideal Magnetohydrodynamics Equation

In this paper, we present hybrid weighted essentially non-oscillatory (WENO) schemes with several discontinuity detectors for solving the compressible ideal magnetohydrodynamics (MHD) equation. Li and Qiu (J Comput Phys 229:8105–8129, 2010) examined effectiveness and efficiency of several different troubled-cell indicators in hybrid WENO methods for Euler gasdynamics. Later, Li et al. (J Sci Comput 51:527–559, 2012) extended the hybrid methods for solving the shallow water equations with four better indicators. Hybrid WENO schemes reduce the computational costs, maintain non-oscillatory properties and keep sharp transitions for problems. The numerical results of hybrid WENO-JS/WENO-M schemes are presented to compare the ability of several troubled-cell indicators with a variety of test problems. The focus of this paper, we propose optimal and reliable indicators for performance comparison of hybrid method using troubled-cell indicators for efficient numerical method of ideal MHD equations. We propose a modified ATV indicator that uses a second derivative. It is advantageous for differential discontinuity detection such as jump discontinuity and kink. A detailed numerical study of one-dimensional and two-dimensional cases is conducted to address efficiency (CPU time reduction and more accurate numerical solution) and non-oscillatory property problems. We demonstrate that the hybrid WENO-M scheme preserves the advantages of WENO-M and the ratio of computational costs of hybrid WENO-M and hybrid WENO-JS is smaller than that of WENO-M and WENO-JS.     

An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows

This paper presents the latest developments of a discontinuous Galerkin (DG) method for incompressible flows introduced in [Bassi F, Crivellini A, Di Pietro DA, Rebay S. An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier–Stokes equations. J Comput Phys 2006;218(2):794–815] for the steady Navier–Stokes equations and extended in [Bassi F, Crivellini A. A high-order discontinuous Galerkin method for natural convection problems. In: Wesseling P, Oñate E, Periaux J, editors. Electronic proceedings of the ECCOMAS CFD 2006 conference, Egmond aan Zee, The Netherlands, September 5–8; 2006. TU Delft] to the coupled Navier–Stokes and energy equations governing natural convection flows.

The method is fully implicit and applies to the governing equations in primitive variable form. Its distinguishing feature is the formulation of the inviscid interface flux, which is based on the solution of local Riemann problems associated with the artificial compressibility perturbation of the Euler equations. The tight coupling between pressure and velocity so introduced stabilizes the method and allows using equal-order approximation spaces for both pressure and velocity. Since, independently of the amount of artificial compressibility added, the interface flux reduces to the physical one for vanishing interface jumps, the resulting method is strongly consistent.

In this paper, we present a review of the method together with two recently developed issues: (i) the high-order DG discretization of the incompressible Euler equations; (ii) the high-order implicit time integration of unsteady flows. The accuracy and versatility of the method are demonstrated by a suite of computations of steady and unsteady, inviscid and viscous incompressible flows.     

Stabilized discontinuous Galerkin method for hyperbolic equations

In this work a new stabilization technique is proposed and studied for the discontinuous Galerkin method applied to hyperbolic equations. In order to avoid the use of slope limiters, a streamline diffusion-like term is added to control oscillations for arbitrary element orders. Thus, the scheme combines ideas from both the Runge-Kutta discontinuous Galerkin method [J. Scient. Comput. 16 (2001) 173] and the streamline diffusion method [Comput. Methods Appl. Mech. Engrg. 32 (1982)]. To increase the stability range of the method, the diffusion term is treated implicitly. The result is a scheme with higher order in space with the same stability range as the finite volume method. An optimal relation between the time step and the size of the diffusion coefficient is analyzed for numerical precision. The scheme is implemented using the object oriented programming philosophy based on the environment described in [Comput. Methods Appl. Mech. Engrg. 150 (1997)]. Accuracy and shock capturing abilities of the method are analyzed in terms of two bidimensional model problems: the rotating cone and the backward facing step problem for the Euler equations of gas dynamics.     

Spatially adaptive techniques for level set methods and incompressible flow

Since the seminal work of [Sussman, M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] on coupling the level set method of [Osher S, Sethian J. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49] to the equations for two-phase incompressible flow, there has been a great deal of interest in this area. That work demonstrated the most powerful aspects of the level set method, i.e. automatic handling of topological changes such as merging and pinching, as well as robust geometric information such as normals and curvature. Interestingly, this work also demonstrated the largest weakness of the level set method, i.e. mass or information loss characteristic of most Eulerian capturing techniques. In fact, [Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] introduced a partial differential equation for battling this weakness, without which their work would not have been possible. In this paper, we discuss both historical and most recent works focused on improving the computational accuracy of the level set method focusing in part on applications related to incompressible flow due to both of its popularity and stringent accuracy requirements. Thus, we discuss higher order accurate numerical methods such as Hamilton–Jacobi WENO [Jiang G-S, Peng D. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM J Sci Comput 2000;21:2126–43], methods for maintaining a signed distance function, hybrid methods such as the particle level set method [Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183:83–116] and the coupled level set volume of fluid method [Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. J Comput Phys 2000;162:301–37], and adaptive gridding techniques such as the octree approach to free surface flows proposed in [Losasso F, Gibou F, Fedkiw R. Simulating water and smoke with an octree data structure, ACM Trans Graph (SIGGRAPH Proc) 2004;23:457–62].     

A new 3D parallel SPH scheme for free surface flows

We propose a new robust and accurate SPH scheme, able to track correctly complex three-dimensional non-hydrostatic free surface flows and, even more important, also able to compute an accurate and little oscillatory pressure field. It uses the explicit third order TVD Runge-Kutta scheme in time, following Shu and Osher [Shu C-W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys 1988;89:439-71], together with the new key idea of introducing a monotone upwind flux for the density equation, thus removing any artificial viscosity term. For the discretization of the velocity equation, the non-diffusive central flux has been used. A new flexible approach to impose the boundary conditions at solid walls is also proposed. It can handle any moving rigid body with arbitrarily irregular geometry. It does neither produce oscillations in the fluid pressure in proximity of the interfaces, nor does it have a restrictive impact on the stability condition of the explicit time stepping method, unlike the repellent boundary forces of Monaghan [Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399-406]. To asses the accuracy of the new SPH scheme, a 3D mesh-convergence study is performed for the strongly deforming free surface in a 3D dam-break and impact-wave test problem providing very good results.Moreover, the parallelization of the new 3D SPH scheme has been carried out using the message passing interface (MPI) standard, together with a dynamic load balancing strategy to improve the computational efficiency of the scheme. Thus, simulations involving millions of particles can be run on modern massively parallel supercomputers, obtaining a very good performance, as confirmed by a speed-up analysis. The 3D applications consist of environmental flow problems, such as dam-break flows and impact flows against a wall. The numerical solutions obtained with our new 3D SPH code have been compared with either experimental results or with other numerical reference solutions, obtaining in all cases a very satisfactory agreement.     

On improving mass-conservation properties of the hybrid particle-level-set method

For the computation of multi-phase flows level-set methods are an attractive alternative to volume-of-fluid or front-tracking approaches. For improving their accuracy and efficiency the hybrid particle-level-set modification was proposed by Enright et al. [Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle-level-set method for improved interface capturing. J Comput Phys 2002;183:83-116]. In actual applications the overall properties of a level-set method, such as mass conservation, are strongly affected by discretization schemes and algorithmic details. In this paper we address these issues with the objective of determining the optimum alternatives for the purpose of direct numerical simulation of dispersed-droplet flows. We evaluate different discretization schemes for curvature and unit normal vector at the interface. Another issue is the particular formulation of the reinitialization of the level-set function which significantly affects the quality of computational results. Different approaches employing higher-order schemes for discretization, supplemented either by a correction step using marker particles (Enright et al., 2002) or by additional constraints [Sussman M, Almgren AS, Bell JB, Colella P, Howell LH, Welcome ML. An adaptive level set approach for incompressible two-phase flows. J Comput Phys 1999;148:81-124] are analyzed. Different parameter choices for the hybrid particle-level-set method are evaluated with the purpose of increasing the efficiency of the method. Aiming at large-scale computations we find that in comparison with pure level-set methods the hybrid particle-level-set method exhibits better mass-conservation properties, especially in the case of marginally resolved interfaces.     

Fast algorithms for high-order numerical methods for space-fractional diffusion equations

In this paper, fast numerical methods for solving space-fractional diffusion equations are studied in two stages. Firstly, a fast direct solver for an implicit finite difference scheme proposed by Hao et al. [A fourth-order approximation of fractional derivatives with its applications, J. Comput. Phys. 281 (2015), pp. 787–805], which is fourth-order accurate in space and second-order accurate in time, is developed based on a circulant-and-skew-circulant (CS) representation of Toeplitz matrix inversion. Secondly, boundary value method with spatial discretization of Hao et al. [A fourth-order approximation of fractional derivatives with its applications, J. Comput. Phys. 281 (2015), pp. 787–805] is adopted to produce a numerical solution with higher order accuracy in time. Particularly, a method with fourth-order accuracy in both space and time can be achieved. GMRES method is employed for solving the discretized linear system with two preconditioners. Based on the CS representation of Toeplitz matrix inversion, the two preconditioners can be applied efficiently, and the convergence rate of the preconditioned GMRES method is proven to be fast. Numerical examples are given to support the theoretical analysis.     

Numerical study of a nonlinear heat equation for plasma physics

This paper is devoted to the numerical approximation of a nonlinear parabolic balance equation, which describes the heat evolution of a magnetically confined plasma in the edge region of a tokamak. The nonlinearity implies some numerical difficulties, in particular for the long-time behaviour approximation, when solved with standard methods. An efficient numerical scheme is presented in this paper, based on a combination of a directional splitting scheme and the implicit–explicit scheme introduced in Filbet and Jin [A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources, J. Comput. Phys. 229 (2010), pp. 7625–7648].     

All Mach Number Second Order Semi-implicit Scheme for the Euler Equations of Gas Dynamics

This paper presents an asymptotic preserving (AP) all Mach number finite volume shock capturing method for the numerical solution of compressible Euler equations of gas dynamics. Both isentropic and full Euler equations are considered. The equations are discretized on a staggered grid. This simplifies flux computation and guarantees a natural central discretization in the low Mach limit, thus dramatically reducing the excessive numerical diffusion of upwind discretizations. Furthermore, second order accuracy in space is automatically guaranteed. For the time discretization we adopt an Semi-IMplicit/EXplicit (S-IMEX) discretization getting an elliptic equation for the pressure in the isentropic case and for the energy in the full Euler case. Such equations can be solved linearly so that we do not need any iterative solver thus reducing computational cost. Second order in time is obtained by a suitable S-IMEX strategy taken from Boscarino et al. (J Sci Comput 68:975–1001, 2016). Moreover, the CFL stability condition is independent of the Mach number and depends essentially on the fluid velocity. Numerical tests are displayed in one and two dimensions to demonstrate performance of our scheme in both compressible and incompressible regimes.     

A Conservative Method for the Simulation of?the?Isothermal Euler System with the van-der-Waals Equation of State

In this article, we are interested in the simulation of phase transition in compressible flows, with the isothermal Euler system, closed by the van-der-Waals model. We formulate the problem as an hyperbolic system, with a source term located at the interface between liquid and vapour. The numerical scheme is based on (Abgrall and Saurel, J. Comput. Phys. 186(2):361?C396, 2003; Le Métayer et al., J. Comput. Phys. 205(2):567?C610, 2005). Compared with previous discretizations of the van-der-Waals system, the novelty of this algorithm is that it is fully conservative. Its Godunov-type formulation allows an easy implementation on multi-dimensional unstructured meshes.     

Uniformly Effective Numerical Methods for Hyperbolic Systems

A unified method to compute compressible and incompressible flows is presented. Accuracy and efficiency do not degrade as the Mach number tends to zero. A staggered scheme solved with a pressure correction method is used. The equation of state is arbitrary. A Riemann problem for the barotropic Euler equations with nonconvex equation of state is solved exactly and numericaly. A hydrodynamic flow with cavitation in which the Mach number varies between 10⁻³ and 20 is computed. Unified methods for compressible and incompressible flows are further discussed for the flow of a perfect gas. The staggered scheme with pressure correction is found to have Mach-uniform accuracy and efficiency, and for the fully compressible case the accuracy is comparable with that of established schemes for compressible flows. Received October 20, 1999; revised May 26, 2000     

Maximum-Principle-Satisfying and Positivity-Preserving High Order Discontinuous Galerkin Schemes for?Conservation Laws on Triangular Meshes

In Zhang and Shu (J. Comput. Phys. 229:3091–3120, 2010), two of the authors constructed uniformly high order accurate finite volume and discontinuous Galerkin (DG) schemes satisfying a strict maximum principle for scalar conservation laws on rectangular meshes. The technique is generalized to positivity preserving (of density and pressure) high order DG or finite volume schemes for compressible Euler equations in Zhang and Shu (J. Comput. Phys. 229:8918–8934, 2010). The extension of these schemes to triangular meshes is conceptually plausible but highly nontrivial. In this paper, we first introduce a special quadrature rule which is exact for two-variable polynomials over a triangle of a given degree and satisfy a few other conditions, by which we can construct high order maximum principle satisfying finite volume schemes (e.g. essentially non-oscillatory (ENO) or weighted ENO (WENO) schemes) or DG method solving two dimensional scalar conservation laws on triangular meshes. The same method can preserve the maximum principle for DG or finite volume schemes solving two-dimensional incompressible Euler equations in the vorticity stream-function formulation, or any passive convection equation with an incompressible velocity field. We also obtain positivity preserving (for density and pressure) high order DG or finite volume schemes solving compressible Euler equations on triangular meshes. Numerical tests for the third order Runge-Kutta DG (RKDG) method on unstructured meshes are reported.     

The BR1 Scheme is Stable for the Compressible Navier–Stokes Equations

In this work we prove that the original (Bassi and Rebay in J Comput Phys 131:267–279, 1997) scheme (BR1) for the discretization of second order viscous terms within the discontinuous Galerkin collocation spectral element method (DGSEM) with Gauss Lobatto nodes is stable. More precisely, we prove in the first part that the BR1 scheme preserves energy stability of the skew-symmetric advection term DGSEM discretization for the linearized compressible Navier–Stokes equations (NSE). In the second part, we prove that the BR1 scheme preserves the entropy stability of the recently developed entropy stable compressible Euler DGSEM discretization of Carpenter et al. (SIAM J Sci Comput 36:B835–B867, 2014) for the non-linear compressible NSE, provided that the auxiliary gradient equations use the entropy variables. Both parts are presented for fully three-dimensional, unstructured curvilinear hexahedral grids. Although the focus of this work is on the BR1 scheme, we show that the proof naturally includes the Local DG scheme of Cockburn and Shu.     

Flux-based level-set method for two-phase flows on unstructured grids

In this paper we deal with the application of the flux-based level set method to moving interface computations on unstructured grids. The focus lies on the overcoming of the known difficulties of level set methods, e.g. accurate computations of important geometric properties, reliable and precise reinitialization of the level set function and the adaption of standard discretization methods to the moving boundary case. The basic building block of our approach is the high-resolution flux-based level set method for general advection equation (Frolkovi? and Mikula in SIAM J Sci Comput 29(2):579–597, 2007, Frolkovi? and Wehner in Comput Vis Sci 12(6):626–650, 2009). We extend this method for the problem of reinitialization of the level set function on unstructured grids by using quadratic interpolation to compute distances for nodes close to the interface. To realize numerical simulation for some applications with moving boundaries, we adapt the approach of ghost fluid method (Gibou and Fedkiw in J Comput Phys 202:577–601, 2005) for unstructured grids. The idea is to describe the development of the moving boundary with a level set formulation while the computational grid remains fixed and the boundary conditions are enforced using some extrapolation. Our main motivation is the numerical solution of two-phase incompressible flow problems. Additionally to previously mentioned steps, we introduce further numerical schemes in the framework of finite volume discretization for the flow. Possible jumps of the pressure and the directional derivative of velocity at the interface are modeled directly within the method using the approach of extended approximation spaces. Besides that, an algorithm for the computations of curvature is considered that exhibits the second order accuracy for some examples. Numerical experiments are provided for the presented methods.