On the characteristics-based ACM for incompressible flows

In this paper, the revised characteristics-based (CB) method for incompressible flows recently derived by Neofytou [P. Neofytou, Revision of the characteristic-based scheme for incompressible flows, J. Comput. Phys. 222 (2007) 475–484] has been further investigated. We have derived all the formulas for pressure and velocities from this revised CB method, which is based on the artificial compressibility method (ACM) [A.J. Chorin, A numerical solution for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967) 12]. Then we analyze the formulations of the original CB method [D. Drikakis, P.A. Govatsos, D.E. Papatonis, A characteristic based method for incompressible flows, Int. J. Numer. Meth. Fluids 19 (1994) 667–685; E. Shapiro, D. Drikakis, Non-conservative and conservative formulations of characteristics numerical reconstructions for incompressible flows, Int. J. Numer. Meth. Eng. 66 (2006) 1466–1482; D. Drikakis, P.K. Smolarkiewicz, On spurious vortical structures, J. Comput. Phys. 172 (2001) 309–325; F. Mallinger, D. Drikakis, Instability in three-dimensional, unsteady stenotic flows, Int. J. Heat Fluid Flow 23 (2002) 657–663; E. Shapiro, D. Drikakis, Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Parts I. Derivation of different formulations and constant density limit, J. Comput. Phys. 210 (2005) 584–607; Y. Zhao, B. Zhang, A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids, Comput. Meth. Appl. Mech. Eng. 190 (5–7) (2000) 733–756] to investigate their consistency with the governing flow equations after convergence has been achieved. Furthermore we have implemented both formulations in an unstructured-grid finite volume solver [Y. Zhao, B. Zhang, A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids, Comput. Meth. Appl. Mech. Eng. 190 (5–7) (2000) 733–756]. Detailed numerical experiments show that both methods give almost identical solutions and convergence rates. Both can generate solutions which agree well with published results and experimental measurements. We thus conclude that both methods, being upwind schemes designed for the ACM, have the same performances in terms of accuracy and convergence speed, even though the revised method is more complex with less stringent assumptions made, while the original CB method is simpler due to the use of extra simplifying assumptions.     

A direct Eulerian GRP scheme for relativistic hydrodynamics: Two-dimensional case

This paper develops a direct Eulerian generalized Riemann problem (GRP) scheme for two-dimensional (2D) relativistic hydrodynamics (RHD). It is an extension of the GRP scheme for one-dimensional (1D) RHDs [Z.C. Yang, P. He, H.Z. Tang, J. Comput. Phys. 230 (2011) 7964–7987] and the GRP scheme for the non-relativistic hydrodynamics [M. Ben-Artzi, J.Q. Li, G. Warnecke, J. Comput. Phys. 218 (2006) 19–43]. In order to derive the direct Eulerian GRP scheme, the (local) GRP of the split 2D RHD equations in the Eulerian formulation has to be directly resolved by using corresponding Riemann invariants and Rankine–Hugoniot jump conditions so that the crucial and delicate Lagrangian treatment in the original GRP scheme [M. Ben-Artzi, J. Falcovitz, J. Comput. Phys. 55 (1984) 1–32] may be avoided. An important difference of resolving the GRP of the split 2D RHD equations from the GRP of the 1D RHD equations or the non-relativistic hydrodynamical equations is coming from the fact that the flow regions across the shock or rarefaction wave in the GRP of the split 2D RHD equations are nonlinearly coupled through the Lorentz factor which is also built in terms of the tangential velocities. It is a purely multi-dimensional relativistic feature. Several numerical examples are given to demonstrate the accuracy and effectiveness of the proposed 2D GRP scheme.     

Compact finite volume schemes on boundary-fitted grids

