2.

22 Computational Fluid Dynamics

2.4.2 First-Order Wave Equation

Before we proceed further, let us look at the system of first-order equations
which are frequently encountered in a class of fluid flow problems. Consider the
second-order wave equation
∂2u 2
2∂ u
= c (2.55)
∂t2 ∂x2
Here c is the wave speed and u is the wave amplitude. This can be written as
a system of two first-order equations. If v = ∂u/∂t and w = c(∂u/∂x), then we
may write ∂v/∂t = c(∂w/∂x) and ∂w/∂t = c(∂v/∂x).
Rather, the system of equations may be written as
∂U ∂U
+ [A] =0
∂t ∂x
which is a first-order equation.
It is implicit that U = { wv } and A = −c
 0 −c 
0 . The eigenvalues λ of the [A]
matrix are found by
det [A − λI] = 0, or λ2 − c2 = 0
Roots of the characteristic equation are λ1 = +c and λ2 = −c, representing two
travelling waves with speeds given by
   
dx dx
= c and = −c
dt 1 dt 2
The system of equations in this example is hyperbolic and it has also been seen
that the eigenvalues of the A matrix represent the characteristic differential
representation of the wave equation. Euler’s equation may be treated as a
system of first-order wave equations. For Euler’s equations, in two dimensions,
we can write a system of first order as
∂E ∂E ∂E
+ [A] + [B] = [S] (2.56)
∂t ∂x ∂y
where
   
u 0 u
E= , A= ,
v u 0
1 ∂p
 
  −
v 0
B= and S =  ρ1 ∂x
 
0 v ∂p 
−
ρ ∂y

2.4.3 Stability of Hyperbolic and Elliptic Equations

Let us examine the characteristics of the first-order wave equation given by
∂u ∂u
+c =0 (2.57)
∂t ∂x
 Finite Difference Method 2.23

Here we shall represent the spatial derivative by the central difference form

∂u un − uni−1
= i+1 (2.58)
∂x 2∆x
We shall replace the time derivative with a first-order difference , where u(t) is
represented by an average value between grid points (i + 1) and (i − 1), i.e
1 n
u(t) = (u + uni−1 )
2 i+1
Then
∂u un+1 − 21 (uni+1 + uni−1 )
= i (2.59)
∂t ∆t
Substituting Eqns. (2.58) and (2.59) into (2.57), we have

uni+1 + uni−1 ∆t (uni+1 − uni−1 )

un+1
i = −c (2.60)
2 ∆x 2
The time derivative is called Lax method of discretization, after the well known
mathematician Peter Lax who first proposed it. If we once again assume an
error of the form
ǫm (x, t) = eat eIkm x (2.61)
as done previously, and substitute this form into Eq. (2.60), following the same
arguments as applied to the analysis of Eq. (2.26), the amplification factor
becomes
ea∆t = cos(km ∆x) − IC sin(km ∆x)
where C = c(∆t/∆x). The stability requirement is |ea∆t | ≤ 1. Finally the
condition culminates in
∆t
C=c ≤1 (2.62)
∆x
In Eq. (2.62), C is the Courant number. This equation restricts ∆t ≤ ∆x/c for
the solution of Eq. (2.62) to be stable. The condition posed by Eq. (2.62) is
called the Courant-Friedrichs-Lewy condition, generally referred to as the CFL
condition.

Physical Example of Unstable Calculation

Let us take the heat conduction once again,

∂u ∂2u
=α 2 (2.63)
∂t ∂x
Applying FTCS discretization scheme depicts simple explicit representation
as
un+1
 n
− uni ui+1 − 2uni + uni−1

i
=α (2.64)
∆t (∆x2 )
2.24 Computational Fluid Dynamics

or
un+1
i = r (uni+1 + uni−1 ) + (1 − 2r) uni , where r = α∆t/(∆x2 ) (2.65)
This is stable only if r ≤ 1/2.
Let us consider a case when r > 1/2. For r = 1 (which is greater than the
stability restriction), we get un+1
i = 1 · (100 + 100) + (1 − 2) · 0 = 200o C, (which
is impossible). The values of u are shown in Fig. 2.6.

Next, an example demonstrating the application of Von Neumann method

0
200 C
n+1 t+ ∆ t

0 0 0
100 C 0 C 100 C
n
i−1 i i+1

