堆芯设计

clc clear clcclear allNe=10^9; %电站电功率nt=0.3333; %效率Fu=0.974; %燃料释热占堆芯总释热的份额P=1.55*10^7; %冷却剂压力x=0; %质量含汽率Hsat=1629850; %饱和水焓Hin=1289420; %进口冷却剂比焓Hout=1459050; %出口冷却剂比焓a=0.040; %格架混流扩散系数Fc=1.05; %热流密度不均匀修正因子MDNBR=1.34; %最小烧毁比F=2.5752; %总热流密度热点因子ds=0.0095; %燃料元件外径HEU=2.1; %水铀比Oc=0.00057; %包壳厚度Og=0.00017; %燃料与包壳间隙n=264; %组件内的元件数LED=1.15; %高径比Ea=0.04; %旁流系数G=100; %假设冷却剂质量流速Gm=200; %冷却剂质量流速Nt=Ne/nt; %反应堆热功率nn=1;%反应堆稳态热工计算p1=sqrt(pi/4*HEU*(ds-2*Oc-2*Og)^2+pi*ds^2/4); %正方形栅元栅距De=4*(p1^2-pi*ds^2/4)/(pi*ds); %当量直径while abs(G-Gm)/Gm>0.005if abs(G-Gm)/Gm>0.005&&nn>1G=Gm;endnn=nn+1;Fs=1.0+0.03*G*(a/0.019)^0.35/(4.882*10^6); %定位格架搅混修正qDNBREU=3.154*10^6*((2.022-6.238*10^(-8)*P)+(0.1722-1.427*10^(-8)*P)*exp((18.177-5.987* 10^(-7)*P)*x))*(1.157-0.869*x)*((0.1484-1.596*x+0.1729*x*abs(x))*G/(10^6)*0.2048+1.037)*(0 .2664+0.8357*exp(-124*De))*(0.8258+0.341*10^(-6)*(Hsat-Hin))*Fs; %临界热流密度qDNBn=qDNBREU/Fc; %非均匀加热时的临界热流密度qmax=qDNBn/MDNBR; %最大允许热流密度q0=qmax/F; %平均热流密度At=Nt*Fu/q0; %燃料元件总传热面积l=17*p1; %组件的边长Def=(4*l^2*Nt*Fu/(n*pi^2*ds*q0*LED))^(1/3); %堆芯等效直径L=Def*LED; %活性段高度N1=Nt*Fu/(pi*ds*q0*L); %燃料元件总根数N=n*round(N1/n); %取整后的燃料元件根数Weff=Nt/(Hout-Hin); %反应堆有效流量Wt=3600*Weff/(1000*(1-Ea)); %总流量Ab=p1^2-pi*ds^2/4; %一根燃料元件栅元的冷却剂通流截面Gm=3600*Weff/(N*Ab); %冷却剂质量流速AL=pi*ds; %一根燃料元件单位长度上的外表面积Endz=sym('z'); %使得z成为一个自变量Hfmz=Hin+vpa(3600*q0*AL/(G*Ab),8)*z; %平均管冷却剂的焓场pp0=713.75; %平均密度pp2=677.54; %出口密度pp1=744.31; %进口密度uf=8.616*10^(-5); %平均温度下的流体粘度us=7.615*10^(-5); %壁面温度下的流体的粘度g=9.8; %重力加速度A=0.1166; %燃料元件正方形排列的系数Ke=8; %一组燃料组件轴向定位格架数Kgd=0.1107; %定位架形阻系数oo=0.05; %下腔室流量系数m=0.2; %指数系数Pel=pp0*g*L; %提升压降Pa=(Gm/3600)^2*(1/pp2-1/pp1); %加速压降V1=Gm/3600/pp1; %截面1冷却剂流速V2=Gm/3600/pp2; %截面2冷却剂流速V=Gm/3600/pp0; %冷却剂平均流速Re=De*V*pp0/uf; %雷诺数fa=0.3164/(Re^0.25); %等温流动圆形通道摩擦系数B=740.3*10^(-6)*(p1/ds)^3*(1.273*(p1/ds)^2-1)^(3/4)/(1.122*(p1/ds)-1)^(9/2);%燃料元件正方形排列的系数fb=fa*(1+A*(p1/ds)^(4/3))*(0.58+0.42*exp(-B));%湍流等温流动状态棒束的摩擦阻力系数f=fb*(us/uf)^0.6; %飞等温流动是的阻力系数Pf=f*L*pp0*V^2/(2*De); %摩擦压降Pgd=Kgd*Ke*pp0*V^2/2; %定位格架型阻压降Kfh=(1-oo)^(2-m); %热管摩擦压降的下腔室修正因子Kah=(1-oo)^m; %热管各型阻压降及加速压降的下腔室修正因子Phe=Kfh*Pf+Kah*(Pa+Pgd)+Pel; %热管有效驱动压头tin=291; %冷却剂进口温度FHN=1.12; %焓升核热点因子FHE=1.02; %焓升工程热管因子Cp=5618.9; %按平均温度计算冷却剂比热J/kgsGh=1.05*Gm; %热管冷却剂的质量流速du=ds-2*(Oc+Og); %燃料芯块直径Prf=0.87884; %按平均温度计算冷却剂的普朗特常数Wf=0.5484; %按平均温度计算的冷却剂的导热系数Wu=2.32; %燃料芯块导热系数Wc=14.3; %包壳导热系数Wg=0.29; %气隙导热系数Hfhz=Hin+vpa(3600*q0*FHN*FHE*AL/(Gh*Ab),8)*z; %热管冷却剂焓场Au=pi*du^2/4; %燃料芯块面积Le=1.0664*L; %堆芯外推高度XS=Le*(2*sin(pi*L/2/Le))/(pi*L); %轴向平均功率计算的系数Ar=5.92*ds*q0/du^2; %轴向功率表达式系数Br=pi/Le; %热管轴向功率表达式系数Firz=vpa(Ar,8)*cos(vpa(Br,8)*z); %轴向功率分布表达式Atf=3600*Ar*Au*FHN*FHE/(Br*Ab*Gh*Cp); %热管冷却剂温度表达式系数thz=vpa(Atf,6)*sin(vpa(Br,5)*z)+vpa(tin+Atf*sin(Br*L/2),7);%热管冷却剂温度表达式thfhL=tin+Atf*(sin(Br*L/2)+sin(Br*L/2)); %热管出口温度C0=0.042*p1/ds-0.024; %冷却剂对流换热系数的系数h=C0*Wf*Re^0.8*Prf^0.4/De; %冷却剂对流换热系数Atc=tin+Atf*sin(Br*L/2); %燃料包壳温度表达式系数Btc=Atf; %燃料包壳温度表达式系数Ctc=Au*Ar/(pi*h*ds); %燃料包壳温度表达式系数thcz=vpa(Atc,7)+vpa(Btc,5)*sin(vpa(Br,4)*z)+vpa(Ctc,5)*cos(vpa(Br,4)*z); %包壳表面温度分布函数Zc=atan(Btc/Ctc)/Br; %燃料包壳最高温度位置（相对于堆芯中分面）thcm=Atc+Btc*sin(Br*Zc)+Ctc*cos(Br*Zc); %燃料包壳最高温度Aus=Au*((log(ds/(ds-2*Oc)))/(2*pi*Wc)+Og/(pi*du*Wg)); %系数Ctu=Aus*Ar+Ctc+Ar*du^2/(16*Wu); %燃料中心温度表达式系数thucz=vpa(Atc,7)+vpa(Btc,7)*sin(vpa(Br,4)*z)+vpa(Ctu,7)*cos(vpa(Br,4)*z);%燃料中心温度的分布函数Zu=atan(Btc/Ctu)/Br; %燃料中心最高温度位置（相对于堆芯中分面）thum=Atc+Btc*sin(Br*Zu)+Ctu*cos(Br*Zu); %燃料中心最高温度%反应堆稳态工况水力计算Sloop=3; %环路数neff=0.75; %主泵的总效率Pin1=0.0275*10^6; %反应堆入口管嘴Pin2=0.201*10^6; %反应堆入口段Pex1=0.037*10^6; %反应堆出口段Pex2=0.0275*10^6; %反应堆出口管嘴PR=Pf+Pa+Pgd+Pel; %堆芯平均管压降P1=Pin1+Pin2+Pex1+Pex2+PR; %反应堆总压降P2=0.175*10^6; %一回路管道、主泵及蒸汽发生器一次侧、阀门压降Pt=P1+P2; %一回路总压降Wul=Wt/(3.6*pp0); %冷却剂的总质量流量Npe=Wul*Pt/(1000*Sloop); %主泵的有效功率，KW Npz=Npe/neff; %主泵的轴功率%自然循环能力计算L0=12; %热源与冷源之间的高度Kzb=20; %主泵阻力系数dzb=0.8; %主泵管道的直径Pnd=(pp1-pp2)*g*L0; %自然循环的驱动压头Vn=0; %自然循环流量假定值Pnt=100; %自然循环时一回路总压降假定值while abs(Pnd-Pnt)/Pnd>0.0001 Vn=Vn+0.00001;Pnin1=Pin1*(Vn/V)^2; %自然循环时反应堆入口管嘴压降Pnin2=Pin2*(Vn/V)^2; %自然循环时反应堆入口段压降Pnex1=Pex1*(Vn/V)^2; %自然循环时反应堆出口段压降Pnex2=Pex2*(Vn/V)^2; %自然循环时反应堆出口管嘴压降Pnel=0; %自然循环时提升压降Pna=Pa*(Vn/V); %自然循环时加速压降Pnf=Pf*(Vn/V)^(1.75); %自然循环时摩擦压降Pngd=Pgd*(Vn/V)^2; %自然循环时定位格架型阻压降PnR=Pnf+Pna+Pngd+Pnel; %自然循环时堆芯平均管压降Pn1=Pnin1+Pnin2+Pnex1+Pnex2+PnR; %自然循环时反应堆总压降Gn=pp0*Vn; %自然循环质量流速Wneff=Gn*N*Ab; %自然循环时的有效质量流量Wnt=Wneff/(1-Ea); %自然循环时的总质量流量Vzb=4*Wnt/(pi*Sloop*pp1*dzb^2); %自然循环时主泵管道中的冷却剂流速Pn2=P2*(Vn/V)^2+0.5*Kzb*pp1*Sloop*Vzb^2;%自然循环时一回路管道、主泵及蒸汽发生器一次侧、阀门压降Pnt=Pn1+Pn2; %自然循环时一回路总压降EndNnt=Wneff*(Hout-Hin); %自然循环时的热功率Ablity=Nnt/Nt; %自然循环能力%压力容器强度计算%Di=sqrt((3.49+0.176)^2+4*Wt/(7*pi*3.6*pp1))Di=4.3; %圆柱壳内径Sm=184; %材料许用应力Pdes=1.11*P*10^(-6); %设计压力Fr=1; %焊缝系数c=0.1; %壁厚附加量tcmin=Pdes*Di/(2*Sm*Fr-Pdes)+c; %圆柱壳最小壁厚tsmin=Pdes*Di/(4*Sm*Fr-2*Pdes)+c; %球壳最小壁厚temin=Pdes*Di/(2*Sm*Fr-Pdes)+c; %椭球壳最小壁厚Do=4.95; %圆柱壳外径Dgi1=0.8; %进口管道内径ttmin1=(Pdes*Dgi1/(2*Sm*Fr+Pdes)+c)/(1-2*Pdes/(2*Sm*Fr+Pdes)); %进口管道最小壁厚Dgi2=0.838; %出口管道内径ttmin2=(Pdes*Dgi2/(2*Sm*Fr+Pdes)+c)/(1-2*Pdes/(2*Sm*Fr+Pdes)); %出口管道最小壁厚%结果输出disp('反应堆热功率W')Ntdisp('燃料元件总传热面积㎡')Atdisp('均匀加热时的临界热流密度W/㎡')qDNBREUdisp('冷却剂质量流速kg/(m2*h)')Gdisp('当量直径m')Dedisp('正方形栅元栅距m')p1disp('定位格架修正因子')Fsdisp('非均匀加热时的修正因子')qDNBndisp('最大允许热流密度W/㎡')qmaxdisp('平均热流密度W/㎡')q0disp('燃料元件的总根数根')Ndisp('堆芯高度m')Ldisp('燃料组件边长m')ldisp('堆芯等效直径m')Defdisp('组件数量个')N/ndisp('反应堆有效流量Kg/s')Weffdisp('堆芯总流量T/h')Wtdisp('平均管冷却剂质量流速Kg/(㎡*h)')Gmdisp('栅元冷却剂流通截面㎡')disp('平均管冷却剂的焓场表达式')Hfmzdisp('燃料元件单位长度外表面积㎡/m')ALdisp('提升压降Pa')Peldisp('加速压降Pa')Padisp('截面1冷却剂流速m/s')V1disp('截面2冷却剂流速m/s')V2disp('摩擦压降Pa')Pfdisp('等温流动圆形通道摩擦系数')fadisp('雷诺数')Redisp('等温棒束摩擦阻力系数')fbdisp('燃料元件排列系数')Bdisp('摩擦系数')fdisp('定位格架型阻系数')Pgddisp('棒束中的平均流速m/s')Vdisp('热管有效驱动压头Pa')Phedisp('热管摩擦压降下腔室修正因子')Kfhdisp('热管各型阻压降及加速压降下腔室修正因子') Kahdisp('热管冷却剂焓场')Hfhzdisp('热管冷却剂的质量流速kg/(㎡h)')Ghdisp('堆芯外推高度m')Ledisp('轴向功率分布表达式')Firzdisp('热管轴向功率表达式系数W/m3')disp('热管轴向功率表达式系数/m')Brdisp('燃料棒芯块截面积㎡')Audisp('热管冷却剂温度表达式')thzdisp('热管冷却剂温度表达式系数℃')Atfdisp('热管出口温度℃')thfhLdisp('冷却剂对流换热系数W/(㎡K)')hdisp('燃料包壳温度表达式系数℃')Atcdisp('燃料包壳温度表达式系数℃')Btcdisp('燃料包壳温度表达式系数℃')Ctcdisp('包壳表面温度分布函数')thczdisp('燃料包壳最高温度位置（相对于堆芯中分面）m') Zcdisp('燃料包壳最高温度℃')thcmdisp('系数(m3K)/W')Ausdisp('燃料中心温度表达式系数℃')Ctudisp('燃料中心温度的分布函数')thuczdisp('燃料中心最高温度位置（相对于堆芯中分面）m') Zudisp('燃料中心最高温度℃')thumdisp('反应堆总压降Pa')P1disp('一回路总压降Pa')Ptdisp('主泵有效功率KW')Npedisp('冷却剂的总体积流量m3/s')Wuldisp('主泵轴功率KW')disp('自然循环能力')Ablitydisp('圆柱壳最小壁厚m') tcmindisp('球壳最小壁厚m') tsmindisp('椭球壳最小壁厚m') temindisp('进口管道最小壁厚m') ttmin1disp('出口管道最小壁厚m') ttmin2disp('压力容器内径m')Didisp('压力容器外径m')Dodisp('进口管道内径m')Dgi1disp('出口管道内径m')Dgi2。