蒙特卡罗方法及其在中子输运问题中得应用

蒙特卡罗方法及其在中子输运问题中得应用目录蒙特卡罗方法及其在中子输运问题中得应用 (1)1蒙特卡罗方法简介 (3)1.1蒙特卡罗方法的基本原理 (3)1.2 蒙特卡罗方法的误差 (4)2 随机变量的抽样方法 (4)2.1 直接抽样方法 (5)2.1.1 离散型随机变量的抽样 (5)2.1.2 连续型随机变量的抽样 (5)2.2 挑选抽样法 (5)2.3 复合抽样法 (6)3 蒙特卡罗方法模拟中子输运过程 (6)3.1 源抽样 (6)3.2 输运距离的抽样 (7)3.3 碰撞核素的抽样值 (7)3.4 反应类型的抽样值 (7)3.5 反应后中子状态的确定 (7)3.5.1 弹性散射 (7)3.5.2 非弹性散射 (8)3.5.3 裂变反应 (8)4 蒙特卡罗方法的减方差技巧 (8)4.1 权 (8)4.2 统计估计法 (9)4.3 权窗 (10)5 蒙特卡罗方法求解通量 (10)5.1 通量的定义 (10)5.2 点通量的计算 (11)5.3 面通量的计算 (11)5.3.1 统计估计法 (11)5.3.2 加权法 (12)5.4 体通量的计算 (12)5.4.1 统计估计法 (12)5.4.2 径迹长度法 (13)5.4.3 碰撞密度法 (13)5.4.4 几种体通量计算方法的比较 (14)5.5 最终结果的统计 (14)6 蒙特卡罗方法求解k eff (15)6.1 有效增值因子k eff的定义 (15)6.2 蒙特卡罗方法求解k eff (15)6.2.1 吸收估计法 (15)6.2.2 碰撞估计法 (15)6.2.3 径迹长度估计法 (16)1蒙特卡罗方法简介1.1蒙特卡罗方法的基本原理蒙特卡罗方法（Mento Carlo Method ）也叫统计模拟方法，是二十世纪四十年代由于计算机科学与技术发展和电子计算机的发明而提出来的一种基于概率论与数理统计的方法。

蒙特卡罗方法广泛应用与金融工程、经济学、粒子输运模拟、热力学与统计物理学等领域。

为了说明蒙特卡罗方法的基本原理，先看两个例子。

例1 用蒲丰方法求解π1977年，法国数学家蒲丰提出了一种计算圆周率的方法——随机投针法，及著名的蒲丰投针问题。

