蒙特卡罗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轴建立柱坐标系，闪烁体位于0z4且r2的区域。Cs电源位于z=-2，r=0点。 二、编程思想一个粒子的状态用表示。其中E为粒子能量，为粒子的位置坐标，为粒子的运动方向。其中为粒子运动方向与 z 轴的夹角的余弦，为粒子运动方向在 x y 平面上投影的方位角。 对于一个粒子经历以下过程：1 源抽样：由于是点源，能量和位置的分布均为函数，抽样得到 。抽样得：。2到达闪烁体：如果1/2，无法到达闪烁体，丢弃，返回源抽样重新产生粒子。粒子进入闪烁体的瞬间，状态为：其他量不变。 3输运过程：抽样得到到下次碰撞的距离，根据当前粒子状态中的算出下次碰撞的坐标，如果不在闪烁体区域（0z45且r 碰撞有两种可能：光电效应和康普顿散射。根据粒子当前的能量，（由NaI(Tl)闪烁体宏观界面数据）线性插值确定它的光电效应截面和康普顿散射界面。抽样得到本次反应的类型。如果光电效应，E=0，输运过程结束。如果康普顿效应，抽样获取碰撞后的能量和运动方向（康普顿散射的能量分布密度函数知道，具体抽样方法参考讲义。）如果E1KeV输运过程结束，反之，重复本过程直到输运过程结束。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 miu1%一个粒子在闪烁体中的输运过程sigmae=interp1(Edata,sigmaedata,E,linear);sigmac=interp1(Edata,sigmacdata,E,linear);sigmat=sigmae+sigmac;%线性插值得到截面数据L=-log(rand(1)/sigmat;%下次碰撞的距离%计算下次碰撞位置坐标rnew=sqrt(r2+L2*(1-miu2)+2*r*L*sqrt(1-miu2)*cos(fai-theta);z=z+L*miu;cdth=(rnew2+r2-L2*(1-miu2)/2/r/rnew;sdth=L*sqrt(1-miu2)*sin(fai-theta)/rnew;dtheta=asin(sdth); if cdth2)|(z=4)|(z0)%判断是否在闪烁体内 break; else collidetime=collidetime+1; endif 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)/x3 flag=1; end end end E=E/x; alphat=alpha/x; miuL=1-1/alphat+1/alpha;%散射后的方向 a=miuL; b=sqrt(1-a2); randangle=2*pi*rand(1); miunew=a*miu+b*sqrt(1-miu2)*cos(randangle); sdf=b*sin(randangle)/sqrt(1-miunew2); cdf=(a-miu*miunew)/sqrt(1-miu2)/sqrt(1-miunew2); sfn=sdf*cos(fai)+cdf*sin(fai); cfn=cdf*cos(fai)-sdf*sin(fai); fainew=asin(sfn); if(cfn0nget=nget+1;%探测到的计数Edown=E0-E;%沉积能量FWHM=0.01+0.05*sqrt(Edown+0.4*Edown2);delta=0.4247*FWHM;Eget=Edown+delta*randn(1);%记录能量num=ceil(Eget/5); %寻道址 if num=0 channel(num)=channel(num)+1; endendend%for countplot(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富密度U235U238O161.80.018754260.98174595622.40.024356890.97656461723.10.031401840.9685918273 2锆合金包壳材料为：Sn 百分之一点五，Fe百分之二 ，氧百分之一 3 冷却剂材料：含硼水，硼浓度为1000ppm，原子密度比为H:2，O:1 B10:0.00033，B11:0.00133，温度定义为290，质量密度为0.7457燃料棒及包壳尺寸已知，但栅距随水铀比的不同而改变，其具体数值如下：水铀比栅距水铀比栅距0.510.31730315.68380011.591103.516.549601.512.73820417.372301.76213.3518.91060213.79020620.332902.514.76740721.661905：材料描述： 冷却剂温度为290C,密度为0.745克每立方米厘米，燃料富密集固定为百分之3.1，硼浓度及相应的B10，B11的原子密度列表如下：B(ppm)B10(H:2.O:1)B11(H:2,O:1)00.0000000.0000002000.0000650.0002673000.0000980.0001074000.0001310.0005335000.0001640.0006687000.0002310.00093510000.0004290.00133513000.0006600.00173820000.0007850.0026766MCNP程序：1 1 -10.4 -1 -9 10 imp:n=12 0 1 -2 -9 10 imp:n=13 2 -6.55 2 -3 -9 10 imnp:n=14 3 -0.745121551 3 4 -5 6 -7 -9 10 imp:n=15 0 -4:5 -6:7:9:-10 imp:n=01 cz 0.42152 cz 0.43003 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.0m1 9225 0.0314081743 92238 0.9685918257 8016 2m2 50000 -0.015 2600 -0.002 2400 -0.001 8016 -0.001 40000 -0.981m3 1001 2 8016 1 5010 0.000657421 5011 0.0002671395 $B-200ppmkcode 5000 1.0 50 100 ksrc 0 0 0print 7:结果分析：将结果画成表，可知：水铀比栅题硼浓度(ppm)0200300400500700100013000.310.31.071.061.0661.0651.0551.0551.18531.0432111.51.251.231.0061.2221.2711.2001.2281.1631.512.71.331.301.2301.2881.2801.2561.2331.1981.76213.31.361.331.3031.3071.3021.2691.2321.199213.71.381.341.3211.3161.2971.2741.2211.1952.314.71.401.361.3311.3241.3051.2721.221.176315.61.411.361.3431.3211.2981.2581.2101.1493.516.51.411.341.3311.3081.2811.2361.1721.117417.31.411.311.321.2911.2161.2121.1141.085518.91.391.321.2851.2511.2181.1631.0871.0188：目的意义分析： 经过软件模拟得到数据，可