蒙特卡罗方法与MCNP程序入门.docx

蒙特卡罗方法与MCNP程序入门.docx

《蒙特卡罗方法与MCNP程序入门.docx》由会员分享，可在线阅读，更多相关《蒙特卡罗方法与MCNP程序入门.docx（167页珍藏版）》请在冰豆网上搜索。

蒙特卡罗方法与MCNP程序入门目录

蒙特卡罗方法

与MCNP程序

入门

-i-

蒙特卡罗方法与MCNP程序入门目录

目录

第1章蒙特卡罗方法简介 1

§1.1　蒙特卡罗方法的基本原理 1

§1.2蒙特卡罗方法的解题手续和特点 6

§1.3用蒙特卡罗方法模拟粒子输运 7

第2章MCNP程序入门 9

§2.1MCNP简介 9

§2.2MCNP程序的组成及特点 11

§2.3MCNP4C程序安装、运行与源程序编译 13

§2.4MCNP输入文件 16

第3章MCNP程序中的几何构建 23

§3.1基础知识 23

§3.2几何描述卡 26

§3.3有效地构建几何 30

第4章MCNP程序的数理基础 33

§4.1物理 33

§4.2记数 41

§4.3减小方差技巧 47

第5章MCNP程序中的数据卡 56

§5.1问题类型卡 56

§5.2栅元参数和曲面参数卡 56

§5.3源的描述 65

§5.4记数方式的指定 74

§5.5材料的指定 88

§5.6能量和热处理方式的指定 89

§5.7问题截断条件 93

§5.8外围卡 95

§5.9MCNP输入文件综述 97

第6章经验 101

§6.1一般应用步骤 101

§6.2需注意的问题 101

第7章应用实例 103

§7.1医学物理中的应用 103

§7.2反应堆物理计算中的应用 129

附录 157

连续能量中子截面库ENDL851数据目录 157

中子热截面库BMCC1数据目录 158

离散中子截面库D91数据目录 159

光子截面库MCPLIB1数据目录 160

特殊材料S（a,b）热截面库TMCC1数据目录 162

参考文献 164

-I-

蒙特卡罗方法与MCNP程序入门

第1章蒙特卡罗方法简介

蒙特卡罗方法，又称随机抽样方法，是一种与一般数值计算方法有本质区别的计算方法，属于试验数学的一个分支，起源于早期的用几率近似概率的数学思想，它利用随机数进行统计试验，以求得的统计特征值（如均值、概率等）作为待解问题的数值解。

随着现代计算机技术的飞速发展，蒙特卡罗方法已经在原子弹工程的科学研究中发挥了极其重要的作用，并正在日益广泛地应用于物理工程的各个方面，如气体放电中的粒子输运过程等。

§1.1　蒙特卡罗方法的基本原理

就数学特性而言,蒙特卡罗方法的发展可以追溯到18世纪著名的蒲丰问题。

1777年，法国科学家蒲丰（Buffon）提出用投针试验计算圆周率π值的问题。

这里我们用蒲丰问题来初步说明蒙特卡罗方法的基本原理和解决问题的基本手续。

蒲丰问题是这样一个古典概率问题:

在平面上有彼此相距为2a的平行线，向此平面任意投一长度为2l的针，假定l

在这里，任投一针的概率含义有以下三点:

（1）针的中点Ml在平行线之间等概率落入，即Ml距平行线的距离x均匀分布在区间[0,a]之内；

（2）针与线的夹角θ均匀分布在区间之内；（3）x与θ互相独立。

如图1.1所示，建立与平行线垂直且原点在某一条平行线上的x轴，不失一般性，假定针的中心处于图示中的x轴上。

由于对称性，我们只需分析针中心处在范围的情况即可。

令探针中心的坐标值为x,显然，只有x≤l时才可能发生相交的事件。

我们来分析在条件x≤l满足时,针与线相交的概率：

只有当时才能相交，且相交的概率为

（1.1）

下面再来分析针中心位置在轴上的分布,显然,这是一个均匀分布,即针中心处于区间（x,x+dx）内的概率为

图1.1　蒲丰问题的概率分析

（1.2）

这样，一次投掷，针中心落入（x,x+dx）且与线相交的概率为

（1.3）

则一次投掷，针与线相交的总概率为

（1.4）

即:

（1.5）

从（5）式可见,可利用投针试验计算π值：

设投针N次，其中n次针与线相交，则可用频率值n/N作为概率P的估计值，从而求得π的估计值为

（1.6）

这就是早期的用频率值作为概率近似值的方法的应用实例，表1是在历史上一些有名的用投针试验计算π值的结果，其中针长以a为单位。

