Analysis of central and upwind compact schemes

https://doi.org/10.1016/j.jcp.2003.07.015Get rights and content

Abstract

Central and upwind compact schemes for spatial discretization have been analyzed with respect to accuracy in spectral space, numerical stability and dispersion relation preservation. A von Neumann matrix spectral analysis is developed here to analyze spatial discretization schemes for any explicit and implicit schemes to investigate the full domain simultaneously. This allows one to evaluate various boundary closures and their effects on the domain interior. The same method can be used for stability analysis performed for the semi-discrete initial boundary value problems (IBVP). This analysis tells one about the stability for every resolved length scale. Some well-known compact schemes that were found to be G-K-S and time stable are shown here to be unstable for selective length scales by this analysis. This is attributed to boundary closure and we suggest special boundary treatment to remove this shortcoming. To demonstrate the asymptotic stability of the resultant schemes, numerical solution of the wave equation is compared with analytical solution. Furthermore, some of these schemes are used to solve two-dimensional Navier–Stokes equation and a computational acoustic problem to check their ability to solve problems for long time. It is found that those schemes, that were found unstable for the wave equation, are unsuitable for solving incompressible Navier–Stokes equation. In contrast, the proposed compact schemes with improved boundary closure and an explicit higher-order upwind scheme produced correct results. The numerical solution for the acoustic problem is compared with the exact solution and the quality of the match shows that the used compact scheme has the requisite DRP property.

Introduction

With the advent of powerful computers it is now common to solve differential equations numerically in various disciplines of science and engineering. This includes the solution of partial differential equations with stringent requirements of resolving wide range of spatial and temporal scales. For example, in CFD it is now common to solve the governing Navier–Stokes equation resolving all the scales of a turbulent flow in DNS for moderate Reynolds numbers. Similarly in many wave propagation problems one solves governing hyperbolic partial differential equations and such solutions are required to be accurate in the far field and for long time periods. These requirements demand that the adopted numerical method be highly accurate and dispersion error free. Lighthill [1] and Taflove [2] have discussed, respectively, the problems of computational aero-acoustics (CAA) and computational electromagnetics (CEM) for numerical solutions with respect to these issues.

Such requirements are best satisfied by spectral methods [3], [4]. Spectral methods have been used mostly for problems involving simple geometries and boundary conditions, though applications to complex domains are available [26], [27]. Alternatives to spectral method are the (i) higher-order explicit upwind methods as in [5], [6] and (ii) methods based on Padé approximation as was originally given in [7]. The latter methods offer higher-order approximations to differential operators using compact stencils that relate various derivatives with function values at discrete nodes. All higher-order compact schemes can be expressed by the following set of linear algebraic equation[A]u=[B]u.Each row of this represents an implicit relation between the derivatives and function values for computational nodes. For numerical convenience, it is desired that A and B be band-limited. The most often used structure of A matrix is tridiagonal. This equation can also be written down in an equivalent explicit form byu=[A]−1[B]u=[C]u,where C is not necessarily compact. The early work on compact differencing scheme can be seen in [8], [9], [10]. Such methods have been optimized to solve the time-domain Maxwell equations in [11] and for problems relating to acoustics in [12], [13], [28]. In [12] the important concept of dispersion relation preservation (DRP) is discussed with respect to high accuracy schemes for acoustics problems.

For problems with periodic boundary conditions, A and B are periodic symmetric matrices. However, for many practical problems periodic boundary conditions are not applicable and one-sided stencils are needed near boundaries, making A and B matrices non-symmetric. Computationally, symmetric B matrix corresponds to non-dissipative central schemes and non-symmetric B matrix arises from upwind schemes. Schemes given in [9], [10], [14] are typical examples of central non-dissipative method of spatial discretization. In [14] a periodic problem is solved and hence the symmetric stencil is used at all points. However, in [9] a non-periodic problem is solved that requires taking asymmetric stencils at and near-boundary points while the inner stencil is still symmetric. Such forced upwinding near boundaries can cause the overall method to become unstable. One of the aim here is to use a Fourier spectral frame-work to analyze schemes including inner and boundary closure schemes simultaneously. For upwind compact schemes, B matrix is asymmetric for all the nodes. Typical examples are as in [15], [16].