The paper focuses on the development of a framework for high-order compact finite volume discretization of the three-dimensional scalar advection–diffusion equation. In order to deal with irregular domains, a coordinate transformation is applied between a curvilinear, non-orthogonal grid in the physical space and the computational space. Advective fluxes are computed by the fifth-order upwind scheme introduced by Pirozzoli [S. Pirozzoli, Conservative hybrid compact-WENO schemes for shock–turbulence interaction, J. Comp. Phys. 178 (2002) 81] while the Coupled Derivative scheme [M.H. Kobayashi, On a class of Padé finite volume methods, J. Comp. Phys. 156 (1999) 137] is used for the discretization of the diffusive fluxes.Numerical tests include unsteady diffusion over a distorted grid, linear free-surface gravity waves in a irregular domain and the advection of a scalar field. The proposed methodology attains high-order formal accuracy and shows very favorable resolution characteristics for the simulation of problems with a wide range of length scales.     

A high-order accurate hybrid scheme using a central flux scheme and a WENO scheme for compressible flowfield analysis

A high-order accurate hybrid central-WENO scheme is proposed. The fifth order WENO scheme [G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228] is divided into two parts, a central flux part and a numerical dissipation part, and is coupled with a central flux scheme. Two sub-schemes, the WENO scheme and the central flux scheme, are hybridized by means of a weighting function that indicates the local smoothness of the flowfields. The derived hybrid central-WENO scheme is written as a combination of the central flux scheme and the numerical dissipation of the fifth order WENO scheme, which is controlled adaptively by a weighting function. The structure of the proposed hybrid central-WENO scheme is similar to that of the YSD-type filter scheme [H.C. Yee, N.D. Sandham, M.J. Djomehri, Low-dissipative high-order shock-capturing methods using characteristic-based filters, J. Comput. Phys. 150 (1999) 199–238]. Therefore, the proposed hybrid scheme has also certain merits that the YSD-type filter scheme has. The accuracy and efficiency of the developed hybrid central-WENO scheme are investigated through numerical experiments on inviscid and viscous problems. Numerical results show that the proposed hybrid central-WENO scheme can resolve flow features extremely well.     

A fourth-order divergence-free method for MHD flows

This paper extends our previous third-order method [S. Li, High order central scheme on overlapping cells for magneto-hydrodynamic flows with and without constrained transport method, J. Comput. Phys. 227 (2008) 7368–7393] to the fourth-order. Central finite-volume schemes on overlapping grid are used for both the volume-averaged variables and the face-averaged magnetic field. The magnetic field at the cell boundaries falls within the dual grid and is naturally continuous so that our method eliminates the instability triggered by the discontinuity in the normal component of the magnetic field. Our fourth-order scheme has much smaller numerical dissipation than the third-order scheme. The divergence-free condition of the magnetic field is preserved by our fourth-order divergence-free reconstruction and the constrained transport method. Numerical examples show that the divergence-free condition is essential to the accuracy of the method when a limiter is used in the reconstruction. The high-order, low-dissipation, and divergence-free properties of this method make it an ideal tool for direct magneto-hydrodynamic turbulence simulations.     

A systematic methodology for constructing high-order energy stable WENO schemes

A third-order Energy Stable Weighted Essentially Non-Oscillatory (ESWENO) finite difference scheme developed by the authors of this paper [N.K. Yamaleev, M.H. Carpenter, Third-order energy stable WENO scheme, J. Comput. Phys. 228 (2009) 3025–3047] was proven to be stable in the energy norm for both continuous and discontinuous solutions of systems of linear hyperbolic equations. Herein, a systematic approach is presented that enables “energy stable” modifications for existing WENO schemes of any order. The technique is demonstrated by developing a one-parameter family of fifth-order upwind-biased ESWENO schemes including one sixth-order central scheme; ESWENO schemes up to eighth order are presented in the Appendix. We also develop new weight functions and derive constraints on their parameters, which provide consistency, much faster convergence of the high-order ESWENO schemes to their underlying linear schemes for smooth solutions with arbitrary number of vanishing derivatives, and better resolution near strong discontinuities than the conventional counterparts.     

