核工程与核技术毕业设计（论文）不同中子源下补偿中子测井仪关键参数的MCNP模拟.doc

《核工程与核技术毕业设计（论文）不同中子源下补偿中子测井仪关键参数的MCNP模拟.doc》由会员分享，可在线阅读，更多相关《核工程与核技术毕业设计（论文）不同中子源下补偿中子测井仪关键参数的MCNP模拟.doc（36页珍藏版）》请在三一办公上搜索。

1、 毕业设计（论文）题 目 不同中子源下补偿中子测井仪关键参数的MCNP模拟 英文题 Monte-carlo comparison of Compensated Neutron Logger pivotal parameters Under Various neutron source学生姓名 专 业 核工程与核技术班 级 指导教师 二零一一年六月摘 要补偿中子测井是利用中子与钻井周围岩石和井内物质相互作用，探测地层岩石及其孔隙流体的化学元素和含量的测井技术。当地层空隙度小，其单位体积含氢量就小，利用补偿中子测井时含氢指数识别油气层则不敏感，可以通过调整合适的屏蔽体厚度、近、远探测器间距r和中子

2、源距r等参数，以提高补偿中子测井对含氢指数识别油气层的灵敏度。本论文便是利用蒙特卡罗方法对补偿中子测井仪模型关键参数进行数值模拟计算，得出如下主要结论：本次模拟得到252Cf和Am-Be中子源下各参数的变化规律为：铅屏蔽体厚度取6cm时远探测器碳氧比值取得最大值。中子源为252Cf、探测器间距大于25cm时，远探测器碳氧比值随探测器间距的增加而不变；中子源为Am-Be、探测器间距大于20cm时，远探测器碳氧比值随探测器间距的增加而不变。252Cf和Am-Be中子源下，近探测器热中子相对注量率和远探测器相对注量率都随源距的增加而减小。Am-Be中子源下远探测器碳氧比值在空隙度为25%时取得最大值

3、，252Cf中子源下远探测器碳氧比值随空隙度的增大而减小。这一结果，在补偿中子测井仪研制过程中，对仪器选取合适的屏蔽体厚度、近、远探测器间距r和中子源距r等参数，减少测试时间，降低成本，人员安全等方面，均有重要的意义。关键词： 蒙特卡罗； MCNP； Logging； 中子注量率 ABSTRACT Compensated Neutron Logger(CNL) is one of the logging technology, Which use the interaction between neutron and rock around drilling to detect the con

4、cent of chemcial elements of pore fluid and rock .Compensated Neutron Logger(CNL) becomes less sensitive to Hydrogen Index(HI) in identification of oil/gas bearing layers when formation porosity is small due to the low hydrogen quantity per unit volume. This paper aims to investigate the relationshi

5、ps between the CNL response and its design parameters, such as lead shield thickness, near-to-far detector spacing, and source-to-detector spacing, using the Monte Carlo method. The main study results as follow: The Simulation shows the change conclusion of each parameter under 252Cf and Am-Be neutr

6、on source: ,Far-detector Carbon-oxygen ratio Obtain maximum while the thickness of lead shielding is 6cm. For 252Cf neutron source and Detector spacing is More than 25cm,far-detector carbon-oxygen ratio is unchanged With the increase of detector spacing; For Am-Be neutron source and detector spacing

7、 is more than 20cm,far-detector carbon-oxygen ratio is also unchanged With the increase of detector spacing.Under Am-Be neutron source, near-detectors thermal neutron relative flux rates and far-detector relative flux rates are reduced as an increase in source distance. Under Am-Be neutron, far-dete

8、ctor carbon-oxygen ratio Obtain maximum when interstice degree up to 25%.Under 252Cf neutron, far-detector carbon-oxygen ratio reduce with the interstice degree accretion. This result has important reference value for the selection of lead shield thickness, near-to-far detector spacing, and source-t

