A Residual-Based Compact Scheme for the Compressible Navier–Stokes Equations

A simple and efficient time-dependent method is presented for solving the steady compressible Euler and Navier–Stokes equations with third-order accuracy. Owing to its residual-based structure, the numerical scheme is compact without requiring any linear algebra, and it uses a simple numerical dissipation built on the residual. The method contains no tuning parameter. Accuracy and efficiency are demonstrated for 2-D inviscid and viscous model problems. Navier–Stokes calculations are presented for a shock/boundary layer interaction, a separated laminar flow, and a transonic turbulent flow over an airfoil.     

Wave Simulation in Frozen Porous Media

We propose a numerical algorithm for simulation of wave propagation in frozen porous media, where the pore space is filled with ice and water. The model, based on a Biot-type three-phase theory, predicts three compressional waves and two shear waves and models the attenuation level observed in rocks. Attenuation is modeled with exponential relaxation functions which allow a differential formulation based on memory variables. The wavefield is obtained using a grid method based on the Fourier differential operator and a Runge–Kutta time-integration algorithm. Since the presence of slow quasistatic modes makes the differential equations stiff, a time-splitting integration algorithm is used to solve the stiff part analytically. The modeling is second-order accurate in the time discretization and has spectral accuracy in the calculation of the spatial derivatives.     

A High-Order WENO Finite Difference Scheme for the Equations of Ideal Magnetohydrodynamics

We present a high-order accurate weighted essentially non-oscillatory (WENO) finite difference scheme for solving the equations of ideal magnetohydrodynamics (MHD). This scheme is a direct extension of a WENO scheme, which has been successfully applied to hydrodynamic problems. The WENO scheme follows the same idea of an essentially non-oscillatory (ENO) scheme with an advantage of achieving higher-order accuracy with fewer computations. Both ENO and WENO can be easily applied to two and three spatial dimensions by evaluating the fluxes dimension-by-dimension. Details of the WENO scheme as well as the construction of a suitable eigen-system, which can properly decompose various families of MHD waves and handle the degenerate situations, are presented. Numerical results are shown to perform well for the one-dimensional Brio–Wu Riemann problems, the two-dimensional Kelvin–Helmholtz instability problems, and the two-dimensional Orszag–Tang MHD vortex system. They also demonstrate the importance of maintaining the divergence free condition for the magnetic field in achieving numerical stability. The tests also show the advantages of using the higher-order scheme. The new 5th-order WENO MHD code can attain an accuracy comparable with that of the second-order schemes with many fewer grid points.     

Divergence- and Curl-Preserving Prolongation and Restriction Formulas

We present new second-order prolongation and restriction formulas which preserve the divergence and, in some cases, the curl of a discretized vector field. The formulas are suitable for adaptive and hierarchical mesh algorithms with a factor-of-2 linear resolution change. We examine both staggered and collocated discretizations for the vector field on two- and three-dimensional Cartesian grids. The new formulas can be used in combination with numerical schemes that require a divergence-free solution in some discrete sense, such as the constrained transport schemes of computational magnetohydrodynamics. We also obtain divergence-preserving interpolation functions which may be used for streamline or field line tracing.     

On Performance of Methods with Third- and Fifth-Order Compact Upwind Differencing

The difference schemes for fluid dynamics type of equations based on third- and fifth-order Compact Upwind Differencing (CUD) are considered. To validate their properties following from a linear analysis, calculations were carried out using the inviscid and viscous Burgers' equation as well as the compressible Navier–Stokes equation written in the conservative form for curvilinear coordinates. In the latter case, transonic cascade flow was chosen as a representative example. The performance of the CUD methods was estimated by investigating mesh convergence of the solutions and comparing with the results of second-order schemes. It is demonstrated that the oscillation-free steep gradients solutions obtained without using smoothing techniques can provide considerable increase of accuracy even when exploiting coarse meshes.     

Monotonicity Preserving Weighted Essentially Non-oscillatory Schemes with Increasingly High Order of Accuracy

