3.7粒子輸运问题的蒙特卡洛模拟 棋拟法 真换棋拟油是基于粒子贛运过程的隨机鴕计特性的考建,认 为物理上的可巩测量就是大量粒子的行为失同货的统计结。 因此,该方法就是考一个一个粒子的传翁,棋拟它们在物质中 随机运动的历史。记录其在运动中对感兴越的物理棋拟量的货 就。在对单个粒子运动历史进行大量的量复模拟之后,就可以对 物理模拟量进行计平均,;到所曇的物狸結果。 中子与物作用后,一部分会被吸收。另一部分过多次散 射后会穿還物孜层射出去。中子与物孜层作用后可能会产生次 級粒子,仉不考虑这些次级粒子的迁移。如是中子和物质材料 中第m种原子核作用的全微面为 "(E)=o"(E)+aa(E) dm(E)和am(E)分别衰示中子被一个原子核m的散射和吸收微面。 如果单位体积第m种原子核的数量记为pn,则中子作用在单位体
3. 7 粒子输运问题的蒙特卡洛模拟 一、直接模拟法 直接模拟法是基于粒子输运过程的随机统计特性的考虑，认 为物理上的可观测量就是大量粒子的行为共同贡献的统计结果。 因此，该方法就是考虑一个一个粒子的传输，模拟它们在物质中 随机运动的历史，记录其在运动中对感兴趣的物理模拟量的贡 献。在对单个粒子运动历史进行大量的重复模拟之后，就可以对 物理模拟量进行统计平均，得到所需要的物理结果。 中子与物质作用后，一部分会被吸收，另一部分经过多次散 射后会穿透物质层透射出去。中子与物质层作用后可能会产生次 级粒子，我们不考虑这些次级粒子的迁移。如果中子和物质材料 中第m 种原子核作用的全截面为 ( ) E (E) (E) m a m s m σ t = σ +σ (E) m σ s 和 分别表示中子被一个原子核 的散射和吸收截面。 如果单位体积第 种原子核的数量记为 (E) m σ a m m ρ m，则中子作用在单位体
积内第m种元上的总氰面为 ∑"(E)=pno"(E) 假如村料中有多种元景,该中子与村料作用的总微面为 ∑,(E)=∑Em 假定中子与某一原子被散射后的分布表示为dG(E)4,当射 角分布对方位角是各向同性时,方角口可以被积萍,脊到微 分散射薇面do(E)/d(os0)。无论微分戴窗do(E)或者 dσ(E)d(os),我们部可以得到相蠃的理论公式。 设在O点有一个能量为E的中子兽直入射到物质层中。我仉 记录这时该中子的状世形为50=(x0=0, Eo. cos e=1),经过第一次 碰撞后散射到状庵世形S1=(x1,E1,cos1),再经过第二次磁撞后射 到状形s2=(x2,E2,cosB2),…,如此我们依次记下该中子在物 败层中动历史上的形点的孰遼 S0→)S1→S2…→SM 直到在米庵。该中子被物质屈吸收。射背射出来。或者在 s处该中子的能量E氐于某一,则程序就停止跟腺。 親在我仉讨论程序如何具体模拟暇麝这个运动历翟 物始状忞位形如前所逃已給出,假定为 s0=(x0=0,E0,cos0=1) 现在要由位形31=(x-1,E1osO-)确突下一个状庵位形 s,=(x,E;,cos)。 我们采用如下步骤来确定状乏的各个參数:
积内第m 种元素上的总截面为 (cos Ms G = 0 Ei , (E) (E m m t m Σt = ρ σ ). 假如材料中有多种元素，该中子与材料作用的总截面为: Σ ( ) = ∑Σ ( ) m m t E t E . 假定中子与某一原子核散射后的角分布表示为dσ (E)/ dΩ，当散射 角分布对方位角ϕ 是各向同性时，方位角ϕ 可以被积掉，得到微 分散射截面 dσ ( ) E / d(cosθ ) 。无论微分截面 dσ ( ) E / dΩ 或 者 dσ ( ) E / d θ )，我们都可以得到相应的理论公式。 设在O点有一个能量为 的中子垂直入射到物质层中。我们 记录这时该中子的状态位形为 E0 ( 0, ,cos 1) s0 = x0 = E0 θ 0 = G ，经过第一次 碰撞后散射到状态位形 ( ) 1 1 ,cos 1 1 s = x , E θ G 2 ，再经过第二次碰撞后散射 到状态位形 ( 2 2 2 s = x , E ,cosθ ) G ，…，如此我们依次记下该中子在物 质层中运动历史上的位形点的轨迹： M s s s s G G G G 0 → 1 → 2 ⋅⋅⋅ → . 直到在 状态，该中子被物质层吸收，透射或背射出来，或者在 sM 处该中子的能量 低于某一阈值，则程序就停止跟踪。 G EM 现在我们讨论程序如何具体模拟跟踪这个运动历程。 初始状态位形如前所述已经给出，假定为 s0 = ( ) x0 , E0 ,cosθ 0 = 1 G 。 现在要由位形 ( ) 1 1 1 1 , ,cos i− = i− Ei− i− s x θ G 确定下一个状态位形 ( ) i i i s = x , cosθ G 。 我们采用如下步骤来确定状态 i s G 的各个参数：
(1)首先确定坐标参数x。中子到达状痕以前,经历过第 i-1次碰贛后敵訇遠直織运动,其运动的自由翟y演足分布度 函数f(y)=∑,(E-)exp,(E)。我们可以采用直接抽样法得到自 由程y的抽样值 则x由下式给出 x=x+ycos6-1=x压cose (2)确定碰撞的原子核种类。中子与物质层中m原子核碰 的几率为 pm=m(E)、(E) 则可以由高散型分布随机变量的直接抽禅法很容易确定发生磁 撞的是何种原子核。 (3)确定碰撞的性质是吸收还是散射。中子与第m原子核发 生觉射的几事为 pm)=om(E-)σ"(E) 同样可以采用高散型随机变量的直接法抽取。若抽绪是为吸 临,则停止跟隙回到状,,开始下一个中子的跟腺;若抽桦绪 最为散封,则进入第(4)乡。 (4)确定中子散射角日和能量E。由于理论上一般给出的是质 心系中的散射微分微面公式d(E)dcos-,因此我们要首先 换照质心系的微分截面抽取散射角汆弦c0s,cosb足的分布
（1） 首先确定坐标参数 xi。中子到达 i s G 状态点以前，经历过第 次碰撞后做匀速直线运动，其运动的自由程 满足分布密度 函数 i −1 y ( ) [ ( ) 1 1 ( ) = Σt i− ⋅ exp − Σt f y E y y ] Ei− 。我们可以采用直接抽样法得到自 由程 的抽样值， ( )lnξ 1 Σ −1 = − t Ei y . 则xi由下式给出 ( ) 1 1 1 1 1 cos ln cos − − − − − Σ = + = − i t i i i i i E x x y x θ ξ θ （2） 确定碰撞的原子核种类。中子与物质层中第m 原子核碰撞 的几率为 ( ) ( ) 1 1 ( 1) − − − = Σ i Σt i m t i pm E E . 则可以由离散型分布随机变量的直接抽样法，很容易确定发生碰 撞的是何种原子核。 （3） 确定碰撞的性质是吸收还是散射。中子与第 原子核发 生散射的几率为 m ( ) ( ) 1 1 ( 1) , − − − = i m i t m s i pm s σ E σ E . 同样可以采用离散型随机变量的直接法抽取。若抽样结果为吸 收，则停止跟踪回到 状态，开始下一个中子的跟踪；若抽样结 果为散射，则进入第（4）步。 0 s G （4） 确定中子散射角θ i和能量 。由于理论上一般给出的是质 心系中的散射微分截面公式 Ei ( ) 1 cos dσ Ei−1 d θ i− ，因此我们需要首先 按照质心系的微分截面抽取散射角余弦 θ c cos ， θ c cos 满足的分布
粤度函数为 d cos e d 0c/J-1 d cos 0 理论上。射后的中子能量E由下弌计扩得到 E[(1+r)+(1-r)cos] 2 其中r=1-,A是原子核质量与中子质量之比。质心系散射 A+1 0可以用下面公式换算为对应的实验宜系的散射角日2 cosL=(1+ Acose\n+A+2Acos8c 开根据下面的歌面三角公式,通过验蜜系散射角来确定日。 b 6. cos 6, + sin 6, sin 6, cos 其中φ为方位角。由守们考的中子散射过程是各向同性的, 方肩φ是通过抽样φ=2n确定抽祥。 按照上面的计算步暻,我们就完成了从S到§状虍的取瞭。 量复上述中宁眾计算过程,真则中宁在物质暴中运动历程的终 点。如果们一失棋拟了N个中子的动过程,接下来,我们就 要对感兴慜的棋拟物颯量选行统计平均。 以计算投射率來计犷可观测的物理量 我们定义在拟过程中邛n个中子对造射率的贡就n为 ≥ x≤a,或被吸收 其中下标M为该中子在物败层中碰的次数。我们得到穿造物质 层的中子数N为
密度函数为 ( ) ( ) ( ) ∫− − − = 1 1 1 1 cos cos cos cos c c i c i c d d d E d d E f θ θ σ θ σ θ . 理论上，散射后的中子能量Ei由下式计算得到 [ ] ( ) ( ) i i c E E 1 r 1 r cosθ 2 1 = −1 + + − . 其中 2 1 1 + − = A A r ，A是原子核质量与中子质量之比。质心系散射角 θ c可以用下面公式换算为对应的实验室系的散射角θ L。 ( ) θ L A θ c A A θ c cos 1 cos 1 2 cos 2 = + + + . 再根据下面的球面三角公式，通过实验室系散射角θ L来确定θ i。 cosθ i = cosθ i−1 cosθ L + sinθ i−1 sinθ L cosϕ 其中ϕ 为方位角。由于我们考虑的中子散射过程是各向同性的， 方位角ϕ 是通过抽样ϕ = 2πξ 确定抽样值。 按照上面的计算步骤，我们就完成了从 i−1 s G 到 i s G 状态的跟踪。 重复上述中子跟踪计算过程，直到中子在物质层中运动历程的终 点。如果我们一共模拟了 个中子的运动过程，接下来，我们就 要对感兴趣的模拟物理量进行统计平均。 N 以计算投射率来计算可观测的物理量： 我们定义在模拟过程中第n个中子对透射率的贡献η n为 ≤ ≥ = 0 ,或被吸收 1 x a x a M M η n 其中下标M 为该中子在物质层中碰撞的次数。我们得到穿透物质 层的中子数N1为
由此孙到射亭的一个佑计值为 P n 在1-a量信水平下,P的误鑾佔计为 E1>E2>…>Em, 0=b0<61<62<…<bn=/2 将适射中子的能量E和极扇0记入国中对应区间,计落入吝个 能量区间或肩度区间的中子数,并画出直方图。这桦觉们就到 相应的射中子的能量分布域角分布图形。 〖接拟法的优痕:模拟过程置了物理过程的机侧,模拟 思丸朴录、蒿阜。 直接模拟法的峡痕:当物赝层毓厚时,适射事会很小,导政 误大。 二、杈量法 实称上在模拟大量粒子穿过校尾物放层时,只有很少教的
∑ . = = N n N n 1 1 η 由此得到透射率的一个估计值为 ∑= = = N n n N N N P 1 1 1 η 在1−α 置信水平下，P 的误差估计为 P − P Emin > > > ⋅⋅⋅ ， 0 = θ 0 < θ1 < θ 2 < ⋅⋅⋅ < θ M = π / 2 将透射中子的能量E 和极角θ 记入图中对应区间，统计落入各个 能量区间或角度区间的中子数，并画出直方图。这样我们就得到 相应的散射中子的能量分布或角分布图形。 直接模拟法的优点：模拟过程重复了物理过程的机制，模拟 思想朴素、简单。 直接模拟法的缺点：当物质层较厚时，透射率会很小，导致 误差较大。 二、权重法 实际上，在模拟大量粒子穿过较厚物质层时，只有很少数的
粒子对射率有獻,这时造射率佑计值的误孌涨痦也比較大。 为了克服这个缺京,可以采用对散射过程加权置的办法。 方法:记录到达某一状恋的中子一定以一个概率权量ν驶放 射,不开判斷中子是否会被吸收。当中子高开物质或能量小于某 值时,就停止包的坛动历程的跟腺。将微射的概率权置w加 到状恋的位形参教中,此时中子的状描写为 s=(x, E, cos 0, w) 其中w即是微射权量因子。类似在直接模拟法中,在跟腺一个中 子前,先给出它的初始状夷形S。=(x0,E0, cos e,w0),并取0=1。 在世形31=(x1,E1 确定后,下一个米庵位形 s=(x,E,cos日,m,)中的參数x,E,cos的确定方法与前面迷的直 接模拟法相同。W由下面公式确定 上式中上m仍指的是第m种原子核。这时第n个中子对适射 率的货觥为 x>a 0其它 假定我们一共跟了N个中子,则灒射率P的佑计值为 ∑δn 它的方为 ∑2-(P)
粒子对透射率有贡献，这时透射率估计值的误差涨落也比较大。 为了克服这个缺点，可以采用对散射过程加权重的办法。 方法：记录到达某一状态的中子一定以一个概率权重 被散 射，不再判断中子是否会被吸收。当中子离开物质或能量小于某 一阈值时，就停止它的运动历程的跟踪。将散射的概率权重 加 到状态的位形参数中，此时中子的状态描写为 wi wi s = (x, E,cosθ ,w) G . 其中 即是散射权重因子。类似在直接模拟法中，在跟踪一个中 子前，先给出它的初始状态位形 w ( ) 0 0 0 0 0 s = x , E ,cosθ ,w G ，并取w 。 在位形 0 = 1 ( ) 1 1 1 1 1 , ,cos , i− = i− Ei− i− wi− s x θ G ( ) i wi cos , 确定后，下一个状态位形 i Ei s , i = x , θ G 中的参数 i Ei i x , ,cosθ wi 的确定方法与前面描述的直 接模拟法相同。 由下面公式确定 ( ) ( ) 1 1 1 − − − Σ Σ = ⋅ i m t i m s i i E E w w . 上式中上标 仍然指的是第 种原子核。这时第 个中子对透射 率的贡献为 m m n > = − 0 其它 wi 1 x a δ n . 假定我们一共跟踪了N 个中子，则透射率P′的估计值为 ∑= ′ = N n n N P 1 1 δ . 它的方差为 ∑ ( ) = = − ′ N n n P N 1 2 2 1 2 σ δ δ
两种方法的方的剔 O-O≈ 由于δ≤1,所以存在不式 这个不式说明杈权量法的方塾小于真接棋拟法的方鎏。 三、筑计佑计法 第n个被跟踪的中子在第i个状由世形3=(x1,E,cos日,w,)描 迷。处在逭一状忞的中子真接穿物质层的几率显为 P w,exp,(E cos0>0 其它 在这个中子的坛动历程中,每个碰京部有可能以几率P没射出 去。我们可以充分利用这一信丸,利用下式计扩这个中子对射 卒的货 =∑P 式中M是被跟赓邛n中子在物败层中的碰京数。如景爱们一类 取驗了N个中子,则后得到灒射率的计值为这N个中子觥 的兔加。 F=1Σ两=1∑8可 它的方接为
两种方法的方差的差别 ∑( ) = − ≈ − N n n n N 1 2 2 1 2 σ η σ δ δ δ . 由于δ n ≤ 1，所以存在不等式 . 2 2 σ η > σ δ 这个不等式说明权重法的方差小于直接模拟法的方差。 三、统计估计法 第 n 个被跟踪的中子在第i 个状态由位形 ( ) i i Ei i wi s = x , ,cosθ , G 描 述。处在这一状态的中子直接穿透物质层的几率显然为 ( ) > − − Σ = 0 其它 , cos 0 cos exp i i i i t i i n a x w E P θ θ . 在这个中子的运动历程中，每个碰撞点都有可能以几率 i Pn 透射出 去。我们可以充分利用这一信息，利用下式计算这个中子对透射 率的贡献： ∑ − = = 1 0 M i i Pn Pn 式中M 是被跟踪第 中子在物质层中的碰撞点数。如果我们一共 跟踪了 个中子，则最后得到透射率的估计值为这 个中子贡献 的叠加。 n N N ∑ ∑ ∑= − = = ′′ ≈ = N n M i i n N n n P N P N P 1 1 1 0 1 1 . 它的方差为
a2=∑(2)-(P) 这种计算射率的方法就叫统计佑讣法。 除了上面介绍的直接棋拟法和在此基础上晨態来的权量 法和计佔计法外,还有其它许多墩展出來的棋拟方法,如:碰 撞痕积分法、半解析方法檄拟方法。这些才法歇昃的物哀就是 要有效地降氐萇拟计犷的方,节约计算时间
∑( ) ( ) = ≈ − ′′ N n w Pn P N 1 2 2 1 2 σ . 这种计算透射率的方法就叫统计估计法。 除了上面介绍的直接模拟法和在此基础上发展起来的权重 法和统计估计法外，还有其它许多发展出来的模拟方法，如：碰 撞点积分法、半解析方法等模拟方法。这些方法发展的初衷就是 要有效地降低模拟计算的方差，节约计算时间