An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws

In this article we develop an improved version of the classical fifth-order weighted essentially non-oscillatory finite difference scheme of [G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228] (WENO-JS) for hyperbolic conservation laws. Through the novel use of a linear combination of the low order smoothness indicators already present in the framework of WENO-JS, a new smoothness indicator of higher order is devised and new non-oscillatory weights are built, providing a new WENO scheme (WENO-Z) with less dissipation and higher resolution than the classical WENO. This new scheme generates solutions that are sharp as the ones of the mapped WENO scheme (WENO-M) of Henrick et al. [A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567], however with a 25% reduction in CPU costs, since no mapping is necessary. We also provide a detailed analysis of the convergence of the WENO-Z scheme at critical points of smooth solutions and show that the solution enhancements of WENO-Z and WENO-M at problems with shocks comes from their ability to assign substantially larger weights to discontinuous stencils than the WENO-JS scheme, not from their superior order of convergence at critical points. Numerical solutions of the linear advection of discontinuous functions and nonlinear hyperbolic conservation laws as the one dimensional Euler equations with Riemann initial value problems, the Mach 3 shock–density wave interaction and the blastwave problems are compared with the ones generated by the WENO-JS and WENO-M schemes. The good performance of the WENO-Z scheme is also demonstrated in the simulation of two dimensional problems as the shock–vortex interaction and a Mach 4.46 Richtmyer–Meshkov Instability (RMI) modeled via the two dimensional Euler equations.     

A direct Eulerian GRP scheme for relativistic hydrodynamics: One-dimensional case

The paper proposes a direct Eulerian generalized Riemann problem (GRP) scheme for one-dimensional relativistic hydrodynamics. It is an extension of the Eulerian GRP scheme for compressible non-relativistic hydrodynamics proposed in [M. Ben-Artzi, J.Q. Li, G. Warnecke, A direct Eulerian GRP scheme for compressible fluid flows, J. Comput. Phys. 218 (2006) 19–43]. Two main ingredients, the Riemann invariant and the Rankine–Hugoniot jump condition, are directly used to resolve the local GRP in the Eulerian formulation, and thus the crucial and delicate Lagrangian treatment in the original GRP scheme [3] can be avoided. Several numerical examples are given to demonstrate the accuracy and effectiveness of the proposed GRP scheme.     

Solvers for the high-order Riemann problem for hyperbolic balance laws

We study three methods for solving the Cauchy problem for a system of non-linear hyperbolic balance laws with initial condition consisting of two smooth vectors, with a discontinuity at the origin, a high-order Riemann problem. Two of the methods are new; one of the them results from a re-interpretation of the high-order numerical methods proposed by Harten et al. [A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high order accuracy essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303] and the other is a modification of the solver in [E.F. Toro, V.A. Titarev, Solution of the generalised Riemann problem for advection-reaction equations, Proc. Roy. Soc. London A 458 (2002) 271–281]. A systematic assessment of all three solvers is carried out and their relative merits are discussed. We also implement the solvers, locally, in the context of high-order finite volume numerical methods of the ADER type, on unstructured meshes. Schemes of up to fifth order of accuracy in space and time for the two-dimensional compressible Euler equations and the shallow water equations with source terms are constructed. Empirically obtained convergence rates are studied systematically and, for the tests considered, these correspond to the theoretically expected orders of accuracy. We also address the question of balance between flux gradients and source terms, for steady flow. We find that the ADER schemes may be termed asymptotically well-balanced, in the sense that the well-balanced property is attained as the order of the method increases, and this without introducing any ad-hoc fixes to the schemes or the equations.     

An improved particle correction procedure for the particle level set method

