Python热工水力程序开发（基于实验报告）

此程序为热工水力计算程序，按照热工水力设计说明书进行开发（借鉴了“不想穿靴子的猫”大佬的代码并进行改动优化）
通过Python计算基础参数、冷却剂出口温度、六段控制体冷却剂出口温度、膜温差、包壳内外温度；六段控制体的qDNB及DNBR，各段总压降及压力损失，程序代码如下:

O_Base.py —>基本参数计算模块

print('\033[1;35m----------------------基本参数--------------------------------\033[0m')
import math
a=math.pi
def csjc():
    fnr,fnz,feh,l = 1.35,1.53,1.08,3.8
    dcs=9.6e-3
    Fa,Nt,plxs =0.974,1840000000,0.05
    W=32600/3.6

    S=a*121*265*dcs*l
    q_ =Nt*Fa/S
    Af=121*265*((12.5e-3)*(12.5e-3)-a/4*dcs*dcs)+121*4*17*(12.5e-3)*(0.4e-3)
    Ab=(12.5e-3)*(12.5e-3)-a/4*dcs*dcs
    Wu=W*((1-plxs)/Af)*Ab
    De=4*Ab/a/dcs
    q_max=q_*fnr*feh*fnz
    ql_=q_*a*dcs
    qlmax=ql_*fnr*fnz*feh
    hmax=Nt*feh*fnr/W/(1-0.05)/1000
    Cm=W/Af*3600  
    Fs=1+0.03*(Cm/4.882/1000000)*(0.04/0.019)**0.35
    b=12.5/dcs/1000#栅距比包壳外径
    hl=0.042*b-0.024#水纵向流过平行棒束时的传热系数

    print("燃料棒传热面积："+str(S))
    print("燃料棒表面热流密度："+str(q_))
    print("燃料棒最大热流密度："+str(q_max))
    print("燃料棒平均线功率："+str(ql_))
    print("燃料棒最大线功率："+str(qlmax))
    print("单元流道面积："+str(Ab))
    print("堆芯燃料元件周围的冷却剂总有效流通面积："+str(Af))
    print("单元通道流量："+str(Wu))
    print("一个栅元中冷却剂的当量直径："+str(De))
    print("堆芯最大焓差："+str(hmax))
    print("冷却剂质量流速："+str(Cm))    
    print("定位格架搅混因素修正系数："+str(Fs))   
    print("栅距比包壳外径："+str(b))   
    print("水纵向流过平行棒束时的传热系数："+str(hl))  

    print('\033[1;35m-----------------冷却剂出口温度-------------------------------\033[0m')
print(csjc())

A_Out_tem.py —>冷却剂出口温度及平均管法压降计算模块

