1D Linear Diffusion

The Heat Equation

  • Understand the Problem

The physics of diffusion are:

  • An expotentially damped wave in time
  • Isotropic in space - the same in all spatial directions - it does not distinguish between upstream and downstream
  • The phenomenon of diffusion is isotropic - so the finite difference formula that represents that physics is central differencing CD, because CD takes values from upstream and downstream equally.

Design Algorithm

---write equations---

Implement Algorithm in Python

In [8]:
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm
%matplotlib inline
In [9]:
def diffusion(nt, nx, tmax, xmax, nu):
    dt = tmax/(nt-1)
    dx = xmax/(nx-1)
    
    # Initilisation of data structures
    u = np.zeros((nx, nt))
    x = np.zeros(nx)
    
    # BCs
    u[0, :] = u[nx-1, :] = 1
    
    # ICs
    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(0,nx-1):
                u[i,n+1] = u[i,n] + nu*(dt/dx**2.0)*(u[i+1,n]-2.0*u[i,n]+u[i-1,n])
           
    # X Loop
    for i in range(0,nx):
           x[i] = i*dx
           
    return u, x

def plot(u, x, nt, title):
    """
    Plots the Velocity Field Results in 1D
    """
    plt.figure()
    plt.show()
    colour=iter(cm.rainbow(np.linspace(0,10,nt)))
    for i in range(0, nt, 10):
        c=next(colour)
        plt.plot(x, u[:,i], c=c)
        plt.xlabel('x [m]')
        plt.ylabel('u [m/s]')
        plt.ylim([0, 2.2])
        plt.title(title)
        plt.grid(True, linestyle='-.', linewidth='0.5', color='black')
In [11]:
u,x = diffusion(151, 51, 0.5, 2.0, 0.1)
plot(u,x,151,'Figure 1: nu=0.1, nt=151, nx=51, tmax=0.5s')

<Figure size 432x288 with 0 Axes>
In [13]:
u,x = diffusion(151, 51, 0.5, 2.0, 0.242)
plot(u,x,151,'Figure 1b: nu=0.242, nt=151, nx=51, tmax=0.5s')

<Figure size 432x288 with 0 Axes>
In [14]:
u,x = diffusion(151, 79, 0.5, 2.0, 0.1)
plot(u,x,151,'Figure 2: nu=0.1, nt=151, nx=79, tmax=0.5s')

<Figure size 432x288 with 0 Axes>
In [15]:
u,x = diffusion(151, 51, 1.217, 2.0, 0.1)
plot(u,x,151,'Figure 3: nu=0.1, nt=151, nx=51, tmax=1.217s')

<Figure size 432x288 with 0 Axes>

Conclusions

  • Why isn’t the square wave maintained?

    • The square wave isn’t maintained because the system is attempting to reach equilibrium - the rate of change of velocity being equal to the shear force per unit mass. There are no external forces and no convective acceleration terms.

  • Why does increasing the viscosity, spatial points and time period cause instability?

    • If the viscosity is too large, or if the number of spatial points is too large or if the timestep is too large, then the central differencing method becomes unstable. This is due to the ratio, r. If r is too large, the method becomes unstable: $$r = \nu\frac{\Delta{t}}{\Delta{x^2}}$$
In [ ]: