Courant Number in CFD

The Courant number is given by the following formula:

where \(\vec U\) indicates the velocity vector, \(\Delta t\) represents time step and \(\Delta x\) is the characteristic size of the mesh (in 1D case). The Courant number indicates how much the information travels across a numerical grid.

In order to fully understand the Courant number definition, let us first consider the 1D case, where we have a cell of size (length) \(\Delta x\) and a velocity field \(\vec U\) – Figure 1.

Over the time step \(\Delta t\), the fluid moves \(\vec U \Delta t\) across the cell, which is marked in Figure 2. It is obvious that the bigger the time step \(\Delta t\), the further the fluid moves across the cell.

Based on that we can explain that physically the Courant number corresponds to the ratio of \(\vec U \Delta t\) and \(\Delta x\). Or we can think about it as the fraction of the cell that the flow moves across in the given time step \(\Delta t\).

We can now easily interpret the value of Courant number. Let us take as an example two values \(Co = 0.5\) and \(Co = 1.5\), and visualize them in Figure 3.

Therefore, the Courant number tells us how fast the fluid (information) travels across the cell within a given \(\Delta t\). Based on the equation above it can be also stated that the Courant number changes with the velocity and cell size.

The 1D example above explains in an easy way the physical meaning of the Courant number. However, OpenFOAM supports any shapes of cell (hexahedron, wedge, prism, pyramid, tetrahedron, tetrahedral wedge, polyhedra). Apart from that, those cells are most likely not perfect elements. The question appears how we can evaluate the Courant number on an arbitrary element with some imperfections of the shape.

Based on the previous equation of Co, we know that for any arbitrary element shape, the distance across the cell \(\Delta x\) and the velocity U must be specified. To define the distance \(\Delta x\), OpenFOAM divides the cell volume V by the total surface area A according to:

\(\Delta x = \frac{V}{A}\)

The total surface area is simply equal to the sum of each face area that closes the volume. It should be noted that such a measure of element size \(\Delta x\) is good for regular elements, because we can get a proper representation of the element edge. For non-regular elements, highly elongated this measure will give us incorrect element size estimation.

To calculate the Courant number we also need the velocity normal to the face, from each of the faces that make up the cell. The velocity normal to the face can be easily calculated by taking the dot product of the face velocity \(\vec U_f\) and unit vector \(\vec n\) normal to the face, pointing towards the cell centroid (\(\vec U_f \cdot \vec n\)) – Figure 4.

At this point we need to note that if we calculate the velocity normal to each face and sum up those terms, based on the continuity equation we will get zero value. It is so, because the mass flow in and out to each cell must give 0. To overcome this we will take to sum of the magnitude of (\(\vec U_f \cdot \vec n\)) for each face rather than sum of normal velocity components. This will reverse half of the flux so that the total is not zero. But in order to get the true flux, we need to take only a half of it – Figure 5.Taking the last piece of the information into account, it finally brings us to the following form of the Courant number:

A unique feature of OpenFOAM is the so-called mesh Courant number. The feature is used together with dynamic mesh capabilities and an example of a solver which uses the mesh Courant number is pimpleFoam.

The mesh Courant number refers to handling of dynamic meshes. When the control volume is changing (mesh motion), instead of using just the velocity of the fluid inside the given cell, we should take into account that the faces of the cell are moving. So, for the Courant number calculations, we rather take the relative velocity between the fluid and face. This concept comes from the so-called control volume equation or Reynolds’ transport theorem and the detailed description can be found in [Reynolds’ transport theorem].

Control volume equation assumes that conservation laws are expressed in terms of intensive properties – properties that do not depend on the amount, like density (in contradiction to extensive properties like mass, volume, energy, momentum). Let’s define any conserved intensive property \(\phi\) (for mass conservation \(\phi\)=1, for momentum conservation \(\phi\)=v, for conservation of a scalar, \(\phi\) represents the conserved property per unit mass), then the corresponding extensive property is given by:

\(\Phi = \int_{\Omega_CM} \rho \phi d\Omega\)

where \(\Omega_(CM)\) stands for volume occupied by the control mass. Therefore, the rate of change of mass for each control volume can be written as:

\(\frac{d}{dt} \int_{\Omega_(CM)} \rho \phi d\Omega = \frac{d}{dt} \int_{\Omega_{(CV)}} \rho \phi d\Omega + \int_{S_(CV)} \rho \phi (\vec v – \vec v_b) \cdot \vec n dS\)

where \(\Omega_{(CV)}\) is the control volume, \(S_(CV)\) is the surface enclosing the control volume, \(\vec n\) is the unit normal vector to \(S_CV\), \(\vec v\) is the fluid velocity and \(\vec v_b\) is the velocity of the moving surface of the control volume. The equation above describes the rate of change of the property within the control volume plus the net flux of it through the boundary faces of the control volume.