The particle level set method [D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys. 183 (2002) 83–116.] can substantially improve the mass conservation property of the level set method by using Lagrangian marker particles to correct the level set function in the under-resolved regions. In this study, the limitations of the particle level set method due to the errors introduced in the particle correction process are analyzed, and an improved particle correction procedure is developed based on a new interface reconstruction scheme. Moreover, the zero level set is “anchored” as the level set functions are reinitialized; hence the additional particle correction after the level set reinitialization is avoided. With this new scheme, a well-defined zero level set can be obtained and the disturbances to the interface are significantly reduced. Consequently, the particle reseeding operation will barely result in the loss of interface characteristics and can be applied as frequently as necessary. To demonstrate the accuracy and robustness of the proposed method, two extreme particle reseeding strategies, one without reseeding and the other with reseeding every time step, are applied in several benchmark advection tests and the results are compared with each other. Three interfacial flow cases, a 2D surface tension driven oscillating droplet, a 2D gas bubble rising in a quiescent liquid, and a 3D drop impact onto a liquid pool are simulated to illustrate the advantages of the current method over the level set and the original particle level set methods with regard to the smoothness of geometric properties and mass conservation in real physical applications.     

A conservative phase field method for solving incompressible two-phase flows

In this paper a conservative phase-field method based on the work of Sun and Beckermann [Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, J. Comput. Phys. 220 (2007) 626–653] for solving the two- and three-dimensional two-phase incompressible Navier–Stokes equations is proposed. The present method can preserve the total mass as the Cahn–Hilliard equation, but the calculation and implementation are much simpler than that. The dispersion-relation-preserving schemes are utilized for the advection terms while the Helmholtz smoother is applied to compute the surface-tension force term. To verify the proposed method, several benchmarks are examined and shown to have good agreements with previous results. It also shows that the satisfactions of mass conservations are guaranteed.     

A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids

A balanced force refined level set grid method for two-phase flows on structured and unstructured flow solver grids is presented. To accurately track the phase interface location, an auxiliary, high-resolution equidistant Cartesian grid is introduced. In conjunction with a dual-layer narrow band approach, this refined level set grid method allows for parallel, efficient grid convergence and error estimation studies of the interface tracking method. The Navier–Stokes equations are solved on an unstructured flow solver grid with a novel balanced force algorithm for level set methods based on the recently proposed method by Francois et al. [M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173] for volume of fluid methods on structured grids. To minimize spurious currents, a second order converging curvature evaluation technique for level set methods is presented. The results of several different test cases demonstrate the effectiveness of the proposed method, showing good mass conservation properties and second order converging spurious current magnitudes.     

Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws

In this article, we propose a new class of finite volume schemes of arbitrary accuracy in space and time for systems of hyperbolic balance laws with stiff source terms. The new class of schemes is based on a three stage procedure. First a high-order WENO reconstruction procedure is applied to the cell averages at the current time level. Second, the temporal evolution of the reconstruction polynomials is computed locally inside each cell using the governing equations. In the original ENO scheme of Harten et al. and in the ADER schemes of Titarev and Toro, this time evolution is achieved via a Taylor series expansion where the time derivatives are computed by repeated differentiation of the governing PDE with respect to space and time, i.e. by applying the so-called Cauchy–Kovalewski procedure. However, this approach is not able to handle stiff source terms. Therefore, we present a new strategy that only replaces the Cauchy–Kovalewski procedure compared to the previously mentioned schemes. For the time-evolution part of the algorithm, we introduce a local space–time discontinuous Galerkin (DG) finite element scheme that is able to handle also stiff source terms. This step is the only part of the algorithm which is locally implicit. The third and last step of the proposed ADER finite volume schemes consists of the standard explicit space–time integration over each control volume, using the local space–time DG solutions at the Gaussian integration points for the intercell fluxes and for the space–time integral over the source term. We will show numerical convergence studies for nonlinear systems in one space dimension with both non-stiff and with very stiff source terms up to sixth order of accuracy in space and time. The application of the new method to a large set of different test cases is shown, in particular the stiff scalar model problem of LeVeque and Yee [R.J. LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with stiff source terms, Journal of Computational Physics 86 (1) (1990) 187–210], the relaxation system of Jin and Xin [S. Jin, Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Communications on Pure and Applied Mathematics 48 (1995) 235–277] and the full compressible Euler equations with stiff friction source terms.     

