Origin of Numerical Dissipation

Discretization of a partial differential equation into a finite difference equation, in general, leads to numerical dissipation in addition to numerical dispersion. To illustrate this, again consider solving the convective wave equation (4.11) using a large-stencil finite difference scheme. Suppose the spatial derivative is approximated by a finite difference quotient and solve time exactly. The exact solution of the semidiscretized problem is given by Eq. (4.21) and, upon evaluating the m-integral, the solution is given by Eq. (4.23).

TO

u(x, t) = J Ф(a)ei(ax-ca(a)t]da. (4.33)

— TO

Now, whether the numerical solution is damped or not depends critically on whether a central difference stencil or an unsymmetric difference stencil is used. If a central difference stencil is used, then a(a) is a real function for real a. In this case, Eq. (4.33) represents a dispersive wave packet without being damped in time. If an unsymmetric stencil is used, then a(a) is complex for real a. In this case, the solution is damped in
time provided that Im(ca) is negative for all a. The time rate of damping for the wave component with wave number a is given by c Im[a(a)]. For many popular numerical schemes c Im(a) is negative. Such schemes have built-in numerical damping. The damping rate can be calculated precisely once the stencil size and stencil coefficients are known.

It is worthwhile to point out that it is possible that c and Im(a) have the same sign. In this case, the numerical solution will grow exponentially in time leading to numerical instability. Therefore, a necessary condition for numerical stability is

c Im [a(a)] < 0, —n < aAx < n.

Eq. (4.34) is the upwinding requirement. It is easy to show numerically that only upwind unsymmetric stencils could satisfy condition (4.34). The upwinding require­ment is extremely important when solving the full Euler equations. Euler equations support several modes of waves (acoustic, vorticity, and entropy waves). The upwind – ing condition must be satisfied by each wave mode if the numerical solution is to remain stable.

Recently, a number of large-stencil upwind schemes were proposed by Lockard et al. (1995), Li (1997), Zhuang et al. (1998, 2002), Zingg (2000), and others for use in CAA. Each of these authors offered examples to demonstrate the numerical damping of his/her upwind scheme that could, indeed, damp out spurious waves to produce an acceptable solution. However, strictly speaking, the spatial derivative terms of the Euler equations are wave propagation terms. They are not damping terms. To use the discretized form of the derivative terms to perform two functions, namely, to support wave propagation and to damp out spurious short waves at the same time, might be imposing too many constraints on the stencil coefficients. This would make it less likely that the stencil design is optimal. A more natural strategy is to use a central difference scheme and incorporate needed numerical damping by the addition of artificial selective damping terms.

Temporal discretization also introduces numerical damping. The damping rate can be determined quantitatively from the dispersion relation. Because ш(ш) is complex, the root ш = ш(а) of Ш(ш) = ca(a) has an imaginary part even when a(a) is real for real a. For the convective wave equation, the exact solution of the fully discretized system is given by Eq. (4.19). The integrals can be evaluated by the Residue Theorem. This leads to the replacement of ш by ш = ш(а), the physical root of ш(ш) = ca(a). Since ш(а) is now complex, the solution is damped in time. The time rate of damping for wave component with wave number a is Im[®(a)].

A simple way to reduce numerical damping due to temporal discretization is to use a smaller At. It is not too difficult to see why reducing time step size would reduce numerical damping. Suppose a spatial discretization scheme is chosen so that the problem would be solved accurately by using N mesh points per wavelength. In other words, if X is the shortest wavelength, the scheme is designed to resolve, then X/Ax > N. This yields, in terms of wave number, 2n/N > alargestAx. Now ш is related to a by the dispersion relation (4.20), thus,

By choosing At /Ax small, «largestAt would be confined to a range that is closer to zero. For most temporal marching schemes, including the Runge-Kutta method, the DRP scheme, the «At versus «At relation is such that Im(«At) and, hence, Im(«At) become smaller when «At is closer to zero. It follows that the root «(a) of Eq. (4.20) would be nearly real when At is reduced. This means that the temporal damping rate Im(«(a)) would be reduced.