9、o-detector spacing，reduce test time, lower cost worker safety etc. on the development process in the CNL.Key words: Monte Carlo; MCNP; Logging; neutron fluence rate目 录绪 论21. 补偿中子测井仪基本原理41.1 伽马射线与物质的相互作用41.1.1 光电效应51.1.2 康普顿效应51.1.3 电子对效应61.1.4 伽马射线的吸收61.2 中子与物质的相互作用71.2.1 弹性散射71.2.2 非弹性散射71.2.3 辐射俘获

10、71.3 补偿中子测井仪的物理基础71.3.1 碰撞机率71.3.2 减速指数81.3.3 热扩散91.4 探测中子的基本方法111.4.1 俘获伽马射线探测111.4.2 热中子探测111.4.3 超热中子探测112. 蒙特卡罗方法与MCNP程序122.1 蒙特卡罗方法122.2 蒙特卡罗方法的解题步骤132.3 基于蒙特卡罗方法的MCNP程序142.3.1 MCNP程序的发展142.3.2 MCNP程序的特点152.3.3 MCNP的应用状况163. 不同中子源下补偿中子测井仪模型的建立与程序模拟183.1 模拟模型183.2 模拟计算及分析193.2.1 屏蔽体厚度模拟193.2.2 探

11、测器间距的模拟203.2.3 源距的模拟203.2.4 孔隙度可变的水砂岩和油砂岩的模拟22结 论24致 谢25参考文献26附录127附录230绪 论当前，核测井已成为石油天然气和煤炭探测与开发的一种重要手段。长期以来，中子测井设计主要是依赖建造大量不同岩性和不同孔隙度的模型井，通过实验测得其响应值，以供制作图版曲线应用,从而用来设计补偿中子测井仪。近年来，随着核能技术和计算机技术的发展，国外开始把核反应堆物理和屏蔽计算中的核数据库和计算软件移植应用到中子测井的数值模拟中来，取得了良好的效果1。目前采用的数值模拟方法基本上可以分为两类：一类是确定性方法，它是从数学上直接对中子输运方程数值求解；

12、另一类为蒙特卡罗(MC)方法2，它是利用统计试验方法，以随机抽样为基础，在计算机中模拟中子在介质中的游动过程，最后获得统计结果。MC方法以其本身固有的一些特点，例如不受复杂几何条件的限制等，使其在核测井的数值模拟中占有重要的地位。近年来，在核能领域中已研制了一些用于核辐射问题的通用MC程序，并建立了一些完整的核数据库，例如美国Los Alams国家实验室的MCNP3 、英国原子能委员会的MCBEND和美国橡树岭国家实验室的SAMCE等程序。这些程序目前均已开始用于测井数值模拟并取得了良好的结果1。补偿中子测井是利用中子与钻井周围岩石和井内物质相互作用，探测地层岩石及其孔隙流体的化学元素和含量，

13、在孔隙地层中，骨架物质的含氢量极少，充满孔隙的水和油却含氢量很高，在物理特性上，氢是最有效的中子慢化剂之一。补偿中子测井仪是根据中子与氢的相互作用来探测地层中含氢量，进而推求地层孔隙度的一种仪器。在实际应用中，为对补偿中子测井仪进行优化设计，需要选择合适的中子源距、近、远探测器间距以及探测器灵敏区长度等参数来提高补偿中子测井对含氢指数识别油气层的灵敏度4,5。但由于现场利用中子源对各种补偿中子测井仪进行优化设计，不仅费时费力，而且也不安全。而MCNP程序是由美国Los Alamos实验室在20世纪70年代开发的基于蒙特卡罗方法的粒子输运通用程序。中子和光子在物质中的输运的宏观表现是大量粒子与原

