蒙特卡羅方法解粒子輸運(yùn)問題

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

2、勻 , 核反應(yīng)截面與能量、位置有關(guān)時，難以用數(shù)值方法求解，用蒙特卡羅方法能夠得到滿意的結(jié)果。粒子的輸運(yùn)問題帶有明顯的隨機(jī)性質(zhì)，粒粒子的輸運(yùn)問題帶有明顯的隨機(jī)性質(zhì)，粒子的輸運(yùn)過程是一個隨機(jī)過程子的輸運(yùn)過程是一個隨機(jī)過程。粒子的運(yùn)動規(guī)律是根據(jù)大量粒子的運(yùn)動狀況總結(jié)出來的，是一種統(tǒng)計規(guī)律。 蒙特卡羅模擬，實(shí)際上就是模擬相當(dāng)數(shù)量蒙特卡羅模擬，實(shí)際上就是模擬相當(dāng)數(shù)量的粒子在介質(zhì)中運(yùn)動的狀況，使粒子運(yùn)動的統(tǒng)的粒子在介質(zhì)中運(yùn)動的狀況，使粒子運(yùn)動的統(tǒng)計規(guī)律得以重現(xiàn)計規(guī)律得以重現(xiàn)。 不過，這種模擬不是用實(shí)驗(yàn)方法，而是利用數(shù)值方法和技巧，即利用隨機(jī)數(shù)來實(shí)現(xiàn)的。事件注釋事件注釋1.中子散射 光子產(chǎn)生2.裂變 光子產(chǎn)

3、生3.中子俘獲4.中子逃逸5.光子散射6.光子逃逸7.光子俘獲 為方便起見，選用平板屏蔽模型，在厚度為 a，長、寬無限的平板左側(cè)放置一個強(qiáng)度已知，具有已知能量、方向分布的輻射源 S 。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。 穿透率就是由源發(fā)出的平均一個粒子穿透屏蔽的數(shù)目穿透率就是由源發(fā)出的平均一個粒子穿透屏蔽的數(shù)目。同時，假定粒子在兩次碰撞之間按直線運(yùn)動 , 且粒子之間的相互作用可以忽略。1) 狀態(tài)參數(shù)與狀態(tài)序列2) 模擬運(yùn)動過程3) 記錄結(jié)果 直接模擬方法就是直接從物理問題出發(fā)，模擬粒子的真實(shí)物理過程。粒子在介質(zhì)中的運(yùn)動的狀態(tài)，可用一組參數(shù)來描述，稱之為狀態(tài)參數(shù)狀態(tài)參數(shù)。它通常包括

4、：粒子的空間位置 r， 能量 E 和運(yùn)動方向，以 S( r , E , ) 表示。有時還需要其他的參數(shù)，如粒子的 時間 t 和附帶的權(quán)重W ，這時狀態(tài)參數(shù)狀態(tài)參數(shù) 為 S( r , E , , t ,W ) 。狀態(tài)參數(shù)狀態(tài)參數(shù) 通常要根據(jù)所求問題的類型和所用的方法來確定。對于無限平板幾何，取 S( z , E , cos)其中 z 為粒子的位置坐標(biāo)，為粒子的運(yùn)動方向與 Z 軸的夾角。對于球?qū)ΨQ幾何 , 取 S( r , E , cos)其中 r 表示粒子所在位置到球心的距離，為粒子的運(yùn)動方向與其所在位置的徑向夾角。粒子第 m 次碰撞后的狀態(tài)參數(shù)為或它表示一個由源發(fā)出的粒子，在介質(zhì)中經(jīng)過它表示一

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

6、動過程可用以下碰撞點(diǎn)的狀態(tài)序列狀態(tài)序列 描述： S0 ，S1 ，SM-1 ，SM或者更詳細(xì)些 , 用來描述。這里 S0 為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，SM 為粒子的終止?fàn)顟B(tài)。M 稱為粒子運(yùn)動的鏈長。這樣的序列稱為粒子隨機(jī)運(yùn)動的歷史，模擬一個粒子的運(yùn)動過程，就變成確定狀態(tài)序列的問題。MMMMMMEEEErrrr,110110110為簡單起見，這里以中子穿透均勻平板的模型來說明，這時狀態(tài)參數(shù)狀態(tài)參數(shù) 取 S( z , E , cos)。模擬的步驟如下：(1) 確定初始狀態(tài)確定初始狀態(tài) S0 ：確定粒子的初始狀態(tài)，實(shí)際上就是要從中子源的空間位置、能量和方向分布中抽樣。設(shè)源分布為則分別從各自的分布中

