蒙特卡罗mncp.doc

蒙特卡罗方法在实验核物理中的应用
结课论文
学院：核工程技术学院
专业：核技术
指导老师：吴和喜
2011年12月9号

作业一

1. HPGe探测器，35cm*35cm,面源r=20cm,两中心相距4.5cm，求面源射线为1.0Mev，3.0Mev，0.7Mev（各自的几率为0.25，0.15，0.6）时，HPGe探测能谱？
一、建立坐标系
以闪烁体中轴线为z轴建立柱坐标系，闪烁体位于0<z<4且r<2的区域。Cs电源位于z=-2，r=0点。
二、编程思想
一个粒子的状态用表示。其中E为粒子能量，为粒子的位置坐标，为粒子的运动方向。其中为粒子运动方向与 z 轴的夹角的余弦，为粒子运动方向在 x y 平面上投影的方位角。
对于一个粒子经历以下过程：
1． 源抽样：
由于是点源，能量和位置的分布均为δ函数，抽样得到 。
抽样得：。
2．到达闪烁体：
如果1/2，无法到达闪烁体，丢弃，返回源抽样重新产生粒子。粒子进入闪烁体的瞬间，状态为：其他量不变。
3．输运过程：
抽样得到到下次碰撞的距离，根据当前粒子状态中的算出下次碰撞的坐标，如果不在闪烁体区域（0<z<45且r<20），输运过程结束。如果在，继续－>
碰撞有两种可能：光电效应和康普顿散射。
根据粒子当前的能量，（由NaI(Tl)闪烁体宏观界面数据）线性插值确定它的光电效应截面和康普顿散射界面。抽样得到本次反应的类型。
如果光电效应，E=0，输运过程结束。如果康普顿效应，抽样获取碰撞后的能量和运动方向（康普顿散射的能量分布密度函数知道，具体抽样方法参考讲义。）
如果E<1KeV输运过程结束，反之，重复本过程直到输运过程结束。
4．记录与统计：
记录末态能量，计算沉积能量，考虑到测量系统分辨率，多道记录能量为沉积能量的高斯展宽。

记录能量。其中，。由标准正态分布抽样得到。
五：程序如下：
count=input('input the count:');%输入模拟粒子数

sigmaedata=[28370,13845,6908,2555,1223,2602,1925,905.3,479.7,164.5,74.24,23.86,66.60,36.62,22.29,9.978,5.298,1.668,0.7378,0.2361,0.1099,0.06211,0.03939,0.02030];

sigmacdata=[0.0220,0.0393,0.0568,0.0904,0.1209,0.1479,0.1722,0.2136,0.2480,0.3092,0.3486,0.3932,0.4153,0.4268,0.4319,0.4291,0.4215,0.3969,0.3691,0.3269,0.2944,0.2709,0.2512,0.2209];

Edata=[1,1.5,2,3,4,5,6,8,10,15,20,30,40,50,60,80,100,150,200,300,400,500,600,800];%截面数据

channel=zeros(1,ceil(662/5)+10);%多道数组

nget=0;%探测到的总计数

ntotal=0;%进入探测器的总计数

for ii=1:count%count个粒子循环

collidetime=0;%当前粒子碰撞次数
E0=input('input the energ0:');%输入模拟面源的一个能量
E1=input('input the energy1:');%输入模拟面源的另一个能量
E2= input('input the energy2:') %输入模拟面源的另另一个能量

%粒子状态初始化

E0=622;

E=E0;

z=-2;

r=0;

theta=2*pi*rand(1);% 源抽样，z,r,theta坐标

miu=2*rand(1)-1;

fai=2*pi*rand(1);%方向角抽样

if miu<cos(pi/4) %是否能进入探测器

continue;

else

z=0;

r=2*sqrt(1-miu^2)/miu;

theta=fai;

end

while E>1%一个粒子在闪烁体中的输运过程

sigmae=interp1(Edata,sigmaedata,E,'linear');

sigmac=interp1(Edata,sigmacdata,E,'linear');

sigmat=sigmae+sigmac;%线性插值得到截面数据

L=-log(rand(1))/sigmat;%下次碰撞的距离

%计算下次碰撞位置坐标

rnew=sqrt(r^2+L^2*(1-miu^2)+2*r*L*sqrt(1-miu^2)*cos(fai-theta));

