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

AlexYoung

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

 
 
 

日志

 
 

二维圆柱(R,Theta)几何的单群SN解法  

2011-06-18 01:30:55|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

!               二维圆柱(R,Theta)几何的单群SN计算程序
!                          ——《核反应堆数值分析与模拟》期末大作业报告
!
!                                        /杨万奎  2011.6.15
!
!有参考解的计算网格为R方向5个网格,Theta方向4个网格,故在调试阶段做相同设
!置;采用S6计算,飞行方向坐标系的每个象限6个离散方向;在单群调试成功后试图
!加上扩散综合加速(DSA),这需要在增加网格点数的基础上才能体现出来加速效
!果。相关参数符号定义如下:
!
!                 R1=116cm---------裂变区半径;
!                 R2=145cm---------反射层外半径;
!                 NI、NJ-----------分别为r和Theta方向的离散网格数;
!                 NM---------------飞行方向坐标系一个坐标上的离散点数;
program SN2DCylinder
integer,parameter::NI=5,  NJ=4,  NM=6
real,parameter::pi=3.14159265, R1=116,  R2=145,  delta_theta=2*pi/NJ,  &
    P_weight=0.1666667
real,dimension(NI,NJ)::sigma_s,  sigma_t,  vsigma_f, phi,phiold=10,keff=1
real:: mu(4*NM),eta(4*NM), Psi_c(NI,NJ,4*NM),Psi_br(NI+1,NJ,4*NM), &
  Psi_btheta(NI,NJ+1,4*NM),Psi_bangle(NI,NJ,4*NM+1), &
  S(NI,NJ,4*NM),E,F,G,r_1,r_2,tempkeff, totalkeff,guiyi, &
  alpha(4*NM+1), V(NI,NJ), A(NI+1), B(NI),totalkeffold=1,totalkeff2=1, &
  fissionsold,fissions,tol=1,time_begin, time_end
integer::i,j,k,ds,ds1, m,countd=1,ct
integer::b1(3)=[1,2,4],   b2(3)=[3,5,6],   b3(3)=[8,10,11], b4(3)=[7,9,12], &
   b5(3)=[13,14,16],b6(3)=[15,17,18],b7(3)=[20,22,23],b8(3)=[19,21,24]
!=======================================================================
call cpu_time(time_begin)
!Parameters Initinization!
mu(1:24)=[-0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846, &
  -0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846]
eta(1:24)=[-0.2561429,-0.2663443,-0.6815646,-0.2561429,-0.6815646,-0.9320846, &
   -0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846, &
   0.2561429,0.6815646,0.2663443,0.9320846,0.6815646,0.2561429]
!Region I CrossSection
sigma_s(1:NI-1,:)=0.8804
sigma_t(1:NI-1,:)=1.0386
vsigma_f(1:NI-1,:)=0.2902
!Region II CrossSection
sigma_s(NI,:)=1.6245
sigma_t(NI,:)=1.6445
vsigma_f(NI,:)=0.0
!alpha
alpha(1)=0.0
alpha(4*NM+1)=0.0
do m=2,4*NM
 alpha(m)=alpha(m-1)-mu(m-1)*P_weight
end do
!A,B,V
do j=1,NJ
 do i=1,NI-1
 r_1=(i-1)*R1/(NI-1)
 r_2=i*R1/(NI-1)
 A(i)=r_1*delta_theta
 B(i)=r_2-r_1
 !注意教材公式(2.98)是错误的,应用theta角度间隔!
 V(i,j)=delta_theta*(r_2**2-r_1**2)/2
 end do
 V(NI,j)=delta_theta*(R2**2-R1**2)/2
end do
 A(NI)=R1*delta_theta
 A(NI+1)=R2*delta_theta
 B(NI)=R2-R1
!Vaccum boundary condition (exact and correct value)
Psi_br(NI+1,:,1:NM*(NM+2)/8)=0.0
Psi_br(NI+1,:,NM*(NM+2)/4+1 : 3*NM*(NM+2)/8)=0.0
!first iteration setup (source & psi_btheta())
!S=1
psi_c=1.0
do j=1,NJ
do i=1,NI
 phi(i,j)=2.0*4.0*pi*sum(psi_c(i,j,:))*P_weight/8.0
 do m=1,4*NM
 S(i,j,m)=(sigma_s(i,j)+vsigma_f(i,j)/totalkeff2)*phi(i,j)/(4.0*pi)
 end do
