蒙特卡羅方法教學課件第四章蒙特卡羅方法解粒子輸運問題兩份文件

第四章蒙特卡羅方法解粒子輸運問題屏蔽問題模型直接模擬方法簡單加權(quán)法統(tǒng)計估計法指數(shù)變換法蒙特卡羅方法的效率作業(yè)第四章蒙特卡羅方法解輻射屏蔽問題 輻射（光子和中子）屏蔽問題是蒙特卡羅方法最早廣泛應(yīng)用的領(lǐng)域之一。本章主要從物理直觀出發(fā)，說明蒙特卡羅方法解決這類粒子輸運問題的基本方法和技巧。而這些方法和技巧對于諸如輻射傳播、多次散射和通量計算等一般粒子輸運問題都是適用的。屏蔽問題模型 在反應(yīng)堆工程和輻射的測量與應(yīng)用中，常常要用一些吸收材料做成屏蔽物擋住光子或中子。我們所關(guān)心的是經(jīng)過屏蔽后射線的強度及其能量分布，這就是屏蔽問題。 當屏蔽物的形狀復(fù)雜，散射各向異性，材料介質(zhì)不均勻,核反應(yīng)截面與能量、位置有關(guān)時，難以用數(shù)值方法求解，用蒙特卡羅方法能夠得到滿意的結(jié)果。 粒子的輸運問題帶有明顯的隨機性質(zhì)，粒子的輸運過程是一個隨機過程。粒子的運動規(guī)律是根據(jù)大量粒子的運動狀況總結(jié)出來的，是一種統(tǒng)計規(guī)律。蒙特卡羅模擬，實際上就是模擬相當數(shù)量的粒子在介質(zhì)中運動的狀況，使粒子運動的統(tǒng)計規(guī)律得以重現(xiàn)。不過，這種模擬不是用實驗方法，而是利用數(shù)值方法和技巧，即利用隨機數(shù)來實現(xiàn)的。 為方便起見，選用平板屏蔽模型，在厚度為a，長、寬無限的平板左側(cè)放置一個強度已知，具有已知能量、方向分布的輻射源S。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源發(fā)出的平均一個粒子穿透屏蔽的數(shù)目。 同時，假定粒子在兩次碰撞之間按直線運動,且粒子之間的相互作用可以忽略。直接模擬方法 直接模擬方法就是直接從物理問題出發(fā)，模擬粒子的真實物理過程。狀態(tài)參數(shù)與狀態(tài)序列模擬運動過程記錄結(jié)果 粒子在介質(zhì)中的運動的狀態(tài)，可用一組參數(shù)來描述，稱之為狀態(tài)參數(shù)。它通常包括：粒子的空間位置r，

能量E和運動方向Ω，以S＝(r,E,Ω)表示。 有時還需要其他的參數(shù)，如粒子的時間t和附帶的權(quán)重W，這時狀態(tài)參數(shù)為S'＝(r,E,Ω,t,W)。

狀態(tài)參數(shù)通常要根據(jù)所求問題的類型和所用的方法來確定。 對于無限平板幾何，取S＝(z,E,cosα) 其中z為粒子的位置坐標，α為粒子的運動方向與Z軸的夾角。 對于球?qū)ΨQ幾何,取S＝(r,E,cosθ) 其中r表示粒子所在位置到球心的距離，θ為粒子的運動方向與其所在位置的徑向夾角。狀態(tài)參數(shù)與狀態(tài)序列 粒子第m次碰撞后的狀態(tài)參數(shù)為 或 它表示一個由源發(fā)出的粒子，在介質(zhì)中經(jīng)過m次碰撞后的狀態(tài)，其中

rm：粒子在第m次碰撞點的位置

Em：粒子第m次碰撞后的能量

Ωm：粒子第m次碰撞后的運動方向

tm：粒子到第m次碰撞時所經(jīng)歷的時間

Wm：粒子第m次碰撞后的權(quán)重 有時，也可選為粒子進入第m次碰撞時的狀態(tài)參數(shù)。 一個由源發(fā)出的粒子在介質(zhì)中運動，經(jīng)過若干次碰撞后，直到其運動歷史結(jié)束（如逃出系統(tǒng)或被吸收等）。假定粒子在兩次碰撞之間按直線運動，其運動方向與能量均不改變，則粒子在介質(zhì)中的運動過程可用以下碰撞點的狀態(tài)序列描述：

S0，S1，…，SM-1，SM 或者更詳細些,用 來描述。這里S0為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，SM為粒子的終止狀態(tài)。M稱為粒子運動的鏈長。 這樣的序列稱為粒子隨機運動的歷史，模擬一個粒子的運動過程，就變成確定狀態(tài)序列的問題。 為簡單起見，這里以中子穿透均勻平板的模型來說明，這時狀態(tài)參數(shù)取S＝(z,E,cosα)。 模擬的步驟如下： (1)確定初始狀態(tài)

