I am solving the steady-state heat equation in two dimensions:
$$\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) = 0$$
I have verified my code against analytical solutions. The temperature field solves without a problem, but when I use the temperatures to calculate the magnitude of the gradient I get results that don't make sense (at least not to me).
Here is the temperature field result (which appears normal to me):
Basically, there is a high temperature boundary condition on the top of the right-hand side, and low temperature boundary condition on the left-hand side separated by a low conductivity barrier with a zero-flux boundary above it. But when I try to plot the calculated magnitude of the heat flux gradient, I get this:
My understanding is that the heat flux contours should be perpendicular to the temperature contours in the first plot. Instead, I'm seeing these localized hot-spots in areas where the heat flux is changing directions.
I have tried finer meshes and uniform meshes, they all show the same oddity. This is how I am calculating the magnitude of the heat flux:
$$|q| = \sqrt{q_x^2 + q_y^2}$$
where
$$q_x = -k\frac{dT}{dx}=-k\frac{T_{i+1,j}-T_{i-1,j}}{2\Delta x}$$
and
$$q_y = -k\frac{dT}{dy}=-k\frac{T_{i,j+1}-T_{i,j-1}}{2\Delta y}$$
Am I missing something here? I've spent days pouring over my code and making sure my T's are crossed and my I's are dotted. So I'm reaching out to the mathematics community in hopes that you all can identify something that I missed.