end do
end do
fissionsold=sum(vsigma_f(:,:)*phi(:,:)*V(:,:))
!=======================================================================
do while(countd>0)
!do while(tol>=1e-6)
!do ds1=1,3000
 !----------------------第III象限--------------------
 do ds=1,10  !让其迭代赋值10次
 do j=NJ,1,-1
 do i=NI,1,-1
  !求m=1,2,4方向
  do k=1,3
  m=b1(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2*abs(eta(m))*B(i)
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j+1,m)+V(i,j)* &
   S(i,j,m))/(E+F+V(i,j)*sigma_t(i,j)) !G=0
  Psi_bangle(i,j,m+1)=max(Psi_c(i,j,m), 0.0)
  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
     end do
  !求m=3,5,6方向
  do k=1,3
  m=b2(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j+1,m)+G* &
   Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
  end do
 end do
 end do
  Psi_btheta(:,NJ+1,1:NM*(NM+2)/8)=Psi_btheta(:,1,1:NM*(NM+2)/8)
  !周期边界条件
 end do

 !传给第IV象限 Psi_bangle
 do j=NJ,1,-1
 do i=NI,1,-1
  Psi_bangle(i,j,10)=Psi_bangle(i,j,7) !6->10
  Psi_bangle(i,j,7)=Psi_bangle(i,j,2)  !1->7
  Psi_bangle(i,j,8)=Psi_bangle(i,j,4)  !3->8
 end do
 end do

 !第IV象限计算初值 Psi_br,由圆柱中心点处的反射边界条件给出!(必须要的)
 !看教材编号图上各点的对称关系!
 do j=1,NJ
  Psi_br(1,j,7)=Psi_br(1,j,1) !1->7
  Psi_br(1,j,9)=Psi_br(1,j,2) !2->9
  Psi_br(1,j,8)=Psi_br(1,j,3) !3->8
  Psi_br(1,j,12)=Psi_br(1,j,4)!4->12
  Psi_br(1,j,11)=Psi_br(1,j,5)!5->11
  Psi_br(1,j,10)=Psi_br(1,j,6)!6->10
 end do

 !m=7,8,9,10,11,12离散方向
 do ds=1,10
 do j=NJ,1,-1
 do i=1,NI                                                                                                                                                                                                                                                                 
  !先求离散方向m=8,10,11(注必需先求10,后求11)
  do k=1,3
  m=b3(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j+1,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
  end do
     !后求离散方向m=7,9,12
  do k=1,3
  m=b4(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j+1,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  !Psi_bangle(i,j,m+1)=max(2*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
  end do
 end do
 end do
  Psi_btheta(:,NJ+1,NM*(NM+2)/8+1:NM*(NM+2)/4)= &
   Psi_btheta(:,1,NM*(NM+2)/8+1:NM*(NM+2)/4)
 end do


 !-----------------------第II象限--------------------
 !求13,14,15,16,17,18离散方向
 do ds=1,10
 do j=1,NJ
 do i=NI,1,-1
  !求m=13,14,16方向的
  do k=1,3
  m=b5(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j,m)+ &
   V(i,j)*S(i,j,m))/(E+F+V(i,j)*sigma_t(i,j)) !G=0
  Psi_bangle(i,j,m+1)=max(Psi_c(i,j,m), 0.0)
  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m), 0.0)
  end do
     !求m=15,17,18方向的
  do k=1,3
  m=b6(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m), 0.0)
  end do
 end do
 end do
  Psi_btheta(:,1,NM*(NM+2)/4+1:3*NM*(NM+2)/8)= &
   Psi_btheta(:,NJ+1,NM*(NM+2)/4+1:3*NM*(NM+2)/8)
   !注这与III、IV的不同
 end do


 !传给第I卦限 Psi_bangle
 do j=NJ,1,-1
 do i=NI,1,-1
  Psi_bangle(i,j,22)=Psi_bangle(i,j,19) !18->22
  Psi_bangle(i,j,19)=Psi_bangle(i,j,14) !13->19
  Psi_bangle(i,j,20)=Psi_bangle(i,j,16) !15->20
 end do
 end do


 !传给第I卦限 Psi_br
 do j=1,NJ
  Psi_br(1,j,19)=Psi_br(1,j,13) !13->19
  Psi_br(1,j,21)=Psi_br(1,j,14) !14->21
  Psi_br(1,j,20)=Psi_br(1,j,15) !15->20
  Psi_br(1,j,24)=Psi_br(1,j,16) !16->24
  Psi_br(1,j,23)=Psi_br(1,j,17) !17->23
  Psi_br(1,j,22)=Psi_br(1,j,18) !18->22
 end do

 !-----------------------第I象限---------------------
 !求m=19,20,21,22,23,24的离散方向
 do ds=1,10
 do j=1,NJ
 do i=1,NI                                                                                                                                                                                                                                                                 
 !先求离散方向m=20,22,23(注必需先求22,后求23)!
  do k=1,3
  m=b7(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m), 0.0)
  end do
  !后求离散方向m=19,21,24
  do k=1,3
  m=b8(k)
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  !Psi_bangle(i,j,m+1)=max(2*Psi_c(i,j,m)-Psi_bangle(i,j,m),0.0)
  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m),0.0)
  end do
 end do
 end do
  Psi_btheta(:,1,3*NM*(NM+2)/8+1:NM*(NM+2)/2)= &
   Psi_btheta(:,NJ+1,3*NM*(NM+2)/8+1:NM*(NM+2)/2)
 end do