Figure 2.6: Physical violations resulting from r=1.

to multidimensional elliptic problems is taken up. Let us take the vorticity

transport equation:
 2
∂ ω ∂2ω

∂ω ∂ω ∂ω
+u +v =ν + (2.66)
∂t ∂x ∂y ∂x2 ∂y 2
We shall extend the Von Neumann stability analysis for this equation, assum-
ing u and v as constant coefficients (within the framework of linear stability
analysis). Using FTCS scheme
n+1 n  n n  n n
ωi,j − ωi,j ωi+1,j − ωi−1,j ωi,j+1 − ωi,j−1
 
=−u −v
∆t 2∆x 2∆y
 n n n
ωi+1,j − 2ωi,j + ωi−1,j

+ν
(∆x2 )
 n n n
ωi,j+1 − 2ωi,j + ωi,j−1

+ν (2.67)
(∆y 2 )
Let us consider N = D + ǫ with
X
ǫ(x, y, t) = eat e(Ikm x+Ikm y) (2.68)
m

where N is the numerical solution obtained from computer, D the exact solution
of the FDE and ǫ error. Substituting Eq. (2.68) into Eq. (2.67) and using the
trignometric identities, we finally obtain
ǫn+1
i,j ea(t+∆t) eIkm (x+y)
n = = ea∆t = G
ǫi,j eat eIkm (x+y)
 Finite Difference Method 2.25

where

G =1 − 2(dx + dy ) + 2dx cos(km ∆x) + 2dy cos(km ∆y)

− I[Cx sin(km ∆x) + Cy sin(km ∆y)]

where
ν∆t ν∆t u∆t v∆t
dx = 2, dy = 2, Cx = , Cy =
(∆x) (∆y) ∆x ∆y
The obvious stability condition |G| ≤ 1, finally leads to