S0： 確定粒子的初始狀態(tài)，實際上就是要從中子源的空間位置、能量和方向分布中抽樣。設(shè)源分布為 則分別從各自的分布中抽樣確定初始狀態(tài)。 對于平板情況， 抽樣得到z0＝0。模擬運動過程(2)確定下一個碰撞點： 已知狀態(tài)Sm-1，要確定狀態(tài)Sm，首先要確定下一個碰撞點的位置zm。在相鄰兩次碰撞之間，中子的輸運長度l服從如下分布： 對于平板模型，l服從分布： 其中，Σt為介質(zhì)的中子宏觀總截面， 積分稱為粒子輸運的自由程數(shù)， 系統(tǒng)的大小通常就是用系統(tǒng)的自由程數(shù)表示的。

顯然，粒子輸運的自由程數(shù)服從指數(shù)分布， 因此從f(l)中抽樣確定l，就是要從積分方程 中解出l。 對于單一介質(zhì) 則下一個碰撞點的位置 如果zm≥a，則中子穿透屏蔽，若zm≤0,則中子被反射出屏蔽。這兩種情況，均視為中子歷史終止。(3)確定被碰撞的原子核： 通常介質(zhì)由幾種原子核組成，中子與核碰撞時，要確定與哪一種核碰撞。設(shè)介質(zhì)由A、B、C三種原子核組成，其核密度分別為NA、NB、NC，則介質(zhì)的宏觀總截面為： 其中分別為核A、B、C的宏觀總截面。其定義如下： 分別表示(·)核的宏觀總截面、核密度和微觀總截面。 由于中子截面表示中子與核碰撞可能性的大小，因此，很自然地，中子與A、B、C核發(fā)生碰撞的幾率分別為： 利用離散型隨機變量的抽樣方法，確定碰撞核種類：＞≤＞≤(4)確定碰撞類型： 確定了碰撞的核(比如B核)后，就要進一步確定碰撞類型。中子與核的反應(yīng)類型有彈性散射、非彈性散射、(n,2n)反應(yīng)，裂變和俘獲等，它們的微觀截面分別為 則有 各種反應(yīng)發(fā)生的幾率分別為 利用離散型隨機變量的抽樣方法，確定反應(yīng)類型。 在屏蔽問題中，中子與核反應(yīng)常只有彈性散射和吸收兩種類型，吸收截面為： 這時，總截面為： 發(fā)生彈性散射的幾率為： 若，則為彈性散射；否則為吸收，發(fā)生吸收反應(yīng)意味著中子的歷史終止。(5)確定碰撞后的能量與運動方向： 如果中子被碰撞核吸收，則其輸運歷史結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運動方向。中子能量Em為：

A是碰撞核的質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量；θC為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角?？蓮馁|(zhì)心系中彈性散射角分布

fC(μC)中抽樣產(chǎn)生。實驗室系散射角θL的余弦μL為： 如果給出實驗室系散射角余弦分布fL(μL)，可直接從fL(μL)中抽取μL，此時能量Em與μL的關(guān)系式為： 確定了實驗室系散射角θL后， 再使用球面三角公式 確定cosαm： 其中χ為在［0,2π］上均勻分布的方位角。 至此，由Sm-1完全可以確定Sm。因此，當中子由源出發(fā)后，即S0確定后，重復(fù)步驟(2)～(5)，直到中子游動歷史終止。于是得到了一個中子的隨機游動歷史S0，S1，…，SM-1，SM，即 也就是模擬了一個由源發(fā)出的中子的運動過程。 以上模擬過程可分為兩大步：第一步確定粒子的初始狀態(tài)S0，第二步由狀態(tài)Sm-1來確定狀態(tài)Sm。這第二步又分為兩個過程：第一個過程是確定碰撞點位置zm，稱為輸運過程；第二個過程是確定碰撞后粒子的能量及運動方向，稱為碰撞過程。對于中子而言，碰撞過程是先確定散射角，進而確定能量和運動方向；而對于光子，碰撞過程是先確定能量，再確定散射角以及運動方向。重復(fù)這兩個過程，直至粒子的歷史終止。 這種模擬過程，是解任何類型的粒子輸運問題所共有的，它是蒙特卡羅方法解題的基本手段。 在獲得中子的隨機游動歷史后，我們要對所要計算的物理量進行估計。對于屏蔽問題，我們要計算中子的穿透率?？疾烀總€中子的隨機游動歷史，它可能穿透屏蔽（zM≥a），可能被屏蔽發(fā)射回來（zM≤0），或者被吸收。設(shè)第n個中子對穿透的貢獻為ηn，則 如果我們共跟蹤了N個中子，則穿透屏蔽的中子數(shù)為：記錄結(jié)果 則穿透屏蔽概率的近似值為： 它是穿透率的一個無偏估計。 我們稱這種直觀地模擬過程和估計方法為直接模 擬方法。在置信水平1－α＝0.95時，的誤差為： 其中

為ηn的均方差，由于ηn是一個服從二項分布的隨機變量，所以 或 為得到中子穿透屏蔽的能量、角分布，將能量、角度范圍分成若干個間隔： 其中Emax，Emin分別表示能量的上、下限，對于穿透屏蔽的中子按其能量、方向分間隔記錄。設(shè)一穿透屏蔽的中子能量為EM，其運動方向與Z軸夾角為αM，若能量EM屬于第i個能量間隔ΔEi，角度αM屬于第j個角度間隔Δαj，則分別在第i個能量計數(shù)器及第j個角度計數(shù)器中加"1"。 跟蹤N個中子后，則 分別為穿透中子的能量分布和角分布。其中N1,i和N2,i分別為第i個能量和第j個角度間隔的穿透中子數(shù)。歸一后分別為：簡單加權(quán)法

從模擬物理過程來說，直接模擬法是最簡單、也是最基本的方法。但是，在直接模擬法中，不管中子在屏蔽中經(jīng)過多少次碰撞，只要在介質(zhì)中被吸收，對穿透的貢獻就為零；因此在所跟蹤的粒子中絕大部分都對穿透沒有貢獻。而在許多屏蔽問題中，穿透率的數(shù)量級在10-6到10-8。進一步，如果我們要求穿透率達相對誤差小于1％，即 那么，N要大到驚人的數(shù)量級1010到1012。顯然，這時用直接模擬法計算不是很有效。 屏蔽物一般是由吸收強的介質(zhì)組成，因此在每次碰撞時，粒子很有可能被吸收而停止跟蹤。現(xiàn)在改變 模擬方法，在判斷碰撞類型時，可以認為粒子 的部分是彈性散射，而其余部 分被吸收，即人為地把中子分成兩部分，一部分彈性散射，一部分吸收。彈性散射這部分繼續(xù)跟蹤；吸收部分則停止跟蹤。也就是說，我們利用中子權(quán)重的變化來反應(yīng)繼續(xù)彈性散射的部分。這就是簡單加權(quán)法的基本思想。簡單加權(quán)法 顯然，在加權(quán)法中中子的權(quán)重W已成為中子狀態(tài)參數(shù)的組成部分。這時，中子歷史成為： 對源中子，取W0=1。經(jīng)過碰撞中子權(quán)重的變化為： 因子稱為尚存因子。 這時，第n個中子對穿透的貢獻為： 如果我們共跟蹤了N個中子，則穿透率P的無偏估計為： 類似地，可以得到穿透中子的能量分布和角分布。只不過在對各計數(shù)器進行的加"1"操作改為加WM。 簡單加權(quán)法的方差估計為： 與直接模擬法相比，有 注意到ζn≤1，有 這表明簡單加權(quán)法的方差小于直接模擬法的方差。這是因為加權(quán)法比直接模擬法減少了一次隨機抽樣。簡單加權(quán)法的方差 加權(quán)法的思想在蒙特卡羅方法中用途很廣泛。例如，對于具有中子增殖反應(yīng)，如裂變，(n,2n)，(n,3n)反應(yīng)的中子輸運問題，一個中子與核發(fā)生碰撞后，根據(jù)反應(yīng)的類型會產(chǎn)生不同數(shù)量的次級中子，每個次級中子又會產(chǎn)生新的次級中子，這樣鏈鎖反應(yīng)下去，使得用直接模擬法模擬每一個中子是非常困難的。這種情況可以利用加權(quán)法來處理。權(quán)重方法的其它應(yīng)用 中子與核發(fā)生碰撞后，產(chǎn)生的次級中子平均數(shù)為：這里νf為裂變次級中子數(shù)。于是，碰撞后的權(quán)重為： 而決定碰撞類型的幾率分別為： 其中 加權(quán)法的思想，還可以應(yīng)用到連續(xù)分布情況和偏倚抽樣的問題統(tǒng)計估計法

加權(quán)法雖然改進了直接模擬法，但它同樣只關(guān)心中子是否穿透屏蔽這一信息，因此對每一個中子歷史的信息利用得很不充分。統(tǒng)計估計法能夠較多地利用中子的歷史信息，因而能得到更好的結(jié)果。 一個中子，可能在介質(zhì)內(nèi)不發(fā)生碰撞而直接穿透屏蔽，也可能在介質(zhì)內(nèi)發(fā)生一次碰撞后再穿透屏蔽，或經(jīng)過二次碰撞穿 透屏蔽，等等， 這些事件是互不相 容的，因此穿透概 率P可表示為： 其中Pm是中子恰好經(jīng)過m次碰撞而穿透屏蔽的概率。這表明，可以用求Pm(m=0,1,…)的方法得到P。這樣，中子對穿透概率的貢獻就不只限于末次碰撞了。α0α1S0S1SmαmP1P0PmZa0 設(shè)中子的歷史為： 根據(jù)該中子的歷史，我們可以估計出中子恰好經(jīng) 過m次碰撞后，穿透屏蔽的部分 顯然，具有初態(tài)S0＝(0,E0,cosα0,W0)的中子，未經(jīng)碰撞直接穿透的部分是： 類似地，在經(jīng)過了第m次碰撞后的中子具有狀態(tài)Sm＝(zm,Em,cosαm,Wm)，其可能穿透的部分，正好是一個中子恰好經(jīng)過m次碰撞穿透的部分： 這里的這種估計技巧，由于是對每次碰撞后的狀態(tài)，求其后未經(jīng)碰撞直接穿透的貢獻，因此該方法也稱為最后自由飛行估計。 于是得到該中子對穿透的貢獻： 如果我們共跟蹤了N個中子，則穿透率P的估計為： 其方差估計為： 在直接模擬方法中，相對誤差為 其中為與置信水平1－α相應(yīng)的量。 如果構(gòu)造一個新的概率模型，使得該模型的穿透率P*與原模型的穿透率P之間存在關(guān)系： 使用直接模擬方法,相對誤差為指數(shù)變換法 如果令ε*＝ε，即 這意味著，達到同樣的相對誤差，跟蹤粒子的數(shù)目縮小K倍，從而減少K倍的計算量。指數(shù)變換法就是構(gòu)造一個新的概率模型的一個有效方法。 構(gòu)造如下偽過程：宏觀總截面為 散射截面仍為Σel(E)。其中Emin、Emax分別為能量的下限和上限，α為粒子的運動方向與Z軸的夾角?？梢宰C明這個偽過程的穿透概率P*與原過程的穿透概率P之間有如下關(guān)系: 顯然，。因偽過程與原過程的結(jié)果相差e指數(shù)，所以該方法稱為指數(shù)變換法。 分析一下偽過程的定義,可以明顯看出P*增大的原因。當cosα=1時，粒子運動方向與Z軸方向一致，其截面最小，粒子沿Z軸方向輸運的距離較遠；而當cosα=－1時，粒子運動方向背向Z軸方向，這時其截面最大，粒子向后輸運的距離較短。因此，截面變換的結(jié)果是加強了粒子向前運動的能力，因而使穿透概率增大。 偽過程的構(gòu)造與幾何形狀及所考慮的問題有關(guān)。比如，對球形幾何，使用指數(shù)變換法求穿透概率時,所構(gòu)造的宏觀總截面與平板屏蔽的情況不同，粒子的模擬方法也較復(fù)雜。蒙特卡羅方法的效率 衡量一種蒙特卡羅技巧的好壞，除了看其方差大小外，還要看其所需費用（計算時間）多少，即從該技巧的效率Ef（方差與費用乘積的倒數(shù)）全面考慮： 其中σ2為方差，T為所需費用。Ef大時，所用方法的效率高；否則，效率低。 在一般情況下，有些方法雖然減小了方差，卻增加了費用。例如，加權(quán)法、統(tǒng)計估計法雖然較直接模擬方法減小了方差，卻使每個粒子的運動鏈長增加，或記錄貢獻的計算時間增加。因此，不能認為方差小的方法一定好，要從方法的效率全面考慮。在有些情況下，直接模擬方法仍然是一個被廣泛使用的方法。作業(yè)

光子散射后能量分布的抽樣 把光子散射能量分布改寫成如下形式進行抽樣： 在［1,1+2α］上定義如下函數(shù):作業(yè)

給出分布密度函數(shù) 的抽樣方法。＞≤第五章蒙特卡羅方法在計算機上的實現(xiàn)源分布抽樣過程空間、能量和運動方向的隨機游動過程記錄貢獻和分析結(jié)果過程核截面數(shù)據(jù)的引用蒙特卡羅程序結(jié)構(gòu)作業(yè)第五章蒙特卡羅方法在計算機上的實現(xiàn)

蒙特卡羅方法是隨著計算機的出現(xiàn)和發(fā)展而逐步發(fā)展起來的。在計算機上能夠產(chǎn)生符合要求的隨機數(shù)，實現(xiàn)對已知分布的抽樣，奠定了蒙特卡羅方法在計算機上得以實現(xiàn)的基礎(chǔ)。在計算機上使用蒙特卡羅方法解粒子輸運問題大致包括三個過程：源分布抽樣過程，空間、能量和運動方向的隨機游動過程以及記錄、分析結(jié)果過程。源分布抽樣過程

源分布抽樣的目的是產(chǎn)生粒子的初始狀態(tài) 。下面我們介紹一些常見的特定 類型的源分布抽樣方法。源粒子的位置常見分布的隨機抽樣圓內(nèi)均勻分布 設(shè)圓半徑為R0，粒子在圓內(nèi)均勻分布時，從發(fā)射點到中心的距離r

的分布密度函數(shù)為：

r

的抽樣方法為：圓環(huán)內(nèi)均勻分布 設(shè)圓環(huán)的內(nèi)半徑為R0，外半徑為R1，則粒子在該圓環(huán)內(nèi)均勻分布時，從發(fā)射點到中心的距離r

的分布密度函數(shù)為：

r

的抽樣方法為：≤>球內(nèi)均勻分布 設(shè)球的半徑為R，粒子在球內(nèi)均勻分布時，從發(fā)射點到中心的距離r

的分布密度函數(shù)為：

r

的抽樣方法為： 在直角坐標系下，抽樣方法為：≤>球殼內(nèi)均勻分布 設(shè)球殼的內(nèi)半徑為R0，外半徑為R1，在均勻分布時，從發(fā)射點到中心的距離r

的分布密度函數(shù)為：

r

的抽樣方法為：≤≤>>

在直角坐標系下，球殼內(nèi)點的坐標為： 其中，r

由前面的抽樣方法確定，θ、φ服從各向同性分布，其抽樣方法為：＞≤圓柱內(nèi)均勻分布 圓柱內(nèi)均勻分布是指粒子發(fā)射點均勻地分布在底半徑為R，高為2H

的圓柱內(nèi)。若固定圓柱的中心為原點，圓柱的軸向為z

軸，則分布密度函數(shù)為： 抽樣方法為：≤>點源分布 點源分布是指粒子由一固定點發(fā)射，其分布密度函數(shù)為： 其中，為狄拉克δ函數(shù)，源粒子的抽樣方法為： 在球坐標系中，粒子發(fā)射點到球心的距離r

的分布密度函數(shù)為： 其中，為點源到球心的距離。源粒子的位置抽樣為：球外平行束源分布 球外平行束源分布是指粒子平行入射到半徑為R

的球面上，或球外點源距離球很遠，可以近似地看作平行束源。設(shè)r

為粒子發(fā)射點到球心的距離,其分布密度函數(shù)為：

r

的抽樣方法為： 在直角坐標系中，抽樣方法為：≤>源粒子的能量常見分布的隨機抽樣單能源分布 單能源分布是指粒子的發(fā)射能量為一固定值E0

，其分布密度函數(shù)為： 源粒子的能量為：裂變中子譜分布 裂變中子譜分布的一般形式為： 其中A，B，C，Emin，Emax

均為與元素有關(guān)的量。 對于鈾-235，

A=0.965，B=2.29，C=0.453，Emin=0，Emax=∞。

采用近似修正抽樣，抽樣方法為：

其中，m≈0.8746，M1≈0.2678，λ≈0.5543。

此外，裂變譜分布也有以數(shù)值曲線形式給出的，此時，用數(shù)值曲線抽樣方法抽取E

?！堋埽荆钧溈怂鬼f(Maxwell)譜分布 麥克斯韋譜分布的一般形式為： 該分布的抽樣方法為＞≤源粒子運動方向常見分布的隨機抽樣各向同性分布 各向同性分布密度函數(shù)為： 其中，μ＝cosθ，θ為運動方向與z

軸的夾角，φ為方位角。

在直角坐標系下，各方向余弦u，v，w

為： 其抽樣方法為：＞≤半面各向同性分布 不妨設(shè)在x≥0的半面方向上各向同性發(fā)射粒子，則在前述各向同性分布的抽樣方法中，用ξ2代替η2就能得到所需分布的抽樣。對于其它方向的情況，可用類似的方法處理。球外平行束源分布 令μ＝cosθ，θ為粒子運動方向的徑向夾角，則μ分布密度函數(shù)為：

μ的抽樣方法為：球外各向同性點源分布 設(shè)球外點源S

到球心的距離為D0。點源S

到球的最大張角為θ*， 則球外各向同性點源分布的抽樣方法是： 先抽樣確定，再轉(zhuǎn)換成θ。

在直角坐標系下，取OS

為z

軸，抽樣方法為：次級粒子的源分布

在有關(guān)次級粒子（如裂變中子，中子生成光子，光子生成中子）的輸運過程中，次級粒子源分布的抽樣方法，主要可分為以下兩種：直接生成法 可將生成的次級粒子的位置、能量、方向、權(quán)重等參數(shù)直接作為源分布的抽樣結(jié)果。也就是直接對生成的次級粒子進行跟蹤。這種方法比較簡單、直觀。離散分布法 將生成的次級粒子的權(quán)重，按空間位置、能量、方向分別記錄，得到次級粒子的空間、能量、運動方向的離散的近似分布。再根據(jù)該分布，利用各種抽樣技巧，得到源分布的抽樣，對抽樣的源粒子進行跟蹤、記錄。 當一個問題需要用兩個以上的蒙卡程序處理時，可采用這種方法。空間、能量和運動方向的隨機游動過程

粒子由狀態(tài)Sm到狀態(tài)Sm+1時，需要確定粒子的空間位置rm+1，能量Em+1和運動方向Ωm+1。碰撞點位置的計算公式