14、子核微观作用的平均结果，蒙特卡罗方法通过逐一模拟和记录单个粒子的历程来解决输运问题。通过跟踪每一个光子（中子或电子）与物质的反应过程，记录每一个感兴趣的结果。在跟踪大量的粒子后，平均结果就可以反映出宏观结果。因此，本文通过MCNP程序来对不同中子源下补偿中子测井仪测井过程进行模拟。通过利用MCNP程序对补偿中子测井仪模型分别在252Cf和Am-Be中子源下进行数值模拟计算，得到仪器灵敏度随屏蔽体厚度、近、远探测器间距r、中子源距r等参数变化而变化的规律。本论文中第一章对补偿中子测井仪的基本结构与原理做了详细的阐述，使本论文的研究工作有了充分的理论支持；第二章对蒙特卡罗方法和MCNP程序进行了详

15、细的介绍，使本论文的研究工作有了基本的模拟方法和程序技术支持；第三章对模型的建立和模拟过程进行了详细的介绍，同时对模拟结果进行了分析与研究，使本论文的研究具有一定的实际意义；第四章是本论文最后的结论部分，是对本论文模拟结果的总结，同时对在本课题模拟过程中出现的问题不足以及未来的改进方向进行了介绍。1. 补偿中子测井仪基本原理补偿中子测井6,7(Compensated Neutron Logger如图1所示)是中子测井方法之一。补偿中子测井仪是一种孔隙度测井仪，该仪器响应于周围存在的氢的数量。在孔隙地层中，骨架物质的含氢量极少，可以忽略不计，而充满孔隙的水和油却含有氢，并且水和油的含氢量基本相等

16、。在物理特性上，氢元素的原子核质量与中子的质量最相近，因此氢是最有效的中子减速剂。补偿中子测井仪器就是根据中子与氢的相互作用来探测地层中的含氢量，从而推求地层孔隙度的。图1 补偿中子测井仪示意图CNL还可用作天然气指示，假如孔隙空间含天然气，天然气的含氢量低，所测的视孔隙度将过低。因而如果CNL孔隙度比由其它方法确定的孔隙贬值低得多(例如补偿地层密度测井(PGT)孔隙度)，则该地层可能含天然气。与其它测井方法组合，CNL还可用于帮助确定许多其它的地层参数，例如岩性、粘土含量、含水饱和度等等。补偿中子测井仪应用一个源和两个不同源距的探测器(图1)。 这里远探测器是起响应地层物质放出的射线，从而得

17、出C/O比值，近探测器是探测热中子的，通过探测热中子得出地层中的含氢量，结合远近探测器的结果最后得出地层中的油和水的含量。1.1 伽马射线与物质的相互作用射线、韧致辐射、湮没辐射和特征射线等，虽然它们的起源不一、能量大小不等，但都属于电磁辐射。电磁辐射与物质相互作用的机制、与这些电磁辐射的起源是无关的，只与它们的能量有关8。射线与物质的相互作用和带电粒子与物质的相互作用有着显著的不同。光子不带电，它不像带电粒子那样直接与靶物质原子、电子发生库伦碰撞而使之电离或激发，或者与靶原子核发生碰撞导致弹性碰撞能量损失或辐射损失，因而不能像带电粒子那样用阻止本领dE/dx和射程来描述光子在物质中的行为。带

18、电粒子主要通过连续的与物质原子的核外电子的许多次非弹性碰撞逐渐损失能量的，每一次碰撞中所转移的能量很小。而光子与物质相互作用时，发生一次相互作用就导致损失其大部分或全部能量。光子不是完全消失就是大角度散射掉。光子可以通过与物质的相互作用被间接探测到。这些作用过程产生带电的次级粒子，随后在探测器的灵敏体积内通过电离过程被记录下来。射线与物质相互作用，可以有多种方式。当射线能量在30MeV以下时，在所有的相互作用中最主要的有三种方式：光电效应、康普顿效应和电子对效应。还有一些其它的相互作用方式，如：相干散射、光致核反应和核共振反应等。1.1.1 光电效应射线与靶物质原子的束缚电子作用时，光子把全部