z=z+L*miu;

cdth=(rnew^2+r^2-L^2*(1-miu^2))/2/r/rnew;

sdth=L*sqrt(1-miu^2)*sin(fai-theta)/rnew;

dtheta=asin(sdth);

if cdth<0

dtheta=pi-dtheta;

end

theta=theta+dtheta;

r=rnew;

if(r>2)|(z>=4)|(z<0)%判断是否在闪烁体内

break;

else

collidetime=collidetime+1;

end

if rand(1)<(sigmae/sigmat) %光电效应

E=0;

else%康普顿散射

alpha=E/511;

flag=0;

while flag==0

if rand(1)<=27/(4*alpha+29)

x=(1+2*alpha)/(1+2*alpha*rand(1));

if rand(1)<=0.5*((alpha+1-x/alpha)^2+1)

flag=1;

end

else

x=1+2*alpha*rand(1);

if rand(1)<=27/4*((x-1)^2)/x^3

flag=1;

end

end

end

E=E/x;

alphat=alpha/x;

miuL=1-1/alphat+1/alpha;%散射后的方向

a=miuL;

b=sqrt(1-a^2);

randangle=2*pi*rand(1);

miunew=a*miu+b*sqrt(1-miu^2)*cos(randangle);

sdf=b*sin(randangle)/sqrt(1-miunew^2);

cdf=(a-miu*miunew)/sqrt(1-miu^2)/sqrt(1-miunew^2);

sfn=sdf*cos(fai)+cdf*sin(fai);

cfn=cdf*cos(fai)-sdf*sin(fai);

fainew=asin(sfn);

if(cfn<0)

fainew=pi-fainew;

end

fai=fainew;

miu=miunew;

end%conputon/electtron

%记录结果

end%while

ntotal=ntotal+1;%进入探测器的总计数

if collidetime>0

nget=nget+1;%探测到的计数

Edown=E0-E;%沉积能量

FWHM=0.01+0.05*sqrt(Edown+0.4*Edown^2);

delta=0.4247*FWHM;

Eget=Edown+delta*randn(1);%记录能量

num=ceil(Eget/5); %寻道址

if num~=0

channel(num)=channel(num)+1;

end

end

end%for count

plot(channel);

fprintf(1,'探测效率为%1.5f',nget/ntotal);

五、模拟结果
在计算了四百万个粒子后，得到HGPE探测能谱。
1) 探测效率与峰总比
N=486597――进入闪烁体的总粒子数
Nm=279832――探测到的总粒子数
NP=89752――全能峰的总计数（全能峰半高宽以内道的计数之和）
探测效率：
峰总比：
2:能谱图：

图2 全谱图

3：线性插值求得半高宽

图3 全能峰局部展宽

作业二
热态不同水铀比.硼浓度的栅元计算
1实验目的：在一已知栅元中，在不同的水铀比.硼浓度下，冷却剂吸收效应主导地位不同，所处区域不同（慢化区，过慢化区），为了弄清在不同水铀比.硼浓度下的冷却剂吸收效应主导地位.慢化转折点以及中子的利用率各种关系，采用MCNP软件编写代码，模拟在已知材料下的水铀比.硼浓度不同条件情况下冷却剂吸收效应主导地位。
2实验原理：MCNP是一个通用的蒙特卡罗粒子运输程序，可用于中子，光子，电子运输，可以计算临界系统的特征值，可以处理材料的任意三维结构，栅元的边界可以是一阶，二阶，以及四阶椭圆曲面，使用逐点的截面数据，对于中子可以计算截面中的所有反应，热中子自由气体和S两种模型描述，对于光子考虑了相干和非相干散射，光电子吸收后的荧光发射，，以及韧致辐射。对于电子，使用连续的Slowing down模型，包括正电子，韧致辐射，但不考虑外部或自身的感应场，MCNP可以使用强大的通用源，临界源和曲面源;可以空间几何结构和输出数据作图；有丰富的方差衰减技术：方便的记数结构，广泛的截面数据库

3栅元描述：栅距13.3mm，燃料棒外径：8.46mm，
锆合金包壳外径：8.6mm，锆合金包壳内径：10mm，
边界条件：栅元边界全反射条件。
4几何描述：
材料：1 燃料棒材料为UO2，在不同富密度比原子密度比不同。
原子密度比为：

