反应堆热工分析计算程序

%另一种带入方式%P=15.8;Nt=1.82*10^9;W=8916.66667;t_fin=287;L=3.66 ;m=121;n=265;n0=17;d_cs=0.0095 ;d_ci=0.00836 ; %%有问题d_u=0.00819 ; %%有问题s=0.0126 ;kexi=0.05;Fu=0.974;F_rn=1.35;F_zn=1.528;F_ln=1.11;F_h5e=0.95;F_qe=1.03;F_he=1.085;K_in=0.75;K_out=1.0;K_gr=1.05;f1=0.48;f2=1.02;f3=1.5;f4=1.56;f5=0.96;f6=0.48;hg=5678;%t_fout=input('请任意输入一堆芯出口温度（T=325）：t_fout='); t_fout=325;r=1;kk=0;while r>=0.0001kk=kk+1;t_p=0.5*(t_fout+t_fin);if (280<t_p<=300)cp=19.35*t_p-350;elseif (300<t_p<=320)cp=34.43*t_p-4874;elseif (320<t_p<=340)cp=84.7*t_p-20960.4;endt1=t_fin+Fu*Nt/(W*(1-kexi)*cp);r=(t1-t_fout)/t_fout;t_fout=t1;endt_p=0.5*(t_fout+t_fin);% 计算燃料棒表面平均热流密度fz=m*n*pi*d_cs*L ; %计算堆心燃料棒的总传热面积q=Fu*Nt/fz ; %燃料棒表面平均热流密度F_qn=F_rn*F_zn*F_ln; % 有问题！！q_max=q*F_qn*F_qe;ql=q*pi*d_cs ;ql_max=ql*F_qn*F_qe;%计算平均管的流速dert=0.0008;Af=m*n*(s^2-pi*d_cs^2/4)+4*n0*s*dert*m/2; %堆心内流通面积if(300<t_p<=310)v=4.25*10^(-6)*t_p+1.0002*10^(-4);elseif(310<t_p<=320)v=5.112*10^(-6)*t_p-1.672*10^(-4);elseif(320<t_p<=330)v=6.444*10^(-6)*t_p-5.9344*10^(-4);elseif(330<t_p<=340)v=9.0976*10^(-5)*t_p-1.469128*10^(-3);endV=W*(1-kexi)*v/(Af);%fprintf('V=%.4f\n',V);%计算热管流量Ab=s^2-pi*d_cs^2/4;Wh=W*(1-kexi)*Ab/Af;%fprintf('Wh=%.4f\n', Wh);%热管中的计算De=4*(s^2-pi*d_cs^2/4)/(pi*d_cs);%t_fh1=input('请任意输入第一控制体热管流体温度（T=291）：t_fh1=');t_fh1=291;r1=1;kk1=0;while r1>=0.0001kk1=kk1+1;t_fhp1=0.5*(t_fh1+t_fin);if (t_fhp1<=300)cp1=19.35*t_fhp1-350;elseif (t_fhp1<=320)cp1=34.43*t_fhp1-4874;elseif (t_fhp1<=340)cp1=84.7*t_fhp1-20960.4;endt11=t_fin+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp1*6)*f1;r1=(t11-t_fh1)/t_fh1;t_fh1=t11;endt_fhp1=0.5*(t_fh1+t_fin);%fprintf('t_fh1=%.4f\n',t_fh1);ts=346.19;%计算第一控制体出口处的包壳外壁温if(t_fhp1<=300)u1=(-0.40072*t_fhp1+209.504)/10^6;k1=(-1.42504*t_fhp1+994.816)/1000;pr1=(-0.40072*t_fhp1+209.504)*(0.01935*t_fhp1-0.35)/(-1.42504*t_fhp1+994.816); elseif(t_fhp1<=320)u1=(-1.044416*t_fhp1+402.6128)/10^6;k1=(-6.916*t_fhp1+2642.104)/1000;pr1=(-1.044416*t_fhp1+402.6128)*(0.03443*t_fhp1-4.874)/(-6.916*t_fhp1+2642.104); elseif(t_fhp1<=340)u1=(-1.044416*t_fhp1+402.6128)/10^6;k1=(-6.916*t_fhp1+2642.104)/1000;pr1=(-1.044416*t_fhp1+402.6128)*(0.0847*t_fhp1-20.9604)/(-6.916*t_fhp1+2642.104); endre1=Wh*De/(Ab*u1) ;h1=0.023*re1^0.8*pr1^0.4*k1/De ;t_csh11=t_fh1+q*F_rn*f1*F_qe/h1;t_csh12=ts+25*(q*F_rn*f1*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh11<=t_csh12)t_csh1=t_csh11;elseif(t_csh11>=t_csh12)t_csh1=t_csh12;end%fprintf('t_csh1=%.3f\n',t_csh1);%t_csh1=t_csh1;%计算包壳内壁温度%t_cih1=input('请任意输入一包壳内壁温度（T=315）:t_cih1='); t_cih1=315;r12=1;kk11=0;while r12>=0.0001kk11=kk11+1;t_cp1=0.5*(t_cih1+t_csh1);k12=0.00547*(1.8*t_cp1+32)+13.8;t12=t_csh1+ql*F_rn*F_qe*f1*log(d_cs/d_ci)/(2*pi*k12);r12=(t12-t_cih1)/t_cih1;t_cih1=t12;end%fprintf('t_cih1=%.4f\n',t_cih1);%计算芯块表面温度t_uh1=t_cih1+ql*F_rn*F_qe*f1/(pi*0.5*(d_ci+d_u)*hg);%fprintf('t_uh1=%.4f\n',t_uh1);%计算芯块中心温度if(t_uh1<=400) %300 400ku1=0.051*t_uh1+6.02;elseif(t_uh1<=500) %400 500ku1=0.0451*t_uh1+8.38;elseif(t_uh1<=600) %500 600ku1=0.0404*t_uh1+10.73elseif(t_uh1<=700) %600 700ku1=0.0368*t_uh1+12.89;elseif(t_uh1<=800) %700 800ku1=0.0337*t_uh1+15.06;elseif(t_uh1<=900) %800 900ku1=0.0312*t_uh1+17.06;elseif(t_uh1<=1000) %900 1000ku1=0.0292*t_uh1+18.86;elseif(t_uh1<=1100) %1000 1100ku1=0.0255*t_uh1+22.56;elseif(t_uh1<=1200) %1100 1200ku1=0.028*t_uh1+19.81;endkuto1=ku1+ql*F_rn*F_qe*f1/(4*pi*100) ;if(kuto1<=34.93) %500 600to1=22.173*kuto1-185.8 ;elseif(kuto1<=38.65) %600 700to1=26.74*kuto1-335.02;elseif(kuto1<=42.02) %700 800to1=27.17*kuto1-446.88;elseif(kuto1<=45.14) %800 900to1=32.05*kuto1-546.79;elseif(kuto1<=48.06) %900 1000to1=34.25*kuto1-645.89;elseif(kuto1<=50.61) %1000 1100to1=39.22*kuto1-884.71;elseif(kuto1<=53.41)% 1100 1200to1=35.71*kuto1-707.5;elseif(kuto1<=55.84) %1200 1298to1=40.32*kuto1-953.63;elseif(kuto1<=58.4) %1298 1405to1=41.84*kuto1-1037;elseif(kuto1<=61.95) %1405 1560to1=43.668*kuto1-1145;elseif(kuto1<=66.87) %1560 1738to1=36.23*kuto1-682.25;elseif(kuto1<=68.86) %1378 1876to1=69.35*kuto1-2899.2;elseif(kuto1<=71.31) %1876 1990to1=46.53*kuto1-1328.1;elseif(kuto1<=74.88) %1990 2155to1=46.22*kuto1-1305.8;elseif(kuto1<=79.16) %2155 2343to1=43.93*kuto1-1134;end%fprintf('to1=%.3f\n',to1);%热管中的计算%t_fh2=input('请任意输入第二控制体热管流体温度（T=300）：t_fh2='); t_fh2=300;r2=1;kk2=0;while r2>=0.0001kk2=kk2+1;t_fhp2=0.5*(t_fh2+t_fh1);if (t_fhp2<=300)cp2=19.35*t_fhp2-350;elseif (t_fhp2<=320)cp2=34.43*t_fhp2-4874;elseif (t_fhp2<=340)cp2=84.7*t_fhp2-20960.4;endt22=t_fh1+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp2*6)*f2;r2=(t22-t_fh2)/t_fh2;t_fh2=t22;endt_fhp2=0.5*(t_fh2+t_fh1);%fprintf('t_fh2=%.4f\n', t_fh2);%计算第二控制体出口处的包壳外壁温if(t_fhp2<=300)u2=(-0.40072*t_fhp2+209.504)/10^6;k2=(-1.42504*t_fhp2+994.816)/1000;pr2=(-0.40072*t_fhp2+209.504)*(0.01935*t_fhp2-0.35)/(-1.42504*t_fhp2+994.816); elseif(t_fhp2<=320)u2=(-1.044416*t_fhp2+402.6128)/10^6;k2=(-6.916*t_fhp2+2642.104)/1000;pr2=(-1.044416*t_fhp2+402.6128)*(0.03443*t_fhp2-4.874)/(-6.916*t_fhp2+2642.104); elseif(t_fhp2<=340)u2=(-1.044416*t_fhp2+402.6128)/10^6;k2=(-6.916*t_fhp2+2642.104)/1000;pr2=(-1.044416*t_fhp2+402.6128)*(0.0847*t_fhp2-20.9604)/(-6.916*t_fhp2+2642.104); endre2=Wh*De/(Ab*u2);h2=0.023*re2^0.8*pr2^0.4*k2/De;t_csh21=t_fh2+q*F_rn*f2*F_qe/h2;t_csh22=ts+25*(q*F_rn*f2*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh21<=t_csh22)t_csh2=t_csh21;elseif(t_csh21>=t_csh22)t_csh2=t_csh22;end%fprintf('t_csh2=%.4f\n',t_csh2);%计算第二控制体包壳内壁温度%t_cih2=input('请任意输入第二控制体包壳内壁温度（例如T=350）:t_cih1='); t_cih2=350;r22=1;kk22=0;while r22>=0.0001kk22=kk22+1;t_cp2=0.5*(t_cih2+t_csh2);k22=0.00547*(1.8*t_cp2+32)+13.8;t22=t_csh2+ql*F_rn*F_qe*f2*log(d_cs/d_ci)/(2*pi*k22);r22=(t22-t_cih2)/t_cih2;t_cih2=t22;end%fprintf('t_cih2=%.4f\n',t_cih2);%计算芯块表面温度t_uh2=t_cih2+ql*F_rn*F_qe*f2/(pi*0.5*(d_ci+d_u)*hg) ;;%fprintf('t_uh2=%.4f\n',t_uh2);%计算芯块中心温度if(t_uh2<=400) %300 400ku2=0.051*t_uh2+6.02;elseif(t_uh2<=500) %400 500ku2=0.0451*t_uh2+8.38;elseif(t_uh2<=600) %500 600ku2=0.0404*t_uh2+10.73elseif(t_uh2<=700) %600 700ku2=0.0368*t_uh2+12.89;elseif(t_uh2<=800) %700 800ku2=0.0337*t_uh2+15.06;elseif(t_uh2<=900) %800 900ku2=0.0312*t_uh2+17.06;elseif(t_uh2<=1000) %900 1000ku2=0.0292*t_uh2+18.86;elseif(t_uh2<=1100) %1000 1100ku2=0.0255*t_uh2+22.56;elseif(t_uh2<=1200) %1100 1200ku2=0.028*t_uh2+19.81;endkuto2=ku2+ql*F_rn*F_qe*f2/(4*pi*100) ;if(kuto2<=38.65) %600 700to2=26.74*kuto2-335.02;elseif(kuto2<=42.02) %700 800to2=27.17*kuto2-446.88;elseif(kuto2<=45.14) %800 900to2=32.05*kuto2-546.79;elseif(kuto2<=48.06) %900 1000to2=34.25*kuto2-645.89;elseif(kuto2<=50.61) %1000 1100to2=39.22*kuto2-884.71;elseif(kuto2<=53.41)% 1100 1200to2=35.71*kuto2-707.5;elseif(kuto2<=55.84) %1200 1298to2=40.32*kuto2-953.63;elseif(kuto2<=58.4) %1298 1405to2=41.84*kuto2-1037;elseif(kuto2<=61.95) %1405 1560to2=43.668*kuto2-1145;elseif(kuto2<=66.87) %1560 1738to2=36.23*kuto2-682.25;elseif(kuto2<=68.86) %1378 1876to2=69.35*kuto2-2899.2;elseif(kuto2<=71.31) %1876 1990to2=46.53*kuto2-1328.1;elseif(kuto2<=74.88) %1990 2155to2=46.22*kuto2-1305.8;elseif(kuto2<=79.16) %2155 2343to2=43.93*kuto2-1134;end%fprintf('to2=%.4f\n',to2);%热管中的计算%t_fh3=input('请任意输入第三控制体热管流体温度（T=314）：t_fh3='); t_fh3=314;r3=1;kk3=0;while r3>=0.0001kk3=kk3+1;t_fhp3=0.5*(t_fh3+t_fh2);if (t_fhp3<=300)cp3=19.35*t_fhp3-350;elseif (t_fhp3<=320)cp3=34.43*t_fhp3-4874;elseif (t_fhp3<=340)cp3=84.7*t_fhp3-20960.4;endt33=t_fh2+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp3*6)*f3;r3=(t33-t_fh3)/t_fh3;t_fh3=t33;endt_fhp3=0.5*(t_fh3+t_fh2);%fprintf('t_fh3=%.3f\n', t_fh3);%计算第三控制体出口处的包壳外壁温if(t_fhp3<=300)u3=(-0.40072*t_fhp3+209.504)/10^6;k3=(-1.42504*t_fhp3+994.816)/1000;pr3=(-0.40072*t_fhp3+209.504)*(0.01935*t_fhp3-0.35)/(-1.42504*t_fhp3+994.816); elseif(t_fhp3<=320)u3=(-1.044416*t_fhp3+402.6128)/10^6;k3=(-6.916*t_fhp3+2642.104)/1000;pr3=(-1.044416*t_fhp3+402.6128)*(0.03443*t_fhp3-4.874)/(-6.916*t_fhp3+2642.104); elseif(t_fhp3<=340)u3=(-1.044416*t_fhp3+402.6128)/10^6;k3=(-6.916*t_fhp3+2642.104)/1000;pr3=(-1.044416*t_fhp3+402.6128)*(0.0847*t_fhp3-20.9604)/(-6.916*t_fhp3+2642.104); endre3=Wh*De/(Ab*u3);h3=0.023*re3^0.8*pr3^0.4*k3/De;t_csh31=t_fh3+q*F_rn*f3*F_qe/h3;t_csh32=ts+25*(q*F_rn*f3*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh31<=t_csh32)t_csh3=t_csh31;elseif(t_csh31>=t_csh32)t_csh3=t_csh32;end%fprintf('t_csh3=%.3f\n',t_csh3);%计算第三控制体包壳内壁温度%t_cih3=input('请任意输入第三控制体包壳内壁温度（例如T=385）:t_cih3=');t_cih3=385;r33=1;kk33=0;while r33>=0.0001kk33=kk33+1;t_cp3=0.5*(t_cih3+t_csh3);k33=0.00547*(1.8*t_cp3+32)+13.8;t33=t_csh3+ql*F_rn*F_qe*f3*log(d_cs/d_ci)/(2*pi*k33); r33=(t33-t_cih3)/t_cih3;t_cih3=t33;end%fprintf('t_cih3=%.3f\n',t_cih3);%计算芯块表面温度t_uh3=t_cih3+ql*F_rn*F_qe*f3/(pi*0.5*(d_ci+d_u)*hg); %fprintf('t_uh3=%.3f\n',t_uh3);%计算芯块中心温度if(t_uh3<=400) %300 400ku3=0.051*t_uh3+6.02;elseif(t_uh3<=500) %400 500ku3=0.0451*t_uh3+8.38;elseif(t_uh3<=600) %500 600ku3=0.0404*t_uh3+10.73 ;elseif(t_uh3<=700) %600 700ku3=0.0368*t_uh3+12.89;elseif(t_uh3<=800) %700 800ku3=0.0337*t_uh3+15.06;elseif(t_uh3<=900) %800 900ku3=0.0312*t_uh3+17.06;elseif(t_uh3<=1000) %900 1000ku3=0.0292*t_uh3+18.86;elseif(t_uh3<=1100) %1000 1100ku3=0.0255*t_uh3+22.56;elseif(t_uh3<=1200) %1100 1200ku3=0.028*t_uh3+19.81;endkuto3=ku3+ql*F_rn*F_qe*f3/(4*pi*100) ;if kuto3<=38.65 %600 700to3=26.74*kuto3-335.02;elseif kuto3<=42.02 %700 800to3=27.17*kuto3-446.88;elseif(kuto3<=45.14) %800 900to3=32.05*kuto3-546.79;elseif kuto3<=48.06 %900 1000to3=34.25*kuto3-645.89;elseif kuto3<=50.61 %1000 1100to3=39.22*kuto3-884.71;elseif kuto3<=53.41 % 1100 1200to3=35.71*kuto3-707.5;elseif kuto3<=55.84 %1200 1298to3=40.32*kuto3-953.63;elseif kuto3<=58.4 %1298 1405to3=41.84*kuto3-1037;elseif kuto3<=61.95 %1405 1560to3=43.668*kuto3-1145;elseif kuto3<=66.87 %1560 1738to3=36.23*kuto3-682.25;elseif kuto3<=68.86 %1378 1876to3=69.35*kuto3-2899.2;elseif kuto3<=71.31 %1876 1990to3=46.53*kuto3-1328.1;elseif kuto3<=74.88 %1990 2155to3=46.22*kuto3-1305.8;elseif kuto3<=79.16 %2155 2343to3=43.93*kuto3-1134;end%fprintf('to3=%.3f\n',to3);%热管中的计算De=4*(s^2-pi*d_cs^2/4)/(pi*d_cs);%t_fh4=input('请任意输入第四控制体热管流体温度（T=326）：t_fh4='); t_fh4=326;r4=1;kk4=0;while r4>=0.0001kk4=kk4+1;t_fhp4=0.5*(t_fh4+t_fh3);if t_fhp4<=300cp4=19.35*t_fhp4-350;elseif t_fhp4<=320cp4=34.43*t_fhp4-4874;elseif t_fhp4<=340cp4=84.7*t_fhp4-20960.4;endt44=t_fh3+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp4*6)*f4;r4=(t44-t_fh4)/t_fh4;t_fh4=t44;endt_fhp4=0.5*(t_fh4+t_fh3);%fprintf('t_fh4=%.3f\n', t_fh4);ts=346.19;%计算第四控制体出口处的包壳外壁温if(t_fhp4<=300)u4=(-0.40072*t_fhp4+209.504)/10^6;k4=(-1.42504*t_fhp4+994.816)/1000;pr4=(-0.40072*t_fhp4+209.504)*(0.01935*t_fhp4-0.35)/(-1.42504*t_fhp4+994.816); elseif(t_fhp4<=320)u4=(-1.044416*t_fhp4+402.6128)/10^6;k4=(-6.916*t_fhp4+2642.104)/1000;pr4=(-1.044416*t_fhp4+402.6128)*(0.03443*t_fhp4-4.874)/(-6.916*t_fhp4+2642.104); elseif(t_fhp4<=340)u4=(-1.044416*t_fhp4+402.6128)/10^6;k4=(-6.916*t_fhp4+2642.104)/1000;pr4=(-1.044416*t_fhp4+402.6128)*(0.0847*t_fhp4-20.9604)/(-6.916*t_fhp4+2642.104); endre4=Wh*De/(Ab*u4);h4=0.023*re4^0.8*pr4^0.4*k4/De;t_csh41=t_fh4+q*F_rn*f4*F_qe/h4;t_csh42=ts+25*(q*F_rn*f4*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh41<=t_csh42)t_csh4=t_csh41;elseif(t_csh41>=t_csh42)t_csh4=t_csh42;end%fprintf('t_csh4=%.3f\n',t_csh4);%计算第四控制体包壳内壁温度%t_cih4=input('请任意输入第五控制体包壳内壁温度（T=386）:t_cih4=');t_cih4=386;r44=1;kk44=0;while r44>=0.0001kk44=kk44+1;t_cp4=0.5*(t_cih4+t_csh4);k44=0.00547*(1.8*t_cp4+32)+13.8;t44=t_csh4+ql*F_rn*F_qe*f4*log(d_cs/d_ci)/(2*pi*k44); r44=(t44-t_cih4)/t_cih4;t_cih4=t44;end%fprintf('t_cih4=%.4f\n',t_cih4);%计算芯块表面温度t_uh4=t_cih4+ql*F_rn*F_qe*f4/(pi*0.5*(d_ci+d_u)*hg); %fprintf('t_uh4=%.3f\n',t_uh4);%计算芯块中心温度if t_uh4<=400 %300 400ku4=0.051*t_uh4+6.02;elseif t_uh4<=500 %400 500ku4=0.0451*t_uh4+8.38;elseif t_uh4<=600 %500 600ku4=0.0404*t_uh4+10.73elseif t_uh4<=700 %600 700ku4=0.0368*t_uh4+12.89;elseif t_uh4<=800 %700 800ku4=0.0337*t_uh4+15.06;elseif t_uh4<=900 %800 900ku4=0.0312*t_uh4+17.06;elseif t_uh4<=1000 %900 1000ku4=0.0292*t_uh4+18.86;elseif t_uh4<=1100 %1000 1100ku4=0.0255*t_uh4+22.56;elseif t_uh4<=1200 %1100 1200ku4=0.028*t_uh4+19.81;endkuto4=ku4+ql*F_rn*F_qe*f4/(4*pi*100);if kuto4<=38.65 %600 700to4=26.74*kuto4-335.02;elseif kuto4<=42.02 %700 800to4=27.17*kuto4-446.88;elseif kuto4<=45.14 %800 900to4=32.05*kuto4-546.79;elseif kuto4<=48.06 %900 1000to4=34.25*kuto4-645.89;elseif kuto4<=50.61 %1000 1100to4=39.22*kuto4-884.71;elseif kuto4<=53.41 % 1100 1200to4=35.71*kuto4-707.5;elseif kuto4<=55.84 %1200 1298to4=40.32*kuto4-953.63;elseif kuto4<=58.4 %1298 1405to4=41.84*kuto4-1037;elseif kuto4<=61.95 %1405 1560to4=43.668*kuto4-1145;elseif kuto4<=66.87 %1560 1738to4=36.23*kuto4-682.25;elseif kuto4<=68.86 %1738 1876to4=69.35*kuto4-2899.2;elseif kuto4<=71.31 %1876 1990to4=46.53*kuto4-1328.1;elseif kuto4<=74.88 %1990 2155to4=46.22*kuto4-1305.8;elseif kuto4<=79.16 %2155 2343to4=43.93*kuto4-1134;end%fprintf('to4=%.3f\n',to4);%热管中的计算%t_fh5=input('请任意输入第五控制体热管流体温度（T=332）：t_fh5='); t_fh5=332;r5=1;kk5=0;while r5>=0.0001kk5=kk5+1;t_fhp5=0.5*(t_fh5+t_fh4);if t_fhp5<=300cp5=19.35*t_fhp5-350;elseif t_fhp5<=320cp5=34.43*t_fhp5-4874;elseif t_fhp5<=340cp5=84.7*t_fhp5-20960.4;endt55=t_fh4+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp5*6)*f5;r5=(t55-t_fh5)/t_fh5;endt_fhp5=0.5*(t_fh5+t_fh5);%fprintf('t_fh5=%.4f\n', t_fh5);%计算第五控制体出口处的包壳外壁温if(t_fhp5<=300)u5=(-0.40072*t_fhp5+209.504)/10^6;k5=(-1.42504*t_fhp5+994.816)/1000;pr5=(-0.40072*t_fhp5+209.504)*(0.01935*t_fhp5-0.35)/(-1.42504*t_fhp5+994.816); elseif(t_fhp5<=320)u5=(-1.044416*t_fhp5+402.6128)/10^6;k5=(-6.916*t_fhp5+2642.104)/1000;pr5=(-1.044416*t_fhp5+402.6128)*(0.03443*t_fhp5-4.874)/(-6.916*t_fhp5+2642.104); elseif(t_fhp5<=340)u5=(-1.044416*t_fhp5+402.6128)/10^6;k5=(-6.916*t_fhp5+2642.104)/1000;pr5=(-1.044416*t_fhp5+402.6128)*(0.0847*t_fhp5-20.9604)/(-6.916*t_fhp5+2642.104); endre5=Wh*De/(Ab*u5);h5=0.023*re5^0.8*pr5^0.4*k5/De;t_csh51=t_fh5+q*F_rn*f5*F_qe/h5;t_csh52=ts+25*(q*F_rn*f5*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh51<=t_csh52)t_csh5=t_csh51;elseif(t_csh51>=t_csh52)t_csh5=t_csh52;end%fprintf('t_csh5=%.4f\n',t_csh5);%计算第五控制体包壳内壁温度%t_cih5=input('请任意输入第五控制体包壳内壁温度（T=371）:t_cih5=');t_cih5=371;r55=1;kk55=0;while r55>=0.00001kk55=kk55+1;t_cp5=0.5*(t_cih5+t_csh5);k55=0.00547*(1.8*t_cp5+32)+13.8;t55=t_csh5+ql*F_rn*F_qe*f5*log(d_cs/d_ci)/(2*pi*k55);r55=(t55-t_cih5)/t_cih5;end%fprintf('t_cih5=%.4f\n',t_cih5);%计算芯块表面温度t_uh5=t_cih5+ql*F_rn*F_qe*f5/(pi*0.5*(d_ci+d_u)*hg); %fprintf('t_uh5=%.4f\n',t_uh5);%计算芯块中心温度if t_uh5<=400 %300 400ku5=0.051*t_uh5+6.02;elseif t_uh5<=500 %400 500ku5=0.0451*t_uh5+8.38;elseif t_uh5<=600 %500 600ku5=0.0404*t_uh5+10.73;elseif t_uh5<=700 %600 700ku5=0.0368*t_uh5+12.89;elseif t_uh5<=800 %700 800ku5=0.0337*t_uh5+15.06;elseif t_uh5<=900 %800 900ku5=0.0312*t_uh5+17.06;elseif t_uh5<=1000 %900 1000ku5=0.0292*t_uh5+18.86;elseif t_uh5<=1100 %1000 1100ku5=0.0255*t_uh5+22.56;elseif t_uh5<=1200 %1100 1200ku5=0.028*t_uh5+19.81;endkuto5=ku5+ql*F_rn*F_qe*f5/(4*pi*100);if kuto5<=38.65 %600 700to5=26.74*kuto5-335.02;elseif kuto5<=42.02 %700 800to5=27.17*kuto5-446.88;elseif kuto5<=45.14 %800 900to5=32.05*kuto5-546.79;elseif kuto5<=48.06 %900 1000to5=34.25*kuto5-645.89;elseif kuto5<=50.61 %1000 1100to5=39.22*kuto5-884.71;elseif kuto5<=53.41 % 1100 1200to5=35.71*kuto5-707.5;elseif kuto5<=55.84 %1200 1298to5=40.32*kuto5-953.63;elseif kuto5<=58.4 %1298 1405to5=41.84*kuto5-1037;elseif kuto5<=61.95 %1405 1560to5=43.668*kuto5-1145;elseif kuto5<=66.87 %1560 1738to5=36.23*kuto5-682.25;elseif kuto5<=68.86 %1738 1876to5=69.35*kuto5-2899.2;elseif kuto5<=71.31 %1876 1990to5=46.53*kuto5-1328.1;elseif kuto5<=74.88 %1990 2155to5=46.22*kuto5-1305.8;elseif kuto5<=79.16 %2155 2343to5=43.93*kuto5-1134;end%fprintf('to5=%.3f\n',to5);%热管中的计算%t_fh6=input('请任意输入第六控制体热管流体温度（T=336）：t_fh6='); t_fh6=336;r6=1;kk6=0;while r6>=0.0001kk6=kk6+1;t_fhp6=0.5*(t_fh6+t_fh5);if t_fhp6<=300cp6=19.35*t_fhp6-350;elseif t_fhp6<=320cp6=34.43*t_fhp6-4874;elseif t_fhp6<=340cp6=84.7*t_fhp6-20960.4;endt66=t_fh5+q*F_rn*F_he*F_h5e*pi*L*d_cs/(Wh*cp6*6)*f6;r6=(t66-t_fh6)/t_fh6;t_fh6=t66;endt_fhp6=0.5*(t_fh6+t_fh6);%fprintf('t_fh6=%.4f\n', t_fh6);%计算第六控制体出口处的包壳外壁温if(t_fhp6<=300)u6=(-0.40072*t_fhp6+209.504)/10^6;k6=(-1.42504*t_fhp6+994.816)/1000;pr6=(-0.40072*t_fhp6+209.504)*(0.01935*t_fhp6-0.35)/(-1.42504*t_fhp6+994.816); elseif(t_fhp6<=320)u6=(-1.044416*t_fhp6+402.6128)/10^6;k6=(-6.916*t_fhp6+2642.104)/1000;pr6=(-1.044416*t_fhp6+402.6128)*(0.03443*t_fhp6-4.874)/(-6.916*t_fhp6+2642.104); elseif(t_fhp6<=340)u6=(-1.044416*t_fhp6+402.6128)/10^6;k6=(-6.916*t_fhp6+2642.104)/1000;pr6=(-1.044416*t_fhp6+402.6128)*(0.0847*t_fhp6-20.9604)/(-6.916*t_fhp6+2642.104); endre6=Wh*De/(Ab*u6);h6=0.023*re6^0.8*pr6^0.4*k6/De;t_csh61=t_fh6+q*F_rn*f6*F_qe/h6;t_csh62=ts+25*(q*F_rn*f6*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh61<=t_csh62)t_csh6=t_csh61;elseif(t_csh61>=t_csh62)t_csh6=t_csh62;end%fprintf('t_csh6=%.4f\n',t_csh6);%计算第六控制体包壳内壁温度%t_cih6=input('请任意输入第六控制体包壳内壁温度（T=358）:t_cih6=');t_cih6=358;r66=1;kk66=0;while r66>=0.00001kk66=kk66+1;t_cp6=0.5*(t_cih6+t_csh6);k66=0.00547*(1.8*t_cp6+32)+13.8;t66=t_csh6+ql*F_rn*F_qe*f6*log(d_cs/d_ci)/(2*pi*k66);r66=(t66-t_cih6)/t_cih6;t_cih6=t66;end%fprintf('t_cih6=%.4f\n',t_cih6);%计算芯块表面温度t_uh6=t_cih6+ql*F_rn*F_qe*f6/(pi*0.5*(d_ci+d_u)*hg); %fprintf('t_uh6=%.4f\n',t_uh6);%计算芯块中心温度if t_uh6<=400 %300 400ku6=0.051*t_uh6+6.02;elseif t_uh6<=500 %400 500ku6=0.0451*t_uh6+8.38;elseif t_uh6<=600 %500 600ku6=0.0404*t_uh6+10.73elseif t_uh6<=700 %600 700ku6=0.0368*t_uh6+12.89;elseif t_uh6<=800 %700 800ku6=0.0337*t_uh6+15.06;elseif t_uh6<=900 %800 900ku6=0.0312*t_uh6+17.06;elseif t_uh6<=1000 %900 1000ku6=0.0292*t_uh6+18.86;elseif t_uh5<=1100 %1000 1100ku6=0.0255*t_uh6+22.56;elseif t_uh6<=1200 %1100 1200ku6=0.028*t_uh6+19.81;endkuto6=ku6+ql*F_rn*F_qe*f6/(4*pi*100);if kuto6<=38.65 %600 700to6=26.74*kuto6-335.02;elseif kuto6<=42.02 %700 800to6=27.17*kuto6-446.88;elseif kuto6<=45.14 %800 900to6=32.05*kuto6-546.79;elseif kuto6<=48.06 %900 1000to6=34.25*kuto6-645.89;elseif kuto6<=50.61 %1000 1100to6=39.22*kuto6-884.71;elseif kuto6<=53.41 % 1100 1200to6=35.71*kuto6-707.5;elseif kuto6<=55.84 %1200 1298to6=40.32*kuto6-953.63;elseif kuto6<=58.4 %1298 1405to6=41.84*kuto6-1037;elseif kuto6<=61.95 %1405 1560to6=43.668*kuto6-1145;elseif kuto6<=66.87 %1560 1738to6=36.23*kuto6-682.25;elseif kuto6<=68.86 %1738 1876to6=69.35*kuto6-2899.2;elseif kuto6<=71.31 %1876 1990to6=46.53*kuto6-1328.1;elseif kuto6<=74.88 %1990 2155to6=46.22*kuto6-1305.8;elseif kuto6<=79.16 %2155 2343to6=43.93*kuto6-1134;end%fprintf('to6=%.4f\n',to6);%热管中的计算f0=1/L;t_fh0=t_fin;ts=346.19;%计算第一控制体出口处的包壳外壁温if(t_fh0<=300)u0=(-0.40072*t_fh0+209.504)/10^6;k0=(-1.42504*t_fh0+994.816)/1000;pr0=(-0.40072*t_fh0+209.504)*(0.01935*t_fh0-0.35)/(-1.42504*t_fh0+994.816); elseif(t_fh0<=320)u0=(-1.044416*t_fh0+402.6128)/10^6;k0=(-6.916*t_fh0+2642.104)/1000;pr0=(-1.044416*t_fh0+402.6128)*(0.03443*t_fh0-4.874)/(-6.916*t_fh0+2642.104); elseif(t_fh0<=340)u0=(-1.044416*t_fh0+402.6128)/10^6;k0=(-6.916*t_fh0+2642.104)/1000;pr0=(-1.044416*t_fh0+402.6128)*(0.0847*t_fh0-20.9604)/(-6.916*t_fh0+2642.104); endre0=Wh*De/(Ab*u0) ;h0=0.023*re0^0.8*pr0^0.4*k0/De ;t_csh01=t_fh0+q*F_rn*f0*F_qe/h0;t_csh02=ts+25*(q*F_rn*f0*F_qe/10^6)^0.25*exp(-P/6.2);if(t_csh01<=t_csh02)t_csh0=t_csh01;elseif(t_csh01>=t_csh02)t_csh0=t_csh02;end%fprintf('t_csh0=%.3f\n',t_csh0);%计算包壳内壁温度t_cih0=300;r0=1;kk0=0;while r0>=0.0001kk0=kk0+1;t_cp0=0.5*(t_cih0+t_csh0);k02=0.00547*(1.8*t_cp0+32)+13.8;t02=t_csh0+ql*F_rn*F_qe*f0*log(d_cs/d_ci)/(2*pi*k02); r0=(t02-t_cih0)/t_cih0;t_cih0=t02;end%计算芯块表面温度t_uh0=t_cih0+ql*F_rn*F_qe*f0/(pi*0.5*(d_ci+d_u)*hg); %计算芯块中心温度if(t_uh0<=400) %300 400ku0=0.051*t_uh0+6.02;elseif(t_uh0<=500) %400 500ku0=0.0451*t_uh0+8.38;elseif(t_uh0<=600) %500 600ku0=0.0404*t_uh0+10.73elseif(t_uh0<=700) %600 700ku0=0.0368*t_uh0+12.89;elseif(t_uh0<=800) %700 800ku0=0.0337*t_uh0+15.06;elseif(t_uh0<=900) %800 900ku0=0.0312*t_uh0+17.06;elseif(t_uh0<=1000) %900 1000ku0=0.0292*t_uh0+18.86;elseif(t_uh0<=1100) %1000 1100ku0=0.0255*t_uh0+22.56;elseif(t_uh0<=1200) %1100 1200ku0=0.028*t_uh0+19.81;endkuto0=ku0+ql*F_rn*F_qe*f0/(4*pi*100) ;if(kuto0<=26.42) %300 400to0=19.61*kuto0-118.04;elseif(kuto0<=34.93) %500 600to0=22.173*kuto0-185.8 ;elseif(kuto0<=38.65) %600 700to0=26.74*kuto0-335.02;elseif(kuto0<=42.02) %700 800to0=27.17*kuto0-446.88;elseif(kuto0<=45.14) %800 900to0=32.05*kuto0-546.79;elseif(kuto0<=48.06) %900 1000to0=34.25*kuto0-645.89;elseif(kuto0<=50.61) %1000 1100to0=39.22*kuto0-884.71;elseif(kuto0<=53.41)% 1100 1200to0=35.71*kuto0-707.5;elseif(kuto0<=55.84) %1200 1298to0=40.32*kuto0-953.63;elseif(kuto0<=58.4) %1298 1405to0=41.84*kuto0-1037;elseif(kuto0<=61.95) %1405 1560to0=43.668*kuto0-1145;elseif(kuto0<=66.87) %1560 1738to0=36.23*kuto0-682.25;elseif(kuto0<=68.86) %1378 1876to0=69.35*kuto0-2899.2;elseif(kuto0<=71.31) %1876 1990to0=46.53*kuto0-1328.1;elseif(kuto0<=74.88) %1990 2155to0=46.22*kuto0-1305.8;elseif(kuto0<=79.16) %2155 2343to0=43.93*kuto0-1134;end%fprintf('to0=%.3f\n',to0);fprintf('q=%.4f\nq_max=%.4f\nql=%.4f\nql_max=%.4f\n ',q,q_max,ql,ql_max); fprintf('t_fout=%.4f\n',t_fout)fprintf('V=%.4f\n',V);fprintf('Wh=%.4f\n', Wh);fprintf('t_fh0=%.4f\nt_csh0=%.4f\nt_cih0=%.4ft_uh0=%.4f\n5_to0=%.4f\n',t_fh0,t_csh0,t_cih0,t_ uh0,to0);fprintf('t_fh1=%.4f\n',t_fh1);fprintf('t_csh1=%.4f\n',t_csh1);fprintf('t_cih1=%.4f\n',t_cih1);fprintf('t_uh1=%.4f\n',t_uh1);fprintf('to1=%.4f\n',to1);fprintf('t_fh2=%.4f\n',t_fh2);fprintf('t_csh2=%.4f\n',t_csh2);fprintf('t_cih2=%.4f\n',t_cih2);fprintf('t_uh2=%.4f\n',t_uh2);fprintf('to2=%.4f\n',to2);fprintf('t_fh3=%.4f\n',t_fh3);fprintf('t_csh3=%.4f\n',t_csh3);fprintf('t_cih3=%.4f\n',t_cih3);fprintf('t_uh3=%.4f\n',t_uh3);fprintf('to3=%.4f\n',to3);fprintf('t_fh4=%.4f\n',t_fh4);fprintf('t_csh4=%.4f\n',t_csh4);fprintf('t_cih4=%.4f\n',t_cih4);fprintf('t_uh4=%.4f\n',t_uh4);fprintf('to4=%.4f\n',to4);fprintf('t_fh5=%.4f\n',t_fh5);fprintf('t_csh5=%.4f\n',t_csh5);fprintf('t_cih5=%.4f\n',t_cih5);fprintf('t_uh5=%.4f\n',t_uh5);fprintf('to5=%.3f\n',to5);fprintf('t_fh6=%.4f\n',t_fh6);fprintf('t_csh6=%.4f\n',t_csh6);fprintf('t_cih6=%.4f\n',t_cih6);fprintf('t_uh6=%.4f\n',t_uh6);fprintf('to6=%.4f\n',to6);Hfs=1642.8*10^3;Hgs=2589.98*10^3;Hfg=Hgs-Hfs;G=W*(1-kexi)/Af;g=9.8;if (t_fin<=290)H0=(5.146*t_fin-208.22)*10^3;elseif (t_fin<=300)elseif (t_fin<=310)H0=(5.59*t_fin-339.44)*10^3;elseif (t_fin<=320)H0=(5.932*t_fin-445.46)*10^3;elseif (t_fin<=330)H0=(6.422*t_fin-602.26)*10^3;elseif (t_fin<=340)H0=(7.23*t_fin-868.9)*10^3;elseif (t_fin<=350)H0=(104.63*t_fin-33984.9)*10^3;endXe0=(H0-Hfs)/Hfg;q_DNB0=3.154*10^6*((2.022-6.238*10^(-8)*P)+(0.1722-1.43*10^(-8)*P)*exp((18.177-5.987*1 0^(-7)*P)*Xe0))*((0.1484-1.596*Xe0+0.1729*Xe0*abs(Xe0))*0.2049*G/10^6+1.037)*(1.157-0. 869*Xe0)*(0.2664+0.8357*exp(-124*De))*(0.8258+0.341*10^(-6)*(Hfs-H0));DNBR0=q_DNB0/(q*F_rn*F_qe*f0);%H1if (t_fh1<=290)H1=(5.146*t_fh1-208.22)*10^3;elseif (t_fh1<=300)H1=(5.344*t_fh1-265.64)*10^3;elseif (t_fh1<=310)H1=(5.59*t_fh1-339.44)*10^3;elseif (t_fh1<=320)H1=(5.932*t_fh1-445.46)*10^3;elseif (t_fh1<=330)H1=(6.422*t_fh1-602.26)*10^3;elseif (t_fh1<=340)H1=(7.23*t_fh1-868.9)*10^3;elseif (t_fh1<=350)H1=(104.63*t_fh1-33984.9)*10^3;endXe1=(H1-Hfs)/Hfg;q_DNB1=3.154*10^6*((2.022-6.238*10^(-8)*P)+(0.1722-1.43*10^(-8)*P)*exp((18.177-5.987*1 0^(-7)*P)*Xe1))*((0.1484-1.596*Xe1+0.1729*Xe1*abs(Xe1))*0.2049*G/10^6+1.037)*(1.157-0. 869*Xe1)*(0.2664+0.8357*exp(-124*De))*(0.8258+0.341*10^(-6)*(Hfs-H0));DNBR1=q_DNB1/(q*F_rn*F_qe*f1);if (t_fh2<=290)elseif (t_fh2<=300)H2=(5.344*t_fh2-265.64)*10^3;elseif (t_fh2<=310)H2=(5.59*t_fh2-339.44)*10^3;elseif (t_fh2<=320)H2=(5.932*t_fh2-445.46)*10^3;elseif (t_fh2<=330)H2=(6.422*t_fh2-602.26)*10^3;elseif (t_fh2<=340)H2=(7.23*t_fh2-868.9)*10^3;elseif (t_fh2<=350)H2=(104.63*t_fh2-33984.9)*10^3;endXe2=(H2-Hfs)/Hfg;q_DNB2=3.154*10^6*((2.022-6.238*10^(-8)*P)+(0.1722-1.43*10^(-8)*P)*exp((18.177-5.987*1 0^(-7)*P)*Xe2))*((0.1484-1.596*Xe2+0.1729*Xe2*abs(Xe2))*0.2049*G/10^6+1.037)*(1.157-0. 869*Xe2)*(0.2664+0.8357*exp(-124*De))*(0.8258+0.341*10^(-6)*(Hfs-H1));DNBR2=q_DNB2/(q*F_rn*F_qe*f2);if (t_fh3<=290)H3=(5.146*t_fh3-208.22)*10^3;elseif (t_fh3<=300)H3=(5.344*t_fh3-265.64)*10^3;elseif (t_fh3<=310)H3=(5.59*t_fh3-339.44)*10^3;elseif (t_fh3<=320)H3=(5.932*t_fh3-445.46)*10^3;elseif (t_fh3<=330)H3=(6.422*t_fh3-602.26)*10^3;elseif (t_fh3<=340)H3=(7.23*t_fh3-868.9)*10^3;elseif (t_fh3<=350)H3=(104.63*t_fh3-33984.9)*10^3;endXe3=(H3-Hfs)/Hfg;q_DNB3=3.154*10^6*((2.022-6.238*10^(-8)*P)+(0.1722-1.43*10^(-8)*P)*exp((18.177-5.987*1 0^(-7)*P)*Xe3))*((0.1484-1.596*Xe3+0.1729*Xe3*abs(Xe3))*0.2049*G/10^6+1.037)*(1.157-0. 869*Xe3)*(0.2664+0.8357*exp(-124*De))*(0.8258+0.341*10^(-6)*(Hfs-H2));DNBR3=q_DNB3/(q*F_rn*F_qe*f3);if (t_fh4<=290)H4=(5.146*t_fh4-208.22)*10^3;elseif (t_fh4<=300)。