import math
a=math.pi
def TemperOut():
    tfin,Fa,Nt,plxs,GuessTempOut1 = 288,0.974,1840000000,0.05,325
    p_=709.2584467#热管内冷却剂平均密度
    n1=86.12219066#*10^(-6)平均动力黏度
    n2=0.117575332#平均运动黏度
    De,dcs,l=0.011123299881757214,9.6,3.8    
    W=32600/3.6
    Af=2.7303579692033026
    TempOut=[0]
    t_   = [i for i in range(0,1)] 
    v    = [i for i in range(0,1)] #热管内冷却剂平均流速
    Re   = [i for i in range(0,1)]#雷诺数
    r    = [i for i in range(0,1)]#半径比
    f    = [i for i in range(0,1)]#摩擦阻力系数
    pf   = [i for i in range(0,1)]#摩擦压降
    fpel = [i for i in range(0,1)]#单相流体提升压降
    fpout= [i for i in range(0,1)]#计算出口局部压降
    fpin = [i for i in range(0,1)]#计算进口局部压降
    pgj  = [i for i in range(0,1)]#计算定位格架出口压降
    
    for i in range (0,1):
        flag = True
        count = 0
        while (flag):           
            tav = (GuessTempOut1 + tfin) / 2
            cp=1000*(0.04006*(tav-310)+5.7437)                           
            TempOut[i] =float(tfin + Fa*Nt/(W*(1 - plxs)*cp))
            if((math.fabs(TempOut[i] - GuessTempOut1)) <= 0.5):
                flag = False
            GuessTempOut1 = TempOut[i]
            count += 1
            t_[i]=(tfin+TempOut[i])/2
            v[i]=W*(1-plxs)/Af/p_
            Re[i]=p_*v[i]*De/n1*1000000
            r[i]=(4/a)**0.5*(12.5/dcs)
            f[i]=64*pow((pow(r[i],2)-1),3)/(-3*pow(r[i],4)+4*pow(r[i],4)*math.log(r[i])-1+4*pow(r[i],2))/Re[i]
            pf[i]=f[i]*l/De*p_*pow(v[i],2)/2
            fpel[i]=p_*9.8*l
            fpout[i] =0.75*p_*pow(v[i],2)/2
            fpin[i] =1*p_*pow(v[i],2)/2 
            pgj[i]  =1.05*p_*pow(v[i],2)/2                                   
            
        print( "冷却剂出口温度：" + str(TempOut[i])+ " 迭代次数：" + str(count))
        print( "热管平均温度：" + str(t_[i]))
        print( "热管内冷却剂平均密度：" + str(p_)) 
        print( "热管内冷却剂平均流速：" + str(v[i]))
        print( "雷诺数：" + str(Re[i]))
        print( "半径比：" + str(r[i]))
        print( "摩擦阻力系数：" + str(f[i]))
        print( "摩擦压降：" + str(pf[i]))
        print( "单相流体提升压降：" + str(fpel[i]))
        print( "计算出口局部压降：" + str(fpout[i]))
        print( "计算进口局部压降：" + str(fpin[i]))
        print( "定位格架出口压降：" + str(pgj[i]))
        print('\033[3;0m-------------------控制体出口温度-------------------------\033[0m')
        
print(TemperOut()) 

B_Control_body.py—>控制体冷却剂出口温度计算模块

import math
from A_Out_tem import TemperOut
a=math.pi
#from scipy import integrate
def ControlBody():
    tfin = 288
    GuessTempOut1 = 325
    fnr,feh,fehm,flag,l = 1.35,1.08,0.95,True,3.8
    dcs=9.6e-3
    Fa,Nt,plxs =0.974,1840000000,0.05
    W=32600/3.6

    TempOut=[0,0,0,0,0,0]
    wc=[i for i in range(0,6)]
    q_ =float( Nt*Fa/(a*121*265*dcs*l))
    gyh = [0.48,1.02,1.5,1.56,0.96,0.48]
    Af=float(121*265*((12.5e-3)*(12.5e-3)-a/4*dcs*dcs)+121*4*17*(12.5e-3)*(0.8e-3))
    Ab=float((12.5e-3)*(12.5e-3)-a/4*dcs*dcs)
    Wu=float(W*((1-plxs)/Af)*Ab)

    for i in range (0,6):#一到六号控制体
        flag = True
        count = 0
        while (flag):          
           tav = (GuessTempOut1 + tfin) / 2           
           if(280<tav<300):
                cp=(5455-5068)*(tav-280)/20+5068
           if(300<tav<320):
                cp=(6143-5455)*(tav-300)/20+5455
           if(320<tav<340):
                cp=(7837-6143)*(tav-320)/20+6143

           TempOut[i] = float(tfin +q_*fnr*feh*fehm*a*dcs*l/6*gyh[i]/Wu/cp)
           wc[i]=TempOut[i]-GuessTempOut1
           if((math.fabs(TempOut[i] - GuessTempOut1)) <= 0.0000001):
                flag = False
           GuessTempOut1 = TempOut[i]
           count += 1
        b=(TempOut[i]+tfin)/2
        print("第"+str(i+1) + "级控制体出口温度：" + str(TempOut[i])+ " 迭代次数：" + str(count))
        print("第"+str(i+1) + "级控制体出口温度迭代误差：" + str(wc[i]))
        print("第"+str(i+1) + "级冷却剂平均温度：" + str(b))
        print("第"+str(i+1)+"级冷却剂定压比热容："+str(cp/1000))       
        tfin = TempOut[i]
    print('\033[1;35m----------------------包壳外壁温度--------------------------------\033[0m')
    #print(OutTemp)