設(shè)rm

為粒子第m

次碰撞點的位置，Ωm

為碰撞后的運動方向，則粒子第m+1次碰撞點的位置rm+1

為： 即 其中為的方向余弦，L

為兩次碰撞點間的距離。

L

的分布密度函數(shù)為： 由f(L)抽樣確定L

的方法通常有三種：

直接抽樣方法 確定L

的直接抽樣方法是： 首先由自由程分布 中抽取ρ

再由下列關(guān)系式解出L。

對于均勻介質(zhì)，有 對于多層介質(zhì)，如果 則 其中，為粒子由rm

出發(fā)，沿Ωm

方向在順序經(jīng)過的第i

個介質(zhì)區(qū)域內(nèi)走過的距離，為第i

個介質(zhì) 區(qū)域的宏觀總截面(i=1，2，…，Imax)。當 時，意味著粒子穿出系統(tǒng)。最大截面法 對于多層介質(zhì)，或其他介質(zhì)密度與位置有關(guān)的問題，在求(i=1，2，…，Imax)時，如果系統(tǒng)形狀復(fù)雜，計算是非常煩雜的。在這種情況下，使用最大截面法更方便。最大截面抽樣方法為： 其中≤＞限制抽樣法 當介質(zhì)區(qū)域很小時，如使用直接抽樣法抽取輸運長度，粒子很容易穿出介質(zhì)，此時使用限制抽樣法確定自由程個數(shù)ρ較好，ρ的分布密度函數(shù)為： 其中Dm

為粒子由rm

出發(fā)，沿Ωm

方向到達區(qū)域邊界的自由程個數(shù)。ρ的抽樣方法是： 然后用直接抽樣法中根據(jù)ρ計算L的方法計算輸運長度L

。此時，粒子的權(quán)重需乘以糾偏因子。碰撞后能量Em+1的隨機抽樣

粒子在介質(zhì)中發(fā)生碰撞后，首先要確定與哪種原子核發(fā)生何種反應(yīng)。粒子發(fā)生碰撞后（吸收除外）的能量Em+1

一般只與其碰撞前后運動方向的夾角（散射角）有關(guān)。 粒子碰撞后常見的能量分布有下面幾種情況。

裂變中子譜 中子引起原子核裂變反應(yīng)時，裂變中子的能量服從裂變譜分布。其抽樣方法可參考以前的介紹。中子彈性散射后能量的確定 中子彈性散射后，能量與質(zhì)心系散射角θC的關(guān)系是： 能量與實驗室系散射角θL的關(guān)系是： 其中，A

為碰撞核的質(zhì)量，。 或確定后，即可求出Em+1。中子非彈性散射后能量的確定 中子非彈性散射后，能量與質(zhì)心系散射角θC的關(guān)系是： 其中，為第K

個能級的閾能，為第K

個能級的激發(fā)態(tài)能量。 如果確定了實驗室系散射角θL，則根據(jù)下式 確定后，再計算出Em+1。光子康普頓（Compton）散射后能量的確定 光子發(fā)生康普頓散射后，其能量分布密度函數(shù)為： 其中，K(α)為歸一因子。 ，和分別為光子散射前后的能量，以m0c2

為單位，m0為電子靜止質(zhì)量，c

為光速。

光子康普頓散射能量分布的抽樣方法為：

x

的抽樣確定后，散射后的能量為：＞＞＞≤≤≤碰撞后散射角的隨機抽樣

粒子碰撞后運動方向Ωm+1的確定，一般與散射角有關(guān)。由已知分布抽樣確定散射角后，再確定Ωm+1。常見的散射角分布有如下幾種：

質(zhì)心系各向同性分布 散射角在質(zhì)心系服從各向同性分布時，其抽樣方 法為。質(zhì)心系散射角θC抽樣確定后， 需轉(zhuǎn)換成實驗室系散射角θL。

在中子彈性散射情況下，轉(zhuǎn)換公式為： 其中A

為碰撞核質(zhì)量，。 在中子非彈性散射情況下，轉(zhuǎn)換公式為： 其中，為第K

個能級的閾能。中子彈性散射勒讓德(Legendre)多項式分布 中子彈性散射角分布常以勒讓德多項式的展開形式給定。散射角余弦x＝cosθ的分布密度函數(shù)為： 其中Pl(x)為l

階勒讓德多項式。 該分布即為n

階勒讓德近似展開。 勒讓德多項式由以下遞推公式確定：

考慮新的分布： 當選取x0，x1，…xn

為Pn+1(x)＝0的根，且 時，fa(x)依照勒讓德多項式展開的前n

項與f

(x)的展開形式相同。因此，可以用fa(x)作為f

(x)的近似分布。

在實際問題中，由于勒讓德多項式展開項數(shù)不夠，可能出現(xiàn)某個為負值的現(xiàn)象。此時可以采用如下近似分布： 其中： 對于該近似分布，可用加抽樣方法進行抽樣： 此時，由于偏倚抽樣而引起的糾偏因子為wK

