python实现数值传热模拟VIP

题目：流体以均匀速度u在半径为R，长度为L圆管中流动（L/R=4），忽略径向温度变化。假定在进口处流体的温度为TL,出口处温度为T0，管内壁温度也为T0。其中：无量纲温度φ=(t-T0)/(TL-T0)；无量纲坐标X=x/R；边界条件：【X=0时,φ=0;X=1时,φ=1∂R/λ=1,uR/α=1】问题：当节点个数为4个时，求X=1，2，3时，各节点无量纲温度φ值，并绘图表示。并且，依次求解当节点个数n为8，16，32，64时各节点的φ值；n=8时,求X=1,2，3……6，7，8；时8个节点的φ值；n=16时,求X=1,2，3……14，15，16；时16个节点的φ值；n=32时,求X=1,2,3,……31，32，33；时32个节点的φ值；n=64时,求X=1,2,3……62,63,64；时64个节点的φ值；然后比较，网格划分粗细程度对数值解的影响。（其中解析解为∅(X)=eX-e-2Xe4-平衡关系式：①控制体内热量的增量=②由于热对流进入控制体的净热量+③由于导热进入控制体的

净热量+④流体与管壁间的对流换热进入控制体净热量

其中：

①稳态：①=0；

②热对流②=(𝑐×𝜌×𝑢×𝜋𝑅2)[(𝑡i+1+𝑡i)/2-(𝑡i-1+𝑡i)/2]=(𝑐×𝜌×𝑢×𝜋𝑅2)(𝑡i+1+-𝑡i-1)/2;

③导热③=𝜋𝑅2×λ(𝑡i+1-𝑡i/∆𝑥+(𝑡i-1-𝑡i)/∆𝑥]=𝜋𝑅2×λ(𝑡i+1-2𝑡i+𝑡i-1)/∆𝑥;

④对流换热④=2𝜋𝑅×∆𝑥×∂(𝑇0-𝑡i);代入平衡关系式：(𝑐×𝜌×𝑢×𝜋𝑅2)(𝑡i+1-𝑡i-1)/2+𝜋𝑅2×λ(𝑡i+1-2𝑡i+𝑡i-1)/∆𝑥+2𝜋𝑅×∆𝑥×(𝑇0-𝑡i)=0同时除以∆𝑥，得：(𝑐×𝜌×𝑢×𝜋𝑅2)(𝑡i+1-𝑡i-1)/2∆𝑥+𝜋𝑅2×λ(𝑡i+1-2𝑡i+𝑡i-1)/(∆𝑥)2+2𝜋𝑅×∂(𝑇0-𝑡i)=0

无量纲化控制方程得：(𝑐×𝜌×𝑢×𝜋𝑅)(φi+1-φi-1)/2X+𝜋×λ(φi+1-2φi+φi-1)/X2-2𝜋R×∂φi=0

又根据∂R/λ=1,uR/α=1，得：(φi+1-φi-1)/(2X)+(φi+1-2φi+φi-1)/X2-2φi=0

整理得：-(2-X)φi-1+(4+4X2)φi-(X+2)φi+1=0三节点关系式为：-aiφi-1+biφi-ciφi+1=di(1)于是有a1当节点数为n时,X=4/n；ai=2-4/n;bi=4+64/n2;ci=（4/n)+2;di=0；递推关系式：φi=Giφi+1+Qi(2)φi-1=Gi-1φi+Qi-1(3)联立(1)、(3)式，并与得(2)对比得：Gi=ci/（bi-aiGi-1）Qi=（di+aiQi-1）/（bi-aiGi-1）则有，G1=ci/bi,GN=0；Q1=d1/b1；则可依次求得G1，Q1；G2，Q2；……Gn，Qn；由φn=1,代入(3)后可依次求得φn-1，φn-2，φn-3……φ1。当总节点数为n=4，8，16，32，64时，使用python编程软件，写出代码如下:

#导入画图库importmatplotlib.pyplotasplt