#ControlBody()
    return TempOut
print(ControlBody())

C_Out_cladding.py —>外壁温度计算模块

import math
from B_Control_body import ControlBody
def OutCladding():
    q_,fnr,feh=487686.61178152927,1.35,1.08
    z=[746.5683956,733.4723679,710.8081004,681.6173182,655.467717,639.3956914]#第i级冷却剂平均温度下的密度
    u=[0.1239780,0.1225427,0.1204575,0.1183717,0.1169534,0.1162529]#第i级运动黏度
    hin=[1273.464617,1296.903322,1345.753378,1416.461225,1489.037581,1534.144778]#第i级进口处水的焓值
    t2=[292.9597669,297.3550564,306.2566258,318.4363807,329.8977537,336.3793903]#第i级半高处水的温度
    h2_=[1285.137081,1321.099432,1380.516159,1451.924545,1511.176831,1545.471293]#第i级平均温度下的焓值
    h3=[1459.393355,1392.758221,1291.34336,1181.642597,1111.225501,1075.073471] #第i级平均温度下的汽化潜热
    pr = [0.836168658,0.848066237,0.876189853,0.928548604,0.995219388,1.048761284]
    k  =[0.579573966,0.569416219,0.552003362,0.530123058,0.511274057,0.500130883]#第i级冷却剂的热导率

    
    gyhcs = [0.48,1.02,1.5,1.56,0.96,0.48]
    W=32600/3.6
    fnz=1.53
    q_max=1087902.0323655286
    Cm=11762592.057525925
    Af=2.7303579692033026
    De=0.011123299881757214
    Fs=1.0952091584565196#定位格架搅混因素修正系数
    Fc=1.05#第i级热流密度不均匀修正因子
    hmax=311.843461414272#堆芯最大焓差
    hl=0.030687500000000006
     
    Cp    = [i for i in range(0,6)]
    v_    = [i for i in range(0,6)]
    Re    = [i for i in range(0,6)]
    Nu    = [i for i in range(0,6)]
    hd    = [i for i in range(0,6)]#第i级冷却剂的对流换热系数
    Tq    = [i for i in range(0,6)]
    tw    = [i for i in range(0,6)]#第i级包壳外壁温度
    h2    = [i for i in range(0,6)]#第i级半高处水的焓值
    x     = [i for i in range(0,6)]#第i级平衡含汽率
    qDNB  = [i for i in range(0,6)]#第i级临界热流密度
    DNBR  = [i for i in range(0,6)]#烧毁比
    tfout = [292.5172621141466, 301.8778004061037, 314.8246254054336, 327.23480187283803, 333.99151925799544, 337.17923488170857]    
    f     =     [i for i in range(0,6)]
    f1    = [i for i in range(0,6)]
    f2    = [i for i in range(0,6)]
    
    for i in range(0,6):
        v_[i]=float(W/Af/z[i]*(1-0.05))
        if(280<=tfout[i] and tfout[i]<300):                #差分法求定压比热容
            Cp[i]=(tfout[i]-280)*0.387/20+5.068
        if(300<=tfout[i] and tfout[i]<320):
            Cp[i]=(tfout[i]-300)*0.6886/20+5.455
        if(320<=tfout[i] and tfout[i]<=340):
            Cp[i]=(tfout[i]-320)*1.694/20+6.1436      

        Re[i]=v_[i]/u[i]*1000000*De
        Nu[i]=pow(Re[i],0.8)*pow(pr[i],(1/3))*hl
        hd[i]=Nu[i]*k[i]/De
        Tq[i]=q_*fnr*feh*gyhcs[i]/hd[i]
        tw[i]=Tq[i]+tfout[i]
        h2[i]=hin[i]+hmax/12
        x[i]=(h2_[i]-1649.671945)/h3[i]
        a1 = 2.022 - 6.238*10** (-8)*16000000
        a2 = 0.1722 - 1.427*10** (-8)*16000000
        d = 0.2664 + 0.8357*math.exp(-124 * De)
        a3 = math.exp((18.177 - 5.987*10** (-7)*16000000)*x[i])
        a = 3.154 * 10** 6 * (a1 + a2*a3)
        b = (0.1484 - 1.596*x[i] + 0.1729*x[i]*math.fabs(x[i]))*0.2049*Cm/ 1000000 + 1.037
        c = 1.157 - 0.869*x[i]
        e = (0.8258 + 0.341*10**(-6)*(1649.671945-hin[i])*1000)
        qDNB[i] = a*b*c*d*e*Fs/Fc
        DNBR[i]= qDNB[i]/q_/fnr/feh/gyhcs[i]
                                
        f1[i]=q_max*gyhcs[i]/hd[i]/fnz#第i级控制体膜温差
        f2[i]= (347.3565346-t2[i]) +25*(q_max/fnz*gyhcs[i]/1000000)**(0.25)*math.exp(-16/6.2)#第i级控制体径向膜温差
        if (f1[i]>f2[i]):
            f[i] = f2[i]
        else:
            f[i] = f1[i]
        f[i] = f[i] + tfout[i]#第i级控制体燃料元件表面最高温度

        print("第"+str(i+1) +"级半高处冷却剂流速："+str(v_[i]))
        print("第"+str(i+1) +"级动力粘度："+str(u[i]))
        print("第"+str(i+1) +"级雷诺数："+str(Re[i]))
        print("第"+str(i+1) +"级普朗特数："+str(pr[i]))
        print("第"+str(i+1) +"级努塞尔数："+str(Nu[i]))
        print("第"+str(i+1) +"级冷却剂热导率："+str(k[i]))
        print("第"+str(i+1) +"级冷却剂的对流换热系数："+str(hd[i]))
        print("第"+str(i+1) +"级强迫对流放热温差："+str(Tq[i]))
        print("第"+str(i+1) +"级包壳外壁温度："+str(tw[i]))
        print("第"+str(i+1) +"级半高处水的焓值："+str(h2[i]))
        print("第"+str(i+1) +"级平衡含汽率："+str(x[i]))
        print("第"+str(i+1) +"级临界热流密度："+str(qDNB[i]))
        print("第"+str(i+1) +"级烧毁比："+str(DNBR[i]))
        print("第"+str(i+1) +"级燃料元件表面最高温度："+str(f[i]))
               
        print("θf1=" + str(f1[i]) + ", θf2=" + str(f2[i]) + ", f=" +str(f[i]))
        tcs = [f[i] for i in range(0,6)]
    print('\033[1;35m----------------------包壳内壁温度--------------------------------\033[0m')
    #print(tcs)
    return tcs#传递给下一个值
