Why are upwind schemes stable in convection flow calculation?

It is well known that upwind schemes are stable when calculating convection flows with $|\text{Pe}|>2$, $\text{Pe}$ is the Peclet number. Why is that, and why is central difference unstable?

Is there any intrinsic reason there?

• This is not a direct answer, so I'll just post a comment, but one article in this field which I found remarkably insightful is this: H.-G. Roos, "Ten ways to generate the Il'in and related schemes"
– cfh
May 23, 2015 at 22:32
• The truthfulness of the assertion in the question depends on how you discretize in time. For instance, (centered differences + backward Euler) is stable. May 27, 2015 at 18:28

The reasoning for the stability of upwind schemes based on an understanding of the characteristics of the hyperbolic equation(s). Characteristics are essentially the finite speeds at which information in a hyperbolic system travel, and are found via decomposing a hyperbolic system into independent hyperbolic PDEs.

Now characteristics are essentially just pushing along the initial conditions of a given hyperbolic equation (though nonlinear equations can distort them). The fact this speed is finite results in a need to be careful with your numerical stencil.

The typical example to illustrate the need to sample carefully is to imagine an initial condition of:

$$u_0(x) = \begin{cases} 1 & x < 2\\ 0 & 2 \leq x \end{cases}$$

Since the characteristic just pushes this initial condition along the space-time domain, its general shape and the center discontinuity remains. Now imagine some time later, you are aiming to evaluate the derivative at some x* location. There's a chance that the x* location is at the center of the discontinuity, where the derivative is technically undefined.

To counter this, we sample on one side. Now how do we decide which side to sample on? It should be the side that would change value first.. Which is the side the characteristic would touch first. This means that if the characteristic is moving in the positive x direction, we should sample to the left to ensure we capture any possible changes in the solution.

This sampling of data that is going in the opposite direction of the characteristic is known as upwinded sampling. This helps ensure numerical stability by only sampling data we know we have information for.

An intuitive reason is that upwinding can be viewed as providing additional numerical diffusion, which is typically associated with stability in standard Finite Difference/Finite Element schemes. In finite differences, this means that the first order upwind scheme can be rewritten

$$\frac{u_i-u_{i-1}}{h} = \frac{u_{i+1}-u_{i-1}}{2h} - \frac{h}{2} \frac{u_{i+1}-2u_i + u_{i-1}}{h^2}.$$

In other words, the first order upwind difference can be interpreted as adding additional artificial diffusion relative to the 2nd order central difference scheme.

The reason why central differences is unstable is a little more involved. IMO, it's easier to analyze stability in the Finite Element framework. For a variational problem $b(u,v) = l(v)$, the constant $\gamma_h$ in the discrete coercivity condition (where $\|u\|_h$ is some norm on $u$)

$$b(u,u)\geq \gamma_h \|u\|_h^2$$

gives a measure of the stability of the problem (sort of like the smallest singular value tells you the stability of a matrix equation). For standard Galerkin methods (which are related to 2nd order central finite differencing) and Discontinuous Galerkin with central fluxes, $\gamma_h$ is inversely proportional to the Peclet number, and grows more unstable as Peclet number increases. For example, take the 1D Galerkin scheme for constant advection diffusion $u' - \epsilon u = f$:

$$b(u,v) = \int_0^1 u'v + \epsilon \int_0^1 u'v' = \int f v$$

Assuming $0$ Dirichlet boundary conditions, $b(u,u)$ gives