19、能量转移给某个束缚电子，使之发射出去，而光子本身消失掉，这种过程称为光电效应8。光电效应中发射出去的电子叫做光电子。原子吸收了光子的全部能量，其中一部分消耗于光电子脱离原子核束缚所需的能量即电离能，另一部分就作为光电子的动能。所以，释放出来的光电子的能量就是入射光子的能量和该束缚电子所处的电子壳层的结合能之差。因此，要发生光电效应，光子的能量必须大于电子的结合能。光电子可以从原子的各个电子壳层中发射出来，但是自由电子（非束缚电子）却不能吸收入射光子能量而成为光电子。这是因为动量守恒要求，在光电效应过程中，除入射光子和光电子外，还需要有一个第三者参加，这第三者就是原子核，严格地讲师发射光电子之后

20、剩余下来的整个原子。它带走一些反冲能量，但这能量十分小。由于它的参加，动量和能量守恒才能满足。而且，电子在原子中束缚得越紧，就容易使原子核参加上述过程，产生光电效应的概率也就越大。所以在K壳层上打出光电子的概率最大，L层次之，M，N层更次之。如果入射光子的能量超过K层电子结合能，那么，大约80%的光电吸收发生在这K层电子上。发生光电效应时，从原子内壳层上打出电子，在此壳层上就留下空位，并使原子处于激发状态。这种激发状态是不稳定的，退激过程有两种。一种过程是外层电子向内跃迁，来填补这个空位，使原子恢复到比较低的能量状态。两个壳层的结合能之差就是跃迁时释放的能量此能量以射线形式释放出来。另一种是原

21、子的激发能交给外壳层的其它电子，使它从原子中发射出来，这种电子称俄歇电子。因此在发射光电子的同时，还伴随着原子发射的特征X射线或俄歇电子。这些特征X射线和俄歇电子再与靶物质原子发生作用。1.1.2 康普顿效应康普顿效应8是入射光子与原子的核外电子之间发生的非弹性碰撞过程。这一作用过程中，入射光子的一部分能量转移给电子，使它脱离原子成为反冲电子，而光子的运动量发生变化。康普顿效应与光电效应不同。光电效应中光子本身消失，能量完全转移给电子；康普顿效应中光子只是损失掉一部分能量。光电效应发生在束缚得最紧的内层电子上；康普顿效应总是发生在束缚得最松的外层电子上。虽然光子与束缚电子之间的康普顿散射，严格

22、地讲是一种非弹性碰撞过程。但外层电子的结合能是较小的，一般是电子伏数量级，与入射光子的能力相比较，完全可以忽略，所以可以把外层电子看作是“自由电子”。这样康普顿效应就可以认为是光子与处于静止状态的自由电子之间的弹性碰撞。入射光子的能量和动量就由反冲电子和散射光子两者之间进行分配。1.1.3 电子对效应当光子从原子核旁经过时，在原子核的库伦场的作用下，光子转化为一个正电子和一个负电子，这种过程称为电子对效应8。根据能量守恒定律，只有当入射光子能量hv1.02MeV时，才能发生电子对效应。光子的能量除一部分转变为正负电子对的静止能量（1.02MeV）外，其余就作为他们的动能。除了在原子核库仑场中发

23、生电子对效应外，在电子的库仑场中也会产生正-负电子对。不过电子质量小，反冲能量较大，所以产生电子对的最低入射能量至少是4M0c2而且产生电子对的概率要小得多。对于一定能量的入射光子，电子对效应产生的正电子和负电子的动能从0到hv2M0c2都是可能的，电子和正电子之间的能量分配是任意的。电子对过程中产生的快速正电子和电子，在吸收物质中通过电离损失和辐射损失消耗能量。正电子在吸收体中被很快慢化后，将发生湮没，湮没光子在物质中再发生相互作用。1.1.4 伽马射线的吸收当光子穿过物质时，与吸收物质的原子一旦发生光电效应、康普顿效应或电子对效应，原来能量为hv的光子就消失，或散射后能量改变掉，并偏离原来