In this paper we design a class of numerical schemes that are higher-order extensions of the weighted essentially non-oscillatory (WENO) schemes of G.-S. Jiang and C.-W. Shu (1996) and X.-D. Liu, S. Osher, and T. Chan (1994). Used by themselves, the schemes may not always be monotonicity preserving but coupled with the monotonicity preserving bounds of A. Suresh and H. T. Huynh (1997) they perform very well. The resulting monotonicity preserving weighted essentially non-oscillatory (MPWENO) schemes have high phase accuracy and high order of accuracy. The higher-order members of this family are almost spectrally accurate for smooth problems. Nevertheless, they, have robust shock capturing ability. The schemes are stable under normal CFL numbers. They are also efficient and do not have a computational complexity that is substantially greater than that of the lower-order members of this same family of schemes. The higher accuracy that these schemes offer coupled with their relatively low computational complexity makes them viable competitors to lower-order schemes, such as the older total variation diminishing schemes, for problems containing both discontinuities and rich smooth region structure. We describe the MPWENO schemes here as well as show their ability to reach their designed accuracies for smooth flow. We also examine the role of steepening algorithms such as the artificial compression method in the design of very high order schemes. Several test problems in one and two dimensions are presented. For multidimensional problems where the flow is not aligned with any of the grid directions it is shown that the present schemes have a substantial advantage over lower-order schemes. It is argued that the methods designed here have great utility for direct numerical simulations and large eddy simulations of compressible turbulence. The methodology developed here is applicable to other hyperbolic systems, which is demonstrated by showing that the MPWENO schemes also work very well on magnetohydrodynamical test problems.     

Minimization of the Truncation Error by Grid Adaptation

A new grid adaptation strategy, which minimizes the truncation error of a pth-order finite difference approximation, is proposed. The main idea of the method is based on the observation that the global truncation error associated with discretization on nonuniform meshes can be minimized if the interior grid points are redistributed in an optimal sequence. The method does not explicitly require the truncation error estimate, and at the same time, it allows one to increase the design order of approximation globally by one, so that the same finite difference operator reveals superconvergence properties on the optimal grid. Another very important characteristic of the method is that if the differential operator and the metric coefficients are evaluated identically by some hybrid approximation, then the single optimal grid generator can be employed in the entire computational domain independently of points where the hybrid discretization switches from one approximation to another. Generalization of the present method to multiple dimensions is presented. Numerical calculations of several one-dimensional and one two-dimensional test examples demonstrate the performance of the method and corroborate the theoretical results.     

A Higher-Order Boundary Treatment for Cartesian-Grid Methods

The Euler equations describe the flow phenomena of compressible inviscid gas dynamics. We simulate such flows using a higher-order Cartesian-grid method, together with a special treatment for the cells cut by the boundary of an object. A new method for the treatment of the boundary is described where these cut boundary cells are maintained as whole cells rather than as cut cells, thus avoiding stability problems. The method is second-order accurate in one dimension and higher-order accurate in two dimensions but not strictly conservative; however, we show that this error in the conservation does not lead to spurious phenomena on some representative test calculations. The advantages of the new boundary treatment are that it is higher-order accurate, that it is independent of the applied method, and that it is simple.     

An Efficient Algorithm for Truncating Spatial Domain in Modeling Light Scattering by Finite-Difference Technique

The finite-difference time domain technique is one of the most robust and accurate numerical methods for the solution of light scattering by small particles with arbitrary composition and geometry. In practice, this method requires that the spatial domain for the computation of near-field be truncated. An absorbing boundary condition must be imposed in conjunction with this truncation. The performance of this boundary condition is essential to the stability of numerical computations and the reliability of results. In the present study, a new boundary condition, referred to as the mixed T algorithm, has been developed, which is a generalization of the transmitting boundary condition originally developed by Liao and co-workers. The present algorithm does not require spatial interpolation for wave values at interior grid points. In addition, it produces two minima of spurious reflections at small and large incident angles, allowing efficient absorption of the scattered waves at the boundary for large incident angles. When the third-order mixed T algorithm is used, the reflection coefficient of the boundary is less than 1% for incident angles from 0° to about 70°. We find that the numerical instability associated with the transmitting boundary condition is caused by the location-dependent amplitude of outgoing waves in the vicinity of the boundary. For this reason, the mixed T algorithm is stabilized by consistently introducing diffusive coefficients into the boundary equation. When the stabilized algorithm is applied, the near-field within the truncated domain can be computed by using single-precision arithmetic without overflows for more than 10⁵steps in the time-marching iteration. Finally, the new absorbing boundary condition is validated by carrying out numerical experiments involving the propagation of a TM wave excited by a sinusoidal point source, simultaneous simulation of the wave propagation in small and large domains, and the scattering of a TM wave by an infinite circular cylinder.     

A low numerical dissipation immersed interface method for the compressible Navier–Stokes equations

