蒙特卡洛方法在中子输运中的应用.doc

1、.中子输运理论与数值方法课程作业蒙特卡洛方法.目录1.前言 .32.蒙特卡洛方法概述 .42.1蒙特卡洛方法的基本思想 .52.2蒙特卡洛方法的收敛性、误差 .52.2.1蒙特卡洛方法的收敛性 .62.2.2蒙特卡洛方法的误差 .62.3蒙特卡洛方法的特点 .72.4蒙特卡洛方法的主要应用围 .93.随机数 .93.1线性乘同余方法 .113.2伪随机数序列的均匀性和独立性 .123.2.1伪随机数的均匀性 .123.2.2伪随机数的独立性 .12.4.蒙特卡洛方法在粒子输运上的应用.124.1屏蔽问题模型 .134.2直接模拟方法 .144.2.1状态参数与状态序列 .144.2.2模拟运动

2、过程 .154.2.3记录结果 .184.3蒙特卡洛方法的效率 .205.蒙特卡洛方法应用程序 MCNP.205.1 MCNP简述 .205.2 MCNP误差的估计 .225.3 MCNP效率因素 .236.结论.23参考文献 .241. 前言半个多世纪以来，由于科学技术的发展和电子计算机的发明，蒙特卡洛( Monte Carlo ) 方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。 蒙特卡洛方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡洛方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，

3、因而该方法的应用领域日趋广泛。蒙特卡洛模拟计算是解决中子在介质中输运较为成熟、有效的方法，对于原子能、辐射防护、剂量学和辐射生物物理学等研究领域.实际问题的计算，都可以利用蒙特卡洛方法予以实现。粒子输运过程可以用玻耳兹曼方程加以描述，然而，以此基础上发展起来的近似数值方法如扩散近似法、 离散坐标方法在处理截面与能量相关以及散射各向异性介质、复杂几何条件问题时碰到了较大困难。而蒙特卡洛方法在处理这类问题时得心应手，有很强的解题能力，并且近似较少，接近于真实情况。粒子辐射问题计算通常有输运方程法、蒙特卡洛法 (MC法 ) 、实验测量法以及经验法等几种方法。 蒙特卡洛计算法又称随机抽样法或统计试验法

4、，是基于计算机模拟的思想， 抓住物理过程的数量和几何特征，进行数字模拟试验， 该方法是求解辐射输运问题的一种相当成熟和有效的方法，而且它对于各种复杂问题， 具有良好的通用性， 实用性相当广泛， 几乎涉及核科学的各个领域。本文主要介绍蒙特卡洛的概念、原理和应用及研究现状。2. 蒙特卡洛方法概述蒙特卡洛方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡洛方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡洛方法能够比较逼真地描述事物

5、的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。蒙特卡洛方法的主要组成部分有：(1) 概率密度函数 (pdf) 必须给出描述一个物理系统的一组概率密度函数；(2) 随机数产生器能够产生在区间 0,1 上均匀分布的随机数；(3) 抽样规则如何从在区间 0,1 上均匀分布的随机数出发 , 随机抽取服从给.定的 pdf 的随机变量；(4) 模拟结果记录记录一些感兴趣的量的模拟结果；(5) 误差估计必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；(6) 减少方差的技术利用该技术可减少模拟过程中计算的次数；(7) 并行和矢量化可以在先进的并行计算机上运行的有

6、效算法2.1 蒙特卡洛方法的基本思想可以通俗地说，蒙特卡洛方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r) 的随机变量 (r) 的数学期望gg(r ) f (r )dr(0.1)0通过某种试验，得到N 个观察值 r 1，r 2， ， r N（用概率语言来说，从分布密度函数 f(r)中抽取 N个子样 r 1，r 2 ， ，r N，），将相应的 N个随机变量的值 g(r 1) ，g(r1N) ， ， g(r ) 的算术平均值 gNg(ri ) ，作为积分的估计值（近似值） 。2NN i 1为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的

7、试验相当困难， 甚至是不可能的。 因此，蒙特卡洛方法的基本思想虽然早已被人们提出， 却很少被使用。 本世纪四十年代以来， 由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成， 使得蒙特卡洛方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。2.2 蒙特卡洛方法的收敛性、误差蒙特卡洛方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。.2.2.1蒙特卡洛方法的收敛性由前面介绍可知，蒙特卡洛方法是由随机变量X 的简单子样X1， X2， ， XN的算术平均值1NX ，X NXi . 作为所求解的近似值。由大数定律可知，如1N i 1X