表1.1　投针试验计算π值的结果

实验者

时间（年份）

针长

投针次数

相交次数

π的估值

Wolf

1850

0.8

5,000

2,532

3.1596

Smith

1855

0.6

3,204

1,218.5

3.1554

De

1860

1.0

600

382.5

3.137

Fox

1884

0.75

1,030

489

3.1595

Lezzerini

1901

0.83

3,408

1,808

3.1415929

Reina

1925

0.5419

2,520

859

3.1795

需要指出的是，上述由投针试验求得π的近似值的方法，是进行真正的试验,并统计试验结果,要使获得的频率值与概率值偏差小，就要进行大量的试验，这在实际中，往往难以做到。

可以设想，对蒲丰问题这样一个简单的概率问题，若要进行10万次投针试验，以每次投针、作出是否相交判断并累加相交次数用时5秒钟计算，则需用时50万秒，即大约139个小时。

那么，可以设想，对于象上述确定条件下的核裂变、直流气体放电中粒子的输运过程及粒子输运的总效应，若要用多次掷骰子的方法近似求出就是不可能的了。

所以，在现代计算机技术出现之前，用频率近似概率的方法——抑或称为雏形时代的蒙特卡罗方法——并没有得到实质上的应用。

若用数值模拟方法代替上述的真正的投针试验，是利用均匀分布于（0,1）之间的随机数序列，并构造出随机投针的数学模型，然后进行大量的随机统计并求得π的近似值。

如图1.2建立坐标系，平面上一根针的位置可以用针中心Ml的坐标x和针与平行线的夹角θ来决定,在y方向上的位置不影响相交性质。

图1.2用数值模拟方法计算蒲丰问题

任意投针，意味着x与θ都是任意取的。

但θ的范围可限于[0,π],x的范围可限于[0,a]。

在这种情况下,针与平行线相交的数学条件是

（1.7）

其次，怎样模拟投针呢?

亦即如何产生任意的[x,θ]。

x在[0,a]任意取值，意味着x在[0,a]上取哪一点的概率都一样，即x的概率密度函数为

（1.8）

类似的，θ的概率密度函数为

（1.9）

由此，产生任意（x,θ）的过程就变为由f1（x）抽样x,由f2（θ）抽样θ的过程。

容易得到

（1.10）

式中，ξ1，ξ2均为（0,1）上均匀分布的随机数。

只要随机数的均匀性和独立性良好,如此构造的数值模型就很好地模拟了实际试验中的一次投针,并用下式判断是否相交且记录统计结果:

如果投针N次，那么

（1.11）

是相交几率P的估计值。

这样就实现了用数值方法模拟真正的投针试验。

用此方法计算的π的近似值的情况如表1.2所示。

表1.2　用蒙特卡罗方法计算的π的近似值

投针次数

10,000

20,000

100,000

200,000

π的近似值

3.162233

3.137993

3.141179

3.141354

表2中的计算结果表明，随着模拟投针次数的增大，所计算的π的近似值越来越接近于其真值，而要进行这样的数值模拟，就需要很大的计算量，只有利用计算机才能实现。

从蒲丰问题可以看出，用蒙特卡罗方法求解问题时，应建立一个概率模型，使待解问题与此概率模型相联系，然后通过随机试验求得某些统计特征值作为待解问题的近似解。

与此相似，在一些物理问题，如核裂变、直流气体放电等过程中，粒子的输运过程及粒子输运的总效应，也是可以与某些概率过程联系起来，例如，电子与原子、分子、离子的碰撞过程，实际上就是与碰撞截面有关的概率过程，这样，从数学物理特征来说，类似于用随机投针方法计算π的近似值,确定条件下的核裂变、直流气体放电中粒子的输运过程及粒子输运的总效应可以用多次掷骰子的方法近似求出。

随着现代计算机技术的出现和飞速发展，用计算机模拟概率过程，实现多次模拟试验并统计计算结果,进而可获得所求问题的近似结果.计算机的大存储量、高运算速度使得在短时间内,获得精度极高且内容丰富的模拟结果。

在历史上，也正是原子弹工程研究初期阶段的工作，为模拟裂变物质的中子随机扩散，提出了运用大存储量、高运算速度计算机的要求，这也成为当时推动计算机技术发展的重要动力，也就是在第二次世界大战期间，冯.诺依曼和乌拉姆两人把他们所从事的与研制原子弹有关的秘密工作——对裂变物质的中子随机扩散进行直接模拟——以摩纳哥国的世界闻名赌城蒙特卡罗（MonteCarlo）作为秘密代号来称呼。