A numerical method to solve the compressible Navier–Stokes equations around objects of arbitrary shape using Cartesian grids is described. The approach considered here uses an embedded geometry representation of the objects and approximate the governing equations with a low numerical dissipation centered finite-difference discretization. The method is suitable for compressible flows without shocks and can be classified as an immersed interface method. The objects are sharply captured by the Cartesian mesh by appropriately adapting the discretization stencils around the irregular grid nodes, located around the boundary. In contrast with available methods, no jump conditions are used or explicitly derived from the boundary conditions, although a number of elements are adopted from previous immersed interface approaches. A new element in the present approach is the use of the summation-by-parts formalism to develop stable non-stiff first-order derivative approximations at the irregular grid points. Second-order derivative approximations, as those appearing in the transport terms, can be stiff when irregular grid points are located too close to the boundary. This is addressed using a semi-implicit time integration method. Moreover, it is shown that the resulting implicit equations can be solved explicitly in the case of constant transport properties. Convergence studies are performed for a rotating cylinder and vortex shedding behind objects of varying shapes at different Mach and Reynolds numbers.     

Phase-Accuracy Comparisons and Improved Far-Field Estimates for 3-D Edge Elements on Tetrahedral Meshes

Edge-element methods have proved very effective for 3-D electromagnetic computations and are widely used on unstructured meshes. However, the accuracy of standard edge elements can be criticised because of their low order. This paper analyses discrete dispersion relations together with numerical propagation accuracy to determine the effect of tetrahedral shape on the phase accuracy of standard 3-D edge-element approximations in comparison to other methods. Scattering computations for the sphere obtained with edge elements are compared with results obtained with vertex elements, and a new formulation of the far-field integral approximations for use with edge elements is shown to give improved cross sections over conventional formulations.     

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.     

Accuracy preserving limiter for the high-order accurate solution of the Euler equations

Higher-order finite-volume methods have been shown to be more efficient than second-order methods. However, no consensus has been reached on how to eliminate the oscillations caused by solution discontinuities. Essentially non-oscillatory (ENO) schemes provide a solution but are computationally expensive to implement and may not converge well for steady-state problems. This work studies the extension of limiters used for second-order methods to the higher-order case. Requirements for accuracy and efficient convergence are discussed. A new limiting procedure is proposed. Ringleb’s flow problem is used to demonstrate that nearly nominal orders of accuracy for schemes up to fourth-order can be achieved in smooth regions using the new limiter. Results for the fourth-order accurate solution of transonic flow demonstrates good convergence properties and significant qualitative improvement of the solution relative the second-order method. The new limiter can also be successfully applied to reduce the dissipation of second-order schemes with minimal sacrifices in convergence properties relative to existing approaches.     

MPDATA: A Finite-Difference Solver for Geophysical Flows

This article is a review of MPDATA, a class of methods for the numerical simulation of fluid flows based on the sign-preserving properties of upstream differencing. MPDATA was designed originally as an inexpensive alternative to flux-limited schemes for evaluating the advection of nonnegative thermodynamic variables (such as liquid water or water vapour) in atmospheric models. During the last decade, MPDATA has evolved from a simple advection scheme to a general approach for integrating the conservation laws of geophysical fluids on micro-to-planetary scales. The purpose of this paper is to summarize the basic concepts leading to a family of MPDATA schemes, to review existing MPDATA options, and to demonstrate the use of MPDATA to effectively construct two distinct types of models (elastic and anelastic) for complex geophysical flows.     

A Finite Volume Method for the Approximation of Diffusion Operators on Distorted Meshes

A new finite volume method is presented for discretizing general linear or nonlinear elliptic second-order partial-differential equations with mixed boundary conditions. The advantage of this method is that arbitrary distorted meshes can be used without the numerical results being altered. The resulting algorithm has more unknowns than standard methods like finite difference or finite element methods. However, the matrices that need to be inverted are positive definite, so the most powerful linear solvers can be applied. The method has been tested on a few elliptic and parabolic equations, either linear, as in the case of the standard heat diffusion equation, or nonlinear, as in the case of the radiation diffusion equation and the resistive diffusion equation with Hall term.     

New High-Resolution Semi-discrete Central Schemes for Hamilton–Jacobi Equations

We introduce a new high-resolution central scheme for multidimensional Hamilton–Jacobi equations. The scheme retains the simplicity of the non-oscillatory central schemes developed by C.-T. Lin and E. Tadmor (in press, SIAM J. Sci. Comput.), yet it enjoys a smaller amount of numerical viscosity, independent of 1/Δt. By letting Δt↓0 we obtain a new second-order central scheme in the particularly simple semi-discrete form, along the lines of the new semi-discrete central schemes recently introduced by the authors in the context of hyperbolic conservation laws. Fully discrete versions are obtained with appropriate Runge–Kutta solvers. The smaller amount of dissipation enables efficient integration of convection-diffusion equations, where the accumulated error is independent of a small time step dictated by the CFL limitation. The scheme is non-oscillatory thanks to the use of nonlinear limiters. Here we advocate the use of such limiters on second discrete derivatives, which is shown to yield an improved high resolution when compared to the usual limitation of first derivatives. Numerical experiments demonstrate the remarkable resolution obtained by the proposed new central scheme.     