8、2， ， XN 独立同分布，且具有有限期望值，则P lim X NE( X )1。即随机N变量 X 的简单子样的算术平均值X N ，当子样数 N充分大时，以概率1 收敛于它的期望值 E(X) 。2.2.2蒙特卡洛方法的误差蒙特卡洛方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1 ，X2， ， XN 独立同分布，且具有有限非零的方差 2 ，即 02(xE( X ) 2 f ( x)dx。 f (X) 是 X 的分布密度函数。则PNx1xt2 /2dt(0.2)XN E(X)elimN2x当 N充分大时，有如下的近似式P X NE(X)2e t 2 /

9、2 dt1(0.3)N20其中 称为置信度，1 称为置信水平。这表明，不等式 X NE(X )N近似地以概率 1 成立，且误差收敛速度的阶为O ( N 1/2 ) 。通常，蒙特卡洛方法的误差 定义为(0.4)N上式中与置信度 是一一对应的，根据问题的要求确定出置信水平后，.查标准正态分布表，就可以确定出。常用的 与 的对应关系为： =0.5 ，=0.6745 ； =0.05 ，=0.96 ； =0.003 ，=3. 蒙特卡洛方法的误差为概率误差，这与其他数值计算方法是有区别的。误差中的均方差 是未知的，必须使用其估计值1 N21 NX i )2(0.5)?Xi(N i 1N i 1来代替以求出

10、均方差。由式 (0.4) 可知当给定置信度 后，误差 由 和 N决定。要减小 ，或者是增大 N，或者是减小方差2。在 固定的情况下，要把精度提高一个数量级， 试验次数 N需增加两个数量级。 因此，单纯增大 N 不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半， 这相当于 N增大四倍的效果。 因此降低方差的各种技巧，引起了人们的普遍注意。2.3 蒙特卡洛方法的特点作为一种统计试验方法，蒙特卡洛方法因其优点在诸多领域有着广泛，但同时存在一些缺点。蒙特卡洛的主要优点有：(1) 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。蒙特卡洛方法可以部分代替物理实验

11、， 甚至可以得到物理实验难以得到的结果。 用蒙特卡洛方法解决实际问题， 可以直接从实际问题本身出发， 而不从方程或数学表达式出发。它有直观、形象的特点。(2)受几 何条 件限 制小。在 计算 s 维 空间中的 任一 区域D 上的积分sg Lg( x1 , x2 ,L , xs )dx1 dx2 L dxs时，无论区域 Ds 的形状多么特殊，只要能给出Ds描述 D的几何特征的条件， 就可以从 D( i)(i )( i ) ，得到中均匀产生 N 个点 ( x1, x2,L , xsss.N积分的近似值 gNDsg( x1(i ) , x2(i ) ,L , xs(i ) ) ，其中 Ds 为区域

12、Ds 的体积。这是数值N i 1方法难以作到的。(3) 收敛速度与问题的维数无关。由误差定义可知，在给定置信水平情况下，蒙特卡洛方法的收敛速度为 O( N 1/2 ) ，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。 也就是说， 使用蒙特卡洛方法时，抽取的子样总数N与维数 s 无关。维数的增加， 除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡洛方法对多维问题的适应性。(4) 具有同时计算多个方案与多个未知量的能力。对于那些需要计算多个方案的问题，使用蒙特卡洛方法有时不需要像常规方法那样逐个计算， 而可以同时计算所有的方案， 其全部计算量几

13、乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。(5) 误差容易确定。对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情， 而蒙特卡洛方法则不然。 根据蒙特卡洛方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡洛方法计算问题，也是容易确定的。(6) 程序结构简单，易于实现。在计算机上进行蒙特卡洛方法计算时，程序结构简单，分块性强，易于实现。蒙特卡洛的主要缺点有：(1) 收敛速度慢。如前所述，蒙特卡洛方法的收敛速度为O ( N 1/

14、2 ) ，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方.法好。(2) 误差具有概率性。 由于蒙特卡洛方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。(3) 在粒子输运问题中，计算结果与系统大小有关。经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右） ，蒙特卡洛方法计算的结果较为满意。 但对于大系统或小概率事件的计算问题， 计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡洛方法时，可以考虑把蒙特卡洛方法与解析（或数值）方法相结合，取长补短。2.4 蒙特卡洛方法的主要应

15、用围蒙特卡洛方法所特有的优点，使得它的应用围越来越广。它的主要应用围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用围将更加广泛。蒙特卡洛方法在粒子输运问题中的应用围主要包括：实验核物理，反应堆物理，高能物理等方面。 蒙特卡洛方法在实验核物理中的应用围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。3. 随机数随机数是蒙特卡洛方法的主要组成部分之一。随机数是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关。 在一个均匀分布的随机

16、数中，每一个体出现的概率是均等的。物理中的很多过程需要随机数确定，比如出射粒子的能量、方向等属性， 粒子与介质的相互作用等等。所模拟的物理过程要求随机数应具有下列特性：.1. 随机数序列应是独立的、互不相关的 (uncorrelated) ：即序列中的任一子序列应与其它的子序列无关；2. 长的周期 (long period) ：实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；3. 均匀分布的随机数应满足均匀性 (Uniformity) ：随机数序列应是均匀的、无偏的，即：如果两个子区间的“面积”相等，则落于这两个子区间的随机数的个数应相等。4. 有

17、效性（ Efficiency): 模拟结果可靠，随机数的产生必须快速、有效，最好能够进行并行计算。为了产生随机数，可以使用随机数表。随机数表是由0，1， ， 9十个数字组成，每个数字以 0.1 的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到 n位有效数字的随机数， 只需将表中每 n个相邻的随机数字合并在一起， 且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为 7634258910 ，要想得到三位有效数字的随机数依次为0.763 ，0.425 ，0.891 。可以使用物理方法产生随机数，用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，

18、另一种是利用计算机的固有噪声。但在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：n k T ( n , n 1,L , n k 1), n 1,2,L(3.1)产生随机数序列。对于给定的初始值。1, 2 ， k，确定 n+k， =1，2 。经常使用的是 k=1的情况。a) 用数学方法产生的随机数有两个特点，即：递推公式和初始值1, 2 ，. k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。b) 由于随机数序列是由递推公式确定的，而在计算机上所能表示的0,1 上的数又是有限的， 因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的 n，n(

19、n n) ，使得 n i n iiL,k 成1,2,立随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复， 其后面的随机数全部重复， 这与随机数的要不相符的。由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。3.1 线性乘同余方法线性乘同余方法是由 Lehmer在1951年提出来的，是一种最常用的产生伪随机数的方法。乘同余方法中采用的递推公式为I n 1( aI n c) mod m(3.2)其中：I 0为初始值， a为乘法器， c为增值， m为模数， mod为取模运算， (aI nc) 除以 m后的余数。 a、c、 m皆为整数。实型随机序列：I n0,1),