#定义函数，以所需要的节点数作为变量，函数返回Q_S(无量纲温度)数组

defNumerical_heat(n):

#定义5个变量的空列表，以下将公式计算出的结果依次添加到列表中a,b,c,d,G；Q_S=[],[],[],[],[]

#构造a列表

foriinrange(1,n):

ifi==1:

aa=0

else:

aa=-(n/8)+(n*n/16)

a.append(aa)

#构造b列表

foriinrange(1,n):

3

ifi==1:

bb=2

else:

bb=(n*n/8+2)

b.append(bb)

#构造c列表

foriinrange(1,n):

ifi==1:

cc=1

elifi==n-1:

cc=0

else:

cc=((n*n/16)+(n/8))

c.append(cc)

#构造d列表

foriinrange(1,n):

ifi==n-1:

dd=((n*n/16)+(n/8))

else:

dd=0

d.append(dd)

#构造G列表

foriinrange(1,n):

ifi==1:

g=0

else:

g=c[i-1]/(b[i-1]-(a[i-1]*G[i-2]))

G.append(g)

#求Q_S的数值

foriinrange(n,0,-1):

ifi==n:

qs=1

elifi==0:

qs=0

else:

q_s=G[i-1]*Q_S[n-i-1]+Q[i-1]

Q_S.append(qs)

#以上遍历是从大到小排列，现将Q_S翻转：

new_Q_S=list(reversed(Q_S))

returnnew_Q_S

#计算各个节点的数值

#定义K列表，将对应NUM列表内的节点返回的new_Q_S装入列表

K=[]

NUM=[4，8，16，32，64]

foriinNUM:

k=Numerical_heat(i)

K.append(k)

print(K)

#计算各个节点的横坐标

X=[]

foriinrange(6):

x=range(1,len(K[i])+1)

X.append(x)

#定义画图空间

fig=plt.figure(figsize=(10,5))

foriinrange(1,6):

ax=fig.add_subplot(2,3,i)

ax.set_title("Temperatureprofile")

plt.xlabel("Nodecoordinates")

plt.ylabel("dimensionlesstemperature")

ax.plot(X[i-1],K[i-1],'-*',c="red",markersize=4)

plt.grid()

plt.tight_layout()

plt.show()

代码运行结果整理如下:1.当计算4个节点时，计算结果：节点位置X01234无量纲温度00.05820.15520.394412.当计算8个节点时，计算结果：节点位置X012345678无量纲温度00.02510.05020.08530.14040.22970.37520.612613.当计算16个节点时，计算结果：节点位置X01234567无量纲温度00.01260.02390.03530.04800.06330.08220.1061节点位置X89101112131415无量纲温度0.13640.17510.22470.28830.36970.47420.60810.7798节点位置X16无量纲温度1.00004.当计算32个节点时，计算结果：节点位置X01234567无量纲温度00.00650.01250.01810.02360.02910.03480.0409节点位置X89101112131415无量纲温度0.04750.05470.06260.07150.08140.09260.10520.1193节点位置X1617181920212223无量纲温度0.13540.15350.17400.19720.22340.25320.28690.3251节点位置X2425262728293031无量纲温度0.36830.41730.47280.53570.60690.68760.77900.8826节点位置X32无量纲温度15.当计算64个节点时，计算结果：节点位置X01234567无量纲温度00.00330.00650.00950.01240.01520.01800.0208节点位置X89101112131415无量纲温度0.02350.02620.02900.03180.03470.03770.04080.0440节点位置X1617181920212223无量纲温度0.04740.05090.05450.05840.06250.06680.07130.0761节点位置X2425262728293031无量纲温度0.08120.08660.09240.09850.10490.11180.11910.1268节点位置X3233343536373839无量纲温度0.13510.14390.15320.16310.17370.18490.19690.2096节点位置X4041424344454647无量纲温度0.22310.23750.25290.26920.28650.30500.32470.3457节点位置X4849505152535455无量纲温度0