When working with dynamic meshes in OpenFOAM, what solver does is first to calculate the absolute flux from the mapped surface velocity. Then the relative flux to the mesh motion is calculated. Finally, the Courant number is calculated but based on the relative flux, and not the fluid velocity.

It is time to answer the question of how to set up correctly the time step in CFD analysis. Let us consider two scenarios that can be found in SimFlow – link do opisu panelu RUN.

Before discussing the time step size we should also understand what happens with the information on the numerical mesh when the Courant number is defined below and above unity. It is explained in Figure 6.

Previously, we stated that the Courant number tells us how fast the information travels across the numerical cell. If the CFL condition is set above unity it means that the information propagates in such a way that within a single time step we can jump across a few cells – therefore we are losing the information and can create the numerical errors, which might lead to instabilities.

The choice of the Courant number is therefore very important as it leads directly to the size of the time step. The time step size has two main implications: the first one is the stability of the time discretization scheme, the second one is related to the physical meaning of the time step.

From the stability point of view, it should be noted that time integration schemes can be divided into explicit, implicit and semi-implicit schemes. Explicit schemes find the solution at the next time step using the known values, which can be written as \(y(t_{n+1}) = y(t_n) + \Delta t \cdot y'(t_n)\). Implicit schemes calculate solution at the next time step using unknown values \(y(t_{n+1}) = y(t_n) + \Delta t \cdot y'(t_{n+1})\). Semi-implicit schemes use both known and unknown values from next and current time steps.

Explicit schemes require such a time step that results in the Courant number smaller than 1. We might also find some recommendations to keep the Courant number below 0.7. For many applications, such a stringent condition leads to long and costly CFD simulations. To overcome this difficulty we can work with semi-implicit and implicit time integration schemes which allow the Courant number to be above 1. However, we already know that the Courant number above 1 means that the information can jump across few cells within a single time step, and we are going to lose some information. Therefore, even with the implicit time integration schemes we still should be careful about the time step. At this point we should also mention the quality of the numerical mesh as it is strictly related to the Courant number. Very small cells or cells with high skewness can create unnecessary stringent time step requirements.

Stability of the time integration scheme is one thing, but from the physical point of view time step defines the accuracy of numerical solution. If a certain flow structure varies faster than the time step, we will not be able to represent it correctly in the simulation. The characteristic time scale in a given problem to be solved depends on geometry size, boundary conditions and flow characteristics (for example laminar vs. turbulent flows).

A search for the representative physical time step can be challenging. However, some problems have an intrinsic time varying quantity that might be very helpful. For example, let us consider the rotation of a propeller with rotation rate of \(\Omega = 600\)rpm. The blade moves by 1 degree in about 2.7e-04s, which could be used as a reference time for selecting the time step (subdivision of the reference time should be applied).

For the propeller, alternatively we can calculate the linear velocity of the tip to get the distance covered within the reference time. If we assume that the radius of the propeller is 0.1m, then within the reference time the tip moves about 1.7mm (for \(\Omega = 600\)rpm). We need to take into consideration the mesh size around the tip of the propeller and assess the CFL condition to adjust the time step accordingly.

Another example of the intrinsic time reference is the flow around the airfoil. Using the velocity from the inlet boundary condition and characteristic length (chord or thickness of the airfoil) we can calculate the reference time scale. Applying the subdivision of the reference time scale will give us the proper time step.

In general, it is recommended to keep the Courant number below 1 in order to assure the stability of the numerical schemes and accuracy of the solution. Some specific problems might require even more stringent conditions, for example Large Eddy Simulations (LES) may require Co around 0.5. Very often in can be found in the literature to keep the Courant number equal or smaller than 0.7.

On the other hand, the Reader might find that for some applications, like turbo-machinery, the Courant number might be significantly larger than one – Co >> 1. It is a specific case when transient simulation is run to obtain a steady-state solution. If the time evolution is of no interest a pseudo-transient approach can be used to solve in time with large time step to drive the solution to convergence.

The PIMPLE algorithm in OpenFOAM is often used for transient problems with large time step (Co >> 1). It is the combination of PISO and SIMPLE algorithms. PISO and PIMPLE are used for transient problems, whereas SIMPLE is used for steady-state problems.

The principle of the PIMPLE algorithm is the following: within a single time step, we look for the steady-state solution while an under-relaxation strategy is applied. When the steady-state solution is found we proceed to the next time step.

The Courant number can be displayed as a field. Every cell in the mesh has a certain velocity, size and shape, therefore also a Courant number, which can be plotted as any other variable – Figure 6.

Plotting the Courant number field helps to understand where we can expect the highest velocity or where we have very small elements that might impose unnecessary stringent time step size. Also, in case the model is unstable, such a field can help to identify the critical regions of the flow and introduce proper modifications to our mesh or model set-up.

The Courant number will be larger if the velocity in the given cell is bigger or the size of the cell smaller. It should be also noted that the solver will use the maximum Courant number for time step calculation.