3. 7 粒子输运问题的蒙特卡洛模拟 一、直接模拟法 直接模拟法是基于粒子输运过程的随机统计特性的考虑，认 为物理上的可观测量就是大量粒子的行为共同贡献的统计结果。 因此，该方法就是考虑一个一个粒子的传输，模拟它们在物质中 随机运动的历史，记录其在运动中对感兴趣的物理模拟量的贡 献。在对单个粒子运动历史进行大量的重复模拟之后，就可以对 物理模拟量进行统计平均，得到所需要的物理结果。 中子与物质作用后，一部分会被吸收，另一部分经过多次散 射后会穿透物质层透射出去。中子与物质层作用后可能会产生次 级粒子，我们不考虑这些次级粒子的迁移。如果中子和物质材料 中第m种原子核作用的全截面为 () ( ) ( )EEE m a m s m t σσσ += ()E m s σ和分别表示中子被一个原子核的散射和吸收截面。 如果单位体积第种原子核的数量记为 ()E m a σ m m m ρ，则中子作用在单位体 积内第m种元素上的总截面为 (cos M s G 0= i E , ( )(EE m tm m t σρ=Σ ). 假如材料中有多种元素，该中子与材料作用的总截面为: () () ∑ Σ=Σ m m tt EE . 假定中子与某一原子核散射后的角分布表示为( ) ?dE /dσ，当散射 角分布对方位角?是各向同性时，方位角?可以被积掉，得到微 分散射截面() ( )θσ cos/ dEd。无论微分截面d ()E / ?dσ或者 () )θσ / dEd，我们都可以得到相应的理论公式。 设在O点有一个能量为的中子垂直入射到物质层中。我们 记录这时该中子的状态位形为 0 E ( )1cos,,0 0000 === θExs G ，经过第一次 碰撞后散射到状态位形( ) 11 cos, 11 , θExs = G 2 ，再经过第二次碰撞后散射 到状态位形( 222 cos,, )θExs = G ，…，如此我们依次记下该中子在物 质层中运动历史上的位形点的轨迹： M ssss GGGG →???→→ 210 . 直到在状态，该中子被物质层吸收，透射或背射出来，或者在 处该中子的能量低于某一阈值，则程序就停止跟踪。 M s G M E 现在我们讨论程序如何具体模拟跟踪这个运动历程。 初始状态位形如前所述已经给出，假定为 ()1cos,, 0000 == θExs G 。 现在要由位形( ) 1111 cos,, ???? = iiii Exs θ G 确定下一个状态位形 () iii xs θcos,= G 。 我们采用如下步骤来确定状态 i s G 的各个参数： （1） 首先确定坐标参数。中子到达 i x i s G 状态点以前，经历过第 次碰撞后做匀速直线运动，其运动的自由程满足分布密度 函数 1?i y () ( )[ 11 exp)( ? Σ??Σ= tit yEyf y ] ?i E。我们可以采用直接抽样法得到自 由程的抽样值， () ξln 1 1? Σ ?= it E y . 则由下式给出 i x () 1 1 111 cos ln cos ? ? ??? Σ ?=+= i it iiii E xyxx θ ξ θ （2） 确定碰撞的原子核种类。中子与物质层中第m原子核碰撞 的几率为 ()( ) 11 )1( ?? ? ΣΣ= iti m t i m EEp . 则可以由离散型分布随机变量的直接抽样法，很容易确定发生碰 撞的是何种原子核。 （3） 确定碰撞的性质是吸收还是散射。中子与第原子核发 生散射的几率为 m ()( ) 11 )1( , ?? ? = i m ti m s i sm EEp σσ . 同样可以采用离散型随机变量的直接法抽取。若抽样结果为吸 收，则停止跟踪回到状态，开始下一个中子的跟踪；若抽样结 果为散射，则进入第（4）步。 0 s G （4） 确定中子散射角 i θ和能量。由于理论上一般给出的是质 心系中的散射微分截面公式 i E ( ) 1 cos ?i dE 1?i d θσ，因此我们需要首先 按照质心系的微分截面抽取散射角余弦 c θcos， c θcos满足的分布 密度函数为 () ( ) ( ) ∫ ? ?? = 1 1 11 cos coscos cos c c i c i c d d Ed d Ed f θ θ σ θ σ θ . 理论上，散射后的中子能量由下式计算得到 i E ()()[] cii rrEE θcos11 2 1 1 ?++= ? . 其中 2 1 1 ? ? ? ? ? ? + ? = A A r，A是原子核质量与中子质量之比。质心系散射角 c θ可以用下面公式换算为对应的实验室系的散射角 L θ。 () ccL AAA θθθ cos21cos1cos 2 +++= . 再根据下面的球面三角公式，通过实验室系散射角 L θ来确定 i θ。 ?θθθθθ cossinsincoscoscos 11 LiLii ?? += 其中?为方位角。由于我们考虑的中子散射过程是各向同性的， 方位角?是通过抽样πξ? 2=确定抽样值。 按照上面的计算步骤，我们就完成了从 1?i s G 到 i s G 状态的跟踪。 重复上述中子跟踪计算过程，直到中子在物质层中运动历程的终 点。如果我们一共模拟了个中子的运动过程，接下来，我们就 要对感兴趣的模拟物理量进行统计平均。 N 以计算投射率来计算可观测的物理量： 我们定义在模拟过程中第个中子对透射率的贡献n n η为 ? ? ? ≤ ≥ = 或被吸收,0 1 ax ax M M n η 其中下标M为该中子在物质层中碰撞的次数。我们得到穿透物质 层的中子数为 1 N . ∑ = = N n n N 1 1 η 由此得到透射率的一个估计值为 ∑ = == N n n NN N P 1 1 1 η 在α?1置信水平下，P的误差估计为 NtPP / ηα σ<? . η σ是 n η的均方差。由于 n η是一个二项式分布的随机变量，所以 ()( )PPPP ?≈?= 11 2 η σ 通过对大量中子运动过程的跟踪，我们也很容易求出透射中子的 能量和角分布。只要将能量E和极角θ分成若干个小区间，如： min210 EEEE >???>>>， 2/0 210 πθθθθ =<???<<<= M 将透射中子的能量E和极角θ记入图中对应区间，统计落入各个 能量区间或角度区间的中子数，并画出直方图。这样我们就得到 相应的散射中子的能量分布或角分布图形。 直接模拟法的优点：模拟过程重复了物理过程的机制，模拟 思想朴素、简单。 直接模拟法的缺点：当物质层较厚时，透射率会很小，导致 误差较大。 二、权重法 实际上，在模拟大量粒子穿过较厚物质层时，只有很少数的 粒子对透射率有贡献，这时透射率估计值的误差涨落也比较大。 为了克服这个缺点，可以采用对散射过程加权重的办法。 方法：记录到达某一状态的中子一定以一个概率权重被散 射，不再判断中子是否会被吸收。当中子离开物质或能量小于某 一阈值时，就停止它的运动历程的跟踪。将散射的概率权重加 到状态的位形参数中，此时中子的状态描写为 i w i w ( )wExs ,cos,, θ= G . 其中即是散射权重因子。类似在直接模拟法中，在跟踪一个中 子前，先给出它的初始状态位形 w ( ) 00000 ,cos,, wExs θ= G ，并取w。 在位形 1 0 = ( ) 11111 ,cos,, ????? = iiiii wExs θ G () ii w,cos 确定后，下一个状态位形 ii Es , i x , θ= G 中的参数 iii Ex θcos,, i w 的确定方法与前面描述的直 接模拟法相同。由下面公式确定 ( ) () 1 1 1 ? ? ? Σ Σ ?= i m t i m s ii E E ww . 上式中上标仍然指的是第种原子核。这时第个中子对透射 率的贡献为 m m n ? ? ? > = ? 其它0 1 axw i n δ . 假定我们一共跟踪了个中子，则透射率N P′的估计值为 ∑ = =′ N n n N P 1 1 δ . 它的方差为 () ∑ = ′?= N n n P N 1 2 22 1 δσ δ . 两种方法的方差的差别 () ∑ = ?≈? N n nn N 1 222 1 δδσσ δη . 由于1≤ n δ，所以存在不等式 . 22 δη σσ > 这个不等式说明权重法的方差小于直接模拟法的方差。 三、统计估计法 第个被跟踪的中子在第i个状态由位形n ( ) iiiii wExs ,cos,, θ= G 描 述。处在这一状态的中子直接穿透物质层的几率显然为 () ? ? ? ? ? > ? ? ? ? ? ? ? Σ? = 其它0 0cos, cos exp i i i iti i n xa Ew P θ θ . 在这个中子的运动历程中，每个碰撞点都有可能以几率 i n P透射出 去。我们可以充分利用这一信息，利用下式计算这个中子对透射 率的贡献： ∑ ? = = 1 0 M i i nn PP 式中M是被跟踪第中子在物质层中的碰撞点数。如果我们一共 跟踪了个中子，则最后得到透射率的估计值为这个中子贡献 的叠加。 n N N ∑∑∑ = ? == =≈′′ N n M i i n N n n P N P N P 1 1 01 11 . 它的方差为 ()() ∑ = ′′?≈ N n nw PP N 1 2 22 1 σ . 这种计算透射率的方法就叫统计估计法。 除了上面介绍的直接模拟法和在此基础上发展起来的权重 法和统计估计法外，还有其它许多发展出来的模拟方法，如：碰 撞点积分法、半解析方法等模拟方法。这些方法发展的初衷就是 要有效地降低模拟计算的方差，节约计算时间。