Python循环问题：有限体积/差分迭代

Question

我无法理解为什么这个 for 循环不起作用。完整代码也在下面。

我想使用有限体积法根据前一个时间点（j）的温度计算未来时间点（j+1）的温度。我有设置函数来计算节点内的传导输入、输出和内部生成。我有另一个函数将它们加在一起并根据能量平衡方程确定最终的温度。我的问题具体在于下面的循环：

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 是一个基于能量平衡和有限体积方法计算节点在下一个时间步的温度的函数。传入、传出和内部生成是用于计算由于进入或传出特定节点的每种机制而产生的热传递的函数。我设置了 for 循环来随着时间和位置进行迭代，以确定分析期间每个节点的温度演变。

问题在于传导出项不适用于节点 7。我怀疑第二个节点的传导项也出现同样的问题。循环背后的思维过程是我要计算内部节点的温度，并维护两个外部节点的边界条件。因此，我将循环设置在 1 和 i-1 之间（以 0 和 i 为边界）。

位置节点一共有9个（r、re、rw各有9个值）。材料属性还有 9 个值（rho、cp、k）。

我是 Python 新手，怀疑问题是让我对 Python 索引感到困惑，但我不明白每个边界旁边的两个节点的循环到底出了什么问题。预先感谢。

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]) 

我得到的：

我的期望：

Answer 1

首先，“好”消息：（根据您当前的导电通量安排 - 但请参阅本文的结尾）您只需将括号中的单个字符从 + 更改为 -

(conduction_in - conduction_out )

给予

(conduction_in + conduction_out )

认真思考导电通量以何种方式从控制体积中带走热量。

就索引而言，您正在更新控制卷 1 到 7（这是可以的）并保留节点 0 和 8 作为边界条件。好吧，但是边界条件应该应用于控制体的面，而不是相邻节点。

我假设 r_nodes2 有一天可能会保留“幽灵节点”，但目前您没有使用该数组。

您可以在 q[1] 和 q[8] 中调整内部热源。但是，您没有更新 [8] 控制音量，因此后面对 q[] 的调整没有任何作用，并且可能不是您想要做的。

最后，你有大量的冗余。

conduction_in()

和

conduction_out()

基本上做同样的事情 - 您只需以正确的顺序向它们发送正确的参数，以便它们按预期方向提供传导热通量。