A mass-conserved multiphase lattice Boltzmann method based on high-order difference

The Z–S–C multiphase lattice Boltzmann model [Zheng, Shu, and Chew(ZSC), J. Comput. Phys. 218, 353(2006)]is favored due to its good stability, high efficiency, and large density ratio. However, in terms of mass conservation, this model is not satisfactory during the simulation computations. In this paper, a mass correction is introduced into the ZSC model to make up the mass leakage, while a high-order difference is used to calculate the gradient of the order parameter to improve the accuracy. To verify the improved model, several three-dimensional multiphase flow simulations are carried out,including a bubble in a stationary flow, the merging of two bubbles, and the bubble rising under buoyancy. The numerical simulations show that the results from the present model are in good agreement with those from previous experiments and simulations. The present model not only retains the good properties of the original ZSC model, but also achieves the mass conservation and higher accuracy.     

A proof that a discrete delta function is second-order accurate

It is proved that a discrete delta function introduced by Smereka [P. Smereka, The numerical approximation of a delta function with application to level set methods, J. Comput. Phys. 211 (2006) 77–90] gives a second-order accurate quadrature rule for surface integrals using values on a regular background grid. The delta function is found using a technique of Mayo [A. Mayo, The fast solution of Poisson’s and the biharmonic equations on irregular regions, SIAM J. Numer. Anal. 21 (1984) 285–299]. It can be expressed naturally using a level set function.     

A lattice Boltzmann method for incompressible two-phase flows on partial wetting surface with large density ratio

This paper reports a new numerical scheme of the lattice Boltzmann method for calculating liquid droplet behaviour on particle wetting surfaces typically for the system of liquid–gas of a large density ratio. The method combines the existing models of Inamuro et al. [T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible two-phase flows with large density differences, J. Comput. Phys. 198 (2004) 628–644] and Briant et al. [A.J. Briant, P. Papatzacos, J.M. Yeomans, Lattice Boltzmann simulations of contact line motion in a liquid–gas system, Philos. Trans. Roy. Soc. London A 360 (2002) 485–495; A.J. Briant, A.J. Wagner, J.M. Yeomans, Lattice Boltzmann simulations of contact line motion: I. Liquid–gas systems. Phys. Rev. E 69 (2004) 031602; A.J. Briant, J.M. Yeomans, Lattice Boltzmann simulations of contact line motion: II. Binary fluids, Phys. Rev. E 69 (2004) 031603] and has developed novel treatment for partial wetting boundaries which involve droplets spreading on a hydrophobic surface combined with the surface of relative low contact angles and strips of relative high contact angles. The interaction between the fluid–fluid interface and the partial wetting wall has been typically considered. Applying the current method, the dynamics of liquid drops on uniform and heterogeneous wetting walls are simulated numerically. The results of the simulation agree well with those of theoretical prediction and show that the present LBM can be used as a reliable way to study fluidic control on heterogeneous surfaces and other wetting related subjects.     

Adaptive solution techniques for simulating underwater explosions and implosions

Adaptive solution techniques are presented for simulating underwater explosions and implosions. The liquid is assumed to be an adiabatic fluid and the solution in the gas is assumed to be uniform in space. The solution in water is integrated in time using a semi-implicit time discretization of the adiabatic Euler equations. Results are presented either using a non-conservative semi-implicit algorithm or a conservative semi-implicit algorithm. A semi-implicit algorithm allows one to compute with relatively large time steps compared to an explicit method. The interface solver is based on the coupled level set and volume-of-fluid method (CLSVOF) [M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comput. Phys. 187 (2003) 110–136; M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2000) 301–337]. Several underwater explosion and implosion test cases are presented to show the performances of our proposed techniques.     