原子密度比
UO2富密度
U235
U238
O16
1.8
0.01875426
0.981745956
2
2.4
0.02435689
0.976564617
2
3.1
0.03140184
0.968591827
3

2锆合金包壳材料为：Sn 百分之一点五，Fe百分之二 ，氧百分之一
3 冷却剂材料：含硼水，硼浓度为1000ppm，原子密度比为H:2，O:1 B10:0.00033，B11:0.00133，温度定义为290，质量密度为0.7457
燃料棒及包壳尺寸已知，但栅距随水铀比的不同而改变，其具体数值如下：
水铀比
栅距
水铀比
栅距
0.5
10.31730
3
15.68380
0
11.59110
3.5
16.54960
1.5
12.73820
4
17.37230
1.762
13.3
5
18.91060
2
13.79020
6
20.33290
2.5
14.76740
7
21.66190

5：材料描述：
冷却剂温度为290°C,密度为0.745克每立方米厘米，燃料富密集固定为百分之3.1，硼浓度及相应的B10，B11的原子密度列表如下：

B(ppm)
B10(H:2.O:1)
B11(H:2,O:1)
0
0.000000
0.000000
200
0.000065
0.000267
300
0.000098
0.000107
400
0.000131
0.000533
500
0.000164
0.000668
700
0.000231
0.000935
1000
0.000429
0.001335
1300
0.000660
0.001738
2000
0.000785
0.002676

6MCNP程序：

1 1 -10.4 -1 -9 10 imp:n=1
2 0 1 -2 -9 10 imp:n=1
3 2 -6.55 2 -3 -9 10 imnp:n=1
4 3 -0.745121551 3 4 -5 6 -7 -9 10 imp:n=1
5 0 -4:5 -6:7:9:-10 imp:n=0

1 cz 0.4215
2 cz 0.4300
3 cz 0.5000
*4 px -0.5158
*5 px 0.5158
*6 py -0.5158
*7 py 0.5158
*9 pz 1.0
*10 pz -1.0

m1 9225 0.0314081743 92238 0.9685918257 8016 2
m2 50000 -0.015 2600 -0.002 2400 -0.001 8016 -0.001 40000 -0.981
m3 1001 2 8016 1 5010 0.000657421
5011 0.0002671395 $B-200ppm
kcode 5000 1.0 50 100
ksrc 0 0 0
print

7:结果分析：将结果画成表，可知：

水铀比
栅题
硼浓度(ppm)
0
200
300
400
500
700
1000
1300
0.3
10.3
1.07
1.06
1.066
1.065
1.055
1.055
1.1853
1.0432
1
11.5
1.25
1.23
1.006
1.222
1.271
1.200
1.228
1.163
1.5
12.7
1.33
1.30
1.230
1.288
1.280
1.256
1.233
1.198
1.762
13.3
1.36
1.33
1.303
1.307
1.302
1.269
1.232
1.199
2
13.7
1.38
1.34
1.321
1.316
1.297
1.274
1.221
1.195
2.3
14.7
1.40
1.36
1.331
1.324
1.305
1.272
1.22
1.176
3
15.6
1.41
1.36
1.343
1.321
1.298
1.258
1.210
1.149
3.5
16.5
1.41
1.34
1.331
1.308
1.281
1.236
1.172
1.117
4
17.3
1.41
1.31
1.32
1.291
1.216
1.212
1.114
1.085
5
18.9
1.39
1.32
1.285
1.251
1.218
1.163
1.087
1.018

8：目的意义分析：
经过软件模拟得到数据，可得到在已知栅元条件下，不同水铀比.硼浓度对冷却剂吸收效应主导地位所处区域有如下关系：
1在硼浓度为0~1300ppm时，处于慢化区，而高于此浓度时，则会处于过慢化区，为了提高中子的利用率，核电站都将水铀比设置在欠慢化区。
2 由硼浓度为200，300 ，400的数据及曲线可知，硼浓度越大，其微分价值的绝对值越小
3 由图可知，随着硼浓度的增大，欠慢化区与过慢化区的转折点向左移动，因为硼的加入对冷却剂的慢化性能几乎无影响，但却增大了冷却剂的吸收性能，使平衡向慢化区倾斜，故转折点左移。
4 硼浓度越大，栅元德反应性越小，反应出硼的对中子的吸收效应。