Artificial compressibility, characteristics-based schemes for variable-density, incompressible, multispecies flows: Part II. Multigrid implementation and numerical tests

The paper presents an investigation of the accuracy and efficiency of artificial compressibility, characteristics-based (CB) schemes for variable-density incompressible flows. The CB schemes have been implemented in conjunction with a multigrid method for accelerating numerical convergence and a fourth-order, explicit Runge–Kutta method for the integration of the governing equations in time. The implementation of the CB schemes is obtained in conjunction with first-, second- and third-order interpolation formulas for calculating the variables at the cell faces of the computational volume. The accuracy and efficiency of the schemes are examined against analytical and experimental results for diffusion broadening in two- and three-dimensional microfluidic channels, a problem that has motivated the development of the present methods. Moreover, unsteady, inviscid simulations have been performed for variable-density mixing layer. The computations revealed that accuracy and efficiency depend on the CB scheme design. The best multigrid convergence rates were exhibited by the conservative CB scheme, which is obtained by the fully conservative formulation of the variable-density, incompressible equations.     

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.     

L-stability of vertex-based MUSCL finite volume schemes on unstructured grids: Simulation of incompressible flows with high density ratios

This work is devoted to the design of multi-dimensional finite volume schemes for solving transport equations on unstructured grids. In the framework of MUSCL vertex-based methods we construct numerical fluxes such that the local maximum property is guaranteed under an explicit Courant–Friedrichs–Levy condition. The method can be naturally completed by adaptive local mesh refinements and it turns out that the mesh generation is less constrained than when using the competitive cell-centered methods. We illustrate the effectiveness of the scheme by simulating variable density incompressible viscous flows. Numerical simulations underline the theoretical predictions and succeed in the computation of high density ratio phenomena such as a water bubble falling in air.     

A consistent and conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part III: On a staggered mesh

The consistent and conservative scheme developed on a rectangular collocated mesh [M.-J. Ni, R. Munipalli, N.B. Morley, P. Huang, M.A. Abdou, A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: on a rectangular collocated grid system, Journal of Computational Physics 227 (2007) 174–204] and on an arbitrary collocated mesh [M.-J. Ni, R. Munipalli, P. Huang, N.B. Morley, M.A. Abdou, A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part II: on an arbitrary collocated mesh, Journal of Computational Physics 227 (2007) 205–228] has been extended and specially designed for calculation of the Lorentz force on a staggered grid system (Part III) by solving the electrical potential equation for magnetohydrodynamics (MHD) at a low magnetic Reynolds number. In a staggered mesh, pressure (p) and electrical potential (φ) are located in the cell center, while velocities and current fluxes are located on the cell faces of a main control volume. The scheme numerically meets the physical conservation laws, charge conservation law and momentum conservation law. Physically, the Lorentz force conserves the momentum when the magnetic field is constant or spatial coordinate independent. The calculation of current density fluxes on cell faces is conducted using a scheme consistent with the discretization for solution of the electrical potential Poisson equation, which can ensure the calculated current density conserves the charge. A divergence formula of the Lorentz force is used to calculate the Lorentz force at the cell center of a main control volume, which can numerically conserve the momentum at constant or spatial coordinate independent magnetic field. The calculated cell-center Lorentz forces are then interpolated to the cell faces, which are used to obtain the corresponding velocity fluxes by solving the momentum equations. The “conservative” is an important property of the scheme, which can guarantee computational accuracy of MHD flows at high Hartmann number with a strongly non-uniform mesh employed to resolve the Hartmann layers and side layers. 2D fully developed MHD flows with analytical solutions available have been conducted to validate the scheme at a staggered mesh. 3D MHD flows, with the experimental data available, at a constant magnetic field in a rectangular duct with sudden expansion and at a varying magnetic field in a rectangular duct are conducted on a staggered mesh to verify the computational accuracy of the scheme. It is expected that the scheme for the Lorentz force can be employed together with a fully conservative scheme for the convective term and the pressure term [Y. Morinishi, T.S. Lund, O.V. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, Journal of Computational Physics 143 (1998) 90–124] for direct simulation of MHD turbulence and MHD instability with good accuracy at a staggered mesh.     

