2D Heat advection

Goal: Find the steady state heat distribution of a 2D duct with moving air

Inital condition
At t = 0, the temperature of the air through the duct is 0.

Boundary conditions
The boundary conditions of the wall are T_left = 1, T_right = 0 (Constant temperature) and dT_top/dt, dT_bottom/dt = 0 (Adiabatic).
Each iteration in the simulation represents moving forward in time of dt seconds. dt = 0.001.
The air is moving sinusoidally (swirling) throughout the duct. The arrows represent the velocity vectors.

Just after initial state (0.1 seconds). Temperature is almost uniform.

After 1000 iterations (1 second), the temperature starts to pick up.

After 10000 iterations (10 seconds), the temperature gradients have spread through the whole duct.

After 1e5 iterations (100 seconds), temperature gradients are weak and profile seems resolved.

After 1e6 iterations (17 minutes), the temperature has reached steady state.

Heat equation derivation

To get the numerical code, I started with the analytical solution, then used the Finite-Volume method for the discretization in space and Forward Euler for the discretization in time. The Finite-Volume method ensures that conservation is preserved over cells. The Forward Euler is used to ensure stability. Using an upwind scheme, we also ensure that the information flows downstream. For the final code, the Courant number is used for compactness and for stability checks.

Here is the full derivation of the numerical code:

Numerical heat equation code

Using the worked out equation, for loops are used to plot this temperature equation for each grid cell. This following part of the code is used to define the new temperature of each grid cell. Here is the main part of the heat equation code written manually:

for i in range(1, Nx+1): #<-------- Loops over each grid cell

for j in range(1, Ny+1):

uw = (ubig[j, i] + ubig[j, i-1]) / 2 #<-------- Defining the average velocity between cells

ue = (ubig[j, i+1] + ubig[j, i]) / 2

vs = (vbig[j, i] + vbig[j-1, i]) / 2

vn = (vbig[j+1, i] + vbig[j, i]) / 2

Tw = T_old[j, i-1] if uw >= 0 else T_old[j, i] #<-------- Important for upwind scheme!

Te = T_old[j, i] if ue >= 0 else T_old[j, i+1]

Ts = T_old[j-1, i] if vs >= 0 else T_old[j, i]

Tn = T_old[j, i] if vn >= 0 else T_old[j+1, i]

Cow = uw * dt / dx #<-------- Courant number used in heat equation

Coe = ue * dt / dx

Cos = vs * dt / dy

Con = vn * dt / dy

if any((np.abs(Cow) > 1, np.abs(Coe) > 1, np.abs(Cos) > 1, np.abs(Con) > 1)): #<--------- Stability check

print("CLS instability", "Cow:", Cow,"Coe:", Coe,"Cos:", Cos,"Con:", Con)

sys.exit()

T[j, i] = ( #<------- Full heat equation

(Cow * Tw - Coe * Te) / 4 +

(Cos * Ts - Con * Tn) / 4 +

Fox * (T_old[j, i-1] - 2*T_old[j, i] + T_old[j, i+1]) +

Foy * (T_old[j-1, i] - 2*T_old[j, i] + T_old[j+1, i]) +

T_old[j, i]

)

More detail

Additionally, the effect of grid resolution has been investigated with 40x40 grid points, 60x60 gridpoints and 80x80 gridpoints. Also, the centerline temperature has been plotted in both the horizontal and vertical direction to find the effect of grid more visually clear. It can be seen that the highest grid resolution has the lowest average temperature in the middle.

Furthermore, it is interesting to find when a grid resolution is "good enough". Therefore, a non-dimensional number is constructed comparing the root mean square temperature difference between grid resolutions after N_iter iterations and the left wall temperature. My own interpretation would be that 1% (0.01) difference in this non-dimensional number for a N+20 grid resolution increase is sufficient, which is a 80x80 grid.

Lastly, it was interesting to see the actual improvement in steady state after each iteration by finding the maximum gradient of each state. It can be concluded that the number of grid cells does not affect the speed of convergence, and that after 1e5 iterations, the gradient has a linear-log relationship with the number of iterations.