Solving the heat equation in 2D numerically with a forward in time and centered in space scheme (FTCS) and later visualizing the predicted diffusion of heat given different initial distributions and diffusivity values. The spatial domain is a square and the boundary conditions are Dirichlet type but the code can be easily adapted to any problem by just making a change in one line.

How is it solved?

The model describing the temperature of a 2D object (or 3D which is short enough) is :

with $T(x,y,0) = T_0$ (initial temperature distribution) and prescribed boundary conditions (Dirichlet, Neumann, Periodic, Mixed,...)

where $T(x,y,t)$ is the temperature at position $x,y$ after $t$ units of time.The constant $\alpha$ is known as the diffusivity coefficient and it has units of area per time, representing the speed at which the heat diffuses, which may be different across various materials. The algorithm I implemented discretizes the above equation and iteratively solves the following finite difference equation based on the boundary values and initial distribution.

Important note : For a uniform spatial mesh, the stability condition for the heat equation requires that $\ \Delta t < \frac{\Delta x^2}{4D} $.

Model simulation and visualization:

Some snapshots in time of the temperature distribution (seen from above) with diffusivity coefficient $\alpha = 1$ units of area per time with initial distribution $T_0(x,y) = (1+cos(\pi x))(1+cos(\pi y))$

Movies for $\alpha = 1$

Movies for $\alpha =0.1$

As expected, making the diffusivity coefficient 10 times smaller results in a slower diffusion of heat.

In both cases due to the boundary conditions, the distribution decays exponentially in time to $z = 0$ since the ends of the body are not insulated.