蒙特卡罗方法在计算机上的实现

第五章蒙特卡罗方法在计算机上的实现1.源分布抽样过程2.空间、能量和运动方向的随机游动过程3.记录贡献和分析结果过程4.核截面数据的引用5.蒙特卡罗程序结构作业第五章蒙特卡罗方法在计算机上的实现蒙特卡罗方法是随着计算机的出现和发展而逐步发展起来的。在计算机上能够产生符合要求的随机数，实现对已知分布的抽样，奠定了蒙特卡罗方法在计算机上得以实现的基础。在计算机上使用蒙特卡罗方法解粒子输运问题大致包括三个过程：源分布抽样过程，空间、能量和运动方向的随机游动过程以及记录、分析结果过程。1.源分布抽样过程源分布抽样的目的是产生粒子的初始状态。下面我们介绍一些常见的特定类型的源分布抽样方法。),,(0000ΩrSE1)源粒子的位置常见分布的随机抽样(1)圆内均匀分布设圆半径为R0，粒子在圆内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：其它当002)(020RrRrrf),max(210Rr(2)圆环内均匀分布设圆环的内半径为R0，外半径为R1，则粒子在该圆环内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：其它当02)(102021RrRRRrrf02010320101011)(),max()(RRRrRRRrRRRR≤(3)球内均匀分布设球的半径为R，粒子在球内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：在直角坐标系下，抽样方法为：其它当003)(32RrRrrf),,max(321Rr302010232221,,1RzRyRx≤(4)球壳内均匀分布设球壳的内半径为R0，外半径为R1，在均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：其它当03)(1030312RrRRRrrf001432322211020101211020201)(),,max(),max(33RxRRrxxxRRRRRRRRRRR≤≤在直角坐标系下，球壳内点的坐标为：其中，r由前面的抽样方法确定，θ、φ服从各向同性分布，其抽样方法为：cossinsincossin000rzryrx＞≤2322222123222221023222221310232222212101223222221cos2sinsin2cossin)(AAAArrzAAArryAAArrxAA(5)圆柱内均匀分布圆柱内均匀分布是指粒子发射点均匀地分布在底半径为R，高为2H的圆柱内。若固定圆柱的中心为原点，圆柱的轴向为z轴，则分布密度函数为：抽样方法为：其它当0||,21),,(2222HzRyxHRzyxf3020102221,,1HzRyRx≤(6)点源分布点源分布是指粒子由一固定点发射，其分布密度函数为：其中，为狄拉克δ函数，源粒子的抽样方法为：在球坐标系中，粒子发射点到球心的距离r的分布密度函数为：其中，为点源到球心的距离。源粒子的位置抽样为：)()()(),,(*0*0*0zzyyxxzyxf),,(*0*0*0zyx*0*0*0,,zzyyxx*0rr*0r)()(*0rrrf)((7)球外平行束源分布球外平行束源分布是指粒子平行入射到半径为R的球面上，或球外点源距离球很远，可以近似地看作平行束源。设r为粒子发射点到球心的距离,其分布密度函数为：r的抽样方法为：在直角坐标系中，抽样方法为：)()(RrrfRr20202020102221,1zyRxRzRy≤2)源粒子的能量常见分布的随机抽样(1)单能源分布单能源分布是指粒子的发射能量为一固定值E0，其分布密度函数为：源粒子的能量为：)()(0EEEf0EE(2)裂变中子谱分布裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax均为与元素有关的量。对于铀-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=∞。maxmin,sh)(EEEBEeCEfAE采用近似修正抽样，抽样方法为：其中，m≈0.8746，M1≈0.2678，λ≈0.5543。此外，裂变谱分布也有以数值曲线形式给出的，此时，用数值曲线抽样方法抽取E。3321121ln1)ln()(AAEEAEMEHm≤≤＞＞EAmEBECAAEHexpsh1)(21(3)麦克斯韦(Maxwell)谱分布麦克斯韦谱分布的一般形式为：该分布的抽样方法为0,2)(23EeEEfE＞22221ln23lnEe≤3)源粒子运动方向常见分布的随机抽样(1)各向同性分布各向同性分布密度函数为：其中，μ＝cosθ，θ为运动方向与z轴的夹角，φ为方位角。21)(,21)(41)()()(2121fffffΩ在直角坐标系下，各方向余弦u，v，w为：其抽样方法为：cossinsincossinwvu＞≤2322222123222221232222213123222221211223222221cos2sinsin2cossin)(AAAAwAAAvAAAuAA(2)半面各向同性分布不妨设在x≥0的半面方向上各向同性发射粒子，则在前述各向同性分布的抽样方法中，用ξ2代替η2就能得到所需分布的抽样。对于其它方向的情况，可用类似的方法处理。(3)球外平行束源分布令μ＝cosθ，θ为粒子运动方向的径向夹角，则μ分布密度函数为：μ的抽样方法为：01,2)(f),max(21u(4)球外各向同性点源分布设球外点源S到球心的距离为D0。点源S到球的最大张角为θ*，则球外各向同性点源分布的抽样方法是：先抽样确定，再转换成θ。0220*cosDRDRDR2202*sincos)cos1(1cos在直角坐标系下，取OS为z轴，抽样方法为：cos0sinwvu4)次级粒子的源分布在有关次级粒子（如裂变中子，中子生成光子，光子生成中子）的输运过程中，次级粒子源分布的抽样方法，主要可分为以下两种：(1)直接生成法可将生成的次级粒子的位置、能量、方向、权重等参数直接作为源分布的抽样结果。也就是直接对生成的次级粒子进行跟踪。这种方法比较简单、直观。(2)离散分布法将生成的次级粒子的权重，按空间位置、能量、方向分别记录，得到次级粒子的空间、能量、运动方向的离散的近似分布。再根据该分布，利用各种抽样技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、记录。当一个问题需要用两个以上的蒙卡程序处理时，可采用这种方法。2.空间、能量和运动方向的随机游动过程粒子由状态Sm到状态Sm+1时，需要确定粒子的空间位置rm+1，能量Em+1和运动方向Ωm+1。1)碰撞点位置的计算公式设rm为粒子第m次碰撞点的位置，Ωm为碰撞后的运动方向，则粒子第m+1次碰撞点的位置rm+1为：即其中为的方向余弦，L为两次碰撞点间的距离。mmmmmmmmmwLzzvLyyuLxx111mmmLΩrr1mmmmwvuΩ),,(L的分布密度函数为：由f(L)抽样确定L的方法通常有三种：(1)直接抽样方法确定L的直接抽样方法是：首先由自由程分布中抽取ρ再由下列关系式解出L。0,),(exp),()(01LdlElELfLmmmtmmtΩrrLmmmtdlEl0),(Ωref)(ln对于均匀介质，有对于多层介质，如果则其中，为粒子由rm出发，沿Ωm方向在顺序经过的第i个介质区域内走过的距离，为第i个介质区域的宏观总截面(i=1，2，…，Imax)。当时，意味着粒子穿出系统。)(ln)(mtmtEEL)()()()(,10,100,10,mItIimitiIiiIimitiIimitiEELLLELELiL)(,mitEmax0,)(IimitiEL(2)最大截面法对于多层介质，或其他介质密度与位置有关的问题，在求(i=1，2，…，Imax)时，如果系统形状复杂，计算是非常烦杂的。在这种情况下，使用最大截面法更方便。最大截面抽样方法为：其中iL),(max)(max,EEttrr1max,12max,1111)(),()(ln0LLEELELLLmtmmmtmtΩr≤＞(3)限制抽样法当介质区域很小时，如使用直接抽样法抽取输运长度，粒子很容易穿出介质，此时使用限制抽样法确定自由程个数ρ较好，ρ的分布密度函数为：其中Dm为粒子由rm出发，沿Ωm方向到达区域边界的自由程个数。ρ的抽样方法是：然后用直接抽样法中根据ρ计算L的方法计算输运长度L。此时，粒子的权重需乘以纠偏因子。其它当001)(mDDeefm])1(1ln[mDe)1(mDe2)碰撞后能量Em+1的随机抽样粒子在介质中发生碰撞后，首先要确定与哪种原子核发生何种反应。粒子发生碰撞后（吸收除外）的能量Em+1一般只与其碰撞前后运动方向的夹角（散射角）有关。粒子碰撞后常见的能量分布有下面几种情况。(1)裂变中子谱中子引起原子核裂变反应时，裂变中子的能量服从裂变谱分布。其抽样方法可参考以前的介绍。(2)中子弹性散射后能量的确定中子弹性散射后，能量与质心系散射角θC的关系是：能量与实验室系散射角θL的关系是：其中，A为碰撞核的质量，。或确定后，即可求出Em+1。222211)1(LLmmAAEE12)1(221CmmAAAEELLCCcos,cosLC(3)中子非弹性散射后能量的确定中子非弹性散射后，能量与质心系散射角θC的关系是：其中，为第K个能级的阈能，为第K个能级的激发态能量。如果确定了实验室系散射角θL，则根据下式确定后，再计算出Em+1。KKCmKmKmmAAEAEAAEE11121)1(221KKC11111222LLmKLmKCEAEA(4)光子康普顿（Compton）散射后能量的确定光子发生康普顿散射后，其能量分布密度函数为：其中，K(α)为归一因子。，和分别为光子散射前后的能量，以m0c2为单位，m0为电子静止质量，c为光速。211,1111)(1)(322xxxxxxKxf22)21(21421)21ln()1(21)(Kx光子康普顿散射能量分布的抽样方法为：x的抽样确定后，散射后的能量为：232223221121321)1(427212942711212121xxxxxxxxx＞＞＞≤≤≤xEcmxcmEmm202013)碰撞后散射角的随机抽样粒子碰撞后运动方向Ωm+1的确定，一般与散射角有关。由已知分布抽样确定散射角后，再确定Ωm+1。常见的散射角分布有如下几种：(1)质心系各向同性分布散射角在质心系服从各向同性分布时，其抽样方法为。质心系散射角θC抽样确定后，需转换成实验室系散射角θL。12C在中子弹性散射情况下，转换公式为：其中A为碰