Large eddy simulation using a new set of sixth order schemes for compressible viscous terms

A set of conservative sixth order central differencing schemes is suggested for compressible flows with variable viscosity coefficient. This new set of central differencing schemes has the stencil width matching that of the seventh order weighted essentially non-oscillatory scheme (WENO). This feature is important to maintain the compactness of the seventh order WENO scheme and facilitate boundary condition treatment. As an application example, a large eddy simulation (LES) is conducted for a cylinder flow using the seventh order WENO scheme for the convective terms and the new set of sixth order central differencing scheme for the viscous terms. The results are compared with those from other research groups and those obtained using the fifth order WENO scheme and fourth order central differencing.     

Velocity–vorticity–helicity formulation and a solver for the Navier–Stokes equations

For the three-dimensional incompressible Navier–Stokes equations, we present a formulation featuring velocity, vorticity and helical density as independent variables. We find the helical density can be observed as a Lagrange multiplier corresponding to the divergence-free constraint on the vorticity variable, similar to the pressure in the case of the incompressibility condition for velocity. As one possible practical application of this new formulation, we consider a time-splitting numerical scheme based on an alternating procedure between vorticity–helical density and velocity–Bernoulli pressure systems of equations. Results of numerical experiments include a comparison with some well-known schemes based on pressure–velocity formulation and illustrate the competitiveness on the new scheme as well as the soundness of the new formulation.     

High-order non-uniform grid schemes for numerical simulation of hypersonic boundary-layer stability and transition

The direct numerical simulation of receptivity, instability and transition of hypersonic boundary layers requires high-order accurate schemes because lower-order schemes do not have an adequate accuracy level to compute the large range of time and length scales in such flow fields. The main limiting factor in the application of high-order schemes to practical boundary-layer flow problems is the numerical instability of high-order boundary closure schemes on the wall. This paper presents a family of high-order non-uniform grid finite difference schemes with stable boundary closures for the direct numerical simulation of hypersonic boundary-layer transition. By using an appropriate grid stretching, and clustering grid points near the boundary, high-order schemes with stable boundary closures can be obtained. The order of the schemes ranges from first-order at the lowest, to the global spectral collocation method at the highest. The accuracy and stability of the new high-order numerical schemes is tested by numerical simulations of the linear wave equation and two-dimensional incompressible flat plate boundary layer flows. The high-order non-uniform-grid schemes (up to the 11th-order) are subsequently applied for the simulation of the receptivity of a hypersonic boundary layer to free stream disturbances over a blunt leading edge. The steady and unsteady results show that the new high-order schemes are stable and are able to produce high accuracy for computations of the nonlinear two-dimensional Navier–Stokes equations for the wall bounded supersonic flow.     

Balanced finite volume WENO and central WENO schemes for the shallow water and the open-channel flow equations

The goal of this work is to extend finite volume WENO and central WENO schemes to the hyperbolic balance laws with geometrical source term and spatially variable flux function. In particular, we apply proposed schemes to the shallow water and the open-channel flow equations where the source term depends on the channel geometry. For obtaining stable numerical schemes that are free of spurious oscillations, it becomes crucial to use the decomposed source term evaluation, which maintains the balancing between the flux gradient and the source term. In addition, the open-channel flow equations contain spatially variable flux function. The appropriate definitions of the terms that arise in the source term decomposition, in combination with the Roe approximate Riemann solver that includes the spatial derivative of the flux function, lead to the finite volume WENO scheme that satisfies the exact conservation property – the property of preserving the quiescent flow exactly. When the central WENO schemes are applied, additional reformulations are introduced for the transition from the staggered values to the nonstaggered ones and vice versa by using the WENO reconstruction procedure. The proposed central WENO schemes also preserve the quiescent flow, but only in prismatic channels. In various test problems the obtained balanced schemes show improvements in comparison with the standard versions of the proposed type schemes, as well as with some other first- and second-order numerical schemes.     

