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

1、反应堆热工分析计算程序tfin=279.4;Fa=0.974;Nt=3400000000;Wt=14314;b=0.059;tfout=340;e0=0.01;while e00.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) %燃料元件表面平均热流量F

2、qN=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*(S2-pi/4*dcs2)+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=S2-pi/4*dcs2; %单元流通截面积Wu=W

3、t*(1-b)*Ab/Af; %单元截面流量%第一控制体温度计算e11=0.01;tf1=300;L1=4.2672/6;fai1=0.4;while e110.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 %求出该控制体出口处的温度endDe=4*(S2-pi/4*dcs2)/(pi*dcs); %单元通道当量直径u1=944e-7;Pr1=0.85;k1=575.5e-3;

4、 %查得该温度下的热物性Re1=Wu*De/(Ab*u1);h1=0.023*Re10.8*Pr10.4*k1/De; %该处的对流换热系数dtf11=q_*fai1*FqE/h1; %单相强迫对流放热公式算得的温压ts=346.310791;P=15.51;dtf12=25*(q_*fai1*FqE/106)0.25*exp(-P/6.2)+ts-tf1; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf110.001 t12_=0.5*(tci1+tcs1); kc1=0.0547*(1.8*t12_+32)+13.8; yi=tcs1+ql_*fai1*FqE/(2*pi*kc

5、1)*log(dcs/dci); e12=(yi-tci1)/yi; tci1=yi %采用迭代算法求得包壳内表面温度endhg=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、+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/106+1.037)*. (1.157-0.869*x1)*. (0.2664+0.8357*e

7、xp(-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;while e210.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 %求出该控制体出口处的温度endDe=4*(S2-pi/4

8、*dcs2)/(pi*dcs);u2=857e-7;Pr2=0.90;k2=547e-3; %查得该温度下的热物性Re2=Wu*De/(Ab*u2);h2=0.023*Re20.8*Pr20.4*k2/De; %该处的对流换热系数dtf21=q_*fai2*FqE/h2; %单相强迫对流放热公式算得的温压ts=346.310791;P=15.5;dtf22=25*(q_*fai2*FqE/106)0.25*exp(-P/6.2)+ts-tf2; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf210.001 t22_=0.5*(tci2+tcs2); kc2=0.0547*(1.8*

9、t22_+32)+13.8; zi=tcs2+ql_*fai2*FqE/(2*pi*kc2)*log(dcs/dci); e22=(zi-tci2)/zi; tci2=zi %采用迭代算法求得包壳内表面温度endhg=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)

10、*(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/106

11、+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;while e310.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-

12、tf3)/tf3; tf3=x3i %求出该控制体出口处的温度endDe=4*(S2-pi/4*dcs2)/(pi*dcs);u3=800e-7;Pr3=0.999;k3=512e-3; %查得该温度下的热物性Re3=Wu*De/(Ab*u3);h3=0.023*Re30.8*Pr30.4*k3/De; %该处的对流换热系数dtf31=q_*fai3*FqE/h3; %单相强迫对流放热公式算得的温压ts=347.328;P=15.5;dtf32=25*(q_*fai3*FqE/106)0.25*exp(-P/6.2)+ts-tf3; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf3

13、10.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 %采用迭代算法求得包壳内表面温度endhg=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

14、_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.14

15、84-1.596*x3+0.1729*x3*abs(x3)*0.2049*G/106+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;while e410.001 t41_=0.5*(tf4+tf3); Cp4_=1000*(0.04006*(t41_-310)+5.7437); x4i=tf3+q_*FqN*FDHE*Fq

16、E*pi*dcs*L4*fai4/(Wu*Cp4_); e41=(x4i-tf4)/tf4; tf4=x4i %求出该控制体出口处的温度endDe=4*(S2-pi/4*dcs2)/(pi*dcs);u4=869e-7;Pr4=1.01;k4=533e-3; %查得该温度下的热物性Re4=Wu*De/(Ab*u4);h4=0.023*Re40.8*Pr40.4*k4/De; %该处的对流换热系数dtf41=q_*fai4*FqE/h4; %单相强迫对流放热公式算得的温压ts=347.328;P=15.5;dtf42=25*(q_*fai4*FqE/106)0.25*exp(-P/6.2)+ts

17、-tf4; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf410.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 %采用迭代算法求得包壳内表面温度endhg=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-

18、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)

19、*exp(18.177- 5.987e-7*p)*x4)*. (0.1484-1.596*x4+0.1729*x4*abs(x4)*0.2049*G/106+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;while e510.001 t51_=0.5*(tf5+tf4); Cp5_=1000*(0.04006*(t51_-

20、310)+5.7437); x5i=tf4+q_*FqN*FDHE*FqE*pi*dcs*L5*fai5/(Wu*Cp5_); e51=(x5i-tf5)/tf5; tf5=x5i %求出该控制体出口处的温度endDe=4*(S2-pi/4*dcs2)/(pi*dcs);u5=650e-7;Pr5=1.47;k5=447e-3; %查得该温度下的热物性Re5=Wu*De/(Ab*u5);h5=0.023*Re50.8*Pr50.4*k5/De; %该处的对流换热系数dtf51=q_*fai5*FqE/h5; %单相强迫对流放热公式算得的温压ts=347.328;P=15.5;dtf52=25*

21、(q_*fai5*FqE/106)0.25*exp(-P/6.2)+ts-tf5; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf510.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 %采用迭代算法求得包壳内表面温度endhg=5678;du=8.19e-3;tu5=tci5+ql_*FqE*fai5*2/(pi*(dci+du)*hg) %燃料芯块表面温度d5_ku=ql

22、_*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公式计算出临

23、界热流量 (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/106+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;while e610.001 t61_=0.5*(tf6+tf5); Cp

24、6_=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 %求出该控制体出口处的温度endDe=4*(S2-pi/4*dcs2)/(pi*dcs);u6=600e-7;Pr6=1.263;k6=425e-3; %查得该温度下的热物性Re6=Wu*De/(Ab*u6);h6=0.023*Re60.8*Pr60.4*k6/De; %该处的对流换热系数dtf61=q_*fai6*FqE/h6; %单相强迫对流放热公式算得的温压ts=

25、347.328;P=15.5;dtf62=25*(q_*fai6*FqE/106)0.25*exp(-P/6.2)+ts-tf6; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf610.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 %采用迭代算法求得包壳内表面温度endhg=5678;du=8.19e-3;tu6=tci6+ql_*FqE*fai6*2/(pi*(dci+d

26、u)*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.

27、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/106+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为按主流平均温度取值的流体的粘性系数

28、.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