《中子扩散作业》由会员分享，可在线阅读，更多相关《中子扩散作业（22页珍藏版）》请在金锄头文库上搜索。
1、中子扩散方程数值求解报告小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX1.扩散方程的离散本次求解中子扩散方程采用的方法是有限差分方法，利用直线族把计算平面分成许多矩阵网格。习惯上，我们把网格分界线与不同介质分界重合，称需要求解的未知函数离散点为节点（网点）。取节点为网格交界点，称为角点差分。xi=xi+1-xi，yi=yi+1-yi称为网距。二维(x,y)几何情况下扩散方程式可写成-xDFx-yDFy+SRx,yFx,y=S(x,y)经过推导，可导出(i,j)节点的差分方程：ai,jFi,j-1+bi,jFi-1,j+ci,jFi+1,j+di,jFi,j+1+ei,jFi
2、,j=si,j式中，ai,j=di,j-1=-12D1xi-1+D2xiyi-1bi,j=ci-1,j=-12D1yj-1+D4yixi-1ei,j=14SR,1Dxi-1Dyi-1+SR,2DxiDyi-1+SR,3DxiDyi+SR,4Dxi-1Dyi-ai,j-bi,j-ci,j-di,jsi,j=14S1xi-1yj-1+S2xiyj-1+S3xiyj+S4xi-1yj得到了上诉个参量之后，应用源迭代法来进行中子通量的计算。应用源迭代的一个很大的优点和方便之处在于：它把原来G个联立的多群扩散方程的求解问题在每次迭代中，变换为按顺序解下面G个单群扩散方程问题：-DgFgnr+R,gFgn
3、r=Sgnrg=1,G当按g=1,G的次序对方程求解时，方程右端源项Sgnr便是已知的，是一个含有已知源项的单群扩散方程。略去上下角标，它可以写成下面标准形式-DFr+RFr=S(r)式中S(r)为已知源项。在反应堆的实际计算中，考虑到控制棒的布置、燃料的非均匀装载以及燃耗等问题，对于大型反应堆往往需要计算二位或三位的中子通量密度或功率的分布。2.边界条件对于边界上的点则和前面一样，需根据边界条件来确定，例如对于四周的外推边界上的所有点有Fi,j=0若i=0,N,或j=0,M对于二维问题是五对角矩阵。利用有限次的算术运算来逐次逼近方程组精确解，也就是通过迭代方程来解其近似解。这种利用迭代方法求
4、解差分方程组的过程通常称为内迭代。内迭代过程的收敛通常使用下面的收敛准则来判定的：为给定参数，一般讲当外迭代远未收敛时，不必使内迭代有太高精度，内迭代的要求精度可随着外迭代的收敛精度而提高。常使用雅可比迭代，赛德尔迭代，超松弛迭代，线性逐次超松弛方法，隐式交替方向法等迭代方法求解。3.材料的处理在cross子程序中，根据所给文献将堆芯分成若干区，确定点(i,j) 在不同材料，不同群的条件下的吸收截面，裂变截面，散射截面和移出截面。在main主程序中，内迭代中，根据所给文献划分的区域，首先求出各区第一群的通量，a，b，c，d，e，s的值；然后求出各区第二群的通量，a，b，c，d，e，s的值。如果
5、通量的迭代差值大于规定值，则迭代继续进行下去，直到迭代差值小于规定值。外迭代中，求得的有效增值因数的迭代差值大于规定值，则迭代继续，直到迭代差值小于规定值，迭代终止。将每次迭代的有效增殖因数写入文件。4.程序及计算结果program main implicit none integer,parameter:n=170 real,parameter:detx=1.0 real,parameter:dety=1.0 integer i,j,g,k real(8) keff,keff_old real(8) FL1,FL2,FL_OLD1,FL_OLD2,err1,err2,err real(8)
6、flux(n,n,2),flux_old(n,n,2),S(n,n,2),a(n,n,2),b(n,n,2),c(n,n,2),d(n,n,2),e(n,n,2) real(8) DI(n,n,2) real(8) sigmaA(n,n,2),sigmaF(n,n,2),sigmaS(n,n,1),sigmaR(n,n,1) real(8) Q,Q_OLD call cross(n,DI,sigmaA,sigmaF,sigmaS,sigmaR) err=1.0keff=1.01Q=0.0 do g=1,2 do j=1,170 do i=1,170 flux(i,j,g)=1.0 Q=Q+si
7、gmaF(i,j,g)*flux(i,j,g) end do end do end do do while(err1.0E-07) Q_OLD=Q keff_old=keff ! = !各区第一群的通量 !赋初值 FL1=0.0 FL_OLD1=0.0 err1=1.0 do while(err11.0E-05) do j=1,170 do i=1,170 flux_old(i,j,1)=flux(i,j,1) end do end do !第1行(j=1) !第一个网格(1,1,1)的通量(中心角) c(1,1,1)=-(2*DI(2,1,1)*DI(1,1,1)*dety)/(DI(2,1
8、,1)*detx+DI(1,1,1)*detx) d(1,1,1)=-(2*DI(1,2,1)*DI(1,1,1)*detx)/(DI(1,2,1)*dety+DI(1,1,1)*dety) e(1,1,1)=sigmaR(1,1,1)*detx*dety-c(1,1,1)-d(1,1,1) S(1,1,1)=1/keff*(sigmaF(1,1,1)*flux(1,1,1)+sigmaF(1,1,2)*flux(1,1,2)*detx*dety flux(1,1,1)=-(c(1,1,1)*flux(2,1,1)+d(1,1,1)*flux(1,2,1)-S(1,1,1)/e(1,1,1)
9、!中间各网格的通量(下边界) do i=2,169 b(i,1,1)=-(2*DI(i,1,1)*DI(i-1,1,1)*dety)/(DI(i,1,1)*detx+DI(i-1,1,1)*detx) c(i,1,1)=-(2*DI(i+1,1,1)*DI(i,1,1)*dety)/(DI(i+1,1,1)*detx+DI(i,1,1)*detx) d(i,1,1)=-(2*DI(i,2,1)*DI(i,1,1)*detx)/(DI(i,2,1)*dety+DI(i,1,1)*dety) e(i,1,1)=sigmaR(i,1,1)*detx*dety-b(i,1,1)-c(i,1,1)-d(
10、i,1,1) S(i,1,1)=1/keff*(sigmaF(i,1,1)*flux(i,1,1)+sigmaF(i,1,2)*flux(i,1,2)*detx*dety flux(i,1,1)=-(b(i,1,1)*flux(i-1,1,1)+c(i,1,1)*flux(i+1,1,1)+d(i,1,1)*flux(i,2,1)-S(i,1,1)/e(i,1,1) flux(1,i,1)=flux(i,1,1) end do !第2-70行 !中间各网格的通量 do j=2,69 do i=j,169 a(i,j,1)=-(2*DI(i,j,1)*DI(i,j-1,1)*detx)/(DI(
11、i,j,1)*dety+DI(i,j-1,1)*dety) b(i,j,1)=-(2*DI(i,j,1)*DI(i-1,j,1)*dety)/(DI(i,j,1)*detx+DI(i-1,j,1)*detx) c(i,j,1)=-(2*DI(i+1,j,1)*DI(i,j,1)*dety)/(DI(i+1,j,1)*detx+DI(i,j,1)*detx) d(i,j,1)=-(2*DI(i,j+1,1)*DI(i,j,1)*detx)/(DI(i,j+1,1)*dety+DI(i,j,1)*dety) e(i,j,1)=sigmaR(i,j,1)*detx*dety-a(i,j,1)-b(i
12、,j,1)-c(i,j,1)-d(i,j,1) S(i,j,1)=1/keff*(sigmaF(i,j,1)*flux(i,j,1)+sigmaF(i,j,2)*flux(i,j,2)*detx*dety flux(i,j,1)=-(a(i,j,1)*flux(i,j-1,1)+b(i,j,1)*flux(i-1,j,1)+c(i,j,1)*flux(i+1,j,1)+d(i,j,1)*flux(i,j+1,1)-S(i,j,1)/e(i,j,1) flux(j,i,1)=flux(i,j,1) end do end do do i=70,150 a(i,70,1)=-(2*DI(i,70,1
13、)*DI(i,69,1)*detx)/(DI(i,70,1)*dety+DI(i,69,1)*dety) b(i,70,1)=-(2*DI(i,70,1)*DI(i-1,70,1)*dety)/(DI(i,70,1)*detx+DI(i-1,70,1)*detx) c(i,70,1)=-(2*DI(i+1,70,1)*DI(i,70,1)*dety)/(DI(i+1,70,1)*detx+DI(i,70,1)*detx) d(i,70,1)=-(2*DI(i,71,1)*DI(i,70,1)*detx)/(DI(i,71,1)*dety+DI(i,70,1)*dety) e(i,70,1)=sigmaR(i,70,1)*detx*dety-a(i,70,1)-b(i,70,1)-c(i,70,1)-d(i,70,1) S(i,70,1)=1/keff*(sigmaF(i,70,1)*flux(i,70,1)+sigmaF(i,70,2)*flux(i,70,2)*detx*dety flux(i,70,1)=-(a(i,70,1)*flux(i,69,1)+b(i,70,1)*flux(i-1,70,1)+c(i,70,1)*flux(i+1,70,1)+d(i,70,1)*flux(i,