24、入射方向，即从原来入射的束中移去。没有与物质发生相互作用的光子穿过吸收层，其能量保持不变。射线穿过物质时，强度逐渐减弱，因而没有射程的概念可言。只可用半吸收厚度（强度减至一半时所需厚度）来表示射线对物质的穿透情况8。射线通过吸收物质，其强度的衰减规律遵循指数衰减规律，可以用公式(1-1)来表示： I=I0e-t =ph+c+p (1-1)ph：光电吸收系数c：康普顿系数系数p：电子对吸收系数射线与物质相互作用产生的次级粒子，如光电子、康普顿散射电子、正负电子对、俄歇电子以及康普顿散射光子、湮灭光子和特征射线等，可以继续在物质中发生相互作用，直到全部能量完全耗尽为止。这些级联过程的发生与射线的能

25、量、靶物质的性质以及探测器的几何尺寸有关。这些也是关系到探测器的能量相应的方面。1.2 中子与物质的相互作用中子与原子核发生的相互作用有弹性散射、非弹性散射、辐射俘获。1.2.1 弹性散射弹性散射9在中子的所有能量范围内都可能发生。它可分为共振弹性散射和势散射两种。前者经过复合核的形成过程，后者不经过复合核的形成过程。由于共振现象只对具有特定能量的入射中子才会产生，因此共振弹性散射也只对特定能量的中子才能发生。在弹性散射过程中，由于散射后靶核的内能没有变化，它仍保持在基态，散射前后中子-靶系统的动能和动量是守恒的，所以可以把这一过程看作“弹性球”式的碰撞，根据动能和动量守恒，用经典力学的方法来

26、处理。中子从高能慢化到低能的过程起主要作用的是弹性散射。1.2.2 非弹性散射当发生非弹性散射9时，中子首先被靶核吸收而形成处于激发态的复合核，在这个过程中，入射中子把它的一部分动能（通常为绝大部分）转变成了靶核的内能，使靶核处于激发态，然后靶核通过放出中子并发射射线而返回基态。因此，散射前后中子与靶核系统的动量守恒，但动能不守恒。在发生非弹性散射时，中子能量的损失是可观的，但并不是所有能量的中子都能发生弹性散射，只有当入射中子动能高于靶核的第一激发态的能量时才能使靶核激发，也就是说，只有入射中子的能量高于某一数值时才能发生非弹性散射，由此可知非弹性散射具有阈能的特点。1.2.3 辐射俘获中子

27、和原子核相互作用时，辐射俘获是很主要的作用过程。中子射入靶核后与靶核形成一个复合核，而后复合核通过发射一个或几个特征g光子跃迁到基态。这些特征 g 光子不同于 (n,ng) 的特征 g 光子。由于这些 g 光子的发射与复合核的寿命相关，一般很快，故称为“中子感生瞬发g射线”，同样在核分析技术中有重要的应用。当发生(n,g)反应后，新形成的核素是放射性的，就是常说的“活化”，测量活化核素的放射性可以用来测量中子流的注量率，区分中子的能量范围。1.3 补偿中子测井仪的物理基础1.3.1 碰撞机率中子与给定类型原子核发生碰撞的相对机率7，正比于单位体积内这些原子核的数目。例如，若一中子垂直入射到一个

28、非常薄的地层表面，发生碰撞的机率正比于它所“见”到的原子核的总数目。这一数目又正比于单位体积内原子核的数目（原子核密度）。除原子核密度以外，尤其要考虑原子核的其他性质。不同的原子核对中子的“注意”程度是不同的，一些原子核比另一些要大。用所谓的原子核“截面”来考虑这一事实，截面大的原子核比截面小的原子核受中子碰击的可能性大。这样，中子与给定类型的原子核发生“下一次”碰撞的相对机率不仅与单位体积内原子核的数目有关，而且与它们在中子面前显示的大小有关，即与它们的截面是多大有关。实际上这一机率是由n和这两个量的乘积来决定，此处n是目的层的原子核密度，是原子核的截面。随中子能量而变化。这样中子与给定类型