A conservative level set method for contact line dynamics

A new model for simulating contact line dynamics is proposed. We apply the idea of driving contact-line movement by enforcing the equilibrium contact angle at the boundary, to the conservative level set method for incompressible two-phase flow [E. Olsson, G. Kreiss, A conservative level set method for two phase flow, J. Comput. Phys. 210 (2005) 225–246]. A modified reinitialization procedure provides a diffusive mechanism for contact-line movement, and results in a smooth transition of the interface near the contact line without explicit reconstruction of the interface. We are able to capture contact-line movement without loosing the conservation. Numerical simulations of capillary dominated flows in two space dimensions demonstrate that the model is able to capture contact line dynamics qualitatively correct.     

An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations

Numerical methods for solving the ideal magnetohydrodynamic (MHD) equations in more than one space dimension must either confront the challenge of controlling errors in the discrete divergence of the magnetic field, or else be faced with nonlinear numerical instabilities. One approach for controlling the discrete divergence is through a so-called constrained transport method, which is based on first predicting a magnetic field through a standard finite volume solver, and then correcting this field through the appropriate use of a magnetic vector potential. In this work we develop a constrained transport method for the 3D ideal MHD equations that is based on a high-resolution wave propagation scheme. Our proposed scheme is the 3D extension of the 2D scheme developed by Rossmanith [J.A. Rossmanith, An unstaggered, high-resolution constrained transport method for magnetohydrodynamic flows, SIAM J. Sci. Comput. 28 (2006) 1766], and is based on the high-resolution wave propagation method of Langseth and LeVeque [J.O. Langseth, R.J. LeVeque, A wave propagation method for threedimensional hyperbolic conservation laws, J. Comput. Phys. 165 (2000) 126]. In particular, in our extension we take great care to maintain the three most important properties of the 2D scheme: (1) all quantities, including all components of the magnetic field and magnetic potential, are treated as cell-centered; (2) we develop a high-resolution wave propagation scheme for evolving the magnetic potential; and (3) we develop a wave limiting approach that is applied during the vector potential evolution, which controls unphysical oscillations in the magnetic field. One of the key numerical difficulties that is novel to 3D is that the transport equation that must be solved for the magnetic vector potential is only weakly hyperbolic. In presenting our numerical algorithm we describe how to numerically handle this problem of weak hyperbolicity, as well as how to choose an appropriate gauge condition. The resulting scheme is applied to several numerical test cases.     

Exact boundary condition for time-dependent wave equation based on boundary integral

An exact non-reflecting boundary conditions based on a boundary integral equation or a modified Kirchhoff-type formula is derived for exterior three-dimensional wave equations. The Kirchhoff-type non-reflecting boundary condition is originally proposed by L. Ting and M.J. Miksis [J. Acoust. Soc. Am. 80 (1986) 1825] and numerically tested by D. Givoli and D. Cohen [J. Comput. Phys. 117 (1995) 102] for a spherically symmetric problem. The computational advantage of Ting–Miksis boundary condition is that its temporal non-locality is limited to a fixed amount of past information. However, a long-time instability is exhibited in testing numerical solutions by using a standard non-dissipative finite-difference scheme. The main purpose of this work is to present a new exact boundary condition and to eliminate the long-time instability. The proposed exact boundary condition can be considered as a limit case of Ting–Miksis boundary condition when the two artificial boundaries used in their method approach each other. Our boundary condition is actually a boundary integral equation on a single artificial boundary for wave equations, which is to be solved in conjunction with the interior wave equation. The new boundary condition needs only one artificial boundary, which can be of any shape, i.e., sphere, cubic surface, etc. It keeps all merits of the original Kirchhoff boundary condition such as restricting the temporal non-locality, free of numerical evaluation of any special functions and so on. Numerical approximation to the artificial boundary condition on cubic surface is derived and three-dimensional numerical tests are carried out on the cubic computational domain.