20、 I nm(3.3)rnfloat (m)I n0,1), I nm 1(3.4)rnfloat (m1)上式中，独立性和均匀性取决于参数a和c的选择。 m应尽可能地大，因为序列的周期不可能大于 m。通常将 m取为计算机所能表示的最大的整型量，在32位机上，m2312 109 。1961年， M. Greenberger 证明，用线性乘同余方法产生的随机数序列具有周期 m的条件是： (1)c 和 m互为质数； (2)a-1 是质数 p的倍数，其中 p是 a-1 和m的公约数； (3) 如果 m是 4的倍数， a-1 也是 4的倍数。.3.2 伪随机数序列的均匀性和独立性3.2.1伪随机数的均匀

21、性这里只考虑伪随机数序列1 , 2 ， n全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对于任意的0 x 1，令 Nn(x) 表示伪随机数序列 1, 2 ， n中适合不等式 i x,i=1,2,n的个数，则(n) sup | N n ( x)x |(3.5)0 x 1n将 伪 随 机 数 序 列 1, 2 ， n从 小 至 大 重 新 排列12Ln ， 令00, n 1 1，则由 ( n) 的定义，容易证明(n) max0x1|iii 1i| ，很n |,|n明显，对于固定的 ， ( n) 的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列， 均有不等式(n)

22、1 成立。当 (n) 1 成立时，2n2n所对应的伪随机数序列为最佳分布。3.2.2伪随机数的独立性对于任意0 x, y1，令 N n (x, y) 表示 ( 1, 2),(2, 3),( n, n+1) 中适合不等式 ix, i 1y 。的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令(n)supNn (x, y)Nn ( x) Nn ( y)(3.6)|0 x, y 1nnn则(n) 标志伪随机数序列 1, 2 ， n的独立程度，简称为独立偏度。对于固定的 n, (n) 的值越接近于零，伪随机数序列的独立性越好。4. 蒙特卡洛方法在粒子输运上的应用辐射（光子和中子）屏蔽问题是蒙