返回第五章　反应堆稳态热工设计原理第六章 反应堆瞬态热工分析简介 (1)§6.1 瞬态过程中反应堆功率计算........................................................................1 §6.2 瞬态工况燃料元件温度场计算....................................................................2 §6.3 基本方程组....................................................................................................3 §6.4 反应堆的安全问题........................................................................................4 §6.5 反应堆失流事故............................................................................................6 §6.6 冷却剂丧失事故.. (6)第六章 反应堆瞬态热工分析简介§6.1 瞬态过程中反应堆功率计算0.11101001000100000.010.11.00图6-1 衰变功率裂变产物的衰变功率：对于稳定运行了很长时间的压水堆，停堆后裂变产物的衰变功率在许多人的实验结果上得出曲线图6-2，也可以表示为：其中：A=53.18 ，α=0.3350剩余裂变功率：裂变时瞬间放出的功率大小与堆芯的热中子密度成正比，可由中子动力学方程计算得到。

对于以恒定功率运行了很长时间的压水堆，如果引入的负反应性绝对值大于4%，则在剩余裂变功率其重要作用的期间内，可用下式估算：对于重水堆，中子俘获产物衰变功率：在用天然铀或低浓缩铀作燃料的反应堆中，对中子俘获产物衰变功率贡献最大的是铀-238吸收中子后产生铀-239（T 1/2=23.5分）和由它衰变成的镎-239（ T 1/2=对于停堆前运行了很长时间的压水堆，C=0.6 ，α=0.2 ，由于忽略了其它俘获产物，还要乘1.1的安全系数。