用赌城名比喻随机模拟，风趣又贴切，很快得到广泛接受，此后，人们便把这种计算机随机模拟方法称为蒙特卡罗方法。

需要指出的是，正是由于广泛领域的物理问题中存在着大量的随机过程，如粒子间的碰撞等，使得蒙特卡罗方法在计算机物理和物理工程中得到日益广泛的应用，并成为沟通理论与实验研究的一个桥梁。

需要指出的是，蒙特卡罗方法不仅在处理具有概率性质的问题方面获得广泛的应用，对于具有确定性问题的计算也因其程序简单等优点获得了广泛的应用。

这里以定积分的计算简要说明其处理确定性问题的手续。

对于定积分

通过变量替换，可以转换为下面的形式

（1.12）

其中g（x）∈（0,1）,当x∈（0,1）时即转换为求积分，亦即求边长为1的正方形中一个曲边梯形的面积的问题，如图1.3所示。

图1.3　用蒙特卡罗方法求定积分

我们可以设想这样一种随机投点求定积分的方法：

在一个边长为1的正方形上并以其两边分别为坐标轴画出曲线g（x）,实际上就是图3，然后随机地正方形投掷小球，那么，小球击中g（x）曲线下部分的概率就等于所要求的积分,这样就将确定性的定积分问题转化为一个概率问题，同样可以通过数值模拟方法——蒙特卡罗方法求得其近似解。

用此方法，我们计算了积分,当投球数为1万次时,得到的积分近似值为0.332800,与其真值极为接近。

§1.2蒙特卡罗方法的解题手续和特点

在用蒙特卡罗方法解算问题时，一般需要这样几个过程：

构造或描述概率过程，对于本身就具有随机性质的问题，如粒子输运问题，主要是正确地描述和模拟这个概率过程。

对于本来不是随机性质的确定性问题，比如计算定积分、解线性方程组、偏微分方程边值问题等，要用求蒙特卡罗方法求解，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。

即要将不具有随机

性质的问题，转化为随机性质的问题。

这构成了蒙特卡罗方法研究与应用上的重要问题之一，然后建立各种估计量，使其期望值是所要求解问题的解，最后根据所构造的概率模型编制计算程序并进行计算，获得计算结果。

与其他的数值计算方法相比，蒙特卡罗方法有这样几个优点:

（1）收敛速度与问题维数无关，换句话说，要达到同一精度，用蒙特卡罗方法选取的点数与维数无关；计算时间仅与维数成比例。

但一般数值方法，比如在计算多重积分时，达到同样的误差，点数与维数的幂次成比例，即计算量要随维数的幂次方而增加。

这一特性，决定了对多维问题的适用性；

（2）受问题的条件限制的影响小；

（3）程序结构简单，在电子计算机上实现蒙特卡罗计算时，程序结构清晰简单，便于编制和调试；

（4）对于模拟象粒子输运等物理问题具有其他数值计算方法不能替代的作用。

蒙特卡罗方法的弱点是收敛速度慢，误差大的概率性质，这一情况在解粒子输运问题中仍然存在，除此之外，经验证明,只有当系统的大小与粒子的平均自由程可以相比较时，一般在10个平均自由程左右，这方法算出的结果较为满意。

而对于大系统深穿透问题，算出的结果往往偏低。

对于大系统，其他数值方法往往很适应，能算出较好的结果。

因此，已有人将数值方法与蒙特卡罗方法联合起来使用，克服这种局限性，取得了一定的效果。

自二战以来，随着计算机技术的飞速发展，蒙特卡罗方法以其独特的优点广泛应用于计算物理和物理工程领域，特别是在核物理、辐射物理、数学、电子学等方面得到了广泛的应用。

对蒙特卡罗方法的推广必将使其对物理学学科的发展发挥更大的作用。

§1.3用蒙特卡罗方法模拟粒子输运

中子和光子在物质中输运的宏观表现是大量粒子与原子核微观作用的平均结果，蒙特卡罗方法通过逐一模拟和记录单个粒子的历程来求解输运问题。

要得到比较合理的平均结果需要跟踪大量的粒子，至于单个粒子在其生命中的某一阶段如何度过，可以在已知统计分布规律的前提下通过抽取随机数来决定。

图1.4单个中子随机历程示意

图1.4显示了模拟中一个中子射入物质后的随机历程。

