蒙特卡羅方法解粒子輸運問題

1、1.屏蔽問題模型屏蔽問題模型2.直接模擬方法直接模擬方法3.簡單加權法簡單加權法4.統(tǒng)計估計統(tǒng)計估計法法5.指數變換法指數變換法6.蒙特卡羅方法的效率蒙特卡羅方法的效率作作 業(yè)業(yè)輻射（光子和中子）屏蔽問題是蒙特卡羅方法最早廣泛應用的領域之一。本章主要從物理直觀出發(fā)，說明蒙特卡羅方法解決這類粒子輸運問題的基本方法和技巧。而這些方法和技巧對于諸如輻射傳播、多次散射和通量計算等一般粒子輸運問題都是適用的。 在反應堆工程和輻射的測量與應用中，常常要用一些吸收材料做成屏蔽物擋住光子或中子。我們所關心的是經過屏蔽后射線的強度及其能量分布，這就是屏蔽問題。 當屏蔽物的形狀復雜，散射各向異性，材料介質不均勻

2、, 核反應截面與能量、位置有關時，難以用數值方法求解，用蒙特卡羅方法能夠得到滿意的結果。粒子的輸運問題帶有明顯的隨機性質，粒子的輸運過程是一個隨機過程。粒子的運動規(guī)律是根據大量粒子的運動狀況總結出來的，是一種統(tǒng)計規(guī)律。蒙特卡羅模擬，實際上就是模擬相當數量的粒子在介質中運動的狀況，使粒子運動的統(tǒng)計規(guī)律得以重現。不過，這種模擬不是用實驗方法，而是利用數值方法和技巧，即利用隨機數來實現的。 為方便起見，選用平板屏蔽模型，在厚度為 a，長、寬無限的平板左側放置一個強度已知，具有已知能量、方向分布的輻射源 S 。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源發(fā)出的平均一個粒子穿透屏蔽的數

3、目。同時，假定粒子在兩次碰撞之間按直線運動 , 且粒子之間的相互作用可以忽略。 直接模擬方法就是直接從物理問題出發(fā)，模擬粒子的真實物理過程。1) 狀態(tài)參數與狀態(tài)序列2) 模擬運動過程3) 記錄結果粒子在介質中的運動的狀態(tài)，可用一組參數來描述，稱之為狀態(tài)參數狀態(tài)參數。它通常包括：粒子的空間位置 r， 能量 E 和運動方向，以 S( r , E , ) 表示。有時還需要其他的參數，如粒子的 時間 t 和附帶的權重W ，這時狀態(tài)參數狀態(tài)參數 為 S( r , E , , t ,W ) 。狀態(tài)參數狀態(tài)參數 通常要根據所求問題的類型和所用的方法來確定。對于無限平板幾何，取 S( z , E , cos)

