Implicit heat conduction

For constant diffusivity, $\alpha$, in space,

$$\frac{\partial T}{\partial t} = \alpha \left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}\right) + H$$

For spatially varying diffusivity,

$$\frac{\partial T}{\partial t} = \frac{\partial \left(\alpha \frac{\partial T}{\partial x} \right)}{\partial x} + \frac{\partial \left(\alpha \frac{\partial T}{\partial y} \right)}{\partial y} + H$$

We use finite differences to approximate the second order spatial derivatives. Since we are using a centred differencing scheme, then the diffusivities should be located at half space steps centred on their specific gradients

$$$$\frac{\partial \left(\alpha \frac{\partial T}{\partial x} \right)}{\partial x} = \frac{1}{\Delta x} \left( \frac{ \alpha_{i+1/2,j} (T_{i+1,j}-T_{i,j}) }{\Delta x} - \frac{ \alpha_{i-1/2,j} (T_{i,j}-T_{i-1,j}) }{\Delta x} \right)$$$$

next

$$\frac{\partial \left(\alpha \frac{\partial T}{\partial y} \right)}{\partial y} = \frac{1}{\Delta y} \left( \frac{ \alpha_{i,j+1/2} (T_{i,j+1}-T_{i,j}) }{\Delta y} - \frac{ \alpha_{i,j-1/2} (T_{i,j}-T_{i,j-1}) }{\Delta y} \right)$$

where $\alpha_{i+1/2,j}$ can be averaged by,

$$\alpha_{i+1/2,j} = \frac{\alpha_{i+1,j} + \alpha_{i,j}}{2}$$

Expanding out the difference and collecting like-terms ($$T_{i+1,j}, T_{i,j}, T_{i-1,j}$$) in one dimensions gives,

\begin{align} 2H \Delta x^2 \frac{\partial T}{\partial t} =& [\alpha_{i+1,j}T_{i+1,j} + \alpha_{i,j}T_{i+1,j}] \ & - [\alpha_{i+1,j},T_{i,j} - 2\alpha_{i,j}T_{i,j} - \alpha_{i-1,j}T_{i,j}] \ & + [\alpha_{i-1,j}T_{i-1,j} + \alpha_{i,j}T_{i-1,j}] \end{align}

Boundary conditions

Surface boundary condition a constant temperature is assigned $$T(0,t)=T0$$.

Left and right boundaries are assigned as insulated boundary conditions $$\frac{\partial T}{\partial x}(0,t)=0$$.

$$\frac{D_{N}+D_{N-1}}{2 \Delta x^2} - \frac{H_{N}+H_{N-1}}{2}$$

• Bottom boundary condition is a Neumann (flux) boundary condition $$\frac{\partial T}{\partial x}(0,t)=Q$$.

$$\frac{D_{N}+D_{N-1}}{2 \Delta x^2} - \frac{H_{N}+H_{N-1}}{2} + \frac{Q_{N}}{\Delta x}$$

Matrix notation

This can be expressed in matrix form to solve algebraically, $$Ax = b$$ where $$A$$ is the coefficient matrix. In 1D,

$$\begin{pmatrix} \alpha_{2,1}+\alpha_{1,1} & -\alpha_{2,1}-2\alpha_{1,1}-\alpha_{0,1} & \alpha_{0,1}+\alpha_{1,1} & 0 & 0 \\ 0 & \alpha_{3,1}+\alpha_{2,1} & -\alpha_{3,1}-2\alpha_{2,1}-\alpha_{1,1} & \alpha_{1,1}+\alpha_{2,1} & 0 \\ 0 & 0 & \alpha_{4,1}+\alpha_{3,1} & -\alpha_{4,1}-2\alpha_{3,1}-\alpha_{2,1} & \alpha_{2,1}+\alpha_{3,1} \\ \end{pmatrix} \begin{bmatrix} T_{1,1} \\ T_{2,1} \\ T_{3,1} \end{bmatrix} = \begin{bmatrix} H_{1,1} \\ H_{2,1} \\ H_{3,1} \end{bmatrix}$$
Dr. Ben Mather
Research Fellow in Geophysics

My research interests include Bayesian inversion, Python programming and Geodynamics.