29、原子核发生“下一次”碰撞的相对机率同样与中子能量有关。图2是孔隙度为15%纯砂岩中原子核和中子发生碰撞的机率与中子能量的关系图。图2 中子能量与碰撞机率图最初中子很可能与硅或氧进行碰撞，主要由于氢原子核对快中子显示出相当小的截面（小）。随着中子在其相续碰撞过程中能量的损失，在横坐标数值向左移动，当中子的能量降低到2104eV以下时，中子与氧或氢进行“下一次”碰撞的机率大致相等。而与硅发生碰撞的机率要小的多。1.3.2 减速指数快中子减速过程主要受两种因素影响，一是中子与原子核发生碰撞期间遭受能量损失，另一个是中子与原子核发生碰撞的实际几率。为了衡量给定类型原子核对中子减速的相对有效性，真正有意

30、义的量是两者的结合。将其称为减速指数7。对于给定元素，减速指数定义为n与的乘积。n是中子扩散1厘米时与原子核发生碰撞的几率，是碰撞过程的平均能力损失，因此这两个因子的乘积是1厘米长的中子路径上能量损耗的对数的度量。图3显示了前面讨论的孔隙度为15%纯砂岩的减速指数。图3 孔隙度为15%的砂岩的减速指数 能量为5MeV的中子与下列原子核碰撞后的平均中子：氢2.5MeV、碳4.29MeV、氧4.45MeV、硅4.67MeV、钙4.76MeV。这里清楚地见到了在中子减速过程氢的主要影响。尽管氢与中子发生碰撞（在中子能量小于2104eV时）的相对几率只比氧大一点，但氢原子核显示的减速指数比氧大一个数量

31、级。即使中子能量达到10MeV，但大部分能量都损失在与氢原子核碰撞期间，这就可以看出氢的影响。硅的影响，正如其减速指数表示的那样，比氧的影响小很多。随着含氢密度（通常为孔隙度）的增加，存在的其它元素的重要性变小，像石灰岩（碳酸钙）和砂岩（二氧化硅）之间的差异对中子减速过程的影响几乎微不足道。这样，尽管氢是决定中子减速的主要元素，但原则上地层中存在的其它元素也起作用，并且在某些情况下不可忽略。初始能量高时，碰撞几率小，中子传播相当远的距离。随其能量的减少，两次碰撞之间的平均距离变得更小。下式近似地表达了这些中子的空间分布:中子相对密度 (1-2)式中，r为中子离源的距离；L1为特定地层的减速长度

32、。1.3.3 热扩散当中子减速与地层热平衡后，中子继续进行连续碰撞，但平均能量无任何变化。这就是热扩散7阶段。尽管每一中子始终运动着，且进行随机碰撞，但从总体上看，它们构成了一空间分布，我们可以讨论空间各点的热中子密度，即单位体积内热中子数目。平均来说，中子离开中子源以后，由于重复的碰撞过程连续不断地进行，以及中子不停地运动，因此热中子在空间的分布比能量在热能以上的中子分布更广。两种过程唯一的差别就是热扩散期间平均能量保持恒定。在热扩散阶段，热中子两次碰撞之间迁移的距离不如高能中子远，由于能量恒定，所以连续两次碰撞之间的平均距离恒定。将L2定义为有效的热扩散长度。特有效热扩散长度与有效长度结合

33、起来，得到热中子空间分布参数即迁移长度M，当热中子密度时，其值见（1-3） 式中。 (1-3)这些热中子最后被俘获。n微秒(或更短)后，达到平衡状态，在平衡区域内从源出发的中子数目等于被俘获的中子数目。可以认为在源的周围存在一片“中子云”，中子云的大小取决于孔隙度。孔隙度越低，中子云范围越大。图4是一些理论曲线，表示了热中子密度与离源距离的函数关系，这是对孔隙中充有水的纯硅质砂岩骨架进行的计算。无论孔隙度多大，随着离源距离的增加，中子密度减小。然而每种孔隙度具有其自己的曲线。靠近个子源，热中子密度随孔隙度的增加而增大，这是由于氢原子使中子能量迅速减到热中子能级。图4 热中子密度与离源距离关系在

