Courant Number and Stability Condition

In 1928, a few German mathematicians named Richard Courant, Kurt Friedrichs, and Hans Lewy derived a stability condition for explicit time-marching solution for partial differential equations in their original paper titled "On the partial differential equations of mathematical physics". This is considered a necessary condition for convergence while solving certain partial differential equations (usually hyperbolic PDEs) numerically using explicit time-marching schemes. This condition states that the ratio of the distance traveled by the fluid in a single time step (`v*dt`) to the spatial step size (`dx` - mesh element size) should be below a certain value as dictated by the condition otherwise, the simulation will produce an unphysical result (i.e., solution blows up/diverges). This ratio is what is defined as the Courant number or the CFL number.

This stability condition can be written as 

 `C=v*(Deltat)/(Deltax)le1`

where C is the Courant Number.

This equation says that `Deltatle(Deltax)/v` for the numerical solution to be stable. This equation is called the Courant-Friedrichs-Lewy condition generally written as CFL condition.

 

How is the Courant Number determined by a Solver?

1-D Analysis

 

Let us take a rectangular cell with its flow moving from left to right with the velocity U and length of the cell being Δx.

 

 

Over a time step Δt, the flow moves a distance UΔt across the cell. The Courant number is the ratio of these two lengths or the fraction of the cell that moves across in a given time step. Hence it is a dimensionless quantity.

 Co = `"Fluid Distance"/"Cell Distance"=(UDeltat)/(Deltax)`

  

If Co =0.3, it means that the flow distance covers 30% of the length of the cell, and if Co =1.1, then the flow of fluid covers or surpasses the length of the cell in a single timestep.

 `therefore`, the Courant Number expresses that the distance that any information travels during the timestep length within the mesh must be lower than the distance between mesh elements. In other words, information from a given cell or mesh element must propagate only to its immediate neighbors.

In 3D, the system courant behaves differently keeping the basic definition the same because of the irregular shape of the 3D cell or skewness in the element quality. For cells that are skewed/complex 3D shapes like a tetrahedron, pyramid, prism, etc the procedure extends as follows:

We need to specify the distance across the cell Δx and the velocity U

`implies Deltax=`   `"Cell Volume" /"Total Surface Area"`

Total surface area can be defined as the sum of the individual surface area/length of the face of the cell.

`A=A_1+A_2+A_3+...A_N`

In the case of  3D cells

 

 

 

     

 In 1D cells

 

We need the velocity normal to the face, which happens to be `(U_f*hat n_f)`

The shaded area is the distance moved across the cell in one-time step.

 Co = `"Fluid Distance"/"Cell Distance"=(UDeltat)/(Deltax)`

Substituting for  Δx and the velocity normal to the face, we get the equation as

`C_o=sum_(Faces)(U_f*hat n_f)*Deltat*A_f/V_P`

But if we look at it closely, we can notice that there is a problem here  `sum_(Faces)(U_f*hat n_f)*A_f=0`

We are actually summing the velocity normal to the face and the face area(which has the units of total volume flux) across the faces of the cell is zero due to conservation of mass.

 

To avoid this, we need to take the magnitude and divide it by two to get the actual velocity flux normal to the face.

The flux from the right face is reversed

 

 Initial case   

Example,

         `because -5+5=0`     `(U_f*hat n_f)`

          `5+5=10`                `abs(U_f*hat n_f)`

          `1/2*(5+5)=5`       `1/2*abs(U_f*hat n_f)`

 

Rearranging slightly to get the equation as follows for an arbitrary shaped 3D cell :

`C_o=1/2*Deltat*(sum_fabs(U_f*hat n_f)*A_f)/V_P`

The CFL condition governs the stability of the convective phenomenon (i.e., hyperbolic nature). Using the above relation the Courant number for 1D, 2D, and 3D can be written as

 

 

There is also the diffusion phenomenon that is a part of fluid transport. In explicit CFD simulations, this phenomenon is also governed by a stability condition, generally referred to as the Diffusion stability condition or Diffusion CFL condition. This condition can be written as

[Diffusion CFD condition]

where d = [LHS of D CFL Cond]  is called the Diffusion number.   

This condition can be extended to 2D, and 3D is as shown below;

[Diffusion stability criteria table]

 

 

NOTE :

• The above stability criteria are for EXPLICIT schemes ONLY !!
• In IMPLICIT schemes, stability can be maintained over much larger values of Δt than for a corresponding explicit method. Some IMPLICIT schemes are UNCONDITIONALLY STABLE, ie any CFL can be applied to them !!!

Courant Number Field

In reality, every cell in the mesh has a different Courant Number and hence it is a field that can be plotted like any other variable.

 

 

This is an example of Courant number variation in a backward-facing step.

You can even see this in OpenFOAM simulations, they have a mean value and a maximum value at different timesteps.

 

 

CFD solvers use the maximum Courant Number to set the time step.

 

Courant Number and Numerical Diffusion

