反应堆热工分析计算程序.docx

反应堆热工分析计算程序.docx

《反应堆热工分析计算程序.docx》由会员分享，可在线阅读，更多相关《反应堆热工分析计算程序.docx（13页珍藏版）》请在冰豆网上搜索。

反应堆热工分析计算程序

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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf11

tcs1=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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf21

tcs2=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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf31

tcs3=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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf41

tcs4=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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf51

tcs5=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;%采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压

ifdtf61

tcs6=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