目录一、设计任务 (1)二、课程设计要求 (2)三、计算过程 (2)四、程序设计框图 (8)五、代码说明书 (9)六、热工设计准则和出错矫正 (10)七、重要的核心程序代码 (11)八、计算结果及分析 (17)一、设计任务某压水反应堆的冷却剂及慢化剂都是水，用二氧化铀作燃料，用Zr-4作包壳材料。

燃料组件无盒壁，燃料元件为棒状，正方形排列。

已知下列参数：系统压力 15.8MPa堆芯输出功率 1820MW冷却剂总流量 32100t/h反应堆进口温度287℃堆芯高度 3.66m燃料组件数 121燃料组件形式17×17每个组件燃料棒数 265燃料包壳直径 9.5mm燃料包壳内径 8.36mm燃料包壳厚度 0.57mm燃料芯块直径 8.19mm燃料棒间距（栅距） 12.6mm芯块密度 95%理论密度旁流系数 5%燃料元件发热占总发热的份额 97.4%径向核热管因子 1.35轴向核热管因子 1.528局部峰核热管因子 1.11交混因子 0.95热流量工程热点因子 1.03焓升工程热管因子 1.085堆芯入口局部阻力系数 0.75堆芯出口局部阻力系数 1.0堆芯定位隔架局部阻力系数 1.05若将堆芯自上而下划分为5个控制体，则其轴向归一化功率分布如下表：堆芯轴向归一化功率分布（轴向等分5个控制体）通过计算，得出1. 堆芯出口温度；2. 燃料棒表面平均热流及最大热流密度，平均线功率，最大线功率；3. 热管的焓，包壳表面温度，芯块中心温度随轴向的分布；4. 包壳表面最高温度，芯块中心最高温度；5. DNBR在轴向上的变化；6. 计算堆芯压降；二、课程设计要求1．设计时间为两周；2．独立编制程序计算；3．迭代误差为0.1%；4．计算机绘图；5．设计报告写作认真，条理清楚，页面整洁；6．设计报告中要附源程序。

