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

$$ \begin{equation} \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) \end{equation} $$

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,