Issue of whether to use a central or upwind compact scheme rests on numerical stability of such schemes everywhere in the domain. Numerical stability is usually investigated in two ways. The first is the method given in [17] that is based on normal modal analysis and referred to as G-K-S stability theory. In this theory estimates are developed for inner and boundary schemes to ensure stability. First-order partial differential equations were considered [17] as IBVPs for which the term involving the spatial derivative is multiplied by a square Hermitian matrix. For compact high-order schemes, any problem with non-periodic boundary would require an asymmetric stencil and hence the corresponding matrix is not Hermitian. For example, Carpenter et al. [18] considered two spatial discretizations that lead to the above matrix being asymmetric. The subsequent semi-discrete eigenvalue analysis (by considering a time continuous system) for the explicit and compact fourth-order spatial operators by the authors revealed a spectrum with some eigenvalues crossing over the imaginary axis into the right half plane (Figs. 3 and 4 of [18]) indicating instability. Investigated spatial discretization schemes were, however, stable in the Kreiss sense.

This brings one to the important aspect of asymptotic stability analysis. The strong point of the G-K-S analysis is that this presents a theory that includes interior points and boundaries together. But, the definition of G-K-S stability (also known as Lax stability) might be too weak [18] and hence it will not be a practical option to use those compact schemes for DNS, CAA and CEM that are only G-K-S stable. According to [18], for truly time dependent problems, sufficiently accurate solutions over long time requires excessively large number of grid points. Basically the essential dynamics of the system is then expressed over a small fraction of grid-resolved wave numbers of the very fine grid. However, the fine grid requirement negates the basic advantages of compact schemes with fewer grid points.

The time-stability analysis given in [15], [18] are in the physical plane and do not provide information about the performance of schemes at different length scales. Here we adopt a different approach of investigation for spatial discretization schemes in spectral plane following the method of [19], allowing us to look at interior and boundary points simultaneously. Spatial discretization schemes for interior and boundary points may be found stable when analyzed in isolation. But, overall numerical stability may be lost when these are analyzed together. In this context a few well-known schemes are analyzed here for their accuracy and stability.

Working in the physical plane with a uniform grid of size h, the unknown is related to its bi-lateral Laplace transform byu(xl)=∫U(k)eikxldk,where the integral is performed over the limit −km to km, defined by the Nyquist limit of km=π/h. Hence a first derivative is obtained analytically by u(xl)=∫ikU(k)eikxldk.

If the first derivative is evaluated numerically by discrete method, then the same can be written asu(xl)=∫ikeqU(k)eikxldk.Different numerical schemes have different estimates of keq and it is in general a complex quantity. For numerical stability of any scheme, one must look at the imaginary part of keq. The imaginary part of keq represents numerical dissipation only when it is negative. Any scheme, that produces a positive imaginary part of keq is numerically unstable because a positive imaginary part is equivalent to adding anti-diffusion. This method of analysis for spatial discretization is independent of the actual differential equation that has to be solved and it also indicates the length scales at which instability arises.

In the next section accuracy in spectral space is discussed considering some well-known compact schemes and some schemes that we propose here for spatial discretization. In Section 3, the numerical stability for the solution of linear wave equation is discussed for the spatial schemes used in conjunction with Euler time integration scheme. The DRP property of these space–time discretization schemes are analyzed with respect to linear wave equation in Section 4. In Section 5, numerical solution of linear wave equation is compared with exact solution for discontinuous initial data. We demonstrate the ability of two compact schemes and an explicit higher-order upwind scheme in solving the Navier–Stokes equation for the receptivity of a shear layer to a convecting vortex in the free stream. Finally, we solve an acoustics problem and compare with exact solution to show the effectiveness of the computing schemes.