A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: On a rectangular collocated grid system

A consistent, conservative and accurate scheme has been designed to calculate the current density and the Lorentz force by solving the electrical potential equation for magnetohydrodynamics (MHD) at low magnetic Reynolds numbers and high Hartmann numbers on a finite-volume structured collocated grid. In this collocated grid, velocity (u), pressure (p), and electrical potential (φ) are located in the grid center, while current fluxes are located on the cell faces. The calculation of current fluxes on the cell faces is conducted using a conservative scheme, which is consistent with the discretization scheme for the solution of electrical potential Poisson equation. A conservative interpolation is used to get the current density at the cell center, which is used to conduct the calculation of Lorentz force at the cell center for momentum equations. We will show that both “conservative” and “consistent” are important properties of the scheme to get an accurate result for high Hartmann number MHD flows with a strongly non-uniform mesh employed to resolve the Hartmann layers and side layers of Hunt’s conductive walls and Shercliff’s insulated walls. A general second-order projection method has been developed for the incompressible Navier–Stokes equations with the Lorentz force included. This projection method can accurately balance the pressure term and the Lorentz force for a fully developed core flow. This method can also simplify the pressure boundary conditions for MHD flows.     

The Zero-Mach Limit of Compressible Flows

We prove that compressible Navier-Stokes flows in two and three space dimensions converge to incompressible Navier-Stokes flows in the limit as the Mach number tends to zero. No smallness restrictions are imposed on the external force, the initial velocity, or the time interval. We assume instead that the incompressible flow exists and is reasonably smooth on a given time interval, and prove that compressible flows with compatible initial data converge uniformly on that time interval. Our analysis shows that the essential mechanism in this process is a hyperbolic effect which becomes stronger with smaller Mach number and which ultimately drives the density to a constant. Received: 10 June 1997 / Accepted: 15 July 1997     

Two new upwind difference schemes for a coupled system of convection–diffusion equations arising from the steady MHD duct flow problems

In this paper, we develop two new upwind difference schemes for solving a coupled system of convection–diffusion equations arising from the steady incompressible MHD duct flow problem with a transverse magnetic field at high Hartmann numbers. Such an MHD duct flow is convection-dominated and its solution may exhibit localized phenomena such as boundary layers, namely, narrow boundary regions where the solution changes rapidly. Most conventional numerical schemes cannot efficiently solve the layer problems because they are lacking in either stability or accuracy. In contrast, the newly proposed upwind difference schemes can achieve a reasonable accuracy with a high stability, and they are capable of resolving high gradients near the layer regions without refining the grid. The accuracy of the first new upwind scheme is O(h + k) and the second one improves the accuracy to O(ε²(h + k) + ε(h² + k²) + (h³ + k³)), where 0 < ε ? 1/M ? 1 and M is the high Hartmann number. Numerical examples are provided to illustrate the performance of the newly proposed upwind difference schemes.     

ON EDGE WAVES IN STRATIFIED WATER ALONG A SLOPING BEACH

Building on previous investigations, we show that Gerstner's famous deep water wave and the related edge wave propagating along a sloping beach, found within the context of water of constant density, can both be adapted to provide explicit free surface flows in incompressible fluids with arbitrary density stratification.     

Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model

In this paper we analytically solve the velocity of the lattice Boltzmann BGK equation (LBGK) for several simple flows. The analysis provides a framework to theoretically analyze various boundary conditions. In particular, the analysis is used to derive the slip velocities generated by various schemes for the nonslip boundary condition. We find that the slip velocity is zero as long as fe=0 at boundaries, no matter what combination of distributions is chosen. The schemes proposed by Nobleet al. and by Inamuroet al. yield the correct zeroslip velocity, while some other schemes, such as the bounce-back scheme and the equilibrium distribution scheme, would inevitably generate a nonzero slip velocity. The bounce-back scheme with the wall located halfway between a flow node and a bounce-back node is also studied for the simple flows considered and is shown to produce results of second-order accuracy. The momentum exchange at boundaries seems to be highly related to the slip velocity at boundaries. To be specific, the slip velocity is zero only when the momentum dissipated by boundaries is equal to the stress provided by fluids.     