三、计算过程目前，压水核反应堆的稳态热工设计准则有：（1）燃料元件芯块内最高温度应低于其相应燃耗下的熔化温度。

核反应堆热⼯分析课设⽬录⼀、设计任务 (1)⼆、课程设计要求 (2)三、计算过程 (2)四、程序设计框图 (8)五、代码说明书 (9)六、热⼯设计准则和出错矫正 (10)七、重要的核⼼程序代码 (11)⼋、计算结果及分析 (17)⼀、设计任务某压⽔反应堆的冷却剂及慢化剂都是⽔，⽤⼆氧化铀作燃料，⽤Zr-4作包壳材料。

燃料组件⽆盒壁，燃料元件为棒状，正⽅形排列。

已知下列参数：系统压⼒ 15.8MPa堆芯输出功率 1820MW冷却剂总流量 32100t/h反应堆进⼝温度287℃堆芯⾼度 3.66m燃料组件数 121燃料组件形式17×17每个组件燃料棒数 265燃料包壳直径 9.5mm燃料包壳内径 8.36mm燃料包壳厚度 0.57mm燃料芯块直径 8.19mm燃料棒间距（栅距） 12.6mm芯块密度 95%理论密度旁流系数 5%燃料元件发热占总发热的份额 97.4%径向核热管因⼦ 1.35轴向核热管因⼦ 1.528局部峰核热管因⼦ 1.11交混因⼦ 0.95热流量⼯程热点因⼦ 1.03堆芯出⼝局部阻⼒系数 1.0堆芯定位隔架局部阻⼒系数 1.05若将堆芯⾃上⽽下划分为5个控制体，则其轴向归⼀化功率分布如下表：堆芯轴向归⼀化功率分布（轴向等分5个控制体）通过计算，得出1. 堆芯出⼝温度；2. 燃料棒表⾯平均热流及最⼤热流密度，平均线功率，最⼤线功率；3. 热管的焓，包壳表⾯温度，芯块中⼼温度随轴向的分布；4. 包壳表⾯最⾼温度，芯块中⼼最⾼温度；5. DNBR在轴向上的变化；6. 计算堆芯压降；⼆、课程设计要求1．设计时间为两周；2．独⽴编制程序计算；3．迭代误差为0.1%；4．计算机绘图；5．设计报告写作认真，条理清楚，页⾯整洁；6．设计报告中要附源程序。

