五章蒙特卡罗方法在计算机上的实现.ppt

五章蒙特卡罗方法在计算机上的实现 Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望第五章第五章 蒙特卡罗方法在计算机上的实现蒙特卡罗方法在计算机上的实现 蒙特卡罗方法是随着计算机的出现和发展而逐步发展起来的。在计算机上能够产生符合要求的随机数，实现对已知分布的抽样，奠定了蒙特卡罗方法在计算机上得以实现的基础。在计算机上使用蒙特卡罗方法解粒子输运问题大致包括三个过程：源分布抽样过程，空间、能量和运动方向的随机游动过程以及记录、分析结果过程。1.源分布抽样过程源分布抽样过程 源分布抽样的目的是产生粒子的初始状态 。下面我们介绍一些常见的特定类型的源分布抽样方法。1)源粒子的位置常见分布的随机抽样(1)圆内均匀分布设圆半径为R0，粒子在圆内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：r 的抽样方法为：(2)圆环内均匀分布设圆环的内半径为R0，外半径为R1，则粒子在该圆环内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：r 的抽样方法为：(3)球内均匀分布设球的半径为R，粒子在球内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：r 的抽样方法为：在直角坐标系下，抽样方法为：(4)球壳内均匀分布设球壳的内半径为R0，外半径为R1，在均匀分布时，从发射点到中心的距离 r 的分布密度函数为：r 的抽样方法为：在直角坐标系下，球壳内点的坐标为：其中，r 由前面的抽样方法确定，、服从各向同性分布，其抽样方法为：(5)圆柱内均匀分布圆柱内均匀分布是指粒子发射点均匀地分布在底半径为 R，高为 2H 的圆柱内。若固定圆柱的中心为原点，圆柱的轴向为 z 轴，则分布密度函数为：抽样方法为：(6)点源分布 点源分布是指粒子由一固定点 发射，其分布密度函数为：其中，为狄拉克函数，源粒子的抽样方法为：在球坐标系中，粒子发射点到球心的距离 r 的分布密度函数为：其中，为点源到球心的距离。源粒子的位置抽样为：(7)球外平行束源分布球外平行束源分布是指粒子平行入射到半径为 R 的球面上，或球外点源距离球很远，可以近似地看作平行束源。设 r 为粒子发射点到球心的距离,其分布密度函数为：r 的抽样方法为：在直角坐标系中，抽样方法为：2)源粒子的能量常见分布的随机抽样(1)单能源分布单能源分布是指粒子的发射能量为一固定值 E0，其分布密度函数为：源粒子的能量为：(2)裂变中子谱分布裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax 均为与元素有关的量。对于铀-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。采用近似修正抽样，抽样方法为：其中，m0.8746，M10.2678，0.5543。此外，裂变谱分布也有以数值曲线形式给出的，此时，用数值曲线抽样方法抽取 E。(3)麦克斯韦(Maxwell)谱分布麦克斯韦谱分布的一般形式为：该分布的抽样方法为3)源粒子运动方向常见分布的随机抽样(1)各向同性分布各向同性分布密度函数为：其中，cos，为运动方向与 z 轴的夹角，为方位角。在直角坐标系下，各方向余弦 u，v，w 为：其抽样方法为：(2)半面各向同性分布不妨设在 x0 的半面方向上各向同性发射粒子，则在前述各向同性分布的抽样方法中，用2代替2就能得到所需分布的抽样。对于其它方向的情况，可用类似的方法处理。(3)球外平行束源分布令cos，为粒子运动方向的径向夹角，则分布密度函数为：的抽样方法为：(4)球外各向同性点源分布设球外点源 S 到球心的距离为D0。点源 S 到球的最大张角为*，则球外各向同性点源分布的抽样方法是：先抽样确定 ，再转换成。在直角坐标系下，取 OS 为 z 轴，抽样方法为：4)次级粒子的源分布 在有关次级粒子（如裂变中子，中子生成光子，光子生成中子）的输运过程中，次级粒子源分布的抽样方法，主要可分为以下两种：(1)直接生成法 可将生成的次级粒子的位置、能量、方向、权重等参数直接作为源分布的抽样结果。也就是直接对生成的次级粒子进行跟踪。这种方法比较简单、直观。(2)离散分布法 将生成的次级粒子的权重，按空间位置、能量、方向分别记录，得到次级粒子的空间、能量、运动方向的离散的近似分布。再根据该分布，利用各种抽样技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、记录。当一个问题需要用两个以上的蒙卡程序处理时，可采用这种方法。2.空间、能量和运动方向的随机游动过程空间、能量和运动方向的随机游动过程 粒子由状态Sm到状态Sm+1时，需要确定粒子的空间位置 rm+1，能量 Em+1和运动方向m+1。1)碰撞点位置的计算公式 设 rm 为粒子第 m 次碰撞点的位置，m 为碰撞后的运动方向，则粒子第 m+1 次碰撞点的位置 rm+1 为：即其中 为 的方向余弦，L 为两次碰撞点间的距离。L 的分布密度函数为：由 f(L)抽样确定 L 的方法通常有三种：(1)直接抽样方法确定 L 的直接抽样方法是：首先由自由程分布中抽取再由下列关系式解出 L。对于均匀介质，有对于多层介质，如果则其中，为粒子由 rm 出发，沿m 方向在顺序经过的第 i 个介质区域内走过的距离，为第 i 个介质区域的宏观总截面(i=1，2，Imax)。当时，意味着粒子穿出系统。(2)最大截面法对于多层介质，或其他介质密度与位置有关的问题，在求 (i=1，2，Imax)时，如果系统形状复杂，计算是非常烦杂的。在这种情况下，使用最大截面法更方便。最大截面抽样方法为：其中(3)限制抽样法 当介质区域很小时，如使用直接抽样法抽取输运长度，粒子很容易穿出介质，此时使用限制抽样法确定自由程个数较好，的分布密度函数为：其中 Dm 为粒子由rm 出发，沿m 方向到达区域边界的自由程个数。的抽样方法是：然后用直接抽样法中根据计算 L 的方法计算输运长度 L。此时，粒子的权重需乘以纠偏因子 。2)碰撞后能量Em+1的随机抽样 粒子在介质中发生碰撞后，首先要确定与哪种原子核发生何种反应。粒子发生碰撞后（吸收除外）的能量 Em+1 一般只与其碰撞前后运动方向的夹角（散射角）有关。粒子碰撞后常见的能量分布有下面几种情况。(1)裂变中子谱 中子引起原子核裂变反应时，裂变中子的能量服从裂变谱分布。其抽样方法可参考以前的介绍。(2)中子弹性散射后能量的确定 中子弹性散射后，能量与质心系散射角C的关系是：能量与实验室系散射角L的关系是：其中，A 为碰撞核的质量，。或 确定后，即可求出 Em+1。(3)中子非弹性散射后能量的确定 中子非弹性散射后，能量与质心系散射角C的关系是：其中，为第 K 个能级的阈能，为第 K 个能级的激发态能量。如果确定了实验室系散射角L，则根据下式确定 后，再计算出 Em+1。(4)光子康普顿（Compton）散射后能量的确定光子发生康普顿散射后，其能量分布密度函数为：其中，K()为归一因子。，和 分别为光子散射前后的能量，以 m0c2 为单位，m0为电子静止质量，c 为光速。光子康普顿散射能量分布的抽样方法为：x 的抽样确定后，散射后的能量为：3)碰撞后散射角的随机抽样粒子碰撞后运动方向m+1的确定，一般与散射角有关。由已知分布抽样确定散射角后，再确定m+1。常见的散射角分布有如下几种：(1)质心系各向同性分布散射角在质心系服从各向同性分布时，其抽样方法为 。质心系散射角C抽样确定后，需转换成实验室系散射角L。在中子弹性散射情况下，转换公式为：其中 A 为碰撞核质量，。在中子非弹性散射情况下，转换公式为：其中，为第 K 个能级的阈能。(2)中子弹性散射勒让德(Legendre)多项式分布 中子弹性散射角分布常以勒让德多项式的展开形式给定。散射角余弦 xcos的分布密 度函数为：其中 Pl(x)为 l 阶勒让德多项式。该分布即为 n 阶勒让德近似展开。勒让德多项式由以下递推公式确定：考虑新的分布：当选取 x0，x1，xn 为 Pn+1(x)0 的根，且时，fa(x)依照勒让德多项式展开的前 n 项与 f(x)的展开形式相同。因此，可以用 fa(x)作为 f(x)的近似分布。在实际问题中，由于勒让德多项式展开项数不够，可能出现某个 为负值的现象。此时可以采用如下近似分布：其中：对于该近似分布，可用加抽样方法进行抽样：此时，由于偏倚抽样而引起的纠偏因子为 wK，也就是说，粒子的权重要乘上wK。(3)光子康普顿散射角分布光子的康普顿散射角与其散射前后的能量有关,它的分布密度函数为：抽样方法为：4)碰撞后运动方向m+1的确定实验室系散射角L确定后，依据不同的坐标系的表现形式，有不同的确定方法。(1)确定方向余弦 um+1，vm+1，wm+1其中，方位角 在 0,2 上均匀分布。当 时，不能使用上述公式，可用下面的简单公式：(2)确定m+1的球坐标(m+1，m+1)设m的球坐标分别为(m,m)，其中，为粒子运动方向与 z 轴的夹角，为粒子运动方向在 x y 平面上投影的方位角。则m+1的球坐标(m+1,m+1)分别由下式确定：5)球形几何的随机游动公式 一般几何的随机游动公式可以应用到球形几何，而对球对称问题，使用特殊形式更为方便。(1)下次碰撞点的径向位置 rm+1的确定 两次碰撞点间的距离 L 确定之后，下次碰撞点的径向位置 rm+1的计算公式为：设系统的外半径为R，如 rm+1R，则粒子逃出系统。(2)粒子碰撞后瞬时运动方向的确定 在球对称系统中，粒子运动方向用其与径向夹角余弦来描述。使用球面三角公式，粒子碰撞后瞬时运动方向与径向夹角余弦 cosm+1的计算公式为：其中，为在 0,2 上均匀分布的方位角，为在 rm+1 点进入碰撞前瞬时运动方向与 rm+1 径向之间的夹角。6)点到给定边界面的距离在抽样确定输运距离、判断粒子是否穿透系统时，常遇到求由 rm 出发，沿m 方向到达某个区域表面的距离问题。在记录对结果的贡献时，也常使用类似的量。区域表面通常是平面或二次曲面。求到达区域表面的距离问题，实际上是求直线（或半射线）与平面或二次曲面的交点问题。这是 蒙特卡罗方法解粒子输运的各种实际问题时,所遇到的基本几何问题。(1)点到平面的距离点 沿方向 的直线方程为：该直线到达方程为的平面的距离为：当 与平面平行时，即直线与平面无交点。如果 l 为负值，直线与平面也无交点。这时，粒子的运动方向是背离平面的。(2)点到球面的距离在三维直角坐标系中，设球心为 rc(xc,yc,zc)，球半径为 R，则球面方程为：将直线方程代入该球面方程，得到点 r0沿0方向到达球面的距离 l：其中当 R0R 时，即 r0 点在球内，2，l 只有一个正根：当 R0R 时，即 r0 点在球外，分以下三种情况：a)若0，l 无正实根，直线与球面无交点。b)若0，0，l 无实根，直线与球面无交点。c)若0，0，l 有两个正实根，直线与球面有两个交点。在球坐标系中，不失一般性，设球心为 rc0，则球面方程为 rR。当 r0R 时，即 r0 点在球内，有一个交点：其中0为0与 r0 的径向夹角。当 r0R 时，即 r0 点在球外，令当 cos00 时，直线与球面无交点。当 cos00 时，若 dR，则直线与球面无交点。若 dR，则有两个交点：(3)点到圆柱面的距离设圆柱面的方程为：其中(xc,yc,0)为圆柱的中心，R 为圆柱底半径。点 r0沿0方向到达圆柱面的距离 l 为：其中 当 R0R 时，r0 点在圆柱内，如果 ，则 l 有一个正根：如果 ，即0平行于圆柱的对称轴，直线与圆柱面无交点。当 R0R 时，r0 点在圆柱外，分以下三种情况：a)若0，l 无正实根，直线与圆柱面无交点。b)若0，0，l 无实根，直线与圆柱面无交点。c)若0，0 且 ，l 有两个正实根，直线与圆柱面有两个交点。d)在 的情况下，直线与圆柱面不相交。(4)点到圆锥面的距离 设圆锥顶点在原点，以 z 轴为对称轴，则圆锥面的方程为：点 r0沿0方向到达圆锥面的距离 l 为：其中如果0与锥面某一母线平行，即 ，则(5)空腔处理 在粒子输运问题中，所考虑的系统常有空腔存在，如中空的球壳,平板间的空隙等。粒子输运时，可有两种处理空腔的方法：a)将空腔作为宏观总截面t0 的区域,按通常的方法输运。b)设 分别为由 rm 出发，沿m 方向到空腔区域的近端和远端的交点，当粒子超过 时，以 为新的起点，重新开始输运。显然，这两种方法在统计上是等价的。7)等效的边界条件(1)全反射边界在反应堆活性区中，元件盒常常按正方形或六角形排列。假定元件盒足够多，每个盒结构相同，那么活性区中每个盒所占的栅胞的物理情况，可以代表整个活性区中的状况。进一步假定，元件盒是圆对称的，那么每个栅胞中情况，可以用更小的单位（栅元）来反映。比如对六角形栅胞可取其 1/12 的OAB 来做代表；正方形栅胞可用其 1/8 的OAB 来做代表。这样一来问题就大大简化了。现在的问题是怎样计算直角三角形栅元的物理量（如通量）。用蒙特卡罗方法如何模拟中子在栅元内的运动，反映出整个活性区对它的影响。我们可把OA、OB、AB 作为全反射边界来处理。所谓全反射边界，就是当中子打到该边界上时，按镜面反射的方式，从边界 上全部反射回来，中子的能量与权重均不改变。在这种边界上的反射条件，称之为全反射条件，就是通常的镜面反射条件。在全反射边界条件下，一条通过活性区若干个区域的中子径迹，可以用栅元OAB 中的一条折线轨 道来反映出来。反过来，在直角三角形栅元OAB 中任一条反射成的折线轨道，都代表了中子在活性区内一条直线轨道的作用。由于系统的对称性，在活性区内，凡是与栅元内位置相当的地方，都有相同的物理情况，因此栅元内各处的情况，当然代表了整个活性区的情况。(2)一般曲面全反射条件 对于一般曲面的全反射，设入射方向为，入射点的内法线方向为 n，则反射方向 为：其中设则(3)平面全反射条件 设三角形栅元的横截面OAB 在 X-Y 平面上，OAB。则边界 OA、OB、AB 上的反射都是平面全反射。在任一与 X-Y 平面垂直且与 X 轴成角的平面上，全反射条件为：由此就可得到OA、OB 和 AB 边上的全反射条件，对于 OB 边，=；对于 OA 边，=0；对于 AB 边，=/2。(4)反射层边界条件 对于具有大反射层的系统，如存放，运输和生产裂变物质的仓 库、车厢和车间等，当中子从里面打到四周墙上或反射层时，还要继续对它进行跟踪。这种跟踪常常要花费很大的计算量，并且在结果中引起的方差也比较大。如果在计算这种系统的不同方案中，反 射层条件不变，那么这种大量重复的计算是很不经济的。中子射入反射层后，一部分被介质吸收，只有一部分返回，由于中子的散射慢化，损失一部分能量，因此反射回来的中子有一个能量方向分布。显然，对这种反射层，不能应用全反射条件。不过，我们仍然可以把它当做边界，在边界上按反射层的物理作用来处理。比如，如果反射层是一种平板几何，我们可以用数值方法或蒙特卡罗方法，预先算好在各种不同入射能量 E 下的反照率(E)，反射中子的能量分布 RE(EE)。于是代替在反射层中眼踪中子，我们可在反射层边界上作如下处理：一旦中子打入反射 层，立即返回，反射后权重为其中，E 为射入反射层中子的能量，W 为中子的权重。反射后的能量 E 由反射能谱 RE(EE)中抽样产生。反射后的方向 由半平面各向同性分布或余弦分布中抽样。反射后的中子位置为入射时的位置。计算表明，对于大尺寸的反射层来说，这样的近似，引 起的结果上的误差是可以忽略的，却能带来计算量的大量节省。3.记录贡献与分析结果过程记录贡献与分析结果过程 在粒子输运问题中，除了得到某些量的积分结果外，还需要得到这些量的方差、协方差、以及这些量的空间、能量、方向和时间的分布。这些量可以利用分类记录手续同时得到。1)记录与结果 为了得到所求量的估计值，在粒子输运过程中需进行记录，即求每个粒子对所求量的贡 献。设模拟了 N 个粒子，所求量的估计值为：其中 gi 为第 i 个粒子的总贡献。记录的贡献由所求量决定。对于同一个所求量，又随所用的蒙特卡罗技巧的不同而不同。例如，所求量是粒子穿透屏蔽概率，使用直接模拟法时，如粒子穿透屏蔽，在叠加记录单元加“1”(初始值为零)，否则没有贡献。使用加权法时，如粒子穿透屏蔽，在叠加记录单元加粒子的权重，否则没有贡献。使用统计估计法时，粒子每发生一次碰撞(包括零次碰撞)，都要记录贡献，等等。2)方差和协方差的估计估计量 g 和 g 的方差和协方差为：它们可以用下式估计：因此，要得到 和 的估计，只要对每一个历史记录结果的 和 进行记录，并加以累加即可。方差估计值 确定后，可得到误差其中 为置信限，它随置信水平 而定。在通常情况下，取 。3)位置、能量、方向、时间分布 在前面已经提到，用蒙特卡罗方法求某种量的空间、能量、方向和时间分布，实质上是得到这种分布的阶梯函数近似的估计值。而求这种估计值是很方便的，只要将跟踪过程中所得到的感兴趣量，按其状态的空间、能量、方向、时间特征，分别记录其权 重，最后将这些记录结果适当处理即可。事先，将问题的空间、能量、方向（常按相对于某个方向的夹角余弦）、时间范围，各分为如下不同间隔：再用一批存贮单元 A 记录相应间隔上阶梯函数近似的累计值。4.核截面数据的引用核截面数据的引用 用蒙特卡罗方法解粒子输运问题，需要介质所包含的各种原子核的核数据。以中子核数据为例，需要各种涉及到的核的微观总截面、弹性散射、非弹性散射、n-2n 反应、裂变、俘获等截面；也需要这些反应的相应能量、角度分布、次级粒子数，以及其它关心的粒子数及其能量、方向分布。从输运方程中可以知道，有了这些数据，问题就完全确定了；反映到蒙特卡罗模拟中，有了这些数据，就能决定宏观总截面，决定碰撞核的具体形式，就能实现抽样和跟踪。在蒙特卡罗计算中，引用的核数据有点截面、分段常数截面和群截面三种形式。1)点截面形式 在跟踪粒子时，对粒子的每一种能量，先从截面库中取出需要的核数据，再用插值（或其它方式），求出相应能量的各种截面数据。这种方法是直接的，也是比较精确的。不少通用程序就是这样做的。这样做的条件是要有完备的核截面数据库，计算机有大而快的存贮能力。2)分段常数形式 将问题的能量范围（Emin，Emax）分成许多间隔，截面数据在每个间隔上看作与能量无关，即截面取分段常数形式。这种引用截面数据的好处是，数据量相对地少。问题是它要根据问题的物理特点来划分能量间隔。而为了保证误差较小，所取的间隔数一般是比较多的。3)群截面形式 解中子输运问题，常常采用多群近似方法。在多群近似下，中子能量 E 用群指标 i 代替。为了实现蒙特卡罗跟踪，只需要以下群截面：显然，这种跟踪过程数据量小，程序简单。5.蒙特卡罗程序结构蒙特卡罗程序结构在电子计算机上，蒙特卡罗方法解粒子输运问题的程序，一般都可分为：源抽样，空间输运过程，碰撞过程，记录过程和结果的处理与输出等部分。开始数据预处理,各记录单元清零取一个粒子历史源分布抽样输运过程碰撞过程历史终止否？统计处理做完给定历史数否？结果的处理与输出终止记录过程记录过程记录过程记录过程 至于粒子历史终止条件，根据问题的几何条件、物理假定，处理方法，可归纳为以下几种：1)粒子从系统逃脱；2)粒子经碰撞被吸收；3)经俄国轮盘赌后，历史被终止；4)粒子能量低于给定能量（阈能）;5)粒子位置越过某一界面；6)粒子飞行时间超过给定时间；7)粒子权重小于某个小量。作作 业业 1)