print(OutCladding())

D_In_envelope.py—>包壳内壁温度计算

from C_Out_cladding import OutCladding
import math
a=math.pi
def Ftci():
    fnz,dcs,dci= 1.53,9.6,8.5
    flag = True
    gyhcs = [0.48,1.02,1.5,1.56,0.96,0.48] #归一化参数
    qlmax=32810.35231396752#燃料棒最大线功率
    tcs = OutCladding()   #包壳外壁温度
    #print(tcs)
    tav = 300 #包壳内外假设平均温度
    tci = [i for i in range(0,6)]
    kc  = [i for i in range(0,6)]
    tav1= [i for i in range(0,6)]
    ta  = [i for i in range(0,6)]
    for i in range(0,6):
        flag = True #初始化
        while (flag):
            kc[i] = 0.0547*(1.8*tav + 32) + 13.8
            tci[i] = tcs[i] +qlmax/fnz*gyhcs[i]*(0.55/1000)/a/kc[i]/((dci+dcs)/2000)
            tav1[i]=(tcs[i]+tci[i])/2
            ta[i]= tav1[i]-tav
            if (math.fabs(tav1[i] - tav) <= 0.00001):
                flag = False
            tav = tav1[i]
        print("第"+str(i+1) +"级控制体的包壳热导率" + str(kc[i]))
        print("第"+str(i+1) +"级包壳内壁温度为" + str(tci[i]))  
        print("第"+str(i+1) +"级包壳内壁温度迭代误差为" + str(ta[i])) 
    print('\033[1;31m------------------芯块表面及中心-----------------------------\033[0m')
    tci=[tci[i] for i in range(0,6)]
    return tci