7、抽樣確定初始狀態(tài)。對于平板情況，抽樣得到 z00。)(cos)()()cos,(030201000fEfzfEzf)()(001zzf(2) 確定下一個碰撞點(diǎn)確定下一個碰撞點(diǎn) ：已知狀態(tài)已知狀態(tài)Sm-1，要確定狀態(tài)，要確定狀態(tài)Sm，首先要確定下一個首先要確定下一個碰撞點(diǎn)的位置碰撞點(diǎn)的位置zm。在相鄰兩次碰撞之間，中子的輸運(yùn)長度 l 服從如下分布： lmmmtmmmtl dElEllf0111111),(exp),()(rr對于平板模型對于平板模型，l 服從分布：其中，t 為介質(zhì)的中子宏觀總截面，積分 稱為粒子輸運(yùn)的自由程數(shù)，系統(tǒng)的大小通常就是用系統(tǒng)的自由程數(shù)表示的。 lmmmtmmmtl dE

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

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

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

11、tmBelelEEPEEPEEPEEPEEP利用離散型隨機(jī)變量的抽樣方法離散型隨機(jī)變量的抽樣方法，確定反應(yīng)類型。在屏蔽問題中，中子與核反應(yīng)常只有彈性散射和吸收兩種類型，吸收截面為：這時，總截面為：發(fā)生彈性散射的幾率為：若 ，則為彈性散射；否則為吸收，發(fā)生吸收反應(yīng)意味著中子的歷史終止。)()(11mBtmBelelEEPelP)()()(111mBcmBfmBaEEE)()()(111mBamBelmBtEEE(5) 確定碰撞后的能量與運(yùn)動方向確定碰撞后的能量與運(yùn)動方向：如果中子被碰撞核吸收，則其輸運(yùn)歷史結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運(yùn)動方向。中子能量 Em 為：A是碰撞核的

12、質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量；C 為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角。 可從質(zhì)心系中彈性散射角分布fC(C) 中抽樣產(chǎn)生。實(shí)驗(yàn)室系散射角L的余弦L為： 21)11(cos112AArrrEECmmCCcosCCLAAA2112如果給出實(shí)驗(yàn)室系散射角余弦分布 fL(L)，可直接從 fL(L)中抽取L，此時能量Em與L的關(guān)系式為：確定了實(shí)驗(yàn)室系散射角L后，再使用球面三角公式確定cosm ：其中為在0,2上均勻分布的方位角。cossinsincoscoscos11LmLmm222211) 1(LLmmAAEE至此，由Sm-1完全可以確定Sm。 因此，當(dāng)中子由源出發(fā)后，即S0確定

13、后，重復(fù)步驟 (2)(5)，直到中子游動歷史終止。于是得到了一個中子的隨機(jī)游動歷史 S0 ，S1 ，SM-1 ，SM，即也就是模擬了一個由源發(fā)出的中子的運(yùn)動過程模擬了一個由源發(fā)出的中子的運(yùn)動過程。MMMMMMEEEEzzzzcos,cos,cos,cos,110110110以上模擬過程可分為兩大步以上模擬過程可分為兩大步： 第一步確定粒子的初始狀態(tài)S0，第二步由狀態(tài)Sm-1來確定狀態(tài)Sm。 這第二步又分為兩個過程這第二步又分為兩個過程：第一個過程是確定碰第一個過程是確定碰撞點(diǎn)位置撞點(diǎn)位置zm （自由程抽樣），稱為（自由程抽樣），稱為輸運(yùn)過程輸運(yùn)過程；第二個第二個過程是確定碰撞后粒子的能量及運(yùn)動

14、方向，稱為過程是確定碰撞后粒子的能量及運(yùn)動方向，稱為碰撞碰撞過程過程。 對于中子而言，碰撞過程是先確定散射角，進(jìn)而確定能量和運(yùn)動方向；而對于光子，碰撞過程是先確定能量，再確定散射角以及運(yùn)動方向。重復(fù)這兩個過程，直至粒子的歷史終止。這種模擬過程，是解任何類型的粒子輸運(yùn)問題所這種模擬過程，是解任何類型的粒子輸運(yùn)問題所共有的，它是蒙特卡羅方法解題的基本手段共有的，它是蒙特卡羅方法解題的基本手段。在獲得中子的隨機(jī)游動歷史后，我們要對所要計算的物理量進(jìn)行估計。對于屏蔽問題，我們要計算中子的穿透率穿透率。考察每個中子的隨機(jī)游動歷史，它可能穿透屏蔽（zMa），可能被屏蔽發(fā)射回來（zM0），或者被吸收。 設(shè)第

15、 n 個中子對穿透的貢獻(xiàn)為n ，則如果我們共跟蹤了N 個中子，則穿透屏蔽的中子數(shù)為：，或者被吸收當(dāng)當(dāng)0, 0, 1MMnzazNnnN11則穿透屏蔽概率的近似值為：它是穿透率的一個無偏估計。我們稱這種直觀地模擬過程和估計方法為直接模擬方法。在置信水平 10.95 時， 的誤差為：其中 為n的均方差，由于n是一個服從二項(xiàng)分布的隨機(jī)變量，所以或NnnNNNNP11)1(1)1 (NPNNPPN96. 1)1()1 ()1 ()1()1(22NNPPPP為得到中子穿透屏蔽的能量能量、角分布角分布，將能量、角度范圍分成若干個間隔：其中Emax，Emin分別表示能量的上、下限，對于穿透屏蔽的中子按其能量

16、、方向分間隔記錄。設(shè)一穿透屏蔽的中子能量為EM，其運(yùn)動方向與Z軸夾角為M，若能量EM屬于第 i 個能量間隔Ei，角度M屬于第 j 個角度間隔j，則分別在第 i 個能量計數(shù)器及第 j 個角度計數(shù)器中加 1。2010max01minJIEEEEE跟蹤 N 個中子后，則分別為穿透中子的能量分布和角分布。其中N1,i 和 N2,i 分別為第 i 個能量和第 j 個角度間隔的穿透中 子數(shù)。歸一后分別為：JjNNPIiENNPjjNiiNji, 2 , 1, 2 , 1, 2)1(, 1)1(,2, 1JjNNPPPIiENNPPPjjNNNiiNNNjjii, 2 , 1, 2 , 11, 2)*1()