High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy

In this paper, a finite difference code for Direct and Large Eddy Simulation (DNS/LES) of incompressible flows is presented. This code is an intermediate tool between fully spectral Navier–Stokes solvers (limited to academic geometry through Fourier or Chebyshev representation) and more versatile codes based on standard numerical schemes (typically only second-order accurate). The interest of high-order schemes is discussed in terms of implementation easiness, computational efficiency and accuracy improvement considered through simplified benchmark problems and practical calculations. The equivalence rules between operations in physical and spectral spaces are efficiently used to solve the Poisson equation introduced by the projection method. It is shown that for the pressure treatment, an accurate Fourier representation can be used for more flexible boundary conditions than periodicity or free-slip. Using the concept of the modified wave number, the incompressibility can be enforced up to the machine accuracy. The benefit offered by this alternative method is found to be very satisfactory, even when a formal second-order error is introduced locally by boundary conditions that are neither periodic nor symmetric. The usefulness of high-order schemes combined with an immersed boundary method (IBM) is also demonstrated despite the second-order accuracy introduced by this wall modelling strategy. In particular, the interest of a partially staggered mesh is exhibited in this specific context. Three-dimensional calculations of transitional and turbulent channel flows emphasize the ability of present high-order schemes to reduce the computational cost for a given accuracy. The main conclusion of this paper is that finite difference schemes with quasi-spectral accuracy can be very efficient for DNS/LES of incompressible flows, while allowing flexibility for the boundary conditions and easiness in the code development. Therefore, this compromise fits particularly well for very high-resolution simulations of turbulent flows with relatively complex geometries without requiring heavy numerical developments.     

A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios

We present an efficient time-stepping scheme for simulations of the coupled Navier–Stokes Cahn–Hilliard equations for the phase field approach. The scheme has several attractive characteristics: (i) it is suitable for large density ratios, and numerical experiments with density ratios up to 1000 have been presented; (ii) it involves only constant (time-independent) coefficient matrices for all flow variables, which can be pre-computed during pre-processing, so it effectively overcomes the performance bottleneck induced by variable coefficient matrices associated with the variable density and variable viscosity; (iii) it completely de-couples the computations of the velocity, pressure, and the phase field function. Strategy for spectral-element type spatial discretizations to overcome the difficulty associated with the large spatial order of the Cahn–Hilliard equation is also discussed. Ample numerical simulations demonstrate that the current algorithm, together with the Navier–Stokes Cahn–Hilliard phase field approach, is an efficient and effective method for studying two-phase flows involving large density ratios, moving contact lines, and interfacial topology changes.     

Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows

In this work we present a general strategy for constructing multidimensional HLLE Riemann solvers, with particular attention paid to detailing the two-dimensional HLLE Riemann solver. This is accomplished by introducing a constant resolved state between the states being considered, which introduces sufficient dissipation for systems of conservation laws. Closed form expressions for the resolved fluxes are also provided to facilitate numerical implementation. The Riemann solver is proved to be positively conservative for the density variable; the positivity of the pressure variable has been demonstrated for Euler flows when the divergence in the fluid velocities is suitably restricted so as to prevent the formation of cavitation in the flow.We also focus on the construction of multidimensionally upwinded electric fields for divergence-free magnetohydrodynamical (MHD) flows. A robust and efficient second order accurate numerical scheme for two and three-dimensional Euler and MHD flows is presented. The scheme is built on the current multidimensional Riemann solver and has been implemented in the author’s RIEMANN code. The number of zones updated per second by this scheme on a modern processor is shown to be cost-competitive with schemes that are based on a one-dimensional Riemann solver. However, the present scheme permits larger timesteps.Accuracy analysis for multidimensional Euler and MHD problems shows that the scheme meets its design accuracy. Several stringent test problems involving Euler and MHD flows are also presented and the scheme is shown to perform robustly on all of them.     

