堆芯一维中子扩散方程的数值解法(通用型)

10-11-17 下午1:11
F:\my_matlab\core_1D.m
3 of 3
flux_fast=tdma(a1,b1,c1,d1);%解出快群. %慢群 d2=[]; for i=1:N for j=1:n d2(j+(i-1)*n)=E12(m(i))*flux_fast(j+(i-1)*n)*dx; end end d2(1)=d2(1)+2*f_left*D(m(1),2)/dx; d2(N*n)=d2(N*n)+2*f_right*D(m(N),2)/dx; flux_slow=tdma(a2,b2,c2,d2);%解出慢群. for i=1:N for j=1:n Q_new(j+(i-1)*n)=flux_fast(j+(i-1)*n)*vEf(m(i),1)+flux_slow(j+(i-1)*n)*vEf(m(i),2); end end %更新源项. keff_new=sum(Q_new)*dx/(sum(Q_old)*dx/keff_old);%更新keff. e_keff=abs((keff_new-keff_old)/(keff_new+eps)); keff_old=keff_new; e_Q=max(abs((Q_new-Q_old)./(Q_new+eps))); Q_old=Q_new; cycle=cycle+1;%记录循环次数. if e_keff<error_keff && e_Q<error_Q break; end end %归一化. flux_fast=flux_fast/(sum(flux_fast)); flux_slow=flux_slow/(sum(flux_slow)); Q=Q_new'; Q=N*L*Q/(dx*sum(Q)); keff=keff_new; end
10-11-17 下午1:11
F:\my_matlab\core_1D.m
1 of 3
function [keff flux_fast flux_slow Q]=core_1D(m) %m为燃料排布,L为结块长度,N为结块数目,n为每个结块划分份额. %材料参数. Ea=[1.071240E-02 9.389003E-02 8.915481E-03 9.007146E-02 1.033911E-02 9.685340E-02 9.920272E-03 9.850721E-02 1.033911E-02 9.685340E-02 9.616552E-03 9.824058E-02 3.064667E-03 4.777252E-02 8.603110E-03 7.853440E-02]; vEf=[4.334360E-03 1.029336E-01 6.260639E-03 1.217170E-01 4.769662E-03 1.141428E-01 5.251954E-03 1.232079E-01 4.769662E-03 1.141428E-01 5.555385E-03 1.264411E-01 0.000000E+00 0.000000E+00 6.160544E-03 1.207603E-01]; E12=[1.771023E-02 1.812268E-02 1.770407E-02 1.776033E-02 1.770407E-02 1.785182E-02 2.090091E-02 1.708523E-02]; D=[1.378098E+00 3.437887E-01 1.373819E+00 3.542306E-01 1.376979E+00 3.443704E-01 1.376071E+00 3.457435E-01 1.376979E+00 3.443704E-01 1.375501E+00 3.472169E-01 1.134501E+00 2.587405E-01 1.407447E+00 0.3670093E+0]; L=20;n=20;N=length(m); dx=L/n;%等边距划分,dx为每份宽度. Q_old=10*ones(1,N*n);keff_old=1.7;%设置初始值. error_keff=1e-6;error_Q=1e-5;%方程收敛验收系数. cycle=0;%存储计算循环次数. f_left=0;f_right=0;%边界条件. %三对角矩阵系数的对角系数向量计算-快群. a1=[];b1=[];c1=[]; for i=1:N a1=[a1,(Ea(m(i),1)*dx+E12(m(i))*dx+2*D(m(i),1)/dx)*ones(1,n)]; end for i=1:N-1 a1(i*n)=(Ea(m(i),1)+E12(m(i)))*dx+(D(m(i),1)^2+3*D(m(i),1)*D(m(i+1),1))/((D(m(i),1)+D(m(i+1),1)) *dx); a1(i*n+1)=(Ea(m(i+1),1)+E12(m(i+1)))*dx+(D(m(i+1),1)^2+3*D(m(i),1)*D(m(i+1),1))/((D(m(i),1)+D(m (i+1),1))*dx);
10-11-17 下午1:11
F:\my_maΒιβλιοθήκη Baidulab\core_1D.m
2 of 3
end a1(1)=3*D(m(1),1)/dx+Ea(m(1),1)*dx+E12(m(1))*dx; a1(N*n)=3*D(m(N),1)/dx+Ea(m(N),1)*dx+E12(m(N))*dx; for i=1:N b1=[b1,-(D(m(i),1)/dx)*ones(1,n)]; end for i=1:N-1 b1(i*n)=-2*D(m(i),1)*D(m(i+1),1)/((D(m(i),1)+D(m(i+1),1))*dx); end for i=1:N c1=[c1,(-D(m(i),1)/dx)*ones(1,n)]; end for i=1:N-1 c1(i*n+1)=-2*D(m(i),1)*D(m(i+1),1)/((D(m(i),1)+D(m(i+1),1))*dx); end %三对角矩阵系数的对角系数向量计算-慢群 a2=[];b2=[];c2=[]; for i=1:N a2=[a2,(Ea(m(i),2)*dx+2*D(m(i),2)/dx)*ones(1,n)]; end for i=1:N-1 a2(i*n)=Ea(m(i),2)*dx+(D(m(i),2)^2+3*D(m(i),2)*D(m(i+1),2))/((D(m(i),2)+D(m(i+1),2))*dx); a2(i*n+1)=Ea(m(i+1),2)*dx+(D(m(i+1),2)^2+3*D(m(i),2)*D(m(i+1),2))/((D(m(i),2)+D(m(i+1),2))*dx); end a2(1)=3*D(m(1),2)/dx+Ea(m(1),2)*dx; a2(N*n)=3*D(m(N),2)/dx+Ea(m(N),2)*dx; for i=1:N b2=[b2,(-D(m(i),2)/dx)*ones(1,n)]; end for i=1:N-1 b2(n*i)=-2*D(m(i),2)*D(m(i+1),2)/((D(m(i),2)+D(m(i+1),2))*dx); end for i=1:N c2=[c2,(-D(m(i),2)/dx)*ones(1,n)]; end for i=1:N-1 c2(i*n+1)=-2*D(m(i),2)*D(m(i+1),2)/((D(m(i),2)+D(m(i+1),2))*dx); end %迭代计算 while(1) %快群 d1=[]; d1=Q_old*dx/keff_old; d1(1)=d1(1)+2*f_left*D(m(1),1)/dx; d1(N*n)=d1(N*n)+2*f_right*D(m(1),1)/dx;