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

AlexYoung

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

 
 
 

日志

 
 

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

2011-04-09 22:01:43|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

! 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
!                                               9th April 2011
program SN1DS
 implicit none
 integer:: k, m, i, ct, bad_grid_number=1
 integer,parameter:: N=6, kk=40*N
 real(8),dimension(2*kk+1):: keff=1,phi, phiold=10, src=1e13 !phi--scalar flux
 real(8),dimension(N)     :: mu, omega, l
 real(8),dimension(2*kk+1,N):: psi !psi--angular flux
 real(8),parameter:: deltax=66.0053/kk, 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]
 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(2*kk+1,m)=0 !Right boundary condition of vacuum.
   do k=kk,1,-1
    psi(2*k-1,m)=(2+l(m))/(2-l(m))*psi(2*k+1,m)-2*l(m)/(2-l(m))*src(2*k)/sigma_t
    psi(2*k,m)=(psi(2*k+1,m)+psi(2*k-1,m))/2
   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,kk
    psi(2*k+1,m)=(2-l(m))/(2+l(m))*psi(2*k-1,m)+2*l(m)/(2+l(m))*src(2*k)/sigma_t
    psi(2*k,m)=(psi(2*k+1,m)+psi(2*k-1,m))/2
   enddo
  enddo
  !-------------------------------------------------------------------------------------
  !Upgrade source, including both scattering source and fission source.
  do i=1,kk
   src(2*i)=(sigma_s+nusigma_f)/4*sum(omega(:)*(psi(2*i-1,:)+psi(2*i+1,:)))
  enddo
  !-------------------------------------------------------------------------------------
  !Convergence judgement by counting the number of bad grids not satisfying criteria.
  bad_grid_number=0
  do i=1,2*kk+1
   phi(i)=sum(omega(:)*psi(i,:))
   tempkeff=phi(i)/phiold(i)
   if (abs(tempkeff-keff(i))/keff(i)>=1e-6) then
    bad_grid_number = bad_grid_number + 1
   endif
   keff(i)=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(kk/2+1)/phi(1),phi(kk+1)/phi(1),phi(3*kk/2+1)/phi(1), &
&  phi(2*kk+1)/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'
end program

 

注:本程序为Fortran90/95标准,在CVF6.5中运行通过。相关的试验分析报告另撰文阐述,真正编写完整的程序之后才会真正理解整个SN的方法要点之外的很多细节的问题:哪些量是已知的?那些是需要中间先行求出的?哪些量是需要自己设定的?源如何赋初值及所赋的初值有什么意义?收敛判断标准是什么?在此基础上才谈得上负通量修正、如何加速等后续问题。

  评论这张
 
阅读(420)| 评论(3)
推荐 转载

历史上的今天

评论

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

页脚

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