print(Ftci())

E_Out_core.py—>燃料芯块表面及中心最高温度计算

from D_In_envelope import Ftci
import math
a=math.pi
def CoreBlock():
    tci = Ftci()#包壳外壁温度/℃    
    to =565.273 #假设二氧化铀中心温度/℃       
    fnz,du,dci=1.53,8.35,8.5
    qlmax=32810.35231396752#线功率
    gyhcs = [0.48,1.02,1.5,1.56,0.96,0.48] 
    
    tu = [i for i in range(0,6)] #芯块表面温度
    k  = [i for i in range(0,6)]#误差
    K0 = [i for i in range(0,6)]#第i级燃料芯块表面最高温度下的积分热导率
    K1 = [i for i in range(0,6)]#第i级燃料芯块中心最高温度下的积分热导率
    K2 = [i for i in range(0,6)]#第i级燃料芯块中心假定最高温度下的积分热导率

    for i in range(0,6):
        tu[i] = tci[i] + qlmax/fnz*gyhcs[i]/a/5678/((du+dci)/2000)
        K0[i]=40.5*math.log(464+tu[i])+0.027366*math.exp(2.14*10**(-3)*tu[i])-248.04#第i级燃料芯块表面最高温度下的积分热导率
        print("第"+str(i+1) +"级燃料芯块表面最高温度" + str(tu[i]))
        flag = True #初始化
        step=1.0    #每一次更新的步长(可修改)
        count=0     #反转次数计数
        minus=1     #k2是否小于k1 1说明小于  -1说明大于

        while (flag):   
            to+=step                       
            K1[i]=K0[i]+qlmax/fnz*gyhcs[i]/4/a/100
            K2[i]=40.4*math.log(464+to)+0.027366*math.exp(2.14*10**(-3)*to)-248.04
            k[i]=math.fabs(K2[i]-K1[i])
            #print(to)
            if (k[i] <= 1):
                print('对于第'+str(i+1)+'级燃料芯块中心最高温度t0：'+str(to))
                to=565.273
                flag = False
            elif (K2[i]-K1[i])*minus > 0:
                minus=-minus
                count=count+1
            to=to+minus*step*10**(-count)
            
        print("第"+str(i+1) +"级控制体积分热导率迭代误差" + str(k[i]))
        print("第"+str(i+1) +"级燃料芯块表面最高温度下的积分热导率" + str(K0[i]))
        print("第"+str(i+1) +"级燃料芯块中心最高温度下的积分热导率" + str(K1[i]))
        print("第"+str(i+1) +"级燃料芯块中心假定最高温度下的积分热导率" + str(K2[i]))
    print('\033[1;34m----------------------qDNB及DNBR-------------------------------\033[0m')
    return tu
print(CoreBlock())

F_Critical_heat.py—>控制体DNBR以及qDNBR的计算（更为准确）

