1.2.1.3. 1D Second-order Linear Diffusion - The Heat Equation

1.2.1.3.1. Understand the Problem

  • What is the final temperature profile for 1D diffusion when the initial conditions are a square wave and the boundary conditions are constant?

  • 1D diffusion is described as follows:

\[{\partial u \over \partial t} = \nu {\partial^2 u \over \partial x^2}\]

  • Consider looking for a solution of type:

\[u = \hat u e^{i(kx-\omega t)}\]

  • Represents a wave of amplitude \(\hat u\), \(\omega = 2 \pi f\)

  • Introducing into PDE, we obtain:

\[i \omega = \nu k^2\]

  • Leading to a solution:

\[u = \hat u e^{ikx}e^{- \nu k^2 t}\]

where \(e^{\nu k^2 t}\) is the exponential damping term. So diffusion is an exponentially damped wave.

Note: \(\nu > 0\) for physical diffusion (if \(\nu < 0\) would represent an exponentially growing phenomenon, e.g. an explosion or ‘the rich get richer’ model)

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.

1.2.1.3.2. Formulate the Problem

  • Same as Linear Convection

1.2.1.3.3. Design Algorithm to Solve Problem

1.2.1.3.3.1. Space-time discretisation

  • FD in time

  • CD in space

1.2.1.3.3.2. Numerical scheme

  • i \(\rightarrow\) index of grid in x

  • n \(\rightarrow\) index of grid in t

1.2.1.3.3.3. Discrete equation

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

1.2.1.3.3.4. Transpose

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

1.2.1.3.3.5. Pseudo-code

  • Very similar to Linear Convection

1.2.1.3.5. Conclusions

1.2.1.3.5.1. 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.

1.2.1.3.5.2. 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 {\Delta t \over (\Delta x)^2}\]