The 1-D Linear Convection equation is the simplest, most basic model that can be used to learn something about CFD.

\(\frac{{\partial u}}{{\partial t}} + c\frac{{\partial u}}{{\partial x}} = 0\)

With given initial conditions (understood as a wave), the equation represents the propagation of that initial wave with speed \(c\), without change of shape.

Let the initial condition be \(u(x,0)=u_0(x)\). Then the exact solution of the equation is \(u(x,t)=u_0(x-ct)\).

We discretize this equation, in both space and time, using the Forward Difference scheme for the time derivative and the Backward Difference scheme for the space derivative.

Our discrete equation is:

\(\frac{{u_i^{n + 1} – u_i^n}}{{\Delta t}} + c\frac{{u_i^n – u_{i – 1}^n}}{{\Delta x}} = 0\)

Now let’s try implementing this in Python.

We’ll start by importing a few libraries to help us out.

• numpy is a library that provides a bunch of useful matrix operations akin to MATLAB;
• matplotlib is a 2D plotting library that we will use to plot our results;
• time and sys provide basic timing functions that we’ll use to slow down animations for viewing.

# Remember: comments in python are denoted by the pound sign
import numpy as np                      #here we load numpy
from matplotlib import pyplot as plt      #here we load matplotlib
import time, sys                   #and load some utilities



nx = 41         # spatial discretization: number of grid points
dx = 2 / (nx-1) # the distance between any pair of adjacent grid points
nt = 25         # nt is the number of timesteps we want to calculate
dt = .025       # dt is the amount of time each timestep covers (delta t)
c = 1           # assume wavespeed of c = 1

We have define an evenly spaced grid of points within a spatial domain that is 2 units of length wide, i.e., \(x_i \in [0,2]\). We have define a variable nx, which represent the number of grid points we want and dx is the distance between any pair of adjacent grid points.

We also need to set up our initial conditions. The initial velocity \(u_0\) is given as \(u=2\) in the interval \(0.5 \leq x \leq 1\) and \(u=1\) everywhere else in \((0,2)\) (i.e., a hat function). Here, we use the function ones() defining a numpy array which is nx elements long with every value equal to 1.

u = np.ones(nx)                      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s

x = np.linspace(0,2,nx)
plt.plot(x,u, 'red', label='Initial shape')

plt.show()

Now it’s time to implement the discretization of the convection equation using a finite-difference scheme.

We’ll store the result in a new (temporary) array un, which will be the solution \(u\) for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

We first initialize our placeholder array un to hold the values we calculate for the \(n+1\) timestep, using once again the NumPy function ones().

Then, we may think we have two iterative operations: one in space and one in time (we’ll learn differently later), so we’ll start by nesting one loop inside the other. Note the use of the nifty range() function. When we write: for i in range(1,nx) we will iterate through the u array, but we’ll be skipping the first element (the zero-th element).

un = np.ones(nx) #initialize a temporary array

for n in range(nt):        #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy()          ##copy the existing values of u into un
    for i in range(1, nx): ## you can try commenting this line and...
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
    

plt.plot(x,u, 'green', label='Final shape')

plt.xlabel('space')
plt.ylabel('value')
pplt.title('1-D Linear Convection : nx=41, nt=25, dt=.025 and c=1')


plt.legend()

plt.show()

OK! So our hat function has definitely moved to the right, but it’s no longer a hat. What’s going on?