snnew.rar

%单速一维均匀平板裸堆的SN解法 %算出角通量分布，通量分布和有效增值系数 %变量定义 k=500; %划分区间数 m=4; %离散化的方向数 S=100*ones(1,k); %源项，赋初值为100 H=zeros(4,2*k+1); %中子角通量分布 keff=1; %有效增殖系数，赋初值为1 k_eff=zeros(2*k+1,1); %储存每次迭代的keff值 a=66.0053; %平板半厚度(cm) sigma_t=0.05;% 宏观截面数据(cm^(-1)) sigma_s=0.03; nusigma_f=0.0225; u=[-0.86113 -0.33998 0.33998 0.86113]; w=[0.34785 0.65214 0.65214 0.34785]; OUT=zeros(4,2*k+1); delta=a/k; %差分步长 km=sigma_t*delta./u; tol=0.001; %精度 bad_grid=1; %bad_grid是迭代前后相差超过误差允许值的节点 cycle_count=0; %记录迭代次数 %源迭代 while bad_grid %flag非零，表示精度不够，需要继续迭代 bad_grid=0; cycle_count=cycle_count+1; S1=S; %迭代源必须每次除以keff进行初始化，确保keff收敛 H1=H; %计算u<0的部分 for i=1:2 H(i,2*k+1)=0;%真空边界条件 for j=1:2:2*k-1 H(i,2*k-j)=(2+km(i))/(2-km(i))*H(i,2*k+2-j)... -2*km(i)/(2-km(i))*S((2*k+1-j)/2)/sigma_t; end end %计算u>0的部分 for i=3:4 H(i,1)=H(5-i,1);%对称边界条件 for j=3:2:2*k+1 H(i,j)=(2-km(i))/(2+km(i))*H(i,j-2)... +2*km(i)/(2+km(i))*S((j-1)/2)/sigma_t; end end %重新计算源项 for j=1:k src=0; for i=1:4 src=src+w(i)*(H(i,2*j-1)+H(i,2*j+1)); end S(j)=(sigma_s+nusigma_f)/4*src; %外源项q(x)取0 end keff=sum(S)/sum(S1); %有效增值系数 k_eff(cycle_count)=keff; %计算误差 if max(max(abs(H-H1)))>tol %精度取千分之一 bad_grid=bad_grid+1; end end for i=1:4 for j=1:2:2*k-1 H(i,j+1)=(H(i,j)+H(i,j+2))/2; end; end; %输出结果 for i=1:4 for j=1:2*k+1 OUT(i,j)=H(5-i,j); end end disp(OUT); %中子角通量分布 Phi=w*H; disp(Phi); %中子通量分布 disp(Phi/Phi(1)); %中子通量分布（中心归一化） str=['迭代次数 cycle_count =',num2str(cycle_count)]; disp(str); %迭代次数 str2=['有效增值系数 keff=',num2str(keff)]; disp(str2); %有效增值系数 %作图 x=0:delta/2:a; figure(1); %中子角通量分布图 for i=1:4 plot(x,OUT(i,:)); hold on; end hold off; figure(2); %中子通量分布图 y=Phi/Phi(1); plot(x,y);