34、交叉带热中于密度实际上与孔隙度无关。进步离开源，热中子密度随孔隙度的增大而减小，因为平均来说， 中子在迁移这样的距离之前，巳被减速和俘获。为测量孔隙度，显然探测器与源的距离必须小于10厘米或大于30厘米。1.4 探测中子的基本方法探测中子的基本方法有三种，分别是俘获伽马射线探测，热中子探测和超热中子探测79。1.4.1 俘获伽马射线探测俘获伽马射线密度与中子密度有关，后者又和所存在的氢的数量有关。因为伽马射线探测器制造较简单，其原理应用于最初的中子测井仪。上述方法由于需要对俘获元素作校正，而这种元素的实际数量又可能是未知的，因而相对来说不太精确。井眼影响也引入大的误差。1.4.2 热中子探测该

35、方法用于CNL-A型中子测井仪，应用一个He-3正比计数器直接探测中子。虽然上述方法优于俘获伽马射线法，但由于俘获元素影响中子密度，测量精度任然多少受些影响。1.4.3 超热中子探测这方法只探测超热中子。高于热能级时俘获率很低，可以忽略不计，所以由于俘获元素引起的误差可以忽略。探测器如同用于热中子探测的一样是一个He-3计数器，但有一个镉过滤器在热中子到达探测器的计数管之前就俘获它们。此仪器的缺点是计数率过低，以致统计误差成为一个问题。一般，超热中子密度大约为热中子密度的十分之一。这种探测中子的方法将被引入于新一代的补偿中子测井仪（CNL-G）中。2. 蒙特卡罗方法与MCNP程序MCNP是用来

36、计算中子、光子、电子或者中子光子电子耦合问题的通用蒙特卡罗粒子输运计算程序，它以其灵活、通用的特点以及强大的功能，在诸多领域得到广泛认可和应用11。2.1 蒙特卡罗方法蒙特卡罗方法（Monte Carlo Method）12，也称统计模拟方法，由S.M.乌拉姆和J.冯诺伊曼在20世纪40年代为研制核武器而首先提出，又称随机抽样或统计试验方法，是以概率和统计理论方法为基础的一种计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。蒙特卡罗方法的名字来源于摩纳哥的一个城市蒙特卡罗，该城市以赌博业闻名，而蒙特卡罗方法正是以概率为基础的方法，与它对应的是确定性算法。蒙特卡罗方法在金融

37、工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果。而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。中子和光子在物质中输运的宏观表现是大量粒子与原子核微观作用的平均结果，蒙特卡罗方法通过逐一模拟和记录单个粒子的历程来求解输运问题。要得到比较合理的平均结果需要跟踪大量的粒子，至于单个粒子在其生命中的某一阶段如何度过，可以在已知统计分布规律的前提下通过抽取随机数来决定。图5 中子与物质相互作用过程模拟图图5显示了模拟中一个中子射入物质后的随机历程。

38、首先根据中子与物质作用的物理规律，选取一个随机数决定中子在何处与原子核碰撞，本例中在1点碰撞，然后再用抽取随机数的方法决定中子与原子核发生了哪种反应，这里抽出的是非弹性散射反应。散射中子的能量和向哪个方向飞行也是用抽取随机数的方法从已知分布函数中决定的，碰撞过程中是否产生光子以及光子的能量、飞行方向等参数还是要通过抽取随机数从已知分布中决定，这里产生了一个光子。跟踪光子，确定它在7点与原子核碰撞并被吸收。散射后的中子在2点与原子核发生（n，2n）反应，其中一个出射中子射向探测器，另一个中子在3点被吸收。在2点的碰撞还产生了一个光子，它在5点又与原子核发生了一次散射反应，并离开物质。这一入射中子

