Convergence Criteria

No material is truly incompressible, but this assumption is often a good approximation. When using this assumption in connection with a numerical solution scheme it is necessary to devise some way to impose the physical mechanism that is responsible for the incompressible behavior.

(1)     $latex \displaystyle \frac{d\rho }{dt}+\rho \nabla \cdot \vec{U}=0$

The measure of incompressibility is the divergence of the velocity, Ñ ·U, which is equal to the time rate-of-change of fluid volume per unit volume. The divergence appears in the mathematical statement of mass conservation, where dr/dt is the change in density of a fluid element as it moves with the fluid. If the fluid is incompressible then its density should remain constant and, according to Eq.(1), this implies the velocity divergence should be zero:

$latex \displaystyle \nabla \cdot \overrightarrow{U}=0$

The physical mechanism that leads to incompressible behavior is the rapid propagation of pressure waves, which must move through a fluid much faster than the material speed of the fluid. Most often, the numerical propagation of pressure waves is accomplished by some sort of iteration scheme that couples pressures to velocities. The goal of the iteration is to reduce the magnitude of the velocity divergence below some absolute numerical value, e, called the convergence criteria.

Selecting a Convergence Criteria

A question frequently asked is how should be be chosen? The answer to this is simple, and yet subtle, and depends on the numerical method that is selected for computing a solution.

A discretized version of Eq.(1) indicates the change in density occurring in one time step, dt, is given by:

(2)     $latex \displaystyle \frac{\delta \rho }{\rho }=\delta t\nabla \cdot \overrightarrow{U}$

There are a couple of things to observe about this result. First, the right side has the general order of magnitude of:

$latex \displaystyle \delta t\nabla \cdot \overrightarrow{U}\approx \frac{\delta tU}{\delta x},$

where d x is the length or size of the discrete elements used in the numerical solution. This quantity is often referred to as the flow Courant number. It is well known that the Courant number should be less than one for accurate results. It must also be less than one for numerical stability when using explicit computational methods.

We see from relation (2) that if the Courant limit of one is exceeded, the change in density could be greater than the density itself! This observation is the reason why the Courant limit is so closely tied to issues of accuracy and stability. It also explains why stability conditions become more restrictive as the number of spatial dimensions increases from one to two or three. For stability, the accumulated change resulting from fluxes in all directions must not be allowed to change the density (or other quantity) by more than its current value. Thus, with more dimensions there are more fluxes and a smaller time-step size is needed to ensure that elements are not over emptied.

Returning to relation (2), we see how much change in density to expect from a given convergence criteria e. For example, if e = 0.001 the expected change in density in one unit of time would be one tenth of one percent, certainly small enough to be considered “incompressible” for most applications.

A value of e = 0.001 is the nominal value that is often suggested for incompressible flow calculations. In the case of natural convection or other problems where very small density variations are important, the value of e may have to be smaller, depending on how the governing equations are formulated.

If the value of e is too large, and the pressure solution method employs some amount of over-relaxation, the fluid may behave as though it is compressible. In this case numerical pressure waves might be observed to travel across the computational region. Reducing the value of e, or the amount of over-relaxation, reduces the amplitude of these artificial waves.

One question that needs to be answered is whether or not errors in density generated in one time step might accumulate to something more significant over a large number of time steps. The answer to this question is no, if the proper numerical method is employed.

Selecting a Self-Correcting Numerical Method

When numerical iteration methods are employed, it is generally thought that very fine convergence criteria are required to produce accurate results. This in not true, however, if the equations being solved can be formulated in a “self correcting” way. The general idea can be illustrated in connection with the equations for an incompressible fluid, where the incompressibility condition is a zero velocity divergence, D:

(3)     $latex \displaystyle D=\nabla \cdot \overrightarrow{U}=0.$

The fluid pressure, p, must satisfy an equation that is derived by taking the divergence of the momentum equation (e.g., the Navier-Stokes equation). Retaining all terms, the pressure equation is:

(4)     $latex \displaystyle \frac{\partial D}{\partial t}={{\nabla }^{2}}\left( {p}/{\rho }\; \right)+\upsilon {{\nabla }^{2}}D-Q.$

where:

$latex \displaystyle Q=\nabla \cdot \left[ \nabla \cdot \left( \overline{UU} \right) \right].$

Here r is the fluid density and u its kinematic viscosity. This is a Poisson equation for pressure that is most efficiently solved each time step by an iteration process, for which there must be a convergence criteria.

To minimize the accumulation of iteration convergence errors over many time steps one can employ a high level of convergence at each time step, or use a self-correcting treatment that was first used in connection with the Marker-and-Cell (MAC) method (F.H. Harlow and J.E. Welch, Phys. Fluids 8, 2182 (1965)). In many numerical methods devised to solve Eq.(4) the temptation is to set the left hand side to zero because of Eq.(3). In the self-corrective procedure, however, the left side is replaced by:

$latex \displaystyle \frac{\partial D}{\partial t}\approx \frac{{{D}^{n+1}}-{{D}^{n}}}{\delta t}\to \frac{-{{D}^{n}}}{\delta t},$ where n indicates the n-th time step.

Retaining the $latex \displaystyle {{D}^{n}}$ term is a corrective procedure that tends to remove any error in velocity divergence left from the previous cycle. For instance, if D had a small positive residual in some element at the end of the n-th step, corresponding to a small expansion, then during step n+1 there is a source term added to the equation (i.e., -Dn/d t) to produce a corrective compression by the same amount.

Using this procedure cuts down the accumulation of incompressibility errors, even when a relatively coarse convergence criteria is used for the pressure equation. More details about self-corrective procedures can be found in the paper: C.W. Hirt and F.H. Harlow, J. Comp. Phys. 2, 114 (1967).

One distinct advantage of the self-corrective feature, which is not always appreciated, is its insensitivity to initial conditions. Some incompressible flow solvers require initial conditions that satisfy the incompressibility condition. Because of this, these methods sometimes require the solution of a separate problem to get usable initial conditions. In the MAC method, and other methods that use the self-correcting feature, any non-zero velocity divergence in initial conditions is automatically removed after the first cycle of calculation. The convenience of this feature for users cannot be over emphasized.

FLOW-3D uses the self-corrective procedure as well as an automatic setting of the convergence criteria that adjusts e to whatever is happening during the course of a solution. The latter is another user convenience that can be extremely useful for problems having a variety of flow stages.