23、特卡洛方法最早广泛应用的领域之一。现.主要从物理直观出发， 说明蒙特卡洛方法解决这类粒子输运问题的基本方法和技巧。解决屏蔽问题时可采取多种方法，如直接模拟方法、简单加权法、统计估计法、指数变换法等，这里只对直接模拟方法做介绍。4.1 屏蔽问题模型在反应堆工程和辐射的测量与应用中，常常要用一些吸收材料做成屏蔽物挡住光子或中子。 我们所关心的是经过屏蔽后射线的强度及其能量分布，这就是屏蔽问题。当屏蔽物的形状复杂，散射各向异性，材料介质不均匀,核反应截面与能量、位置有关时， 难以用数值方法求解， 用蒙特卡洛方法能够得到满意的结果。粒子的输运问题带有明显的随机性质，粒子的输运过程是一个随机过程。粒子的

24、运动规律是根据大量粒子的运动状况总结出来的，是一种统计规律。 蒙特卡洛模拟，实际上就是模拟相当数量的粒子在介质中运动的状况，使粒子运动的统计规律得以重现。不过，这种模拟不是用实验方法，而是利用数值方法和技巧，即利用随机数来实现的。为方便起见，选用平板屏蔽模型，在厚度为a，长、宽无限的平板左侧放置一个强度已知，具有已知能量、方向分布的辐射源S，见图 4.1 。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源发出的平均一个粒子穿透屏蔽的数目。同时，假定粒子在两次碰撞之间按直线运动,且粒子之间的相互作用可以忽略。.图4.1屏蔽问题模型4.2 直接模拟方法直接模拟方法就是直接从物理问题出

25、发，模拟粒子的真实物理过程。4.2.1状态参数与状态序列粒子在介质中的运动的状态，可用一组参数来描述，称之为状态参数。它通常包括：粒子的空间位置r ，能量 E和运动方向，以 S(r , E , ) 表示。有时还需要其他的参数， 如粒子的时间 t 和附带的权重 W，这时状态参数为S (r ,E , , t ,W )。状态参数通常要根据所求问题的类型和所用的方法来确定。对于无限平板几何，取 S(z , E , cos) ，其中 z为粒子的位置坐标，为粒子的运动方向与 Z轴的夹角。对于球对称几何，取S (r , E , cos) ，其中 r 表示粒子所在位置到球心的距离，为粒子的运动方向与其所在位置

26、的径向夹角。粒子第 m次碰撞后的状态参数为Sm(rm , Em , m ) ，或 Sm(rm, Em, m ,tm ,Wm) ，它表示一个由源发出的粒子，在介质中经过m 次碰撞后的状态，其中r m：粒子在第m 次碰撞点的位置Em：粒子第m 次碰撞后的能量m：粒子第m 次碰撞后的运动方向t m：粒子到第m 次碰撞时所经历的时间Wm：粒子第m 次碰撞后的权重一个由源发出的粒子在介质中运动，经过若干次碰撞后，直到其运动历史结束（如逃出系统或被吸收等） 。假定粒子在两次碰撞之间按直线运动，其运动方向与能量均不改变，则粒子在介质中的运动过程可用以下碰撞点的状态序列描述，即 S0，S1， ， SM-1，S

27、M。或.r0 ,r1 , L,rM 1 ,rME0, E1, L , EM 1, EM0 ,1, L,M 1,M来描述。这里 S0为粒子由源出发的状态，称为初态，SM 为粒子的终止状态。 M称为粒子运动的链长。4.2.2模拟运动过程这里以中子穿透均匀平板的模型来说明，这时状态参数取S(z,E,cos ) 。模拟的步骤如下。(1) 确定初始状态 S0：确定粒子的初始状态，实际上就是要从中子源的空间位置、能量和方向分布中抽样。设源分布为f ( z0 , E0 ,cos0 )f1(z0 ) f 2 ( E0 ) f3 (cos0 ) ，则分别从各自的分布中抽样确定初始状态。(2) 确定下一个碰撞点：