!求标通量及keff(两类keff:临界keff和真实keff)
countd=0
do j=1,NJ
do i=1,NI
 phi(i,j)=2.0*4.0*pi*sum(Psi_c(i,j,:))*P_weight/8.0
 tempkeff=phi(i,j)/phiold(i,j)
 !这里的收敛标准非常重要,为什么要与1进行比较?因为对裂变源项
 !除以keff之后,即认为系统是临界的,这就是我们可以用临界方程解非
 !临界稳态问题的理论基础......
 if (abs(tempkeff-1)/1>=1e-6) then
  countd=countd+1
 end if
 keff(i,j)=tempkeff
end do
end do
totalkeff=sum(phi)/sum(phiold)
phiold=phi
!totalkeff-----为临界方程中加入keff本征值之后的临界系统的有效增值系数,
!              可用于表征通量收敛。
!totalkeff2----为真实的有效增值系数,可用于表征keff收敛。
fissions=sum(vsigma_f(:,:)*phi(:,:)*V(:,:))
totalkeff2=fissions/fissionsold*totalkeffold
tol2=abs(totalkeff2-totalkeffold)/totalkeff2
totalkeffold=totalkeff2
fissionsold=fissions

ct=ct+1
write(*,*) "iteration times=",ct,"   totalkeff=", totalkeff2,totalkeff
!----------------------------------------------------------------------
!Calculate Source
do i=1,NI
do j=1,NJ
do m=1,4*NM
S(i,j,m)=(sigma_s(i,j)+vsigma_f(i,j)/totalkeff2)*phi(i,j)/(4.0*pi)
end do
end do
end do

end do

!================================输出结果===============================
write(*,'(f15.10)') phi(:,1)
open(unit=11,file='Results.txt')
open(unit=12,file='scalar_flux.txt')
write(12,'(f12.6)') ((phi(i,j),i=1,NI),j=1,NJ)

write(11,*)"scalar flux"
write(11,'(f12.6)')((phi(i,j),i=1,NI),j=1,NJ)

write(11,*)"Psi_c"
write(11,'(f12.6)') (((Psi_c(i,j,k),i=1,NI),j=1,NJ),k=1,4*NM)

write(11,*)"keff"
write(11,'(f12.6)') totalkeff2

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

  评论这张
 
阅读(393)| 评论(0)
推荐 转载

历史上的今天

评论

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

页脚

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