The CFL condition, Von-Neumann stability analysis and related topics 71

The stability conditions so far have been derived using ordinary differential equations only.

If we consider linear partial differential equations instead, more insight can be gained. The conditions obtained in this way are less strict than the SSP conditions, but nevertheless useful in practice. Furthermore, the maximal stable time step ∆t_EE for the explicit Euler method used in the SSP theory can be determined.

First of all, for the hyperbolic inviscid terms, we know that the timestep size has to satisfy the Courant-Friedrichs-Levy (CFL) condition that the domain of dependence of the solution has to contain the domain of dependence of the numerical method [41]. This is illustrated in figure 4.1.

The CFL condition is automatically satisfied by implicit schemes, regardless of the time step. For finite volume methods for one dimensional equations with explicit time integration, we obtain the constraint

∆t < CF L_max· ∆x

k,|n|=1max |λ_k(u,n)|, (4.9)

with the maximal CFL numberCF L_max depending on the method, for example 1.0 for the explicit Euler method. The λ_k are the eigenvalues (2.30) of the Jacobian of the inviscid flux. When comparing this constraint to the one that would come out of applying the eigenvalue analysis as for the linear test equation, it turns out that this constraint is twice as strong, meaning that the time step is twice as small.

Figure 4.1: Illustration of the CFL condition for a linear equation: the exact solution in the point (x_i, t_n+1) can be influenced by the points in the shaded area. The numerical domain of dependence (the grey area) has to contain the shaded area for the scheme to be stable, resulting in a constraint on ∆t/∆x.

It can be shown that the condition (4.9) is the result of a linear stability analysis for an upwind discretization, namely the von Neumann stability analysis. This is, besides the energy method, the standard tool for the stability analysis of discretized partial differential equations. To this end, the equation is considered in one dimension with periodic bound-ary conditions. Then, the discretization is applied and Fourier data is inserted into the approximation. The discretization is stable if no component is amplified, typically leading to conditions on ∆t and ∆x. Unfortunately, this technique becomes extremely complicated when looking at unstructured grids, making it less useful for these. For linear problems, the analysis results in necessary and sufficient conditions. This means that the CFL condition is sharp for linear problems.

For nonlinear equations the interactions between different modes are more than just superposition. Therefore, the von Neumann stability analysis leads to necessary conditions only. Typically, there are additional stability constraints on top of these. For the case of the viscous Burger’s equation, this has been demonstrated in [131]. For general equations, no complete stability analysis is known.

If additional parabolic viscous terms are present, like in the Navier-Stokes equations, the stability constraint changes. For a pure diffusion equation, resp. the linear heat equation, this is sometimes called the DFL condition and for a finite volume method given by

4.2. STABILITY 73

∆t < ∆x²

2 , (4.10)

where is the diffusion coefficient. Here, the dependence on the mesh width is quadratic and thus, this condition can be much more severe for fine grids than the CFL condition.

One way of obtaining a stable time integration method for the Navier-Stokes equations is to require the scheme to satisfy both the CFL and DFL conditions. However, this is too severe, as seen in practice and by a more detailed analysis. For example, the restriction on the time step in theorem 3 is of the form

∆t < α³∆x² α max

k,|n|=1|λ_k(u,n)|∆x+. (4.11)

Here, α is a grid dependent factor that is roughly speaking closer to one the more regular the grid is (see [150] for details). A similar bound has been found in praxis to be useful for the Navier-Stokes equations:

∆t < ∆x²

max

k,|n|=1|λ_k(u,n)|∆x+ 2.

In the case of a DG method, both time step constraints additionally depend on the order of the polynomial basis in that there is an additional factor 1/(2N + 1) for the CFL condition and of 1/N² for the DFL condition. The dependence of the stability constraint on the choice of the viscous flux has been considered in [105].

In more than one space dimension, the situation is, as expected, more difficult. It is not obvious which value to choose in a cell for ∆x and no stability analysis is known. A relation to determine a locally stable time step in a cell Ω that works on both structured and unstructured grids for finite volume schemes is the following:

∆t=σ |Ω|

λ_c,1+λ_c,2+λ_c,3+ 4(λ_v,1 +λ_v,2+λ_v,3). Here, λ_c,i = (|v_i|+c)|s_i| and

λ_v,i = max 4

3ρ,γ p

µ P r

|s_i|

|Ω|,

where s_i is a projection of the control volume on the planes orthogonal to the x_i direction [214]. The parameterσ corresponds to a CFL number and has to be determined in practice for each scheme, due to the absence of theory. Finally, a globally stable time step is obtained by taking the minimum of the local time steps.

4.3 Stiff problems

Obviously, implicit methods are more costly than explicit ones. The reason they are consid-ered nevertheless are stiff problems, which can loosely be defined as problems where implicit methods are more efficient than explicit ones. This can happen since broadly speaking, im-plicit schemes have better stability properties than exim-plicit ones. If the problem is such, that a stability constraint on the time step size leads to time steps much smaller than useful to resolve the physics of the system considered, the additional cost per time step can be justified making implicit schemes the method of choice. To illustrate this phenomenon, consider the system of two equations

x_t = −x, y_t = −1000y.

The fast scale quantity y(t) will decay extremely fast to near zero and thus, large time steps that resolve the evolution of the slow scale quantityx(t) should be possible. However, a typical explicit scheme will be hampered by stability constraints in that it needs to chose the time step according to the fastest scale, even though this has no influence on the solution after an initial transient phase. This example illustrates one possibility of characterising stiffness: we have a large Lipschitz constant of our function, but additionally eigenvalues of small magnitude. In this case, the longterm behavior of the solution is much better characterised by the largest real part of the eigenvalues, instead of the Lipschitz constant itself [47]. This is connected to the logarithmic norm. Another possibility of obtaining stiffness are stiff source terms. Note that if stiffness manifests itself also depends on the initial values chosen, as well as on the maximal time considered, which is why it is better to talk of stiff problems instead of stiff equations.

The Navier-Stokes equations have similar properties as the above system with both slow and fast scales present. First of all for the Euler equations, there are fast moving acoustic terms which correspond to the largest eigenvalue. Then, there are the significantly slower convective terms, which correspond to the convective eigenvalue |v|. The time step size of explicit schemes will be dictated by the acoustic eigenvalue, even if the physical processes of interest usually live on the convective scale. This property becomes extreme for low Mach numbers, when the convective eigenvalue approaches zero and the quotient between these two eigenvalues becomes very large, e.g. stiffness increases.

In the case of the Navier-Stokes equations with a bounding wall, we have additionally the boundary layer. This has to be resolved with extremely fine grid spacing in normal direction, leading to cells which are several orders of magnitude smaller than for the Euler equations. For an explicit scheme, the CFL condition therefore becomes several orders of magnitude more restrictive, independent of the flow physics. Furthermore, large aspect ratios are another source of stiffness, sometimes called geometric stiffness. Therefore, im-plicit methods play an important role in the numerical computation of unsteady viscous