，也就是說，粒子的權(quán)重要乘上wK。光子康普頓散射角分布 光子的康普頓散射角與其散射前后的能量有關(guān),它的分布密度函數(shù)為： 抽樣方法為：碰撞后運動方向Ωm+1的確定

實驗室系散射角θL確定后，依據(jù)不同的坐標系的表現(xiàn)形式，有不同的確定方法。確定方向余弦um+1，vm+1，wm+1

其中， 方位角

在[0,2π]上均勻分布。 當時，不能使用上述公式，可用下面的簡單公式：確定Ωm+1的球坐標(θm+1，φm+1)

設(shè)Ωm的球坐標分別為(θm,φm)，其中，θ為粒子運動方向與z

軸的夾角，φ為粒子運動方向在xy

平面上投影的方位角。則Ωm+1的球坐標(θm+1,φm+1)分別由下式確定：球形幾何的隨機游動公式

一般幾何的隨機游動公式可以應(yīng)用到球形幾何，而對球?qū)ΨQ問題，使用特殊形式更為方便。下次碰撞點的徑向位置rm+1的確定 兩次碰撞點間的距離L

確定之后，下次碰撞點的徑向位置rm+1的計算公式為： 設(shè)系統(tǒng)的外半徑為R，如rm+1≥R，則粒子逃出系統(tǒng)。粒子碰撞后瞬時運動方向的確定 在球?qū)ΨQ系統(tǒng)中，粒子運動方向用其與徑向夾角余弦來描述。使用球面三角公式，粒子碰撞后瞬時運動方向與徑向夾角余弦cosθm+1的計算公式為： 其中，為在[0,2π]上均勻分布的方位角，為在rm+1

點進入碰撞前瞬時運動方向與rm+1

徑向之間的夾角。點到給定邊界面的距離

在抽樣確定輸運距離、判斷粒子是否穿透系統(tǒng)時，常遇到求由rm

出發(fā)，沿Ωm

方向到達某個區(qū)域表面的距離問題。在記錄對結(jié)果的貢獻時，也常使用類似的量。區(qū)域表面通常是平面或二次曲面。求到達區(qū)域表面的距離問題，實際上是求直線（或半射線）與平面或二次曲面的交點問題。這是蒙特卡羅方法解粒子輸運的各種實際問題時,所遇到的基本幾何問題。點到平面的距離 點沿方向的直線方程為： 該直線到達方程為 的平面的距離為： 當與平面平行時，即 直線與平面無交點。 如果l

為負值，直線與平面也無交點。這時，粒子的運動方向是背離平面的。點到球面的距離 在三維直角坐標系中，設(shè)球心為rc＝(xc,yc,zc)，球半徑為R，則球面方程為： 將直線方程代入該球面方程，得到點r0沿Ω0方向到達球面的距離l

： 其中

當R0≤R

時，即r0

點在球內(nèi)，Δ≥δ2，l

只有一個正根： 當R0＞R

時，即r0

點在球外，分以下三種情況：若δ≥0，l

無正實根，直線與球面無交點。若δ＜0，Δ＜0，l

無實根，直線與球面無交點。若δ＜0，Δ≥0，l

有兩個正實根，直線與球面有兩個交點。

在球坐標系中，不失一般性，設(shè)球心為rc＝0，則球面方程為r＝R。 當r0≤R

時，即r0

點在球內(nèi)，有一個交點： 其中θ0為Ω0與r0

的徑向夾角。 當r0＞R

時，即r0

點在球外，令 當cosθ0≥0時，直線與球面無交點。 當cosθ0＜0時，若d≥R，則直線與球面無交點。 若d＜R，則有兩個交點：點到圓柱面的距離 設(shè)圓柱面的方程為： 其中(xc,yc,0)為圓柱的中心，R

為圓柱底半徑。 點r0沿Ω0方向到達圓柱面的距離l

為： 其中

當R0≤R

時，r0

點在圓柱內(nèi)，如果，則l

有一個正根： 如果，即Ω0平行于圓柱的對稱軸，直線與圓柱面無交點。 當R0＞R

時，r0

點在圓柱外，分以下三種情況：若δ≥0，l

無正實根，直線與圓柱面無交點。若δ＜0，Δ＜0，l

無實根，直線與圓柱面無交點。若δ＜0，Δ≥0且，l

有兩個正實根，直線與圓柱面有兩個交點。 在的情況下，直線與圓柱面不相交。點到圓錐面的距離 設(shè)圓錐頂點在原點，以z

軸為對稱軸，則圓錐面的方程為： 點r0沿Ω0方向到達圓錐面的距離l

為： 其中 如果Ω0與錐面某一母線平行，即，則空腔處理 在粒子輸運問題中，所考慮的系統(tǒng)常有空腔存在，如中空的球殼,平板間的空隙等。粒子輸運時，可有兩種處理空腔的方法：將空腔作為宏觀總截面Σt＝0的區(qū)域,按通常的方法輸運。設(shè)分別為由rm