首先根据中子与物质作用的物理规律（分布函数），选取一个随机数决定中子在何处与原子核碰撞，本例中在事件1处碰撞；然后再用抽取随机数的方法决定中子与原子核发生了哪种反应，这里抽出的是非弹性散射反应；散射中子的能量和向哪个方向飞行也是用抽取随机数的方法从已知分布函数中决定的；碰撞过程中是否产生光子以及光子的能量、飞行方向等参数还是要通过抽取随机数从已知分布中决定，这里产生了一个光子，在MCNP程序中该光子暂时被存储起来,散射后的中子在事件2处与原子核发生（n,2n）裂变反应，同时产生两个新的中子和一个光子，裂变产生的光子和一个中子被存储起来，现在开始对另外一个中子进行跟踪。

这个中子在事件3处被俘获，这时取出那个刚才存储起来的中子，通过随机抽样，该中子在事件4处逃逸出平板，进入探测器，MCNP结束对该中子的模拟，重新取出裂变产生那个光子进行跟踪。

这个光子在事件5处与原子核发生了一次散射，在事件6处逃逸出平板。

现在取出最后一个粒子——在事件1处产生的光子进行跟踪，该光子在事件7处与原子核碰撞并被吸收。

至此，入射中子的整个历史也就完成了,其有一个中子到达了探测器，感兴趣的结果被记录下来。

跟踪越来越多的入射粒子历程后，平均结果就能反映出宏观效果。

通过以上描述，读者不难领略蒙特卡罗方法如何通过跟踪粒子历程的方法计算问题，也了解了随机数在蒙特卡罗计算中的独特作用。

需要注意的是，MCNP程序在重新取回以前存储的粒子时，是按照这样的规则进行的：

最后存储的粒子最先被取出，也就是“先进后出”或“后进先出”的堆栈操作。

第2章MCNP程序入门

由于蒙特卡罗（MonteCarlo）方法在解决粒子输运问题方面的优越性，人们开发了很多种MonteCarlo模拟程序，如EGS、MCNP、MORSE和FLUKA等，其中，EGS4（ElectronGammaShower，Version4）和MCNP4C（MonteCarloN—ParticleTransportCode，Version4C）是应用最多并得到广泛验证的两种经典MonteCarlo模拟程序。

EGS4是山美国斯坦福直线加速器中心、日本高能物理国家实验室（KEK）和加拿大国家研究所（NRC）联合推出的——套模拟电子和光子在物质中输运过程的通用MonteCarlo程序系统，是揭示电子和光子在物质中输运规律的有力和较为方便的模拟研究工具，目前KEK发布了EGS5版本。

§2.1MCNP简介

MCNP全名为MonteCarloNeutronandPhotoTransportCode（蒙特卡罗中子-光子输运程序），它是由美国LosAlamos国家实验室应用理论物理部（X部）的MonteCarlo小组（X-6X小组）经过数十年的研究开发的一个基于蒙特卡罗方法的大型的多功能MonteCarlo粒子输运程序。

从1977年开始产生到现在历经十几个版本，解决了核能领域很多关键性问题，功能也越来越强大。

现在MCNP可在微机的的UNIX、LINUX、DOS、WINDOWS98、WindowsXP等操作系统下工作。

现其最新版本的MCNP-51.30具有如下功能：

（1）非带电粒子成相技术。

在用户指定的栅格中，MCNP-5使用多个点探测器来确定某个象素区域的粒子流量。

用户可以根据需求设置尽可能多的探测点以便生成尽可能平滑的图象。

（2）随机几何能力。

该能力可用来分析颗粒燃料，还可用来研究燃料核在石墨矩阵中的随机位置。

（3）可处理复杂三维几何系统的输运问题，几何界面除任意平面和二阶曲面外，也可包括四阶椭环面。

（4）粒子输运方式可以是中子输运、光子输运、电子输运、中子-光子联合输运、光子-电子联合输运、电子-光子联合输运、中子-光子-电子联合输运。

既可用于求解通常的输运方程，也可解多群共轭输运方程。

MCNP-5已经能够处理低能光子相互作用的不连续散射问题。

（5）既可计算穿透问题，也可计算临界特征值问题。

对临界特征值的计算，给出了KEFF、预期寿命和生存时间的计算方法，还可计算各种记数关于介质成分、密度或截面数据的一阶、二阶微扰量。

（6）配备的截面数据覆盖了所有常用的核素和同位素，并可选用点截面方式或多群截面方式。

可处理的中子能量范围为10E-11至20MeV，光子和电子为0.001至1000MeV。

（7）有多种物理量的计算选择，包括点通量、界面通量、任意独立栅格的粒子流及通量、几何体上的通量及能量沉积，可给出按空间、时间、能量的谱（分布）和联合分布，粒子流还可增加角度分布。