39、的历史过程结束了，有一个中子到达了探测器，感兴趣的结果被记录下来。跟踪越来越多的入射粒子历程后，平均结果就能反映出宏观效果。2.2 蒙特卡罗方法的解题步骤当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，可以通过某种试验的方法，得到这种事件出现的频率，或者这个随机变量的平均值，并用它们作为问题的解，这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是一个以概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，得到问题的近似解。蒙特卡罗解题过程可以归结为三个主要步骤11：（1）构造或描述概率

40、过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。（2）实现从已知概率分布抽样构造概率模型以后由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因，最为简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机

41、数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从（0，1）上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本

42、工具。（3）建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计，当然还可以引入其它类型的估计，如最大似然估计，渐进有偏估计等。但是在蒙特卡罗计算中，使用最多的是无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很节省，这些是使用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂，蒙特卡罗方法的应

43、用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学、信息科学、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。2.3 基于蒙特卡罗方法的MCNP程序MCNP程序是基于蒙特卡罗方法的一种粒子输运模拟程序。2.3.1 MCNP程序的发展MCNP程序是由美国Los Alamos国家实验室研制开发的，为其所投入的研究、发展、程序编写及参数制作超过了500人年。20世纪40年代美国Los Alamos实验室的Fermi、Von Neumann和Ulam等人提

44、出用蒙特卡罗方法模拟辐射输运的思想。1947年Fermi发明了第一台用蒙特卡罗方法计算中子链式反应的机器。从50年代开始，Von Neumann领导一个小组研究输运问题的蒙特卡罗处理方法，编写出模拟中子输运的程序MCS。1963年蒙特卡罗方法描述语言标准化。1965年完成的中子输运程序MCN有了很大改进，使用了标准的截面库，并且具有复杂几何描述功能。后来，Los Alamos实验室又开发了模拟光子输运的程序MCG（高能）和MCP（能量低至1keV）。1973年MCN和MCG合并成MCNG，为MCNP的雏形。自那时起，每隔2-3年更新一次，版本不断发展，功能不断增加，适应面也越来越广。已知的MC

45、NP程序研制版本的更新时间表如下11：MCNP-3：1983年写成，为标准的FORTRAN77版本，截面采用ENDF/B2III。MCNP-3A：1986年写成，加进了多种标准源，截面采用ENDF/B2IV。 MCNP-3B：1988年写成，具有阵列几何处理能力(即重复结构描述)，多群截面和计数输出的图形化功能，截面采用ENDF/B2IV和ENDL2851。MCNP-4：1990年7月由LANL写成，截面采用ENDF/B2V。MCNP-4.2：1991年3月由ORNL的RSIC写成，程序有较大改进，增加了基于Sandia国家实验室的ITS(Integrated Tiger Series)2连续

46、能量电子输运包，将其编入MCNP程序，专用于UNIX系统,从此MCNP程序成为中子/光子/电子耦合输运程序。MCNP-4A：1993年诞生，仍为UNIX系统，开始引入PVM并行，适合共享存储并行计算机，截面为ENDF/B2V。MCNP-4B：1997年3月正式推出，有PC版(需要LAHEY编译系统支持)，UNIX版，采用ENDF/B2VI截面库和彩色图形系统，仍采用PVM并行编程。MCNP-4B2：为MCNP4B的升级版，其支持FORTRAN90系统。MCNP4C：2000正式推出,在MCNP4B基础上增加共振自屏、瞬发本征值、微扰和多群伴随中子输运计算等处理，采用F90编译器，工作站版本支持PVM和SMPP并行。 MCNP-5：2003年推出。在这个版本中，完成了从FORTRAN277到FORTRAN290的重新组织，支持以前的MCNP-4C2/4C3全部功能，同时在提高图形显示，易安装性以及更好的在线文档方面有较大改善。 另外，在MCNP系列版本中，出现了MCNPX版本,该版本程序仍为该实验室研制