Python loop issue: finite volume/difference iteration

Question

I'm having trouble understanding why this for loop is not working. Full code is also below.

I want to calculate the temperature at future time points (j+1) based on the temperatures at the previous time point (j) using the finite volume method. I have setup functions to calculate the conduction in, out, and the internal generation within a node. I have another function which sums these together and determines the resultant temperature based on an energy balance equation. My issue though is specifically with the loop below:

for j in range(1, n): # time
    for e in range(1, i-1): # position
        T[e,j] = Tn_plus_1( T[e, j-1] , timestep,rho[e], cp[e], re_nodes[e] , rw_nodes[e],
                
                conduction_in(k[e-1], rw_nodes[e],T[e, j-1], T[e-1, j-1], r_nodes[e], r_nodes[e-1]),
                
                conduction_out(k[e], re_nodes[e], T[e+1, j-1], T[e, j-1], r_nodes[e+1], r_nodes[e] ),
                
                fraction[j]*internal_gen(q_dens[e], timestep, rho[e], cp[e]) )

Tn_plus_1 is a function to calculate the temperature of a node at the next timestep based on an energy balance and finite volume method. Conduction in, out, and internal generation are functions to calculate the heat transfer due to each mechanism into or out of a particular node. I set up for loops to iterate over time and position to determine the temperature evolution in each node for the duration of the analysis.

The issue is that the conduction out term does not apply correctly to node 7. I suspect the same issue is occurring with the conduction in term at the second node. The thought process behind the loop is that I want to calculate the temperatures of the internal nodes, and maintain the boundary condition on the two external nodes. Hence, I set the loop between 1 and i-1 (with 0 and i being the boundaries).

There are 9 position nodes in total (r, re, and rw each have 9 values). There are also 9 values for the material properties (rho, cp, k).

I'm new to Python and suspect the issue is to do me getting muddled up with Python indexing, but I don’t understand where exactly the loop is going wrong for the two nodes next to each boundary. Thanks in advance.

import numpy as np
import matplotlib.pylab as plt

empty = np.zeros(9) # empty array to fill with required properties

#material properties (will vary with position at a later stage)
k = empty + 23
rho = empty + 1600
cp = empty + 1700

#node spacing
dr_internal = 0.15
dr_external = 0.14

r_nodes = np.linspace(0, 1.2, 9) # nodes, m

r_nodes2 = np.insert(r_nodes,9 , r_nodes[-1]+dr_external) #inserting additional point to determine
east boundary on last node
re_nodes = (r_nodes2[1:] + r_nodes2[:-1]) / 2 #east locations
rw_nodes = np.insert(re_nodes, 0, 0)[:-1] #west locations

#q density for the internal generation term
q_dens = empty + 4400000
q_dens[0] = 0
q_dens[1] = q_dens[1]*0.625
q_dens[-1] = q_dens[-1]*0.5

#setting up time increments for the analysis
timestep = 15
steps = 100
total_time = timestep * steps
time = np.arange(0, steps*timestep, timestep)

#defining the number of time points and spacial nodes
i = len(r_nodes)
n = len(time)
T = np.zeros((i,n))

#initialising temperatures based on the boundary condition and initial condition
IC = 873
BC = [873, 873]
T[0,:] = BC[0]
T[-1,:] = BC[1]
T[:,0] = IC

#energy generation is time dependent, so creating a fraction to multiply internal generation term by
fraction=np.zeros(steps)
for j in range(1, steps):
    time2 = j*timestep
    fraction[j] = 0.06*time2**(-0.2)
fraction[0] = 0.06


def conduction_in( kw, rw, Ti, Ti_minus_1, ri, ri_minus_1):
    return (-kw*rw*((Ti-Ti_minus_1)/(ri-ri_minus_1)))

def conduction_out(k, re, Ti_plus_1, Ti, ri_plus_1, ri):
    return (k*re*((Ti_plus_1-Ti)/(ri_plus_1-ri)))

def internal_gen(q_dens, dt, rho, cp):
    return q_dens*dt/(rho*cp)

def Tn_plus_1(Ti, dt, rho, cp, re, rw, conduction_in, conduction_out, internal_gen):    
        return Ti + ((2*dt/(rho*cp*(re**2-rw**2)))*(conduction_in - conduction_out ))+ internal_gen

for j in range(1, n): # time
    for e in range(1, i-1): # position
        T[e,j] = Tn_plus_1( T[e, j-1] , timestep,rho[e], cp[e], re_nodes[e] , rw_nodes[e],
                
                conduction_in(k[e-1], rw_nodes[e],T[e, j-1], T[e-1, j-1], r_nodes[e], r_nodes[e-1]),
                
                conduction_out(k[e], re_nodes[e], T[e+1, j-1], T[e, j-1], r_nodes[e+1], r_nodes[e] ),
                
                fraction[j]*internal_gen(q_dens[e], timestep, rho[e], cp[e]) )
        

plt.plot(T[:,99]) 

What I get:

What I expect:

Answer 1

First, the "good" news: (with your current arrangement for conductive fluxes - but see the end of this post) you only need to correct the sign of your conductive flux

return (k*re*((Ti_plus_1-Ti)/(ri_plus_1-ri)))

should be (noting the minus sign at the start):

return (-k*re*((Ti_plus_1-Ti)/(ri_plus_1-ri)))

Think very hard about which way conductive fluxes are taking heat into or out of your control volume.

Please put plt.show() at the end of your posted code, or it won't show anything.

In terms of indices, you are updating control volumes 1 to 7 (which is OK) and keeping nodes 0 and 8 as boundary conditions. Well, OK, but boundary conditions are supposed to be applied on the FACES of control volumes, not adjacent nodes.

I assume that r_nodes2 may one day hold "ghost nodes", but at the moment you aren't using that array.

You adjust the internal heat source in q[1] and q[8]. However, you aren't updating the [8] control volume, so the latter adjustment to q[] does nothing and probably wasn't what you meant to do.

Finally, you have significant redundancy. conduction_in() and conduction_out() DO FUNDAMENTALLY THE SAME THING - you just have to send them the correct parameters, in the correct order so that they give conductive heat flux in the direction intended. Indeed, if you correct the sign as suggested this will become more obvious, as they would then both give a heat flux in the "forward" (increasing-r) direction.