$$\int_0^1 u'u + \epsilon \int_0^1 (u')^2 = \int_0^1 (u^2/2)' + \epsilon \int_0^1 (u')^2 = [u^2/2]_0^1 + \epsilon\int_0^1 (u')^2 = \epsilon\|u\|_{H^1}^2$$

which implies issues with stability as $\epsilon \rightarrow 0$. Note that this implies this issue with stability is usually worse for boundary layer problems (Dirichlet boundary conditions on both sides).

Upwinded Galerkin (i.e. SUPG, streamline diffusion, etc) and upwind DG restore discrete coercivity with respect to a specific norm $\|u\|_h$, such that $\gamma_h$ does not approach $0$ as Peclet number increases. This is also a part of why adding numerical diffusion tends to stabilizes advective flows - additional numerical diffusion adds a coercive term to the variational form.

• Is there any book or paper talk about this in a similar way like yours? It is very appreciated. May 11, 2015 at 5:56
• I took most of these from some class notes that Tom Hughes wrote up previously - I believe he was planning on making it a book, but I'm not sure if it's out yet. More mathematical discussions may be found in "Robust Methods for Singular Perturbation Problems" by Roos, Stynes and Tobiska springer.com/us/book/9783540344667. May 11, 2015 at 15:01

To address the second question on the instability of the central scheme, it is helpful to consider the scalar advection equation $$u_t + c u_x = 0$$ which is a simplification of the convection equation $$u_t + \nabla \cdot (\boldsymbol{c} u) = 0$$ to one dimension and constant velocity $$c \neq c(x)$$. One can show that the central difference scheme $$u_i^{(n + 1)} - u_i^{(n)} + \frac{c \Delta t}{2 \Delta x} \Big(u_{i+1}^{(n)} - u_{i - 1}^{(n)} \Big) = 0$$ is energy-unstable. Defining the discrete equivalent to the energy $$\int_\Omega \frac{u^2}{2} dx$$ as a Riemann sum $$E^{(n)} = \frac{\Delta x}{2} \sum_{i = 1}^M \Big( u_i^{(n)} \Big)^2$$ one can show that for zero or periodic boundary conditions the following equality holds: $$E^{(n + 1)} = E^{(n)} + \frac{\Delta x}{2} \sum_{i = 1}^M \Big( u_i^{(n + 1)} - u_i^{(n)}\Big)^2$$ Thus, unless the solution stays the same, the energy grows after every timestep for any choice of $$\Delta t, \Delta x$$.

Proof.

Multiply the discretized PDE with $$u_i^{(n)}$$ to obtain $$u_i^{(n)} \Big(u_i^{(n + 1)} - u_i^{(n)} \Big) + \frac{c \Delta t}{2 \Delta x} \Big(u_i^{(n)}u_{i+1}^{(n)} - u_i^{(n)}u_{i - 1}^{(n)} \Big) = 0$$ By means of the identity $$\alpha(\beta - \alpha) = 0.5 \Big(\beta^2 - \alpha^2 - (\beta- \alpha)^2\Big)$$ which holds for $$\alpha, \beta \in \mathbb{R}$$ the scheme can be rewritten as $$\frac{\Big(u_i^{(n + 1)}\Big)^2}{2} - \frac{\Big(u_i^{(n)}\Big)^2}{2} - \frac{\Big(u_i^{(n + 1)} - u_i^{(n)}\Big)^2}{2}+ \frac{c \Delta t}{2\Delta x} \Big(u_i^{(n)}u_{i+1}^{(n)} - u_i^{(n)}u_{i - 1}^{(n)} \Big) = 0$$ Now sum over all grid points $$i = 1, \dots M$$ to obtain the discrete energies: $$E^{(n + 1)} = E^{(n)} + \frac{\Delta x}{2} \sum_{i = 1}^M \Big( u_i^{(n + 1)} - u_i^{(n)}\Big)^2 - \frac{c \Delta t}{2} \sum_{i = 1}^M \Big(u_i^{(n)}u_{i+1}^{(n)} - u_i^{(n)}u_{i - 1}^{(n)} \Big)$$ the second sum on the right-hand-side turns out to be a telescope sum: $$\sum_{i = 1}^M \Big(u_i^{(n)}u_{i+1}^{(n)} - u_i^{(n)}u_{i - 1}^{(n)} \Big) = -u_1^{(n)} u_0^{(n)} + u_{M + 1}^{(n)} u_{M}^{(n)}$$ For zero boundary conditions $$u_0 = u_{M + 1} = 0$$ or periodic boundary conditions $$u_0 = u_M, u_1 = u_{M+1}$$ these two contributions vanish and one obtains the proposition.

Also, here is a pretty extensive discussion with results of various numerical schemes:

http://www.holzmann-cfd.de/index.php/en/numerical-schemes

I hope to answer your question in a general matter without any references to mathematical equations.

To be stable, a numerical scheme needs to capture the physics of the equation you are discretizing. Without going into the details, each term in your equation has eigenvalues which define the direction in which the information of your flow travels.

For example, when discretizing the diffusion term in the Navier-Stokes Equations (NSE), you should use a central scheme since diffusion, by definition, happens in all directions. If you were to use a one-sided numerical scheme, you could make your scheme stabilized by satisfying the CFL condition, but your results may become non-physical, since it can't exactly capture expansion waves. Where we should use Godunov's family schemes. Please note that Godunov's family schemes use upwinding with some constraints. Now, if you are solving the NSE for a supersonic flow, where the flow characteristics travel in ONLY one direction, then you MUST use an upwind scheme for stability. Good reference for upwind schemes in high speed flow is a paper by van Leer

• anyone care to comment on the downvote? My answer is pretty much the same as the accepted answer without the charasteristics background May 28, 2015 at 9:14
• Yes, your answer is the simplest one to understand. And it is correct within my knowledge. Jun 6, 2015 at 9:31