28、已知状态 Sm-1 ，要确定状态 Sm， 首先要确定下一个碰撞点的位置zm。在相邻两次碰撞之间，中子的输运长度l 服从如下分布：f (l )t (rm 1 l m 1 , Em 1 )explt (rm 1 l m 1 , Em 1 )dl(4.1)0对于平板模型， l 服从分布：f (l ) t (zm 1 l cos m 1, Em 1 )explt ( zm 1lcos m 1, Em 1) dl (4.2)0其中， t 为介质的中子宏观总截面。积分lt (rm 1lm 1 , Em 1 )dl 称为粒子输0运的自由程数。 显然，粒子输运的自由程数服从指数分布，因此从 f ( l ) 中

29、抽样确lt (rm 1 l m 1 , Em 1 )dlln中解出 l 。对于单一介定 l ，就是要从积分方程0ln质 l，则下一个碰撞点的位置为t (Em 1).zm zm 1 l cos m 1 zm 1ln(4.3)cosm 1t (Em 1)如果 zma，则中子穿透屏蔽，若zm0,则中子被反射出屏蔽。这两种情况，均视为中子历史终止。(3) 确定被碰撞的原子核：通常介质由几种原子核组成，中子与核碰撞时，要确定与哪一种核碰撞。设介质由 A、B、C 三种原子核组成，其核密度分别为NA、 NB、NC，则介质的宏观总截面为：t (Em 1 )tA (Em 1 ) tB (Em 1 ) tC (E

30、m 1)(4.4)其 中tA ,Bt , Ct分别为核A、B、C的宏观总截面。其定义如下：(t ) ( Em 1 )N( ) t( ) ( Em 1 ) ，t( ) ( Em 1 )、N( )、 t( ) ( Em 1 ) 分别表示 ( ) 核的宏观总截面、核密度和微观总截面。由于中子截面表示中子与核碰撞可能性的大小，因此，很自然地，中子与A、B、C核发生碰撞的几率分别为：PAtA ( Em 1 ) ,PBtB ( Em 1) ,PCtC ( Em 1)(4.5)t ( Em 1 )t ( Em 1 )t ( Em 1)若PA ，则中子与 A核碰撞；若PA PB ，则中子与 B核碰撞；若PA

31、PB，则与 C核碰撞。(4) 确定碰撞类型：确定了碰撞的核 ( 比如 B核) 后，就要进一步确定碰撞类型。中子与核的反应类型有弹性散射、非弹性散射、(n,2n) 反应，裂变和俘获等，它们的微观截面分别为B(Em 1)、B(Em 1)、B、B( Em 1)和B( Em 1)，则有elin( n,2n ) ( Em 1 )fcB( Em 1 )BB( Em 1 )BB( Em 1)B( Em 1 ) (4.6)tel ( Em 1 )in(n ,2 n) ( Em 1)fc.各种反应发生的几率分别为PelPinP(n,2 n)PfPcelB (Em 1 )inB (Em 1 )B(n ,2n )

32、( Em 1 )Bf (Em 1 )cB (Em 1 )B tB tB tB t( Em 1 )( Em 1 )tB (Em 1)(4.7)( Em 1 )( Em 1 )利用离散型随机变量的抽样方法，确定反应类型。 在屏蔽问题中， 中子与核反应常只有弹性散射和吸收两种类型，吸收截面为：B (Em 1)B(E)B(Em 1) 。afm 1c这时，总截面为：tB ( Em 1 )elB ( Em 1 )aB ( Em 1 )(4.8)发生弹性散射的几率为：PelB ( Em 1 ) 若Pel ，则为弹性散射；否则为吸收，eltB ( Em 1 )发生吸收反应意味着中子的历史终止。(5) 确定碰撞

33、后的能量与运动方向：如果中子被碰撞核吸收，则其输运历史结束。如果发生弹性散射，需要确定散射后中子的能量和运动方向。中子能量Em为：EmEm 11 r1 r cos C(4.9)2其中 r( AA11) 2 。 A是碰撞核的质量与中子质量之比，一般就取元素的原子量； C为质心系中中子散射前后方向间的夹角，即偏转角。Ccos C 可从质心系中弹性散 射角 分布 f C( C) 中抽 样 产生 。 实 验室 系散射角 L 的 余弦 L 为：L1A C。A212 A C如果给出实验室系散射角余弦分布f L( L) ，可直接从 f L( L ) 中抽取 L，此时能量 E 与Em 1LA2122。确定了实