1
dx + dy ≤ , Cx + Cy ≤ 1 (2.69)
2
when
1
dx = dy = d (for ∆x = ∆y), d≤
4
which means
ν∆t 1
2 ≤
(∆x) 4
This is twice as restrictive as the one-dimensional diffusive limitation (compare
with Eq. (2.54). Again for the special case (u = v and ∆x = ∆y)

1
Cx = Cy = C, hence C ≤
2
which is also twice as restrictive as one dimensional convective limitation (com-
pare with Eq. (2.62).
Finally, let us look at the stability requirements for the second-order wave
equation given by
∂2u ∂2u
2
= c2 2
∂t ∂x
We replace both the spatial and time derivatives with central difference scheme
(which is second-order accurate)
" #
un+1
i − 2uni + uin−1 n n n
2 ui+1 − 2ui + ui−1
2 =c 2 (2.70)
(∆t) (∆x)

Again assume
N =D+ǫ (2.71)
and
ǫni = eat eIkm x (2.72)
Substituting Eq. (2.72) and (2.71) in (2.70) and dividing both sides by eat eIkm x ,
we get
ea∆t − 2 + e−a∆t = C 2 eIkm ∆x + e−Ikm ∆x − 2
 
(2.73)
2.26 Computational Fluid Dynamics

where
c(∆t)
C, the Courant number = (2.74)
∆x
From Eq. (2.73), using trignometric identities, we get
 
km ∆x
ea∆t + e−a∆t = 2 − 4C 2 sin2 (2.75)
2

and, the amplification factor

ǫn+1
G= i
= |ea∆t | (2.76)
ǫni

However, from Eq. (2.75) we arrive at

  
km ∆x
e2a∆t − 2 1 − 2C 2 sin2 ea∆t + 1 = 0 (2.77)
2

which is a quadratic equation for ea∆t . This equation, quite obviously, has two
roots, and the product of the roots is equal to +1. Thus, it follows that the
magnitude of one of the roots (value of ea∆t ) must exceed 1 unless both the
roots are equal to unity.
But ea∆t is the magnification factor. If its value exceeds 1, the error will grow
exponentially which will lead to an unstable situation. All these possibilities
mean that Eq (2.77) should possess complex roots in order that both have the
values of ea∆t equal to unity. This implies that the discriminant of Eq. (2.77)
should be negative.
  2
km ∆x
1 − 2C 2 sin2 −1<0 (2.78)
2
or
1
C2 < (2.79)
sin2 km ∆x


2

which is always true if C < 1. Hence CFL condition (C < 1), must again be
satisfied for the stability of second-order hyperbolic equations.
In light of the above discussion, we can say that a finite-difference procedure
will be unstable if for that procedure, the solution becomes unbounded, i.e the
error grows exponentially as the calculation progresses in the marching direction.
In order to have a stable calculation, we pose different conditions based on
stability analysis. Here we have discussed the Von Neumann stability analysis
which is indeed a linear stability analysis.
However, situations may arise where the amplification factor is always less
than unity. These conditions are referred to as unconditionally stable. In a
similar way for some procedures, we may get an amplification factor which is
always greater than unity. Such methods are unconditionally unstable.
 Finite Difference Method 2.27

Over and above, it should be realized that such stability analysis are not
really adequate for practical complex problems. In actual fluid flow problems,
the stability restrictions are applied locally. The mesh is scanned for the most
restrictive value of the stability limitations and the resulting minimum ∆t is used
throughout the mesh. For variable coefficients, the Von Neumann condition is
only necessary but not sufficient. As such, stability criterion of a procedure is not
defined by its universal applicability. For nonlinear problems we need numerical
experimentation in order to obtain stable solutions wherein the routine stability
analysis will provide the initial clues to practical stability. In other words, it
will give tutorial guidance only.

2.5 Fundamentals of Fluid Flow Modeling

We have discussed the finite-difference methods with respect to the solution of
linear problems such as heat conduction. The problems of fluid mechanics are
more complex in character. The governing partial differential equations form
a nonlinear system which must be solved for the unknown pressures, densities,
temperatures and velocities.
Before entering into the domain of actual flow modeling, we shall discuss
some subtle points of fluid flow equations with the help of a model equation.
The model equation should have convective, diffusive and time-dependent terms.
Burgers (1948) introduced a simple nonlinear equation which meets the aforesaid
requirements (Burger’s equation).

∂ζ ∂ζ ∂ 2ζ
+u =ν (2.80)
∂t ∂x ∂x2
Here, u is the velocity, ν is the coefficient of diffusivity and ζ is any property
which can be transported and diffused. If the viscous term (diffusive term) on
the right-hand side is neglected, the remaining equation may be viewed as a
simple analog of Euler’s equation.

∂ζ ∂ζ
+u =0 (2.81)
∂t ∂x
Now we shall see the behavior of Burger’s equations for different kinds of dis-
cretization methods. In particular, we shall study their influence on conservative
and transportive property, and artificial viscosity.

2.5.1 Conservative Property

A finite-difference equation possesses conservative property if it preserves in-
tegral conservation relations of the continuum. Let us consider the vorticity
transport equation
∂ω
= −(V · ∇) ω + ν ∇2 ω (2.82)
∂t
 Finite Difference Method 2.27

Over and above, it should be realized that such stability analysis are not
really adequate for practical complex problems. In actual fluid flow problems,
the stability restrictions are applied locally. The mesh is scanned for the most
restrictive value of the stability limitations and the resulting minimum ∆t is used
throughout the mesh. For variable coefficients, the Von Neumann condition is
only necessary but not sufficient. As such, stability criterion of a procedure is not
defined by its universal applicability. For nonlinear problems we need numerical
experimentation in order to obtain stable solutions wherein the routine stability
analysis will provide the initial clues to practical stability. In other words, it
will give tutorial guidance only.

2.5 Fundamentals of Fluid Flow Modeling

We have discussed the finite-difference methods with respect to the solution of
linear problems such as heat conduction. The problems of fluid mechanics are
more complex in character. The governing partial differential equations form
a nonlinear system which must be solved for the unknown pressures, densities,
temperatures and velocities.
Before entering into the domain of actual flow modeling, we shall discuss
some subtle points of fluid flow equations with the help of a model equation.
The model equation should have convective, diffusive and time-dependent terms.
Burgers (1948) introduced a simple nonlinear equation which meets the aforesaid
requirements (Burger’s equation).

∂ζ ∂ζ ∂ 2ζ
+u =ν (2.80)
∂t ∂x ∂x2
Here, u is the velocity, ν is the coefficient of diffusivity and ζ is any property
which can be transported and diffused. If the viscous term (diffusive term) on
the right-hand side is neglected, the remaining equation may be viewed as a
simple analog of Euler’s equation.

∂ζ ∂ζ
+u =0 (2.81)
∂t ∂x
Now we shall see the behavior of Burger’s equations for different kinds of dis-
cretization methods. In particular, we shall study their influence on conservative
and transportive property, and artificial viscosity.

2.5.1 Conservative Property

A finite-difference equation possesses conservative property if it preserves in-
tegral conservation relations of the continuum. Let us consider the vorticity
transport equation
∂ω
= −(V · ∇) ω + ν ∇2 ω (2.82)
∂t
2.28 Computational Fluid Dynamics

where ∇ is nabla or differential operator, V the fluid velocity and ω the

vorticity. If we integrate this over some fixed space region ℜ, we get

∂ω
Z Z Z
dℜ = − (V · ∇)ωdℜ + ν∇2 ωdℜ (2.83)
ℜ ∂t ℜ ℜ

The first term of the Eq. (2.83) can be written as

∂ω ∂
Z Z
dℜ = ω dℜ
ℜ ∂t ∂t ℜ

The second term of the Eq.( (2.83)) may be expressed as

Z Z Z
− (V · ∇)ω dℜ = − ∇ · (V ω) dℜ = − (V ω) · n dA
ℜ ℜ Ao

Ao is the boundary of ℜ, n is unit normal vector and dA is the differential

element of Ao . The remaining term of Eq. (2.83) may be written as
Z Z
ν ∇2 ω dℜ = ν (∇ω) · n dA
ℜ Ao

As because,
Z Z Z
ν (∇ω) · n dA = ν ∇ · (∇ω) dℜ = ν ∇2 ω dℜ
Ao ℜ ℜ

Finally, we can write

∂
Z Z Z
ω dℜ = − (V ω) · n dA + ν (∇ ω) · n dA (2.84)
∂t ℜ Ao Ao

which implies that the time rate of accumulation of ω in ℜ is equal to net

advective flux rate of ω across Ao into ℜ plus net diffusive flux rate of ω across
Ao into ℜ. The concept of conservative property is to maintain this integral
relation in finite difference representation.
For clarity , again let us consider inviscid Burger’s equation ((2.81)). This time
we let ζ = ω = vorticity, which means

∂ω ∂
= − (uω) (2.85)
∂t ∂x

The finite difference analog is given by FTCS method as

ωin+1 − ωin uni+1 ωi+1

n
− uni−1 ωi−1
n
=− (2.86)
∆t 2∆x
 Finite Difference Method 2.29

Let us consider a region ℜ running from i = I1 to i = I2 see (Figure 2.7). We

I2
1 X
evaluate the the integral ω ∆x as
∆t
i=I1
" I i=I
# I2
1 X 2
n+1
X2
n
X (uni+1 ωi+1
n
) − (uni−1 ωi−1
n
)
ωi ∆x − ωi ∆x = −
∆t 2
i=I1 i=I1 i=I1
i=I
1 X2
= [(u ω)ni−1 − (u ω)ni+1 ] (2.87)
2
i=I1

Summation of the right hand side (running i from I1 to I2 ) finally gives

" I I2
#
2
1 X n+1
X
n 1
(u ω)nI1 −1 + (u ω)nI1

ωi ∆x − ωi ∆x =
∆t 2
i=I1 i=I1

1
(u ω)nI2 + (u ω)nI2 +1

−
2
= (u ω)nI1 − 1 − (u ω)nI2 + 1 (2.88)
2 2

Eq. (2.88) states that the rate of accumulation of ωi in ℜ is identically equal

to the net advective flux rate across the boundary of ℜ running from i = I1 to
i = I2 . Thus the FDE analogous to inviscid part of the integral Eq. (2.86) has
preserved the conservative property. As such, conservative property depends on
the form of the continuum equation used. Let us take non-conservative form of
inviscid Burger’s equation (2.81) as
∂ω ∂ω
= −u (2.89)
∂t ∂x
Using FTCS differencing technique as before, we can write

ωin+1 − ωin
 n n 
ωi+1 − ωi−1
= −uni (2.90)
∆t 2∆x

Now, the integration over ℜ running from i = I1 to i = I2 , yields

" I i=I
# I2
1 X 2 X2 (ω n − ωi−1 n
)
−uni i+1
X
n+1 n
ωi ∆x − ωi ∆x =
∆t 2
i=I1 i=I1 i=I1
I2
1X
= [uni ωi−1
n
− uni ωi+1
n
] (2.91)
2
i=I1

While performing the summation of the right-hand side of Eq. (2.91), it can
be observed that terms corresponding to inner cell fluxes do not cancel out.
2.30 Computational Fluid Dynamics

Consequently an expression in terms of fluxes at the inlet and outlet section, as

it was found earlier, could not be obtained. Hence the finite-difference analog
Eq. (2.90) has failed to preserve the integral Gauss-divergence property, i,e., the
conservative property of the continuum.
The quality of preserving the conservative property is of special importance

I1 I2

Figure 2.7: Domine running from i = I1 to i = I2 .

with regards to the methods involving finite-volume approach(a special form

of finite-difference equation). The use of conservative form depicts that the
1
advective flux rate of ω out of a control volume at the interface i = I2 + is
2
exactly equal to flux rate of ω in to the next control volume and so on. The
meaning of calling Eq. (2.85)as “conservative form” is now clearly understood.
However, the conservative form of advective part is of prime importance for
modeling fluid flow and is often referred to as weak conservative form. For the
incompressible flow in Cartesian coordinate system this form is :

∂u ∂u2 ∂uv ∂uw 1 ∂p

+ + + = − + ν∇2 u
∂t ∂x ∂y ∂z ρ ∂x
∂v ∂uv ∂v 2 ∂vw 1 ∂p
+ + + = − + ν∇2 v
∂t ∂x ∂y ∂z ρ ∂y
∂w ∂uw ∂vw ∂w2 1 ∂p
+ + + = − + ν∇2 w (2.92)
∂t ∂x ∂y ∂z ρ ∂z

If all the terms in the flow equation are recast in the form of first- order deriva-
tives of x, y, z and t , the equations are said to be in strong“ conservative
form”. We shall write the strong conservation form of Navier-Stokes equation
in Cartesian coordinate system:

∂u ∂ 2 p ∂u ∂ ∂u ∂ ∂u
+ (u + − ν ) + (uv − ν ) + (uw − ν ) =0
∂t ∂x ρ ∂x ∂y ∂y ∂z ∂z
∂v ∂ ∂u ∂ 2 p ∂v ∂ ∂v
+ (uv − ν ) + (v + − ν ) + (vw − ν ) =0
∂t ∂x ∂x ∂y ρ ∂y ∂z ∂z
∂w ∂ ∂w ∂ ∂w ∂ p ∂w
+ (uw − ν )+ (wv − ν )+ (w2 + − ν ) = 0 (2.93)
∂t ∂x ∂x ∂y ∂y ∂z ρ ∂z
 Finite Difference Method 2.31

2.5.2 The Upwind Scheme

Once again,we shall start with the inviscid Burger’s equation. (2.81). Regarding
discretization, we can think about the following formulations
ζin+1 − ζin ζ n − ζin
+ u i+1 =0 (2.94)
∆t ∆x
ζin+1 − ζin ζ n − ζi−1
n
+ u i+1 =0 (2.95)
∆t 2 ∆x
If Von Neumann’s stability analysis is applied to these schemes, we find that
both are unconditionally unstable.
A well known remedy for the difficulties encountered in such formulations is
the upwind scheme which is described by Gentry, Martin and Daly (1966) and
Runchal and Wolfshtein (1969). Eq. (2.94) can be made stable by substituting
the forward space difference by a backward space difference scheme, provided
that the carrier velocity u is positive. If u is negative, a forward difference
scheme must be used to assure stability. For full Burger’s equation. (2.80), the
formulation of the diffusion term remains unchanged and only the convective
term (in conservative form) is calculated in the following way (Figure 2.8):
ζin+1 − ζin u ζ n − u ζi−1
n
=− i + viscous term, for u > 0 (2.96)
∆t ∆x
ζin+1 − ζin u ζ n − u ζin
= − i+1 + viscous term, for u < 0 (2.97)
∆t ∆x
It is also well known that upwind method of discretization is very much necessary
in convection (advection) dominated flows in order to obtain numerically stable
results. As such, upwind bias retains transprotative property of flow equation.
Let us have a closer look at the transportative property and related upwind
bias.

∆x ∆x
u
u
i−2 i−1 i i+1 i+2

Figure 2.8: The Upwind Scheme.

2.5.3 Transportive Property

A finite-difference formulation of a flow equation possesses the transportive
property if the effect of a perturbation is convected (advected) only in the di-
rection of the velocity.
2.32 Computational Fluid Dynamics

Consider the model Burger’s equation in conservative form

∂ζ ∂u ζ
=− (2.98)
∂t ∂x
Let us examine a method which is central in space. Using FTCS we get
ζin+1 − ζin u ζ n − u ζi−1
n
= − i+1 (2.99)
∆t 2∆x
Consider a perturbation ǫm = δ in ζ. A perturbation will spread in all directions
due to diffusion. We are taking an inviscid model equation and we want the
perturbation to be carried along only in the direction of the velocity. So, for
u > 0, ǫm = δ (perturbation at mth space location), all other ǫ = 0. Therefore,
at a point (m + 1) downstream of the perturbation
n+1 n
ζm+1 − ζm+1 0−uδ uδ
=− =+
∆t 2∆x 2∆x
which is acceptable. However, at the point of perturbation (i = m),
n+1 n
ζm − ζm 0−0
=− =0
∆t 2∆x
which is not very reasonable. But at the upstream station (i = m − 1) we
observe
n+1 n
ζm−1 − ζm−1 uδ −0 uδ
=− =−
∆t 2∆x 2∆x
which indicates that the transportive property is violated.
On the contrary, let us see what happens when an upwind scheme is used.
We know that for u > 0
ζin+1 − ζin u ζ n − u ζi−1
n
=− i (2.100)
∆t ∆x
Then for ǫm = δ at the downstream location (m + 1)
n+1 n
ζm+1 − ζm+1 0−uδ uδ
=− =+
∆t ∆x ∆x
which follows the rationale for the transport property.
At point m of the disturbance
n+1 n
ζm − ζm uδ −0 uδ
=− =−
∆t ∆x ∆x
which means that the perturbation is being transported out of the affected
region.
Finally, at (m − 1) station, we observe that
n+1 n
ζm−1 − ζm−1 0−0
=− =0
∆t ∆x
 Finite Difference Method 2.33

This signifies that no perturbation effect is carried upstream. In other words, the
upwind method maintains unidirectional flow of information. In conclusion, it
can be said that while space centred differences are more accurate than upwind
differences, as indicated by the Taylor series expansion, the whole system is not
more accurate if the criteria for accuracy includes the transportive property as
well.

2.5.4 Upwind Differencing and Artificial Viscosity

Consider the model Burger’s equation. (2.80) and focus the attention on the
inertia terms
∂ζ ∂ζ ∂2ζ
+u =ν 2
∂t ∂x ∂x
As seen, the simple upwind scheme gives

ζin+1 − ζin u ζ n − u ζi−1

n
=− i + ··· for u > 0
∆t ∆x
ζin+1 − ζin u ζ n − u ζin
= − i+1 + ··· for u < 0
∆t ∆x
From Taylor series expansion, we can write
n 2 n
∂ζ (∆t) ∂ 2 ζ
ζin+1 = ζin + ∆t + + ··· (2.101)
∂t i 2 ∂t2 i
n 2 n
n ∂ζ (∆x) ∂ 2 ζ
ζi±1 = ζin ± ∆x + ± ··· (2.102)
∂x i 2 ∂x2 i

Substituting Eqns. (2.101) and (2.102) into (2.96) gives (dropping the subscript
i and superscript n)
" #
2
1 ∂ζ (∆t) ∂ 2 ζ
∆t + + O(∆t)3
∆t ∂t 2 ∂t2
" #
u ∂ζ (∆x)2 ∂ 2 ζ 3
=− ∆x − + O(∆x) + [Diffusive term]
∆x ∂x 2 ∂x2

or  2
∂2ζ
 
∂ζ ∂ζ 1 u ∆t ∂ ζ
= −u + u ∆x 1 − + ν + O(∆x)2
∂t ∂x 2 ∆x ∂x2 ∂x2
which may be rewritten as

∂ζ ∂ζ ∂2ζ ∂2ζ
= −u + ν 2 + νe 2 + higher-order terms (2.103)
∂t ∂x ∂x ∂x
where
1 u ∆t
νe = [u ∆x (1 − C)], C (Courant number) =
2 ∆x