三、计算过程⽬前，压⽔核反应堆的稳态热⼯设计准则有：（1）燃料元件芯块内最⾼温度应低于其相应燃耗下的熔化温度。

⽬前，压⽔堆⼤多采⽤UO2作为燃料。

第六章反应堆动态热工分析6.1 瞬态过程反应堆功率计算衰变功率的衰减(除显热)。

6.2 动态工况下燃料元件温度场的计算6.3 基本方程组6.4 反应堆的安全问题反应堆的事故额外的反应性引入堆芯冷却能力不足正常运行的瞬态过程常见事故出现的可能性很小的事故 极限事故（设计基准事故）反应堆的安全保护:紧急停堆系统，要求接到停堆信号后投入运行的速度要快，在重大事故后有能力连续运行一定的时间，以去除堆芯的衰变热。

还要有足够的停堆深度，以保证在事故过程中堆芯始终处于次临界状态。

工程安全设施；安全壳密封装置；安全壳大气的降压设施；安全壳大气中去除放射性物质的设施；另外，为了保证安全保护系统能够随时投入工作，必须设置应急电源，一般都备有快速启动的柴油发电机组和蓄电池组。

专设安全系统除了用控制系统实现停堆保护外，专门装备的安全系统（即专设安全系统）发挥作用来限制事故的后果。

