反应堆热工分析计算程序.docx
《反应堆热工分析计算程序.docx》由会员分享,可在线阅读,更多相关《反应堆热工分析计算程序.docx(11页珍藏版)》请在冰豆网上搜索。
反应堆热工分析计算程序
tfin=279.4;Fa=0.974;Nt=3400000000;Wt=14314;b=0.059;
tfout=340;e0=0.01;
whilee0>0.001
t0_=0.5*(tfout+tfin);
Cp_=1000*(0.04006*(t0_-310)+5.7437);
xi=tfin+Fa*Nt/(Wt*(1-b)*Cp_);
e0=(tfout-xi)/tfout;
tfout=xi%堆芯出口处温度
end
%热流密度计算
m=157;n=264;dcs=9.5e-3;L=4.2672;
q_=Fa*Nt/(m*n*pi*dcs*L)%燃料元件表面平均热流量
FqN=2.524;FqE=1.03;FDHE=1.085;
qmax=q_*FqN*FqE%最大热流量
ql_=q_*pi*dcs%平均线功率
qlmax=ql_*FqN*FqE%最大线功率
%平均管情况
B=17;S=12.6e-3;dx=0.4e-3;
Af=m*n*(S^2-pi/4*dcs^2)+m*4*B*S*dx;%总的流通截面积
tf_=0.5*(tfout+tfin)%热管平均温度
vf_=5.13e-6*(tf_-310)+0.0014189;
pf_=1/vf_;%平均密度
v=Wt*(1-b)/(Af*pf_);%平均流速
Ab=S^2-pi/4*dcs^2;%单元流通截面积
Wu=Wt*(1-b)*Ab/Af;%单元截面流量
%第一控制体温度计算
e11=0.01;tf1=300;L1=4.2672/6;fai1=0.4;
whilee11>0.001
t11_=0.5*(tf1+tfin);
Cp1_=1000*(0.02155*(t11_-290)+5.2428);
x1i=tfin+q_*FqN*pi*dcs*L1*fai1/(Wu*Cp1_);
e11=(x1i-tf1)/tf1;
tf1=x1i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);%单元通道当量直径
u1=944e-7;Pr1=0.85;k1=575.5e-3;%查得该温度下的热物性
Re1=Wu*De/(Ab*u1);
h1=0.023*Re1^0.8*Pr1^0.4*k1/De;%该处的对流换热系数
dtf11=q_*fai1*FqE/h1;%单相强迫对流放热公式算得的温压
ts=346.310791;P=15.51;
dtf12=25*(q_*fai1*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf1;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf11tcs1=tf1+dtf11
else
tcs1=tf1+dtf12
end
dci=8.93e-3;tci1=349;e12=0.01;
whilee12>0.001
t12_=0.5*(tci1+tcs1);
kc1=0.0547*(1.8*t12_+32)+13.8;
yi=tcs1+ql_*fai1*FqE/(2*pi*kc1)*log(dcs/dci);
e12=(yi-tci1)/yi;
tci1=yi%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu1=tci1+ql_*FqE*fai1*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d1_ku=ql_*FqE*fai1/(4*pi*100);
tu1_ku=(26.42-21.32)/(400-300)*(tu1-300)+21.32;
to1_ku=tu1_ku+d1_ku;
to1=(600-500)/(34.97-30.93)*(to1_ku-30.93)+500%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;
h1=1296.4746e+3;x1=(h1-hfs)/(hgs-hfs);%该点含汽量
qDNB1=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x1))*...
((0.1484-1.596*x1+0.1729*x1*abs(x1))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x1)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR1=qDNB1/(q_*FqE*fai1)%计算烧毁比
%第二控制体温度计算
fai2=0.8;L2=4.2672/6;e21=0.01;tf2=310;
whilee21>0.001
t21_=0.5*(tf1+tf2);
Cp2_=1000*(0.027625*(t21_-300)+5.4583);
x2i=tf1+q_*FqN*FDHE*FqE*pi*dcs*L2*fai2/(Wu*Cp2_);
e21=(x2i-tf2)/tf2;
tf2=x2i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);
u2=857e-7;Pr2=0.90;k2=547e-3;%查得该温度下的热物性
Re2=Wu*De/(Ab*u2);
h2=0.023*Re2^0.8*Pr2^0.4*k2/De;%该处的对流换热系数
dtf21=q_*fai2*FqE/h2;%单相强迫对流放热公式算得的温压
ts=346.310791;P=15.5;
dtf22=25*(q_*fai2*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf2;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf21tcs2=tf2+dtf21
else
tcs2=tf2+dtf22
end
dci=8.93e-3;tci2=349;e22=0.01;
whilee22>0.001
t22_=0.5*(tci2+tcs2);
kc2=0.0547*(1.8*t22_+32)+13.8;
zi=tcs2+ql_*fai2*FqE/(2*pi*kc2)*log(dcs/dci);
e22=(zi-tci2)/zi;
tci2=zi%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu2=tci2+ql_*FqE*fai2*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d2_ku=ql_*FqE*fai2/(4*pi*100);
tu2_ku=(30.93-26.42)/(500-400)*(tu2-400)+26.42;
to2_ku=tu2_ku+d2_ku;
to2=(1000-900)/(48.06-45.14)*(to2_ku-45.14)+900%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;
h2=1341.5988e+3;x2=(h2-hfs)/(hgs-hfs);%该点含汽量
qDNB2=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x2))*...
((0.1484-1.596*x2+0.1729*x2*abs(x2))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x2)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR2=qDNB2/(q_*FqE*fai2)%计算烧毁比
%第三控制体温度计算
fai3=1.20;L3=4.2672/6;e31=0.01;tf3=320;
whilee31>0.001
t31_=0.5*(tf3+tf2);
Cp3_=1000*(0.04006*(t31_-310)+5.7437);
x3i=tf2+q_*FqN*FDHE*FqE*pi*dcs*L3*fai3/(Wu*Cp3_);
e31=(x3i-tf3)/tf3;
tf3=x3i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);
u3=800e-7;Pr3=0.999;k3=512e-3;%查得该温度下的热物性
Re3=Wu*De/(Ab*u3);
h3=0.023*Re3^0.8*Pr3^0.4*k3/De;%该处的对流换热系数
dtf31=q_*fai3*FqE/h3;%单相强迫对流放热公式算得的温压
ts=347.328;P=15.5;
dtf32=25*(q_*fai3*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf3;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf31tcs3=tf3+dtf31
else
tcs3=tf3+dtf32
end
dci=8.60e-3;tci3=349;e32=0.01;
whilee32>0.001
t32_=0.5*(tci3+tcs3);
kc3=0.0547*(1.8*t32_+32)+13.8;
ai=tcs3+ql_*fai3*FqE/(2*pi*kc3)*log(dcs/dci);
e32=(ai-tci3)/ai;
tci3=ai%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu3=tci3+ql_*FqE*fai3*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d3_ku=ql_*FqE*fai3/(4*pi*100);
tu3_ku=(34.97-30.93)/(600-500)*(tu3-500)+30.93;
to3_ku=tu3_ku+d3_ku;
to3=(1560-1405)/(61.95-58.4)*(to3_ku-58.4)+1405%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;
G=pf_*v*3600;
h3=1416.5e+3;x3=(h3-hfs)/(hgs-hfs);%该点含汽量
qDNB3=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x3))*...
((0.1484-1.596*x3+0.1729*x3*abs(x3))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x3)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR3=qDNB3/(q_*FqE*fai3)%计算烧毁比
%第四控制体
fai4=1.20;L4=4.2672/6;e41=0.01;tf4=330;
whilee41>0.001
t41_=0.5*(tf4+tf3);
Cp4_=1000*(0.04006*(t41_-310)+5.7437);
x4i=tf3+q_*FqN*FDHE*FqE*pi*dcs*L4*fai4/(Wu*Cp4_);
e41=(x4i-tf4)/tf4;
tf4=x4i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);
u4=869e-7;Pr4=1.01;k4=533e-3;%查得该温度下的热物性
Re4=Wu*De/(Ab*u4);
h4=0.023*Re4^0.8*Pr4^0.4*k4/De;%该处的对流换热系数
dtf41=q_*fai4*FqE/h4;%单相强迫对流放热公式算得的温压
ts=347.328;P=15.5;
dtf42=25*(q_*fai4*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf4;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf41tcs4=tf4+dtf41
else
tcs4=tf4+dtf42
end
dci=8.93e-3;tci4=349;e42=0.01;
whilee42>0.001
t42_=0.5*(tci4+tcs4);
kc4=0.0547*(1.8*t42_+32)+13.8;
mi=tcs4+ql_*fai4*FqE/(2*pi*kc4)*log(dcs/dci);
e42=(mi-tci4)/mi;
tci4=mi%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu4=tci4+ql_*FqE*fai4*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d4_ku=ql_*FqE*fai4/(4*pi*100);
tu4_ku=(34.97-30.93)/(600-500)*(tu3-500)+30.93;
to4_ku=tu4_ku+d4_ku;
to4=(1560-1405)/(61.95-58.4)*(to4_ku-58.4)+1405%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;
G=pf_*v*3600;
h4=1416.5e+3;x4=(h4-hfs)/(hgs-hfs);%该点含汽量
qDNB4=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x4))*...
((0.1484-1.596*x4+0.1729*x4*abs(x4))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x4)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR4=qDNB4/(q_*FqE*fai4)%计算烧毁比
%第五控制体
fai5=0.80;L5=4.2672/6;e51=0.01;tf5=350;
whilee51>0.001
t51_=0.5*(tf5+tf4);
Cp5_=1000*(0.04006*(t51_-310)+5.7437);
x5i=tf4+q_*FqN*FDHE*FqE*pi*dcs*L5*fai5/(Wu*Cp5_);
e51=(x5i-tf5)/tf5;
tf5=x5i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);
u5=650e-7;Pr5=1.47;k5=447e-3;%查得该温度下的热物性
Re5=Wu*De/(Ab*u5);
h5=0.023*Re5^0.8*Pr5^0.4*k5/De;%该处的对流换热系数
dtf51=q_*fai5*FqE/h5;%单相强迫对流放热公式算得的温压
ts=347.328;P=15.5;
dtf52=25*(q_*fai5*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf5;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf51tcs5=tf5+dtf51
else
tcs5=tf5+dtf52
end
dci=8.93e-3;tci5=350;e52=0.01;
whilee52>0.001
t52_=0.5*(tci5+tcs5);
kc5=0.0547*(1.8*t52_+32)+13.8;
ni=tcs5+ql_*fai5*FqE/(2*pi*kc5)*log(dcs/dci);
e52=(ni-tci5)/ni;
tci5=ni%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu5=tci5+ql_*FqE*fai5*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d5_ku=ql_*FqE*fai5/(4*pi*100);
tu5_ku=(34.97-30.93)/(600-500)*(tu4-500)+30.93;
to5_ku=tu5_ku+d5_ku;
to5=(1560-1405)/(61.95-58.4)*(to5_ku-58.4)+1405%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;
G=pf_*v*3600;
x5=(h5-hfs)/(hgs-hfs);%该点含汽量
qDNB5=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x5))*...
((0.1484-1.596*x5+0.1729*x5*abs(x5))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x5)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR5=qDNB5/(q_*FqE*fai5)%计算烧毁比
%第六控制体
fai6=0.40;L6=4.2672/6;e61=0.01;tf6=360;
whilee61>0.001
t61_=0.5*(tf6+tf5);
Cp6_=1000*(0.04006*(t61_-310)+5.7437);
x6i=tf5+q_*FqN*FDHE*FqE*pi*dcs*L6*fai6/(Wu*Cp6_);
e61=(x6i-tf6)/tf6;
tf6=x6i%求出该控制体出口处的温度
end
De=4*(S^2-pi/4*dcs^2)/(pi*dcs);
u6=600e-7;Pr6=1.263;k6=425e-3;%查得该温度下的热物性
Re6=Wu*De/(Ab*u6);
h6=0.023*Re6^0.8*Pr6^0.4*k6/De;%该处的对流换热系数
dtf61=q_*fai6*FqE/h6;%单相强迫对流放热公式算得的温压
ts=347.328;P=15.5;
dtf62=25*(q_*fai6*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf6;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压
ifdtf61tcs6=tf6+dtf61
else
tcs6=tf6+dtf62
end
dci=8.93e-3;tci6=350;e62=0.01;
whilee62>0.001
t62_=0.5*(tci6+tcs6);
kc6=0.0547*(1.8*t62_+32)+13.8;
ri=tcs6+ql_*fai6*FqE/(2*pi*kc6)*log(dcs/dci);
e62=(ri-tci6)/ri;
tci6=ri%采用迭代算法求得包壳内表面温度
end
hg=5678;du=8.19e-3;
tu6=tci6+ql_*FqE*fai6*2/(pi*(dci+du)*hg)%燃料芯块表面温度
d6_ku=ql_*FqE*fai6/(4*pi*100);
tu6_ku=(34.97-30.93)/(600-500)*(tu5-500)+30.93;
to6_ku=tu6_ku+d6_ku;
to6=(1560-1405)/(61.95-58.4)*(to6_ku-58.4)+1405%根据积分热导率图表查得芯块中心温度
p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;
G=pf_*v*3600;
x6=(h6-hfs)/(hgs-hfs);%该点含汽量
qDNB6=3.154e6*((2.022-6.238e-8*p)+...%根据W-3公式计算出临界热流量
(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x6))*...
((0.1484-1.596*x6+0.1729*x6*abs(x6))*0.2049*G/10^6+1.037)*...
(1.157-0.869*x6)*...
(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))
DNBR6=qDNB6/(q_*FqE*fai6)%计算烧毁比
%热管中的压降
uf=856e-7;L=4.2672;uw=817e-7;%uf为按主流平均温度取值的流体的粘性系数.uw为按照壁面温度取值的流体的粘性系数。
Re_=pf_*v*De/uf;
f=0.3146/Re_^0.25*(uw/uf)^0.6;
%摩擦压降
dPf=f*L*(G/3600)^2*vf_/(2*De)
%单相流体提升压降计算
g=9.8;Kout=1.0;Kin=0.75;Kgr=1.05;vfin=0.0013324;vfout=0.0018959;
dPel=pf_*g*L
%进口局部压降计算
dPin=Kin*(G/3600)^2*vfin/2
%出口局部压降计算
dPout=Kout*(G/3600)^2*vfout/2
%定位搁架出口压降计算
dPgr=Kgr*(G/3600)^2*(vfin+vfout)/2/2
%总的压降计算
dP=dPf+dPel+dPin+dPout+dPgr