Numerical diffusion [what is it?]. This is caused as a result of neglecting the higher-order terms in the Taylor series. The degree of numerical diffusion depends on numerical viscosity which is given by the below-mentioned formula.

`mu_("artificial")=(c*Deltax)/2*(1-C) implies prop(1-c*(Deltat)/(Deltax))`

Also in any explicit CFD simulation, the Courant number must be equal to or smaller than 1 otherwise, the numerical viscosity would be negative. A negative viscosity is unphysical and will cause the solution to blow up. Let us understand this more intuitively using the method of characteristics.

 

CASE-1:  `C_ogt1`

The CFL condition, i.e the Courant number must be less than or at most equal to unity, is also the stability condition that holds for the second-order wave equation shown below:

`(del^2u)/(delt^2)=c^2*(del^2u)/(delx^2)`

There is a connection between the characteristic lines associated with this hyperbolic equation and the CFL condition. The characteristic lines are given by 

 `x={(ct,,,=,"right-running"),(-ct,,,=,"left-running"):}`

 

Let `Deltat_"C=1"`  denote the value of  `Deltat`  when C=1

`=>Deltat_(C=1)=(Deltax)/c`     `because ("C"=c(Deltat)/(Deltax)le1)`

 

Illustration of an unstable case

The numerical domain does not include all the analytical domain. 

C > 1 leads to unstable behavior (where  `Delta"le"" "DEF`   is the

analytical domain and `Delta"le"" "DAC`   is the numerical domain)

 

Information from one point to another can propagate only along lines of characteristics. At point `b`,   we can draw information from i+1 and i-1, since the characteristics at point `b`   go through i+1 and i-1 nodes.

The yellow color circled point  `d` is the one under consideration in all 3 cases

 

CASE-2:  `C_ole1`

 

 

The numerical domain includes all the analytical domains.

(where `Delta"le"" "DPQ`   is the analytical domain and `Delta"le"" "DAC`   is the numerical domain)

CASE-3:  `C_o` << `1`

 

 

Illustration of a stable but inaccurate case.

The numerical domain includes all the analytical domains.

(where `Delta"le"" "DPQ`   is the analytical domain and `Delta"le"" "DAC`   is the numerical domain)

If we apply too small a time step `Deltat_d` represented by point `d`,   then we cannot draw information from points i-1 and i+1 because they lie outside the characteristics running from point `b`. However, in the numerical algorithm (CFD code) we can force the method to still draw information from i-1 and i+1 (represented by the dashed lines) which introduces inaccuracies in predicting the values at point `d`. This inaccuracy is a result of numerical diffusion, which is inversely proportional to the Courant number.

This implies that, if we take `Deltat_d` ie `Deltat_(C<1)` to be very small(`Deltat_d` << `Deltat_(C=1)`), then even though numerical calculations are stable, the results may be quite inaccurate due to the large mismatch between the domain of dependence of point `d` and the location of the actual numerical data used to calculate properties at `d`.

Therefore, while for stability reasons one has to choose a CFL number below the appropriate stability criterion, we should stay as close as possible to this upper limit for accuracy reasons.

 

Types of Time Stepping - Local vs Global Time-Stepping

The CFL number for practical flows can be expressed as:

  `CFL=((u+a)_max*Deltat)/(Deltax)`, where u is the local velocity in the cell

                                                            a is the local speed of sound

                                                           (u+a) is right running wave speed

To satisfy stability criteria, one might set a constant CFL number for every cell. However since both (u, a, Δx) change over the computational domain, clearly each cell will have a different Δt which we call the “local time step”.

`Deltat_("local")=(CFL)_("global")*(Deltax)/((u+a)_max)`

This really means that each cell converges at a different -  fastest possible - rate. So instead of defining a “global time step” valid for all cells, it is much more convenient for steady-state solvers to define a “global CFL number”, which ensures that stability criteria are met in every cell and that overall convergence is maximized by “local time-stepping”.

Note that “local time-stepping” is only applicable to steady-state solvers where time accuracy is not important. For an unsteady solver, the same “real-time step” has to be applied for each cell.

In global time-stepping, we calculate the Δt at all the grid points and then choose the minimum value among them to use for the calculations so that solution does not blow up.

`Deltat="minimum"(Deltat_1^t,Deltat_2^t,...,Deltat_i^t,...Deltat_N^t)`

This time-marching solution is following the actual unsteady flow variations that would exist in nature ie the solution gives a time-accurate solution of the actual transient flow field. This approach is giving us physically meaningful transient variations - which frequently are of intrinsic value by themselves.

General Courant number range:

Different Courant -  Some Recommendations for stability in  Courant Number are :

• Compressible Flows (depends on the speed of sound)
• Explicit Time Stepping  (< 1)
• Free Surface Flows  (< 1)
• Large Eddy Simulation  (0.5 - 1)
• Turbomachinery  (< 80)

In most cases,  the maximum Courant number should be below 1.0.