以压水堆的专设安全系统为例介绍：1.应急堆芯冷却系统（安全注射系统）当一回路系发生冷却剂丧失事故时，把足够的应急冷却水注入堆芯，以防止燃料过热。

分能动和非能动两类。

2. 辅助给水系统在二回路主给水流量丧失的情况下向蒸汽发生器二次侧供水，以维持蒸汽发生器的排热能力，冷却一回路。

3. 安全壳喷淋系统用喷淋水泵把含硼水送到安全壳的顶部，通过喷嘴向壳内空间喷淋，用以抑制一回路或二回路发生大破口事故时安全壳内压力上升过高，防止安全壳超压，喷淋水中可以加氢氧化钠，它有助于除去泄露的冷却剂中的放射性物质(主要是碘)。

4. 其他安全设施放射性去除系统、消氢系统和贯穿件密封装置等6.5 负荷丧失瞬态6.5 失流事故。

tfin=279.4;Fa=0.974;Nt=3400000000;Wt=14314;b=0.059; tfout=340;e0=0.01;while e0>0.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*(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;while e11>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf11<dtf12 %膜温压取两个中较小值，算得包壳外表面温度 tcs1=tf1+dtf11elsetcs1=tf1+dtf12enddci=8.93e-3;tci1=349;e12=0.01;while e12>0.001t12_=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 %采用迭代算法求得包壳内表面温度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;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*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;while e21>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf21<dtf22 %膜温压取两个中较小值，算得包壳外表面温度tcs2=tf2+dtf21elsetcs2=tf2+dtf22enddci=8.93e-3;tci2=349;e22=0.01;while e22>0.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^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;while e31>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf31<dtf32 %膜温压取两个中较小值，算得包壳外表面温度tcs3=tf3+dtf31elsetcs3=tf3+dtf32enddci=8.60e-3;tci3=349;e32=0.01;while e32>0.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^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;while e41>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf41<dtf42 %膜温压取两个中较小值，算得包壳外表面温度tcs4=tf4+dtf41elsetcs4=tf4+dtf42enddci=8.93e-3;tci4=349;e42=0.01;while e42>0.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^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;while e51>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf51<dtf52 %膜温压取两个中较小值，算得包壳外表面温度tcs5=tf5+dtf51elsetcs5=tf5+dtf52enddci=8.93e-3;tci5=350;e52=0.01;while e52>0.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^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;while e61>0.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*(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; %采用詹斯-洛特斯传热方程算得的过冷沸腾膜温压if dtf61<dtf62 %膜温压取两个中较小值，算得包壳外表面温度tcs6=tf6+dtf61elsetcs6=tf6+dtf62enddci=8.93e-3;tci6=350;e62=0.01;while e62>0.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^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为按照壁面温度取值的流体的粘性系数。