17、1()*1(1, 1)*1()1()*1(,2,2, 1, 1從模擬物理過程來說，直接模擬法是最簡單、也是最基本的方法。但是，在直接模擬法中，不管中子在直接模擬法中，不管中子在屏蔽中經(jīng)過多少次碰撞，只要在介質(zhì)中被吸收，對在屏蔽中經(jīng)過多少次碰撞，只要在介質(zhì)中被吸收，對穿透的貢獻(xiàn)就為零穿透的貢獻(xiàn)就為零；因此在所跟蹤的粒子中絕大部分都對穿透沒有貢獻(xiàn)。而在許多屏蔽問題中，穿透率的數(shù)量級在10-6到10-8。進(jìn)一步，如果我們要求穿透率達(dá)相對誤差小于1，即那么，N 要大到驚人的數(shù)量級1010到1012。顯然，這時用直接模擬法計算不是很有效。%1)1 (|PNPNPPPNPPP屏蔽物一般是由吸收強(qiáng)的介質(zhì)組成

18、，因此在每次碰撞時，粒子很有可能被吸收而停止跟蹤。 現(xiàn)在改變模擬方法，在判斷碰撞類型 時，可以認(rèn)為粒子的 部分是彈性散射，而其余部分被吸收，即人為地把中子分成兩部分，人為地把中子分成兩部分，一部分彈性散射，一部分吸收一部分彈性散射，一部分吸收。 彈性散射這部分繼續(xù)跟蹤；吸收部分則停止跟蹤。 也就是說，我們利用中子權(quán)重的變化來反應(yīng)繼續(xù)彈性散射的部分。這就是簡單加權(quán)法的基本思想簡單加權(quán)法的基本思想。elP)()(11mBtmBelelEEP顯然，在加權(quán)法中中子的權(quán)重W 已成為中子狀態(tài)參數(shù)的組成部分。這時，中子歷史成為：對源中子，取W0=1。經(jīng)過碰撞中子權(quán)重的變化為：因子 稱為尚存因子。,cos,c

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

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

21、tnntintel,3,2,)3 ,()2,(ffnnnninelt)3 ,()2,(3211加權(quán)法雖然改進(jì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。這樣，中子對穿透概率的貢獻(xiàn)就不只限于末次碰撞了。0mmP

22、P01S0S1SmmP1P0PmZa0設(shè)中子的歷史為：根據(jù)該中子的歷史，我們可以估計出中子恰好經(jīng)過 m 次碰撞后，穿透屏蔽的部分顯然，具有初態(tài) S0( 0, E0, cos0,W0 ) 的中子，未經(jīng)碰撞直接穿透的部分是：,cos,cos,cos,cos,110110110110MMMMMMMMWWWWEEEEzzzz1, 1 , 0,MmPm0000cos)(expaEWPt類似地，在經(jīng)過了第 m 次碰撞后的中子具有狀態(tài) Sm( zm, Em, cosm,Wm ) ，其可能穿透的部分，正好是一個中子恰好經(jīng)過 m 次碰撞穿透的部分：這里的這種估計技巧，由于是對每次碰撞后的狀態(tài)，求其后未經(jīng)碰撞直接穿透的貢獻(xiàn)求其后未經(jīng)碰撞直接穿透的貢獻(xiàn)，因此該方法也稱為最后自由飛行估計最后自由飛行估計。1, 1 , 000coscos)(expMmzaEWPmmmmtmm其它，于是得到該中子對穿透的貢獻(xiàn)：如果我們共跟蹤了N個中子，則穿透率P的估計為：其方差估計為：10MmmPPNnNnPNP1)3()(12)3(122)(1)(NNnPnPNP在直接模擬方法中，相對誤差為其中 為與置信水平 1相應(yīng)的量。如果構(gòu)造一個新的概率模型，使得該模型的穿透率P*與原模型的穿透率P之間存在關(guān)系：使用直接模擬方法 , 相對誤差為PNPNPP