34、验室系散射的关系式为： EmLmL( A1)2角 L后，再使用球面三角公式确定cos m ：.cosmcos m 1 cos Lsinm 1sinLcos。各角度关系如图 4.2 所示。图 4.2角度关系示意图至此，由 S完全可以确定 S 。因此，当中子由源出发后，即S 确定后，重复步骤m-1m0(2) (5)，直到中子游动历史终止。于是得到了一个中子的随机游动历史S0，1M-1MS ， ，S，S，即z0 ,z1,L ,zM 1 ,zME0,E1,L ,EM 1,EMcos 0 ,cos 1,L ,cos M 1 ,cos M也就是模拟了一个由源发出的中子的运动过程。4.2.3记录结果在获得中

35、子的随机游动历史后，我们要对所要计算的物理量进行估计。对于屏蔽问题，我们要计算中子的穿透率。 考察每个中子的随机游动历史，它可能穿透屏蔽（ zM a），可能被屏蔽发射回来（ zM 0），或者被吸收。设第 n个中子对穿透的贡献为n，则1,当 zMan0,当 zM0N如果我们共跟踪了 N 个中子，则穿透屏蔽的中子数为： N1n 。则穿透屏蔽n1概率的近似值为：.?(1)N11N(4.10)PNNnNn 1我们称这种直观地模拟过程和估计方法为直接模拟方法。在置信水平1?(1)的误差为： 0.95 时， PN?(1)P1.96PN(4.11)NN其中为n 的均方差，由于n是一个服从二项分布的随机变量，

36、所以2P(1 P)(4.12)或2?(1)(1?(1)(4.13)?PNPN为得到中子穿透屏蔽的能量、角分布，将能量、角度围分成若干个间隔：EminEI LE1 E0 Emax ， 001LJ 2 。其中 Emax，Emin分别表示能量的上、下限，对于穿透屏蔽的中子按其能量、方向分间隔记录。设一穿透屏蔽的中子能量为 EM，其运动方向与 Z轴夹角为 M，若能量 EM属于第i个能量间隔Ei ，角度 M属于第 j个角度间隔j，则分别在第i 个能量计数器及第j个角度计数器中加 1。跟踪 N个中子后，则?(1)N1,ii1,2,L, I(4.14)PN1,iNEi?(1)N 2, jj1,2,L, J(

37、4.15)PN 2, jNj分别为穿透中子的能量分布和角分布。其中N1, i 和 N2, i 分别为第 i 个能量和第j个角度间隔的穿透中子数。归一后分别为：.?(1)*?(1)N1,iPN1,ii1,2,L, I(4.16)PN1, i?(1)*N1EiPN?(1)*?(1)N 2, jPN2, jj1,2,L, J(4.17)PN?(1)*N12, jPNj4.3 蒙特卡洛方法的效率衡量蒙特卡洛技巧的好坏，除了看其方差大小外，还要看其所需费用（计算时间）多少，即从该技巧的效率Ef （方差与费用乘积的倒数）全面考虑：Ef1(4.18)2T其中 2为方差， T为所需费用。 Ef 大时，所用方法

38、的效率高；否则，效率低。在一般情况下，直接模拟方法、简单加权法、统计估计法、指数变换法等方法中有些方法虽然减小了方差，却增加了费用。例如，加权法、统计估计法虽然较直接模拟方法减小了方差， 却使每个粒子的运动链长增加，或记录贡献的计算时间增加。因此，不能认为方差小的方法一定好，要从方法的效率全面考虑。在有些情况下，直接模拟方法仍然是一个被广泛使用的方法。5. 蒙特卡洛方法应用程序 MCNP5.1 MCNP简述MCNP(A General Monte Carlo Code for Neutron and Particle Transport)是一套通用的、三维空间中连续能量中子、光子和带电粒子（离

39、子）联合输运过程模拟程序，在军事和工业领域有着广泛应用。是基于蒙特卡洛方法的用于计算三维复杂几何结构中的中子、光子、电子或者耦合中子、光子、电子输运问题的通用软件包， 也具有计算核临界系统 （包括次临界和超临界系统）本征值问题的能力。该软件包通过 FORTRAN语言编程实现。MCNP程序具有超强的几何处理能力，几何系统由几何空间单元(cell)组成，.而几何空间单元的界面 (surface)由平面、二次曲面及特殊的四次椭圆环曲面组成。几何空间单元中的材料由包括同位素在的多种核素组成，使用精确的点截面参数，对特定的评价库 (ENDF/B-IV,V,V,VI库或 ENDL851库) ，考虑了该库给出的所有中子反应类型。 在截面数据文件中收集了多种评