出發(fā)，沿Ωm

方向到空腔區(qū)域的 近端和遠端的交點，當粒子超過時，以為新的起 點，重新開始輸運。 顯然，這兩種方法在統(tǒng)計上是等價的。等效的邊界條件全反射邊界 在反應(yīng)堆活性區(qū)中，元件盒常常按正方形或六角形排列。假定元件盒足夠多，每個盒結(jié)構(gòu)相同，那么活性區(qū)中每個盒所占的柵胞的物理情況，可以代表整個活性區(qū)中的狀況。

進一步假定，元件盒是圓對稱的，那么每個柵胞中情況，可以用更小的單位（柵元）來反映。比如對六角形柵胞可取其

1/12的ΔOAB

來做代表；正方形柵胞可用其

1/8的ΔO'A'B'

來做代表。這樣一來問題就大大簡化了。

現(xiàn)在的問題是怎樣計算直角三角形柵元的物理量（如通量）。用蒙特卡羅方法如何模擬中子在柵元內(nèi)的運動，反映出整個活性區(qū)對它的影響。 我們可把OA'、OB'、A'B'作為全反射邊界來處理。所謂全反射邊界，就是當中子打到該邊界上時，按鏡面反射的方式，從邊界上全部反射回來，中子的能量與權(quán)重均不改變。 在這種邊界上的反射條件，稱之為全反射條件，就是通常的鏡面反射條件。

在全反射邊界條件下，一條通過活性區(qū)若干個區(qū)域的中子徑跡，可以用柵元ΔO'A'B'

中的一條折線軌道來反映出來。 反過來，在直角三角形柵元ΔO'A'B'

中任一條反射成的折線軌道，都代表了中子在活性區(qū)內(nèi)一條直線軌道的作用。由于系統(tǒng)的對稱性，在活性區(qū)內(nèi)，凡是與柵元內(nèi)位置相當?shù)牡胤?，都有相同的物理情況，因此柵元內(nèi)各處的情況，當然代表了整個活性區(qū)的情況。一般曲面全反射條件 對于一般曲面的全反射，設(shè)入射方向為Ω，入射點的內(nèi)法線方向為n

，則反射方向Ω'

為： 其中 設(shè) 則平面全反射條件 設(shè)三角形柵元的橫截面ΔOAB

在X-Y

平面上，∠OAB＝θ。則邊界OA、OB、AB上的反射都是平面全反射。在任一與X-Y

平面垂直且與X

軸成α角的平面上，全反射條件為： 由此就可得到OA、OB

和AB

邊上的全反射條件，對于OB

邊，α=θ；對于OA

邊，α=0；對于AB

邊，α=π/2。反射層邊界條件 對于具有大反射層的系統(tǒng)，如存放，運輸和生產(chǎn)裂變物質(zhì)的倉庫、車廂和車間等，當中子從里面打到四周墻上或反射層時，還要繼續(xù)對它進行跟蹤。這種跟蹤常常要花費很大的計算量，并且在結(jié)果中引起的方差也比較大。如果在計算這種系統(tǒng)的不同方案中，反射層條件不變，那么這種大量重復(fù)的計算是很不經(jīng)濟的。

中子射入反射層后，一部分被介質(zhì)吸收，只有一部分返回，由于中子的散射慢化，損失一部分能量，因此反射回來的中子有一個能量方向分布。顯然，對這種反射層，不能應(yīng)用全反射條件。不過，我們?nèi)匀豢梢园阉斪鲞吔纾谶吔缟习捶瓷鋵拥奈锢碜饔脕硖幚怼?/p>

比如，如果反射層是一種平板幾何，我們可以用數(shù)值方法或蒙特卡羅方法，預(yù)先算好在各種不同入射能量E

下的反照率β(E)，反射中子的能量分布RE(E→E')。于是代替在反射層中眼蹤中子，我們可在反射層邊界上作如下處理： 一旦中子打入反射層，立即返回，反射后權(quán)重為 其中，E

為射入反射層中子的能量，W

為中子的權(quán)重。反射后的能量Eβ

由反射能譜RE(E→Eβ)中抽樣產(chǎn)生。反射后的方向Ωβ

由半平面各向同性分布或余弦分布中抽樣。反射后的中子位置為入射時的位置。 計算表明，對于大尺寸的反射層來說，這樣的近似，引起的結(jié)果上的誤差是可以忽略的，卻能帶來計算量的大量節(jié)省。記錄貢獻與分析結(jié)果過程

在粒子輸運問題中，除了得到某些量的積分結(jié)果外，還需要得到這些量的方差、協(xié)方差、以及這些量的空間、能量、方向和時間的分布。這些量可以利用分類記錄手續(xù)同時得到。記錄與結(jié)果

為了得