A Volume of Fluid Based Method for Fluid Flows with Phase Change

This paper presents a numerical method directed towards the simulation of flows with mass transfer due to changes of phase. We use a volume of fluid (VOF) based interface tracking method in conjunction with a mass transfer model and a model for surface tension. The bulk fluids are viscous, conducting, and incompressible. A one-dimensional test problem is developed with the feature that a thin thermal layer propagates with the moving phase interface. This test problem isolates the ability of a method to accurately calculate the thermal layers responsible for driving the mass transfer in boiling flows. The numerical method is tested on this problem and then is used in simulations of horizontal film boiling.     

Solution of Time-Dependent Diffusion Equations with Variable Coefficients Using Multiwavelets

A new numerical algorithm is developed for the solution of time-dependent differential equations of diffusion type. It allows for an accurate and efficient treatment of multidimensional problems with variable coefficients, nonlinearities, and general boundary conditions. For space discretization we use the multiwavelet bases introduced by Alpert (1993,SIAM J. Math. Anal.24, 246–262), and then applied to the representation of differential operators and functions of operators presented by Alpert, Beylkin, and Vozovoi (Representation of operators in the multiwavelet basis, in preparation). An important advantage of multiwavelet basis functions is the fact that they are supported only on non-overlapping subdomains. Thus multiwavelet bases are attractive for solving problems in finite (non periodic) domains. Boundary conditions are imposed with a penalty technique of Hesthaven and Gottlieb (1996,SIAM J. Sci. Comput., 579–612) which can be used to impose rather general boundary conditions. The penalty approach was extended to a procedure for ensuring the continuity of the solution and its first derivative across interior boundaries between neighboring subdomains while time stepping the solution of a time dependent problem. This penalty procedure on the interfaces allows for a simplification and sparsification of the representation of differential operators by discarding the elements responsible for interactions between neighboring subdomains. Consequently the matrices representing the differential operators (on the finest scale) have block-diagonal structure. For a fixed order of multiwavelets (i.e., a fixed number of vanishing moments) the computational complexity of the present algorithm is proportional to the number of subdomains. The time discretization method of Beylkin, Keiser, and Vozovoi (1998, PAM Report 347) is used in view of its favorable stability properties. Numerical results are presented for evolution equations with variable coefficients in one and two dimensions.     

Multigrid Solution of the Incompressible Navier–Stokes Equations for Density-Stratified Flow past Three-Dimensional Obstacles

The steady incompressible Navier–Stokes equations in three dimensions are solved for neutral and stably stratified flow past three-dimensional obstacles of increasing spanwise width. The continuous equations are approximated using a finite volume discretisation on staggered grids with a flux-limited monotonic scheme for the advective terms. The discrete equations which arise are solved using a nonlinear multigrid algorithm with up to four grid levels using the SIMPLE pressure correction method as smoother. When at its most effective the multigrid algorithm is demonstrated to yield convergence rates which are independent of the grid density. However, it is found that the asymptotic convergence rate depends on the choice of the limiter used for the advective terms of the density equation, and some commonly used schemes are investigated. The variation with obstacle width of the influence of the stratification on the flow field is described and the results of the three-dimensional computations are compared with those of the corresponding computation of flow over a two-dimensional obstacle (of effectively infinite width). Also given are the results of time-dependent computations for three-dimensional flows under conditions of strong static stability when lee-wave propagation is present and the multigrid algorithm is used to compute the flow at each time step.     

A High-Order Discontinuous Galerkin Method for 2D Incompressible Flows

In this paper we introduce a high-order discontinuous Galerkin method for two-dimensional incompressible flow in the vorticity stream-function formulation. The momentum equation is treated explicitly, utilizing the efficiency of the discontinuous Galerkin method. The stream function is obtained by a standard Poisson solver using continuous finite elements. There is a natural matching between these two finite element spaces, since the normal component of the velocity field is continuous across element boundaries. This allows for a correct upwinding gluing in the discontinuous Galerkin framework, while still maintaining total energy conservation with no numerical dissipation and total enstrophy stability. The method is efficient for inviscid or high Reynolds number flows. Optimal error estimates are proved and verified by numerical experiments.