0

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 get

What I expect: What I expect

Improve this question
3
  • 1
    do not post the same question after it is closed (though this appears to be an improvement, it could have been edited into the original question and gone into the re-open queue) stackoverflow.com/questions/78402883/… Commented Apr 29, 2024 at 17:33
  • Are you coming from MATLAB? Looks like you really want to think in terms of 1-based indexing. Embrace the 0-based nature of Python (range(n) instead of range(1,m)) and re-review your code. Commented Apr 29, 2024 at 17:35
  • Also, descriptive variables names will help anyone trying to read your code in the future - including you after a couple weeks away. Commented Apr 29, 2024 at 17:36

1 Answer 1

Reset to default
0

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.

Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

I realise my mistake now. I suspected it was the for loop but was actually a mix up with the signs on the conduction terms. I'm not sure how I would apply the boundary conditions to the faces rather than creating an adjacent node. Have you got an example? Also, apologies - I forgot to mention, r_nodes2 is just created so that I can make an array for the east boundaries, as the final east boundary needs the adjacent node. This is more important in future when I expand the calculation. I will also take your suggestion on board on redundancy, I thought it might be less confusing, maybe not.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.