4、其中 z 為粒子的位置坐標，為粒子的運動方向與 Z 軸的夾角。對于球對稱幾何 , 取 S( r , E , cos)其中 r 表示粒子所在位置到球心的距離，為粒子的運動方向與其所在位置的徑向夾角。粒子第 m 次碰撞后的狀態(tài)參數為或它表示一個由源發(fā)出的粒子，在介質中經過 m 次碰撞后的狀態(tài)，其中 rm ：粒子在第 m 次碰撞點的位置 Em ：粒子第 m 次碰撞后的能量m：粒子第 m 次碰撞后的運動方向 tm ：粒子到第 m 次碰撞時所經歷的時間Wm ：粒子第 m 次碰撞后的權重有時，也可選為粒子進入第 m 次碰撞時的狀態(tài)參數。),(mmmmE rS),(mmmmmmWtErS一個由源發(fā)出的粒子在

5、介質中運動，經過若干次碰撞后，直到其運動歷史結束（如逃出系統(tǒng)或被吸收等）。假定粒子在兩次碰撞之間按直線運動，其運動方向與能量均不改變，則粒子在介質中的運動過程可用以下碰撞點的狀態(tài)序列狀態(tài)序列 描述： S0 ，S1 ，SM-1 ，SM或者更詳細些 , 用來描述。這里 S0 為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，SM 為粒子的終止狀態(tài)。M 稱為粒子運動的鏈長。這樣的序列稱為粒子隨機運動的歷史，模擬一個粒子的運動過程，就變成確定狀態(tài)序列的問題。MMMMMMEEEErrrr,110110110為簡單起見，這里以中子穿透均勻平板的模型來說明，這時狀態(tài)參數狀態(tài)參數 取 S( z , E , cos)。模擬的步驟

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

7、111111),(exp),()(rr lmmmtmmmtl dElzElzlf0111111),cos(exp),cos()( lmmmtl dEl0111),(r顯然，粒子輸運的自由程數服從指數分布，因此從 f ( l ) 中抽樣確定 l，就是要從積分方程 中解出 l。對于單一介質則下一個碰撞點的位置如果 zma，則中子穿透屏蔽，若 zm0, 則中子被反射出屏蔽。這兩種情況，均視為中子歷史終止。ln),(0111 lmmmtl dEl r)(ln1mtEl11111cos)(lncosmmtmmmmEzlzz(3) 確定被碰撞的原子核 ：通常介質由幾種原子核組成，中子與核碰撞時，要確定與哪

8、一種核碰撞。設介質由A、B、C 三種原子核組成，其核密度分別為NA、NB、NC，則介質的宏觀總截面為：其中 分別為核A、B、C 的宏觀總截面。其定義如下： 分別表示()核的宏觀總截面、核密度和微觀總截面。)()()()(1111mCtmBtmAtmtEEEECtBtAt,)()(1)()(1)(mtmtENE)()(1)()(1)(mtmtENE、由于中子截面表示中子與核碰撞可能性的大小，因此，很自然地，中子與A、B、C 核發(fā)生碰撞的幾率分別為：利用離散型隨機變量的抽樣方法，確定碰撞核種類：)()(,)()(,)()(111111mtmCtCmtmBtBmtmAtAEEPEEPEEP核碰撞與核

9、碰撞與核碰撞與CBPPAPBAA(4) 確定碰撞類型 ：確定了碰撞的核(比如B核)后，就要進一步確定碰撞類型。中子與核的反應類型有彈性散射、非彈性散射、(n,2n)反應，裂變和俘獲等，它們的微觀截面分別為則有各種反應發(fā)生的幾率分別為)()()()()(111)2,(11mBcmBfmBnnmBinmBelEEEEE和、)()()()()()(111)2,(111mBcmBfmBnnmBinmBelmBtEEEEEE)()()()()()()()()()(111111)2,()2,(1111mBtmBccmBtmBffmBtmBnnnnmBtmBininmBtmBelelEEPEEPEEPEEP

10、EEP利用離散型隨機變量的抽樣方法，確定反應類型。在屏蔽問題中，中子與核反應常只有彈性散射和吸收兩種類型，吸收截面為：這時，總截面為：發(fā)生彈性散射的幾率為：若 ，則為彈性散射；否則為吸收，發(fā)生吸收反應意味著中子的歷史終止。)()(11mBtmBelelEEPelP)()()(111mBcmBfmBaEEE)()()(111mBamBelmBtEEE(5) 確定碰撞后的能量與運動方向：如果中子被碰撞核吸收，則其輸運歷史結束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運動方向。中子能量 Em 為：A是碰撞核的質量與中子質量之比，一般就取元素的原子量；C 為質心系中中子散射前后方向間的夾角，即偏轉

11、角。 可從質心系中彈性散射角分布fC(C) 中抽樣產生。實驗室系散射角L的余弦L為： 21)11(cos112AArrrEECmmCCcosCCLAAA2112如果給出實驗室系散射角余弦分布 fL(L)，可直接從 fL(L)中抽取L，此時能量Em與L的關系式為：確定了實驗室系散射角L后，再使用球面三角公式確定cosm ：其中為在0,2上均勻分布的方位角。cossinsincoscoscos11LmLmm222211) 1(LLmmAAEE至此，由Sm-1完全可以確定Sm。 因此，當中子由源出發(fā)后，即S0確定后，重復步驟 (2)(5)，直到中子游動歷史終止。于是得到了一個中子的隨機游動歷史 S0

12、 ，S1 ，SM-1 ，SM，即也就是模擬了一個由源發(fā)出的中子的運動過程。MMMMMMEEEEzzzzcos,cos,cos,cos,110110110以上模擬過程可分為兩大步：第一步確定粒子的初始狀態(tài)S0，第二步由狀態(tài)Sm-1來確定狀態(tài)Sm。這第二步又分為兩個過程：第一個過程是確定碰撞點位置zm ，稱為輸運過程；第二個過程是確定碰撞后粒子的能量及運動方向，稱為碰撞過程。對于中子而言，碰撞過程是先確定散射角，進而確定能量和運動方向；而對于光子，碰撞過程是先確定能量，再確定散射角以及運動方向。重復這兩個過程，直至粒子的歷史終止。這種模擬過程，是解任何類型的粒子輸運問題所共有的，它是蒙特卡羅方法解

13、題的基本手段。在獲得中子的隨機游動歷史后，我們要對所要計算的物理量進行估計。對于屏蔽問題，我們要計算中子的穿透率?？疾烀總€中子的隨機游動歷史，它可能穿透屏蔽（zMa），可能被屏蔽發(fā)射回來（zM0），或者被吸收。設第 n 個中子對穿透的貢獻為n ，則如果我們共跟蹤了N 個中子，則穿透屏蔽的中子數為：，或者被吸收當當0, 0, 1MMnzazNnnN11則穿透屏蔽概率的近似值為：它是穿透率的一個無偏估計。我們稱這種直觀地模擬過程和估計方法為直接模擬方法。在置信水平 10.95 時， 的誤差為：其中 為n的均方差，由于n是一個服從二項分布的隨機變量，所以或NnnNNNNP11)1(1)1 (NPNN

14、PPN96. 1)1()1 ()1 ()1()1(22NNPPPP為得到中子穿透屏蔽的能量、角分布，將能量、角度范圍分成若干個間隔：其中Emax，Emin分別表示能量的上、下限，對于穿透屏蔽的中子按其能量、方向分間隔記錄。設一穿透屏蔽的中子能量為EM，其運動方向與Z軸夾角為M，若能量EM屬于第 i 個能量間隔Ei，角度M屬于第 j 個角度間隔j，則分別在第 i 個能量計數器及第 j 個角度計數器中加 1。2010max01minJIEEEEE跟蹤 N 個中子后，則分別為穿透中子的能量分布和角分布。其中N1,i 和 N2,i 分別為第 i 個能量和第 j 個角度間隔的穿透中 子數。歸一后分別為：

15、JjNNPIiENNPjjNiiNji, 2 , 1, 2 , 1, 2)1(, 1)1(,2, 1JjNNPPPIiENNPPPjjNNNiiNNNjjii, 2 , 1, 2 , 11, 2)*1()1()*1(1, 1)*1()1()*1(,2,2, 1, 1從模擬物理過程來說，直接模擬法是最簡單、也是最基本的方法。但是，在直接模擬法中，不管中子在屏蔽中經過多少次碰撞，只要在介質中被吸收，對穿透的貢獻就為零；因此在所跟蹤的粒子中絕大部分都對穿透沒有貢獻。而在許多屏蔽問題中，穿透率的數量級在10-6到10-8。進一步，如果我們要求穿透率達相對誤差小于1，即那么，N 要大到驚人的數量級101

16、0到1012。顯然，這時用直接模擬法計算不是很有效。%1)1 (|PNPNPPPNPPP屏蔽物一般是由吸收強的介質組成，因此在每次碰撞時，粒子很有可能被吸收而停止跟蹤?，F在改變模擬方法，在判斷碰撞類型 時，可以認為粒子的 部分是彈性散射，而其余部分被吸收，即人為地把中子分成兩部分，一部分彈性散射，一部分吸收。彈性散射這部分繼續(xù)跟蹤；吸收部分則停止跟蹤。也就是說，我們利用中子權重的變化來反應繼續(xù)彈性散射的部分。這就是簡單加權法的基本思想。elP)()(11mBtmBelelEEP顯然，在加權法中中子的權重W 已成為中子狀態(tài)參數的組成部分。這時，中子歷史成為：對源中子，取W0=1。經過碰撞中子權重

17、的變化為：因子 稱為尚存因子。,cos,cos,cos,cos,110110110110MMMMMMMMWWWWEEEEzzzz)()(111mBtmBelmmEEWW)()(11mBtmBelEE這時，第 n 個中子對穿透的貢獻為：如果我們共跟蹤了N個中子，則穿透率P的無偏估計為：類似地，可以得到穿透中子的能量分布和角分布。只不過在對各計數器進行的加 1 操作改為加WM。0, 0,MMMnzazW當當NnnNNP1)2(1簡單加權法的方差估計為：與直接模擬法相比，有注意到n1，有這表明簡單加權法的方差小于直接模擬法的方差。這是因為加權法比直接模擬法減少了一次隨機抽樣。2)2(122)(1NN

18、nnPNNnnnN1222)(122加權法的思想在蒙特卡羅方法中用途很廣泛。例如，對于具有中子增殖反應，如裂變，(n,2n)，(n,3n) 反應的中子輸運問題，一個中子與核發(fā)生碰撞后，根據反應的類型會產生不同數量的次級中子，每個次級中子又會產生新的次級中子，這樣鏈鎖反應 下去，使得用直接模擬法模擬每一個中子是非常困難的。這種情況可以利用加權法來處理。中子與核發(fā)生碰撞 后，產生的次級中子平均數為： 這里f 為裂變次級中子數。于是，碰撞后的權重為：而決定碰撞類型的幾率分別為：其中加權法的思想，還可以應用到連續(xù)分布情況和偏倚抽樣的問題tffnnnninel)3 ,()2,(32111mmWWtfft

19、nntnntintel,3,2,)3 ,()2,(ffnnnninelt)3 ,()2,(3211加權法雖然改進了直接模擬法，但它同樣只關心中子是否穿透屏蔽這一信息，因此對每一個中子歷史的信息利用得很不充分。統(tǒng)計估計法能夠較多地利用中子的歷史信息，因而能得到更好的結果。一個中子，可能在介質內不發(fā)生碰撞而直接穿透屏蔽，也可能在介質內發(fā)生一次碰撞后再穿透屏蔽，或經過二次碰撞穿透屏蔽，等等，這些事件是互不相容的，因此穿透概率P 可表示為：其中Pm 是中子恰好經過 m 次碰撞而穿透屏蔽的概率。這表明，可以用求 Pm (m=0,1, ) 的方法得到P。這樣，中子對穿透概率的貢獻就不只限于末次碰撞了。0m

20、mPP01S0S1SmmP1P0PmZa0設中子的歷史為：根據該中子的歷史，我們可以估計出中子恰好經過 m 次碰撞后，穿透屏蔽的部分顯然，具有初態(tài) S0( 0, E0, cos0,W0 ) 的中子，未經碰撞直接穿透的部分是：,cos,cos,cos,cos,110110110110MMMMMMMMWWWWEEEEzzzz1, 1 , 0,MmPm0000cos)(expaEWPt類似地，在經過了第 m 次碰撞后的中子具有狀態(tài) Sm( zm, Em, cosm,Wm ) ，其可能穿透的部分，正好是一個中子恰好經過 m 次碰撞穿透的部分：這里的這種估計技巧，由于是對每次碰撞后的狀態(tài)，求其后未經碰撞

21、直接穿透的貢獻，因此該方法也稱為最后自由飛行估計。1, 1 , 000coscos)(expMmzaEWPmmmmtmm其它，于是得到該中子對穿透的貢獻：如果我們共跟蹤了N個中子，則穿透率P的估計為：其方差估計為：10MmmPPNnNnPNP1)3()(12)3(122)(1)(NNnPnPNP在直接模擬方法中，相對誤差為其中 為與置信水平 1相應的量。如果構造一個新的概率模型，使得該模型的穿透率P*與原模型的穿透率P之間存在關系：使用直接模擬方法 , 相對誤差為PNPNPPPN)1 (PKP*PN 如果令*，即這意味著，達到同樣的相對誤差，跟蹤粒子的數目縮小 K 倍，從而減少 K 倍的計算量。指數變換法就是構造一個新的概率模型的一個有效方法。KNNPKNPNPN*構造如下偽過程：宏觀總截面為散射截面仍為el(E)。其中 Emin、Emax 分別為能量的下限和上限，為粒子的運動方向與 Z 軸的夾角?？梢宰C明這個偽過程的穿透概率P* 與原過程的穿