反应堆热工分析计算程序

tfin=279.4;Fa=0.974;Nt=3400000000;Wt=14314;b=0.059; tfout=340;e0=0.01;while e00.001t0_=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* (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=Wt*(1-b)*Ab/Af; %单元截面流量 %第一控制体温度计算e11=0.01;tf1=300;L1=4.2672/6;fai1=0.4;while e110.001t11_=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; %查得该温度下的热物性 Re1=Wu*De/(Ab*u1);hl=0.023*Re0.8 *Pr0.4 *kl/De; % 该处的对流换热系数 dtf11=q_*fai1*FqE/h1; %单相强迫对流放热公式算得的温压 ts=346.310791;P=15.51;dt fl2=25 *(q_* fail *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-tfl; % 采用詹斯- 洛特斯传热方程算得的过冷沸腾膜温压if dtfll0.00ltl2_=0.5*(tcil+tcsl); kcl=0.0547*(l.8*tl2_+32)+l3.8;yi=tcsl+ql_*fail*FqE/(2*pi*kcl)*log(dcs/dci); el2=(yi-tcil)/yi;tcil=yi %采用迭代算法求得包壳内表面温度endhg=5678;du=8.l9e-3;tul=tcil+ql_*FqE*fail*2/(pi*(dci+du)*hg) %燃料芯块表面温度 dl_ku=ql_*FqE*fail/(4*pi*l00);tul_ku=(26.42-2l.32)/(400-300)*(tul-300)+2l.32; tol_ku=tul_ku+dl_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*3 600;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*xl+0.1729*xl *abs(xl) *0.2049*G/10飞+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;while e210.001t21_=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 *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;dt f22=25*(q_* fai2 *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-t f2; % 采用詹斯-洛 特斯传热方程算得的过冷沸腾膜温压tcs2=tf2+dtf21elsetcs2=tf2+dtf22enddci=8.93e-3;tci2=349;e22=0.01;while e220.001t22_=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% 采用迭代算法求得包壳内表面温度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)*(to2_ku-45.14)+900 %根据积分热导率 图表查得芯块中心温度p=15.8e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3 600;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飞+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.001t31_=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% 求出该控制体出口处的温度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;dt f32=25*(q_* fai3 *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-tf3;% 采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf310.001t32_=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_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飞+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.001t41_=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%求出该控制体出口处的温度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;dt f42=25*(q_* fai4 *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-tf4;% 采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf410.001t42_=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-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飞+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.001t51_=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% 求出该控制体出口处的温度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 *Pr5“0.4 *k5/De;% 该处的对流换热系数dtf51=q_*fai5*FqE/h5;%单相强迫对流放热公式算得的温压ts=347.328;P=15.5;dt f52=25*(q_* fai5 *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-tf5;% 采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf510.001t52_=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_*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飞+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.001t61_=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%求出该控制体出口处的温度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=347.328;P=15.5;dt f62=25*(q_* fai6 *FqE/10飞厂0.25 *exp(-P/6.2)+ ts-tf6;% 采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf610.001t62_=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+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飞+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