High-order finite-volume methods for the shallow-water equations on the sphere

This paper presents a third-order and fourth-order finite-volume method for solving the shallow-water equations on a non-orthogonal equiangular cubed-sphere grid. Such a grid is built upon an inflated cube placed inside a sphere and provides an almost uniform grid point distribution. The numerical schemes are based on a high-order variant of the Monotone Upstream-centered Schemes for Conservation Laws (MUSCL) pioneered by van Leer. In each cell the reconstructed left and right states are either obtained via a dimension-split piecewise-parabolic method or a piecewise-cubic reconstruction. The reconstructed states then serve as input to an approximate Riemann solver that determines the numerical fluxes at two Gaussian quadrature points along the cell boundary. The use of multiple quadrature points renders the resulting flux high-order. Three types of approximate Riemann solvers are compared, including the widely used solver of Rusanov, the solver of Roe and the new AUSM⁺-up solver of Liou that has been designed for low-Mach number flows. Spatial discretizations are paired with either a third-order or fourth-order total-variation-diminishing Runge–Kutta timestepping scheme to match the order of the spatial discretization. The numerical schemes are evaluated with several standard shallow-water test cases that emphasize accuracy and conservation properties. These tests show that the AUSM⁺-up flux provides the best overall accuracy, followed closely by the Roe solver. The Rusanov flux, with its simplicity, provides significantly larger errors by comparison. A brief discussion on extending the method to arbitrary order-of-accuracy is included.     

High-order well-balanced schemes and applications to non-equilibrium flow

The appearance of the source terms in modeling non-equilibrium flow problems containing finite-rate chemistry or combustion poses additional numerical difficulties beyond that for solving non-reacting flows. A well-balanced scheme, which can preserve certain non-trivial steady state solutions exactly, may help minimize some of these difficulties. In this paper, a simple one-dimensional non-equilibrium model with one temperature is considered. We first describe a general strategy to design high-order well-balanced finite-difference schemes and then study the well-balanced properties of the high-order finite-difference weighted essentially non-oscillatory (WENO) scheme, modified balanced WENO schemes and various total variation diminishing (TVD) schemes. The advantages of using a well-balanced scheme in preserving steady states and in resolving small perturbations of such states will be shown. Numerical examples containing both smooth and discontinuous solutions are included to verify the improved accuracy, in addition to the well-balanced behavior.     

Rayleigh-Taylor instability in variable density swirling flows

Using linear instability theory and nonlinear dynamics, the Rayleigh-Taylor instability of variable density swirling flows is studied. It is found that the flow topology could be predicted, when the instability sets in, using a function χ dependent on density and axial and azimuthal velocities. It is shown that even when the inner axial-flow is heavier than the outer one (a favorable case for the development of the Rayleigh-Taylor instability thanks to the centrifugal force) the instability is not necessarily Rayleigh-Taylor-dominated. It is also shown that when the Rayleigh-Taylor instability develops, it is helical.     

Efficient kinetic schemes for steady and unsteady flow simulations on unstructured meshes

This paper presents efficient second-order kinetic schemes on unstructured meshes for both compressible unsteady and incompressible steady flows. For compressible unsteady flows, a time-dependent gas distribution function with a discontinuous particle velocity space at a cell interface is constructed and used for the evaluations of both numerical fluxes and conservative flow variables. As a result, a compact scheme on the unstructured meshes is developed. For incompressible steady flows, a continuous second-order gas-kinetic BGK type scheme is presented, for which the time-dependent gas distribution function with a continuous particle velocity is used on unstructured meshes. The efficiency of the schemes lies in the fact that the slopes of the flow variables inside each cell can be constructed using values of the flow variables within that cell only without involving neighboring cells. Therefore, even with the stencil of a first-order scheme, a high resolution method is constructed. Numerical examples are presented which are compared with the benchmark solutions and the experimental measurements.