Heat conduction in 2D

The heat equation in two dimensions (without advection) can be stated as

(1)\[\rho C_p \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) + H\]

The only addition compared to the one-dimensional case is the second order derivative with respect to \(y\) on the right hand side of the equation.

Exercise

1. Discretize equation (1). Use the implicit discretization, i.e. calculate temperature for time step \(i\) using temperature values from time step \(i\) on the right hand side.

   • How many indices do you need?
   • How many boundary conditions do you need?
   • Draw a stencil to illustrate the spatial grid

2. Reorganize the equation so that on the left hand side you have all the \(T^i\) values and on the right hand side you have all the known temperature values (\(T^{i-1}\)).

   • Further reorganize the equation so that the left hand side is a sum of \(T^i\) values multiplied by some constant coefficients, i.e.

     \[a T_{n-1,k}^i + b T_{n,k-1}^i + c T_{n+1,k}^i + ... = \mathrm{something with T^{i-1}\]

The 2D heat equation can be discretized with forward difference in time and central difference in space (like before):

\[\rho Cp \frac{T_{n,k}^i - T_{n,k}^{i-1}}{\Delta t} = \alpha \left( \frac{T_{n+1,k}^i - 2T_{n,k}^i + T_{n-1,k}^i}{\Delta x^2} + \frac{T_{n,k+1}^i - 2T_{n,k}^i + T_{n,k-1}^i}{\Delta y^2} \right) + H\]

which can be reorganized to find the coefficients of the unknown temperature values:

\[T_{n-1,k}^i \left( \frac{-\alpha}{\Delta x^2} \right) + T_{n+1,k}^i \left( \frac{-\alpha}{\Delta x^2} \right) + T_{n,k}^i \left( \frac{2\alpha}{\Delta x^2} + \frac{2\alpha}{\Delta y^2} + \frac{\rho C_p}{\Delta t} \right) + T_{n,k+1}^i \left( \frac{-\alpha}{\Delta y^2} \right) + T_{n,k-1}^i \left( \frac{-\alpha}{\Delta y^2} \right) = H + T_{n,k}^{i-1} \rho C_p / \Delta t\]

This can again be shortened by assigning \(A = \frac{-\alpha}{\Delta x^2}\), \(C = \frac{-\alpha}{\Delta y^2}\) and \(B = \frac{2\alpha}{\Delta x^2} + \frac{2\alpha}{\Delta y^2} + \frac{\rho C_p}{\Delta t}\), leading to

\[A T_{n-1,k}^i + A T_{n+1,k}^i + B T_{n,k}^i + C T_{n,k+1}^i + C T_{n,k-1}^i = H + T_{n,k}^{i-1} \rho C_p / \Delta t\]

We have one of such equation for each spatial grid point (i.e. if there are \(N=3\) grid points in \(x\) direction and \(K=4\) grid points in \(y\) direction, there are total of twenty of these equation). Some of the equation have to be replaced by boundary conditions.

We can use the method shown in “Finite difference method, part II” to solve these equations. However, this time the indexing of the coefficient array \(\mathbf{M}\) is a bit more complicated.

Since we have \(N \times K\) unknowns, and equally many equations, \(\mathbf{M}\) has \(N \times K\) rows and \(N \times K\) columns. With 3 grid points in \(x\) direction and 4 grid points in \(y\) direction, this means that the size of \(\mathbf{M}\) is \(12 \times 12\). For each unique pair of \(n\) and \(k\), i.e. for each spatial grid point, there is one row in the coefficient matrix.

When calculating the value of \(T_{n,k}^i\), one fills the coefficients of the \(T\) s (as separated above) at row \(n \times K + k\) (assuming we start counting from zero!). This number is called the global index of the grid point and is a unique number for each grid point in both (all) directions.

For example, when filling in coefficients of the equation that calculate \(T\) at \(n = 1\) and \(k = 2\) the row in \(\mathbf{M}\) we are looking at row \(1 \times 4 + 2\). On this row, the coefficients of the uknown temperatures will go to columns determined with the same rule: The coefficient for \(T_{n-1,k}^i\) (which is \(A\)) will go to column \((n-1)\times 4 + k\), i.e. \((1-1)\times 4 + 2\).

The whole right hand side of the discretized expression (that has only temperature values from the previous time steps) will go to the right hand side vector, row \(n \times K + k\).

If one wants to fill in a boundary condition \(T_{1,0}^i\) (i.e. the location \(y=0\)) the coefficient matrix will have a value of one at row \(1\times 4 + 0 = 4\), columnt \(1\times 4 + 0 = 4\).

Exercise

Script heat_diff2d_simple_implicit.py is an example of how to calculate one time step using the implicit finite difference method, in two spatial dimensions.

• Read through the script and make sure you understand it.

• Run the script in IPython (%run scriptname.py) and then print arrays like M or Tnew at your IPython prompt. How do they look like?

• Make the script calculate multiple time steps instead of only one:

  • Modify the script so that you add an for loop around the two existing loops (for this, you need to define nt in your code)
  • Note that the temperature array only records the current known temperature field: it does not store the previous time steps within it like in the previous, 1D cases.
  • At the end of the time loop, you need to copy (np.copy()) the new temperature array Tnew to be the known, “old” temperature array T

–