Numerical Instability
Table of Contents
Many numerical approximations to partial differential equations are unusable because they produce unstable computational results. Computational stability issues have been discussed in two previous articles in the CFD-101 series: Computational Stability and Heuristic Analysis. In this article, a simple mechanical model is described that leads to an understanding of common numerical instabilities associated with approximations of the Navier-Stokes equations. In particular, the simple model applies to instabilities arising from fluid dynamic forces such as viscous stress, surface tension, elasticity and more. The stability conditions obtained with the simple model do not depend on any particular numerical approximations, but instead involve only generic considerations of mass and forces.
The Model System
Imagine a mass M located between two rigid walls and connected to the walls by springs, as shown in Fig. 1. Assuming that the springs satisfy Hook’s law in which a spring force on the block is proportional to the change in length of the spring, an equation of motion for the block, which moves in only the horizontal direction is,
(1) $latex \displaystyle M\frac{\partial U}{\partial t}=-2k\left( X-{{X}^{0}} \right)$
Symbol X0 indicates the initial x position of the block and M is the block mass. Initially X=X0 and the block is at rest, U=0. Now imagine a perturbation given to the block by assigning a velocity of U0 at time t=0. After a small time interval δt the block moves to position X1=X0+U0δt and a simple discretization of Eq. 1 gives the velocity at the end of the time interval as,
(2) $latex \displaystyle M\left( \frac{{{U}^{1}}-{{U}^{0}}}{\delta t} \right)=-2k\left( {{X}^{1}}-{{X}^{0}} \right)$
Replacing X1 by its value X0+U0 δt and rearranging gives an equation for the new velocity U1,
(3) $latex \displaystyle {{U}^{1}}={{U}^{0}}\left( 1-2\frac{k\delta {{t}^{2}}}{M} \right)$
This is a recursion in which for each successive time-step the velocity at the end of the time step is equal to the previous value of the velocity times the bracketed quantity in Eq. 3, so that after n time steps,
(4) $latex \displaystyle {{U}^{n}}={{U}^{n-1}}\left( 1-2\frac{k\delta {{t}^{2}}}{M} \right)={{U}^{0}}{{\left( 1-2\frac{k\delta {{t}^{2}}}{M} \right)}^{n}}$
Note that the superscript n on the bracketed quantity in Eq. 4 is an exponent, not a time level index, although its value in this case is the same as the time level index. From Eq. 4 we see that if the quantity in the bracket has an absolute magnitude larger than 1.0 the velocity Un will increase exponentially with increasing n. Thus, to prevent the exponential growth of the velocity in this model system, the time-step size must be limited to satisfy the following inequality,
(5) $latex \displaystyle \frac{k\delta {{t}^{2}}}{M}\le 1$
When the left hand side of Eq. 5 is greater than one, the velocity Un will oscillate between positive and negative values on consecutive time steps while exponentially increasing in magnitude.
This behavior is characteristic of a classical numerical instability. In this case, Eq. 5 shows that the instability can be prevented by keeping δt small enough to satisfy the inequality. As a general rule when a numerical instability occurs and exhibits the character of increasing plus and minus values on successive time steps it can be cured by reducing the time-step size.
Exploring this simple mechanical model further we can see from Eq. 4 that the instability results from an overreaction to an initial action. That is, when the time-step size large, Eq. 4 predicts a new velocity in the opposite direction and with a larger magnitude. This excessive velocity then becomes the starting condition for the subsequent time step, leading to an exponential increase in the velocity magnitude.
The stability condition in Eq. 5 is based on an explicit formulation, meaning that the current response of the mass is expressed in terms of the previous displacement. An implicit formulation,where the current response of the mass is based on the subsequent position of the mass (i.e., using X2 in Eq. 2 instead of X1) would likely be unconditionally stable, but it requires a knowledge of the unknown final position X2. For most equations implicit methods require an iterative solution, and the additional computational effort required for such solutions is the price that must be paid to eliminate the stability condition.
Application to the Navier-Stokes Equation
To use the above mechanical model of a numerical instability to understand instabilities that may occur in the Navier-Stokes equation imagine two elements of an Eulerian computational grid in which a perturbation is made to a velocity on the boundary separating two elements, as shown in Fig. 2.
For simplicity, think of the elements outlined by solid lines as cubes of equal size, and that the vector represents a velocity in the x direction. The dashed lines in the centers of the elements (i.e., y-z planes) define the extent of the partial volumes of the elements assigned to the u velocity. The total mass of fluid in the partial volumes correlates to the mass M in the mechanical model. When the velocity U moves fluid between elements the elements respond by generating forces that act to counter the velocity, much like the springs in the mechanical model. These forces may arise because of compression or expansion of the fluid, viscous stresses, surface tension (if there is a fluid interface within the fluid mass M) or other forces. By identifying the appropriate stiffness coefficient k in each case we can use Eq. 5 to arrive at a stability criterion for that physical process when using an explicit numerical approximation in a grid like that shown in Fig. 2
Several examples of how this analogy can be applied are given in the following sections. In each case the mass M of the fluid in the partial volumes is given by,
(6) $latex \displaystyle M=\rho \delta x\delta y\delta z,$
where ρ is the density of the fluid and elements have dimensions δx, δy and δz.
Compressible Fluids
To find the stiffness coefficient, k, recall that k is a measure of the force generated to resist an applied perturbation. For a compressible fluid this force is related to the change in fluid pressure because of a change in fluid density according to the thermodynamic relation dp=c2dρ, where c is the speed of sound in the fluid. The fluid mass moved across the boundary between the elements is ρUδt*δyδz and the change in density in the element receiving the mass is this mass change divided by the volume of the element, dρ=ρUδt/(δx). The corresponding change in element pressure is then given by
(7) $latex \displaystyle dp={{c}^{2}}d\rho =\frac{\rho {{c}^{2}}U\delta t}{\delta x}.$
The force responding to the U velocity perturbation in each element is the product of the pressure change from Eq. 7 and the cross sectional area of the element δyδz. The effective stiffness of an element k, is therefore the force divided by the initial displacement Uδt,
(8) $latex \displaystyle k=\frac{\rho {{c}^{2}}U\delta t\delta y\delta z}{\delta xU\delta t}=\frac{\rho {{c}^{2}}\delta y\delta z}{\delta x}$
Substituting this value for k and the definition for M, from Eq. 6, into the stability condition Eq. 5 results in the stability condition for compressible fluids,
(9) $latex \displaystyle \frac{k\delta {{t}^{2}}}{M}={{\left( \frac{c\delta t}{\delta x} \right)}^{2}}\le 1.$
This is the well-known Courant condition than restricts the distance a sound wave travels in one time step to be less than the width of a computational element. An analogy with the simple mechanical model has provided this result without the need to write out an equation for pressure waves in a compressible fluid and then perform a stability analysis on that equation.
Viscous Stresses
The viscous forces that are generated in response to a perturbed velocity U in a fluid of viscosity μ consist of shears in the x, y and z directions. For example, the shear stress on the lower surface of the element for the velocity U is μU/δz, assuming that the velocity in neighboring cells is zero. There is a corresponding stress at the upper surface of the element. Each of these stresses act on a surface of area δxδy (in the current example) to produce a viscous force. Similarly, there are stresses in the x and y direction acting on their corresponding areas. In each direction there are force pairs (similar to the two springs) but because of Eq. 5 it is necessary to use the effective k for a single spring. This is half of the total of all the the viscous forces divided by the initial displacement Uδt,
(10) $latex \displaystyle k=\mu \left( \frac{U\delta y\delta z}{\delta x}+\frac{U\delta x\delta z}{\delta y}+\frac{U\delta x\delta y}{\delta z} \right)\frac{1}{U\delta t}.$
Inserting this value for k and using Eq. 6 for M into the stability condition for the mechanical model given in Eq. 5 yields
(11) $latex \displaystyle \frac{k\delta {{t}^{2}}}{M}=\frac{\mu }{\rho }\left( \frac{1}{\delta {{x}^{2}}}+\frac{1}{\delta {{y}^{2}}}+\frac{1}{\delta {{z}^{2}}} \right)\delta t\le 1.$
Equation 11 is the stability condition for explicit viscous stresses approximated in an Eulerian grid.
Surface Tension
Imagine a fluid interface located between the two elements that is deformed by a U velocity perturbation, as shown in Fig. 2A. In this case the reaction is a surface tension force on each of the segments of the surface illustrated in Fig. 2A. For simplicity, assume a two-dimensional surface (e.g., no variation in the z direction) and a constant surface tension coefficient σ.
Surface tension in each segment acts tangentially along the surface so the force responding to the U velocity is the x component of that force, i.e., the surface tension coefficient times the sine of the angle of the surface segment with respect to the vertical. For a small initial displacement, the x-force from each surface segment can be approximated by σUδtδz/δy giving the resulting stiffness coefficient,
(12) $latex \displaystyle k=\left( \frac{\sigma U\delta t\delta z}{\delta y} \right)\frac{1}{U\delta t}$
Substituting this k and M into Eq. 5 gives the stability condition for surface tension,
(13) $latex \displaystyle \frac{k\delta {{t}^{2}}}{M}=\frac{\sigma }{\rho }\frac{\delta {{t}^{2}}}{\delta x\delta {{y}^{2}}}\le 1.$
This result appears somewhat odd because of the different exponents of the δx and δy factors, but for cubic or square elements this makes no difference. For non-uniform elements, however, we should perform a similar evaluation in the y and z directions and then use the most restrictive of the results. In any case, this is a reasonable and useful result from a very simple model based on action and reaction principles.
Bulk Elasticity
For a fluid with elastic properties the U velocity perturbation is resisted by an elastic stress in a way that closely resembles a spring. If ε is the bulk modulus of the fluid then the stress associated with extension or compression in an element is εUδt/δx. This stress acts over the surface area δyδz, and the stiffness k is this force divided by the displacement Uδt. Substituting this into Eq. 5 provides the stability condition,
(14) $latex \displaystyle \frac{k\delta {{t}^{2}}}{M}=\frac{\varepsilon }{\rho }\frac{\delta {{t}^{2}}}{\delta {{x}^{2}}}\le 1.$
Similar results exist for the y and z directions.
Concluding Remarks
A simple mechanical model has been used to illustrate a common type of numerical instability, that arises from an action-reaction process. Using this simple model it is possible to quickly derive a variety of stability conditions for fluid dynamic forces modeled by explicit finite difference approximations in an Eulerian grid. The stability conditions are derived by using mass and force concepts in fluids that are analogous to the mass and forces in the model mechanical system. Significantly, the stability conditions arrived at are generic and do not depend on specific finite-difference approximations. Additionally, these derivations provide a simple way to understand the mechanisms driving the unstable behavior.
The approach taken here could be extended to other types of physical forces (e.g., electrical, non-inertial, etc.) and even advective processes could be included by using the analogy that a change in momentum resulting from advection could be thought of as the result of an equivalent force.
Furthermore, more refined estimates of the stiffness coefficient, for instance, by including more dimensional effects, could be added to enhance the stability conditions. In any case, the object here is to show that numerical instabilities can often be understood from a simple analysis. It is hoped that the insight this provides might guide the development of more robust and accurate numerical approximations.