Section snippets

Accuracy of compact schemes in spectral space

To evaluate keq for compact schemes for first derivative in Eq. (4), the implicit relation between the first derivatives and the function values at different nodes can be expressed by the equivalent explicit relation as given in Eq. (2). To perform spectral analysis one uses Eq. (3) in Eq. (2) – giving rise to an implicit relation between the derivative at a given point with function values of the whole domain. One can relate the derivative at the lth-node in terms of the function value there

The numerical stability and amplification properties

To evaluate any numerical scheme it is natural to calibrate it with respect to a model equation that mimics the physical processes and preferably possesses exact solution. The linear wave equation serves these requirements. For DNS, physical convection is most important and a method is preferred that solves wave equation most effectively. With this in mind we use different compact schemes to discretized spatial derivative while using standard temporal schemes for the wave equationut+cux=0,

Dispersion relation preservation property of numerical schemes

For many wave propagation problems the group velocity (Vg) is a more meaningful physical quantity than the phase speed – as energy propagates with this velocity. In this context linear wave equation is non-dispersive i.e. the group velocity and phase speed are identical. Thus if we use Eq. (13) as the model equation, then a numerical scheme for space–time dependent problem can be characterized by comparing the numerical group velocity with the phase speed. Spectral method shows equality of

Solution of linear wave equation

The analyses of the previous sections have shown that the Euler time integration scheme has near-neutral amplification factor and good DRP property for low computational efforts. In solving linear wave equation Euler time integration is considered here to highlight certain features of spatial schemes by comparing the computed results with exact solution.

Here a domain of 0⩽x⩽1 is considered with 200 uniformly distributed points. The phase speed is taken as c=0.05 such that the signal given by

Receptivity of shear layer

Vortex-induced eruptions inside a shear layer, as excited by a convecting vortex far outside the shear layer is an important problem with many applications. This has been investigated earlier theoretically [23], [24] and experimentally [25]. This is a typical receptivity problem and can be used to check the efficiency of some of the compact schemes to solve Navier–Stokes equation. Here the Navier–Stokes equation in stream function-vorticity (ψω) formulation has been solved using the K-scheme

Conclusions

Central and upwind compact schemes have been analyzed by a matrix- spectral analysis developed here. It is shown that some of these schemes are unstable due to boundary stencils near one of the boundaries. To cure this problem – arising from the boundary treatment – three new optimal upwind-biased schemes have been proposed (OUCS1, OUCS2 and OUCS3 defined in Eqs. , , ).

The combined space–time discretization methods for linear wave equation show the amplification factor to be unstable near the

References (28)

  • M.T. Nair et al.

    J. Fluids Struct.

    (1997)
  • S.K. Lele

    J. Comput. Phys.

    (1992)
  • D. Gaitonde et al.

    J. Comput. Phys.

    (1997)
  • C.K.W. Tam et al.

    J. Comput. Phys.

    (1993)
  • Z. Haras et al.

    J. Comput. Phys.

    (1994)
  • X. Zhong

    J. Comput. Phys.

    (1998)
  • N.A. Adams et al.

    J. Comput. Phys.

    (1996)
  • M.H. Carpenter et al.

    J. Comput. Phys.

    (1993)
  • S.A. Orszag

    J. Comput. Phys.

    (1980)
  • A. Ditkowski et al.

    J. Comput. Phys.

    (2001)
  • R. Hixon et al.

    J. Comput. Phys.

    (2000)
  • M.J. Lighthill

    Computational Aeroacoustics

    (1993)
  • A. Taflove

    Computational Electrodynamics: The Finite-DifferenceTIme-Domain Method

    (1995)
  • D. Gottlieb et al.

    Numerical Analysis of Spectral Methods

    (1977)
  • Cited by (170)

    • An efficient fourth-order three-point scheme for solving some nonlinear dispersive wave equations

      2023, Communications in Nonlinear Science and Numerical Simulation
    View all citing articles on Scopus
    View full text