1. 1D First-order Linear Convection - The Wave Equation¶

Understand the Problem
Formulate the Problem
- Input Data
- Output Data
Design Algorithm to Solve Problem
Implement Algorithm in Python
Conclusions

1.1. Understand the Problem ¶

What is the final velocity profile for 1D linear convection when the initial conditions are a square wave and the boundary conditions are constant?
1D linear convection is described as follows:

\[{\partial u \over \partial t} + c {\partial u \over \partial x} = 0\]

1.2. Formulate the Problem ¶

1.2.1. Input Data ¶

nt = 51 (number of temporal points)
nx = 21 (number of spatial points)
tmax = 0.5
xmax = 2
c = 1

x	i	t	n	u(x,t)
\(0\)	\(0\)	\(0 \le t \le 0.5\)	\(0 \le i \le 20\)	\(1\)
\(0 < x \le 0.5\)	\(0 < n \le 12.5\)	\(0\)	\(0\)	\(1\)
\(0.5 < x \le 1\)	\(12.5 < n \le 25\)	\(0\)	\(0\)	\(2\)
\(1 < x < 2\)	\(25 < n < 50\)	\(0\)	\(0\)	\(1\)
\(2\)	\(50\)	\(0 \le t \le 0.5\)	\(0 \le i \le 20\)	\(1\)

1.2.2. Output Data ¶

x	t	u(x,t)
\(0 < x < 2\)	\(0.5\)	\(?\)

1.3. Design Algorithm to Solve Problem ¶

1.3.1. Space-time discretisation ¶

i \(\rightarrow\) index of grid in x
n \(\rightarrow\) index of grid in t

1.3.2. Numerical scheme ¶

FD in time
BD in space

1.3.3. Discrete equation ¶

\[{{u_i^{n+1} - u_i^n} \over {\Delta t}} + c {{u_i^n - u_{i-1}^n} \over \Delta x}=0\]

1.3.4. Transpose ¶

\[u_i^{n+1} = u_i^n - c{\Delta t \over \Delta x}(u_i^n - u_{i-1}^n)\]

1.3.5. Pseudo-code ¶

#Constants
nt = 51
tmax = 0.5
dt =  tmax/(nt-1)
nx =  21
xmax = 2
dx = xmax/(nx-1)

#Boundary Conditions
for n between 0 and 20
   u(0,n)=1
   u(50,n)=1

#Initial Conditions
for i between 1 and 49
   if(12.5 < i < 25)
       u(i,0) = 2
   else
       u(i,0) = 1

#Iteration
for n between 1 and 20
   u(i,n+1) = u(i,n)-c*(dt/dx)*(u(i,n)-u(i-1,n)

1.4. Implement Algorithm in Python ¶

def convection(nt, nx, tmax, xmax, c):
   """
   Returns the velocity field and distance for 1D linear convection
   """
   # Increments
   dt = tmax/(nt-1)
   dx = xmax/(nx-1)

   # Initialise data structures
   import numpy as np
   u = np.zeros((nx,nt))
   x = np.zeros(nx)

   # Boundary conditions
   u[0,:] = u[nx-1,:] = 1

   # Initial conditions
   for i in range(1,nx-1):
      if(i > (nx-1)/4 and i < (nx-1)/2):
         u[i,0] = 2
      else:
         u[i,0] = 1

   # Loop
   for n in range(0,nt-1):
      for i in range(1,nx-1):
         u[i,n+1] = u[i,n]-c*(dt/dx)*(u[i,n]-u[i-1,n])

   # X Loop
   for i in range(0,nx):
      x[i] = i*dx

   return u, x

def plot_convection(u,x,nt,title):
   """
   Plots the 1D velocity field
   """

   import matplotlib.pyplot as plt
   plt.figure()
   for i in range(0,nt,10):
      plt.plot(x,u[:,i],'r')
      plt.xlabel('x (m)')
      plt.ylabel('u (m/s)')
      plt.ylim([0,2.2])
      plt.title(title)
      plt.show()

(Source code)

u,x = convection(151, 51, 0.5, 2.0, 0.5)
plot_convection(u,x,151,'Figure 1: c=0.5m/s, nt=151, nx=51, tmax=0.5s')

u,x = convection(151, 1001, 0.5, 2.0, 0.5)
plot_convection(u,x,151,'Figure 2: c=0.5m/s, nt=151, nx=1001, tmax=0.5s')

u,x = convection(151, 51, 2.0, 2.0, 0.5)
plot_convection(u,x,151,'Figure 3: c=0.5m/s, nt=151, nx=51, tmax=2s')

(png, hires.png, pdf)

(png, hires.png, pdf)

(png, hires.png, pdf)

1.5. Conclusions ¶

1.5.1. Why isn’t the square wave maintained?¶

The first order backward differencing scheme in space creates false diffusion.
If the spatial step is reduced, the error reduces - compare Figure 1 and Figure 2

1.5.2. Why does the wave shift to the right?¶

The square wave is being convected by the constant linear wave speed c
For \(c > 0\) profiles shift to the right by \(c \Delta t\) - see Figure 2

1.5.3. What happens at the wall?¶

As there is no viscosity, there is a non-physical change the the profile near the wall.