这一方法如下：1） 取一张白纸，在上面画两条间距为2a 的平行线；2） 取一根长度为2l （l <a ）的针，随即的向平行线间投掷n 次，观察针与平行线中得任一条相交的次数，记为m ；3） 计算针与平行线相交的概率n m p =可以证明这个概率的真实值为)/(2=d πl p ，也即当n 很大时，)/(2≈=d πl n m p 。

由此可以得到π的估计值)/(2ma nl π≈例2 求函数f (x )在[a ,b ]上的积分f (x ) [a ,b ]上的积分值I 就是曲线y =f (x )、y 轴、x =a 以及x =b 所围成的面积。

使用蒙特卡罗方法求解该问题的方法就是：如右图所示，随机在y 轴、y =M 、x =a 以及x =b 围成的矩形均匀投掷n 个点，如果落在y =f (x )下面的有m 个点，则曲线y =f (x )、y 轴、x =a 以及x =b 所围成的面积S矩形S nm S = 即积分值矩形S n m I =Ya由上面的两个例子可以看出，蒙特卡罗方法以一个“概率模型”为基础，将所求解的问题抽象为一个随机过程，使用已知分布抽样的方法求得试验结果的观察值，从而求得问题的近似解。

一、 蒙特卡罗方法及其模拟粒子在物质中的输运过程蒙特卡罗（Monte Carlo ）方法又叫随机模拟法，是一种以概率论与数理统计为基础的方法。

例 用蒙特卡罗方法求解π如下图所示，半径为1的圆外有一外接矩形，则圆与矩形的面积比为π/4。

我们向矩形内均匀地投掷N 个点，记录落在圆内部的点数M ，则有M/N ≈π/4由此可以得到π≈4M/N由上例可以看出，蒙特卡罗方法即是根据所求的问题构造一个概率模型，使得所求问题的解等于该模型的某个变量，根据求解该变量而得到问题的解。

1.2 蒙特卡罗方法的误差加入进行了N 次模拟（对于例1，即投递了N 次针，对于例2，即投掷了N 个点），这N 次模拟的标准差为σ，蒙特卡罗方法的误差为N Const σε×=2 随机变量的抽样方法以下叙述中，f (x ) 表示随机变量的概率密度函数，其分布函数为F (x )，ξ为[0, 1)之间均匀分布的随机数，即ξ~U [0,1)。

2.1 直接抽样方法2.1.1 离散型随机变量的抽样设离散型随机变量x 的分布律如下：∑∑-ij j i j j p p111==<≤ξ2.1.2 连续型随机变量的抽样若F (x )的反函数F -1(x )存在，则随机变量x 的抽样值为)(=1ξF x f -例如，如果x ~U [a ,b ]，即x 的分布函数为()⎪⎪⎩⎪⎪⎨⎧>≤≤--<=bx b x a ab a x a x x F 10 F (x )的反函数x a b a x F )(+=)(1--，所以x 的抽样值为ξa b a x f )(+=-如果随机变量的概率分布函数的反函数不存在或很难求得，就应该采用其它的抽样方法。

2.2 挑选抽样法若分布密度函数f (x )可以分解为一个分布密度函数f 1(x )与一个有界函数h (x )的乘积，即f (x ) = f 1(x )h (x )记M = maxh (x )，则随机变量x 可以采用如下抽样方法：1） 从f 1(x )抽样得到x f ；2） 判断M ξ≤h (x f )是否成立,若成立，则x f 有效并作为x 的抽样值，否则，回到1）直至满足2）为止。

例 标准正太分布的抽样标准正太分布的分布密度函数为())2/ex p(212x x f -π=令()2/)1(122)(---,x x e e x h e x f π==π2))(max (e x h M == 1) 抽取随机数ξ1，由f 1(x ) 抽样得到1ln ξ-=f x2) 抽取随机数ξ2，若)(22f x h e ≤ξπ，则x f 作为x 的抽样值，否则回到1）。

2.3 复合抽样法3 蒙特卡罗方法模拟中子输运过程一个中子的随机历史（中子历史是指中子从产生到消亡的过程）如上图所示。

中子从源产生后，飞行一段距离后进入了裂变材料区域，在该区域的①点，中子与靶核相互作用发生了散射，改变了飞行方向，在②点，中子与靶核作用发生了裂变反应，产生了两个次级中子。

其中一个中子在③处散射后飞出了裂变材料区，而另一个中子经过在④处的散射后被材料吸收，该中子的历史结束。

一个中子的飞行状态可以用(r , E, Ω)，即位置、能量、飞行方向等七个变量描述。

蒙特卡罗方法模拟粒子的输运过程便是用数学的方法计算粒子每行进一步的这七个量。

3.1 源抽样假定源的能量分布函数为f (E )，位置分布函数f (r )，飞行方向分布函数为f (Ω)。

则从源采样出来的中子能量、位置以及飞行方向分别为)()()(ξξr ξf 1-f 1-f F F F E =Ω==1-下面给出几种常见的源的抽样方法。

1 均匀分布的矩形平面源3.2 输运距离的抽样假定材料的宏观总截面为Σt ，则中子在该材料中的自由飞行距离的抽样值为tf l Σ=ξln - 3.3 碰撞核素的抽样值自由飞行了l f 段距离后，中子将于材料的核发生碰撞。

首先应确定中子与材料的那种核素发生碰撞。

假定，该材料又n 种核素组成，每种核素的宏观总截面分别是Σ1,Σ2…Σn 。

确定碰撞核素的方法如同离散型随机变量的直接抽样方法。

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧<≤∑∑-∑∑∑+∑<≤∑∑∑∑<≤=1......................02111ξξξtn t t t t n 种第第二种第一种与中子碰撞的核素 3.4 反应类型的抽样值确定了与哪一种核素碰撞以后，还需要确定与该核素发生了哪一种反应类型。

假定，该核素与中子的反应有弹性散射、非弹性散射、吸收、裂变等，其微观截面分别为σel , σin , σa , σf ，总的微观截面为σt 。

确定反应类型的方法如同离散型随机变量的直接抽样方法。

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧<≤-+<≤<≤=11......................0ξσσσσσξσσσσξtf t in el t el t el 裂变非弹性散射弹性散射反应类型 确定了反应类型后，下一步便是确定反应后中子的状态。

3.5 反应后中子状态的确定如果中子与靶核作用后，中子被吸收，则中子历史结束。

如果中子被散射，则需要计算中子散射后的状态。

3.5.1 弹性散射首先根据核数据库计算得到中子在质心系下色散射角θc 余弦μc ,进而根据下列公式计算在实验室系下的散射角θt 的余弦μt ，其中A 是靶核的原子量。

1212+++=c c t A A A μμμ进而，根据碰撞前中子的能量E in ，可以根据以下公式求得碰撞后的中子能量E out ，22)1(12+++=A A A E E c in outμ 3.5.2 非弹性散射 非弹性散射中子散射后的能量的方向都根据核素的核数据库计算。

3.5.3 裂变反应裂变反应平均产生2~3个中子，每个中子的散射后的能量和方向与非弹性散射一致。

每一次裂变产生的中子数根据公式⎣⎦ξ+v 确定。

ν是平均中子产额。

确定了中子与靶核作用后的状态后，回到3.2节所述步骤继续模拟中子，直至中子被吸收。

4 蒙特卡罗方法的减方差技巧正如1.2节所述，蒙特卡罗方法的误差正比于方差，并随着模拟的次数增大而减小。

实际中，如果只是增大模拟的次数，计算时间也会相应地加长。

所以降低误差更有效的方法是降低方差。

4.1 权权是蒙特卡罗方法中一个极为重要的概念。

在蒙特卡罗方法中，若一个中子的权为ω，则该中子效能相当于ω个实际的中子。

加权法有着重要的应用。

例如，中子与靶核反应后可能发生裂变，产生多于一个中子，也可能发生(n , 2n )，(n , 3n )反应。

这样连锁反应下去，直接模拟会十分麻烦，此时可以使用权来处理。

我们可以认为一个权为ω的中子与靶核反应后，仍旧只产生一个中子，只是权变为t fn n n n in el v σσσσσσωω++++=)3,()2,('32如图所示，如果要模拟一个强的中子吸收材料组成的隔离墙对中子的屏蔽效果。

假如我们模拟了N 个中子，最后有M 个中子穿过墙，则该材料的屏蔽效果为M /N 。

为了模拟得精确，N 就应该大一些，但这样会增加模拟时间。

另外，假如该材料的屏蔽效果为百万分之一，如果我们模拟了一千万个中子，平均会有10个中子穿过屏蔽墙。

但是，某次模拟时，可能会有12个中子穿过，而另一次模拟时有7个中子穿过，两次模拟出现了较大的统计涨落，方差就很大。