另外还有脉冲高度谱的计算。

所有量的计算结果都同时给出了统计误差。

（8）有十余种降低方差技巧可以选用。

（9）具有较好的图形输出功能。

可作几何结构任意方位的二维剖面图，计算结果的二维X-Y曲线图、等值线图、三维曲面图、截面数据图及粒子径迹图。

既可交互方式作图，也可用批命令方式作图，也可在中途暂挂计算进行作图。

MCNP-5可绘制最多有64种颜色的彩图。

（10）既可在串行环境下运行，也可在并行环境下运行。

另外，根据具体应用的需要，MCNP还有一些专门用途的版本。

如MCNP-X和MCNP-PoliMi。

MCNP-X（MonteCarloN-particlecode）是由美国LosAlamos实验室开发的用于多粒子传输模拟计算的蒙特卡罗程序。

它是由MCNP4B版本发展而来,可以跟踪的粒子数达到40余种,可跟踪的中子最高能量达到GeV以上。

在辐射防护研究中也得到了广泛的应用。

MCNP-PoliMi（PolitecnicodiMilano）是由意大利米兰工艺学院核工程学部根据MCNP4C开发用于（时间）相关性测量的蒙特卡罗软件。

它可以很详细的记录粒子-核子的相互作用信息，比如每个粒子关于每一次碰撞的类型、目标核、能量沉积、位置等。

MCNP-PoliMi在处理中子输运问题时，首先抽取样品中子碰装类型,然后产生光子，这与标准MCNP采用的技术相反。

一个与MCNP-PoliMi配套的后加工处理程序也已经开发，用它能够定制模拟特定特征的探测器。

该软件可用于模拟塑料闪烁体对快中子或伽吗射线的探测问题。

MCNP程序的应用范围十分广阔，主要包括：

反应堆设计、核临界安全、辐射屏蔽和核防护、探测器的设计与分析、核测井、个人剂量与物理保健、加速器靶的设计、医学物理与放射性治疗、国家防御、废物处理、射线探伤等。

§2.2MCNP程序的组成及特点

A、组成

MCNP是主程序模块，它根据正在运行的问题之需要分别调用IMCN、MCRUN、XACT和PLOT等主要模块。

各模块之间的关系，如图1所示。

图2.1MCNP程序模块构成简图

IMCN模块总是要被调用的，其任务是读入输入文件INP，并对输入数据进行分析处理。

它根据输入信息定义动态分配存储空间的大小，对相应的动态分配数组建立相对地址索引表，随后将调用所需要的模块如RDPROB、TRFMAT、ITALLY和VOLUM等。

RDPROB模块主要是依序读人输入文件INP的卡片映像，通过分析处理送入相应的存储位置。

TRFMAT模块对非标准形式输入的曲面描述卡，完成坐标变换，给出曲面标准形式的方程系数；对输入的源分布处理成便于抽样的形式。

ITALLY模块对输入的记数及材料说明信息进行加工处理。

VOLUME模块则主要是计算栅元的体积及曲面的面积，但限于具有对称轴的情况。

MCRUN模块是输运计算的实体，是MCNP的核心部分。

它将根据需要调用两个模块TRNSPT和OUTPUT。

TRNSPT模块按用户指定的样本数NPS或时间限CTME完成对粒子历史的模拟。

OUTPUT模块则实现对计算结果的编辑输出。

XACT模块的任务是从截面数据文件中读出问题中所用核素的截面数据，并根据用户给出的信息删去所关心能量范围以外的中子截面数据。

PLOT模块实现在各种图像设备上绘制或显示问题的几何图形。

它分别调用SETUP模块和WRAPUP模块。

SETUP模块根据用户的终端键盘信息绘制几何图形，产生图形文件PIX。

WRAPUP模块在不同的图像设备上实现PIX文件的图形输出。

显然，上述各模块并不是运行任何问题都必须全部调用的，MCNP从分析用户在输入文件INP上指定的信息，决定所需要调用的模块。

B、特点

MCNP程序中物理量的单位规定为：

表2.2MCNP程序中物理量的单位规定

物理量

单位

长度

cm

能量

MeV

时间

shake,=10-8sec

温度

MeV（kT）

原子密度

atoms/barn-cm

质量密度

g/cm3

截面

barns（10-24cm2）

加热量

MeV/collision

此外，原子质量按中子质量为1.0计算，这种单位下阿佛伽德罗常数是0.59703109；程序运行时间以分钟为单位。

MCNP程序具有很强的通用性，主要体现在：

1.可以处理任