注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

AlexYoung

做好大家都必须要做的事,发展自己感兴趣的事,结束圣人所不齿的事...

 
 
 

日志

 
 

一维平几何NTE(中子输运方程)的SN解法(离散纵标法)【new】  

2011-04-10 10:44:25|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

改进之处:

此次改进主要是在下标索引上的改进,即空间网格的划分上的变化:

0=x(1)<x(3/2)<x(2)<...<x(rk-1/2)<x(rk)=a;

推导上与前述完全相同,只是下标的变化,主要是为了使其更便于编程,不至于出现2k+1之类的情形。在推导过程中可以将所有的带二分之一的节点值都用菱形关系式代替,从而显式地不出现在编程中。但是这依然是采用的中心差分格式,精度与前述是一样的,不过个人觉得下标索引上更加赏心悦目了,哈哈。同时,扩展至S16的离散度。

下面附上改进之后的程序源码:

! SN(Discrete Ordinates) Method of NET numerical solver in One-Dimentional Slab Geometry.
! Isotropical scattering and fission neutrons are surposed here.
! Anyone who is interested can edit or correct this version.

!                                               Programmed by Alex Young
!                                               CAEP, Beijing, China
!            10th April 2011
program SN1DS
 implicit none
 integer:: k, m, i, ct, bad_grid_number=1
 integer,parameter:: N=6, rk=40*N+1 !rk--rightmost grid number
 real(8),dimension(N):: mu, omega, l
 real(8),dimension(rk):: keff=1,phi, phiold=10 !phi--scalar flux
 real(8),dimension(rk-1):: src=1e13
 real(8),dimension(rk,N):: psi !psi--angular flux
 real(8),parameter:: deltax=66.0053/(rk-1), sigma_t=0.05, sigma_s=0.03, nusigma_f=0.0225
 real(8):: tempkeff, time_begin, time_end

 call cpu_time(time_begin)
 !generate the quadrature set
 !select number of dicreted directions, which usually tend to be an even integer.
 select case(N)
  case(2)
   mu(1:2)=[-0.57735, 0.57735]
   omega(1:2)=[1.00000, 1.00000]
  case(4)
   mu(1:4)=[-0.86113, -0.33998, 0.33998, 0.86113]
   omega(1:4)=[0.34785, 0.65214, 0.65214, 0.34785]
  case(6)
   mu(1:6)=[-0.93246, -0.66120, -0.23861, 0.23861, 0.66120, 0.93246]
   omega(1:6)=[0.17132, 0.36076, 0.46791, 0.46791, 0.36076, 0.17132]
  case(8)
   mu(1:8)=[-0.96028, -0.79666, -0.52553, -0.18343, 0.18343, 0.52553, 0.79666,0.96028]
   omega(1:8)=[0.10122, 0.22238, 0.31370, 0.36268, 0.36268, 0.31370, 0.22238, 0.10122]

  case(16)
   mu(1:16)=[-0.98940, -0.94457, -0.86563, -0.75540, -0.61787, -0.45801, -0.28160, -0.09501,&
&        0.09501,  0.28160,  0.45801,  0.61787,  0.75540,  0.86563,  0.94457,  0.98940]
   omega(1:16)=[0.02715, 0.06225, 0.09515, 0.12462, 0.14959, 0.16951, 0.18260, 0.18945,&
&       0.18945, 0.18260, 0.16951, 0.14959, 0.12462, 0.09515, 0.06225, 0.02715]
 end select
 
 !Source iteration or power iteration, in which the inner iteration is launched.
 !This kind of convergence criteria is more efficient, which should be recommended.
 do while (bad_grid_number > 0)
  !-----------------------------------------------------------------------------------
  !Start calculation from mu<0, scan grids from right side to left side.
  do m=1,N/2
   l(m)=deltax*sigma_t/mu(m)
   psi(rk,m)=0 !Right boundary condition of vacuum.
   do k=rk-1, 1, -1
    psi(k,m)=(2+l(m))/(2-l(m))*psi(k+1,m) - 2*l(m)/(2-l(m))*src(k)/sigma_t
   enddo
  enddo
  !------------------------------------------------------------------------------------
  !Calculate angular flux with mu>0, scan grids from left to right.
  do m=N/2+1,N
   l(m)=deltax*sigma_t/mu(m)
   psi(1,m)=psi(1,N+1-m) !Left boundary condition of reflection.
   do k=1,rk-1
    psi(k+1,m)=(2-l(m))/(2+l(m))*psi(k,m)+2*l(m)/(2+l(m))*src(k)/sigma_t
   enddo
  enddo
  !-------------------------------------------------------------------------------------
  !Upgrade source, including both scattering source and fission source.
  !src(k)--source in the grid center from k to k+1, ie. the q(k+1/2) in derivation.
  do k=1,rk-1
   src(k)=(sigma_s+nusigma_f)/4*sum(omega(:)*(psi(k,:)+psi(k+1,:)))
  enddo
  !-------------------------------------------------------------------------------------
  !Convergence judgement by counting the number of bad grids not satisfying criteria.
  bad_grid_number=0
  do k=1,rk
   phi(k)=sum(omega(:)*psi(k,:))
   tempkeff=phi(k)/phiold(k)
   if (abs(tempkeff-keff(k))/keff(k)>=1e-6) then
    bad_grid_number = bad_grid_number + 1
   endif
   keff(k)=tempkeff
  enddo
  phiold=phi
  !-------------------------------------------------------------------------------------
  print *, 'bad_grid_number=',bad_grid_number
  ct=ct+1
 enddo
 print *,'cycle times=',ct

 !open(unit=10,file="result.txt")
 !DO m=1,N
 ! write (10,fmt='(201f10.4)') phi(:)/phi(1)
 write (*,fmt='(4f12.8)') phi((rk-1)/4+1)/phi(1),phi((rk-1)/2+1)/phi(1),phi(3*(rk-1)/4+1)/phi(1), &
&  phi(rk)/phi(1),maxval(keff),minval(keff)
 !end do
 write(*,*) "Congratulation! The data have been saved successfully!"

 call cpu_time(time_end)
 print *,'clapsed time is ',time_end - time_begin,'    secondes'
 print *,'Press return carriage key to exit...'
 read(*,*) !Wait for a return carriage key to exit current window.
 stop
end program

  评论这张
 
阅读(372)| 评论(1)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017