import math
from E_Out_core import CoreBlock
def FqDNB():   
    hin=[1273.464617,1296.903322,1345.753378,1416.461225,1489.037581,1534.144778]#第i级进口处水的焓值
    Cm,q_  =11939826.340614386,487686.61178152927
    P=16000000
    Fs,Fc=1.0952091584565196,1.05
    De = 0.011123299881757214
    h2_=[1285.137081,1321.099432,1380.516159,1451.924545,1511.176831,1545.471293]#第i级平均温度下的焓值
    h3=[1459.393355,1392.758221,1291.34336,1181.642597,1111.225501,1075.073471] #第i级平均温度下的汽化潜热
    fnr,feh=1.35,1.08
    gyhcs = [0.48,1.02,1.5,1.56,0.96,0.48]     
        
    a1 = 2.022 - 6.238*pow(10, -8)*P
    a2 = 0.1722 - 1.427*pow(10, -8)*P
    d = 0.2664 + 0.8357*math.exp(-124 * De)
    qDNB = [i for i in range(0,6)]
    DNBR = [i for i in range(0,6)]
    x    = [i for i in range(0,6)]

    for i in range(0,6):
        x[i]=(h2_[i]-1649.671945)/h3[i]
        a3 = math.exp((18.177 - 5.987*pow(10, -7)*P)*x[i])
        a = 3.154 * pow(10, 6) * (a1 + a2*a3)
        b = (0.1484 - 1.596*x[i] + 0.1729*x[i]*math.fabs(x[i]))*0.2048*(Cm/ 1000000) + 1.037
        c = 1.157 - 0.869*x[i]
        e = (0.8258 + 0.341*pow(10, -6)*(1649.671945 - hin[i])*1000)
        qDNB[i] = a*b*c*d*e*Fs/Fc
        DNBR[i]= qDNB[i]/q_/fnr/feh/gyhcs[i]
        print(str(i+1) + "级控制体qDNB：" + str(qDNB[i])+ " 控制体DNBR：" + str(DNBR[i]))
    print('\033[1;35m------------------总压降------------------------------------\033[0m')
print(FqDNB())

G_Pressure（1）.py---->分级算总压降

import math
def press():    
    z=[746.5683956,733.4723679,710.8081004,681.6173182,655.467717,639.3956914]#第i级冷却剂平均温度下的密度
    v_=[4.220360151,4.295713982,4.432683736,4.622516803,4.806930112,4.927758427]#第i级的半高处冷却剂流速
    Re =[378650.3806,389925.297,409323.4066,434374.3015,457181.4023,471497.5239]

    fpout=5248.859896008969 #计算出口局部压降
    r = 1.4692437071556155
    l,De =3.8,0.011123299881757214
    pin=6998.479861345293#计算进口局部压降
    pgj=7348.403854412558#定位格架出口压降
    f=[i for i in range(0,6)]#摩擦阻力系数
    pf=[i for i in range(0,6)]#计算摩擦压降
    fpel=[i for i in range(0,6)]#单相流体提升压降  
    ts,mc=0,0
    
    for i in range(0,6):
        f[i] = 64*pow((pow(r,2)-1),3)/(-3*pow(r,4)+4*pow(r,4)*math.log(r)-1+4*pow(r,2))/Re[i]
        pf[i]=f[i]*l/6/De*z[i]*pow(v_[i],2)/2
        fpel[i]=z[i]*9.8*l/6     
        print("第"+str(i+1) + "级摩擦阻力系数:" + str(f[i]))    
        print("第"+str(i+1) + "级摩擦压降:" + str(pf[i]))    
        print("第"+str(i+1) + "级提升压降:" + str(fpel[i]))    
        mc = pf[0]+pf[1]+pf[2]+pf[3]+pf[4]+pf[5]
        ts = fpel[0]+fpel[1]+fpel[2]+fpel[3]+fpel[4]+fpel[5]
        compress = mc+ts+fpout+pin+pgj       
    print("总摩擦压降:" + str(mc))    
    print("总提升压降:" + str(ts))    
    print("总的压力损失:" + str(compress) )
    print("------------------------END----------------------------------")
print(press())

G_Pressure.py---->平均管算总压降

def press():

    pf=703.3284910230612#计算摩擦压降
    fpel=26574.073121512003#单相流体提升压降
    fpout=5217.002488075618 #计算出口局部压降
    fpin=6956.003317434157#计算进口局部压降
    pgj=7303.8034833058655#定位格架出口压降

    compress =pf+fpel+fpout+fpin+pgj
    print("总的压力损失:" + str(compress) )
    print("------------------------END----------------------------------")
print(press())

使用这些代码的时候你会有很多疑问，有些数值为什么凭空出现，因为那些数值要么是给定的，要么的查表得出的，这里你会发现，我的代码不够‘自动化’：不能自己查对应的数值，例如焓、雷诺数等，这里给以后的学弟学妹门一个建议：你们可以试着找一些库去优化我这些代码：可以找到对应的热力参数，这样的话数值会更准确，代码也更加‘智能’。