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

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

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

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

4、其中 z 為粒子的位置坐標(biāo)，為粒子的運(yùn)動(dòng)方向與 z 軸的夾角。對(duì)于球?qū)ΨQ幾何 , 取 s( r , e , cos)其中 r 表示粒子所在位置到球心的距離，為粒子的運(yùn)動(dòng)方向與其所在位置的徑向夾角。1) 狀態(tài)參數(shù)與狀態(tài)序列粒子第 m 次碰撞后的狀態(tài)參數(shù)為或它表示一個(gè)由源發(fā)出的粒子，在介質(zhì)中經(jīng)過 m 次碰撞后的狀態(tài)，其中 rm ：粒子在第 m 次碰撞點(diǎn)的位置 em ：粒子第 m 次碰撞后的能量m：粒子第 m 次碰撞后的運(yùn)動(dòng)方向 tm ：粒子到第 m 次碰撞時(shí)所經(jīng)歷的時(shí)間wm ：粒子第 m 次碰撞后的權(quán)重有時(shí)，也可選為粒子進(jìn)入第 m 次碰撞時(shí)的狀態(tài)參數(shù)。),(mmmme rs),(mmmmmmwte

5、rs一個(gè)由源發(fā)出的粒子在介質(zhì)中運(yùn)動(dòng)，經(jīng)過若干次碰撞后，直到其運(yùn)動(dòng)歷史結(jié)束（如逃出系統(tǒng)或被吸收等）。假定粒子在兩次碰撞之間按直線運(yùn)動(dòng)，其運(yùn)動(dòng)方向與能量均不改變，則粒子在介質(zhì)中的運(yùn)動(dòng)過程可用以下碰撞點(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)動(dòng)的鏈長。這樣的序列稱為粒子隨機(jī)運(yùn)動(dòng)的歷史，模擬一個(gè)粒子的運(yùn)動(dòng)過程，就變成確定狀態(tài)序列的問題。mmmmmmeeeerrrr,110110110為簡單起見，這里以中子穿透均勻平板的模型來說明，這時(shí)狀態(tài)參數(shù)狀態(tài)參數(shù) 取 s( z , e

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

7、。 lmmmtmmmtl delellf0111111),(exp),()(rr lmmmtmmmtl delzelzlf0111111),cos(exp),cos()( lmmmtl del0111),(r顯然，粒子輸運(yùn)的自由程數(shù)服從指數(shù)分布，因此從 f ( l ) 中抽樣確定 l，就是要從積分方程 中解出 l。對(duì)于單一介質(zhì)則下一個(gè)碰撞點(diǎn)的位置如果 zma，則中子穿透屏蔽，若 zm0, 則中子被反射出屏蔽。這兩種情況，均視為中子歷史終止。ln),(0111 lmmmtl del r)(ln1mtel11111cos)(lncosmmtmmmmezlzz(3) 確定被碰撞的原子核 ：通常介質(zhì)由

8、幾種原子核組成，中子與核碰撞時(shí)，要確定與哪一種核碰撞。設(shè)介質(zhì)由a、b、c 三種原子核組成，其核密度分別為na、nb、nc，則介質(zhì)的宏觀總截面為：其中 分別為核a、b、c 的宏觀總截面。其定義如下： 分別表示()核的宏觀總截面、核密度和微觀總截面。)()()()(1111mctmbtmatmteeeectbtat,)()(1)()(1)(mtmtene)()(1)()(1)(mtmtene、由于中子截面表示中子與核碰撞可能性的大小，因此，很自然地，中子與a、b、c 核發(fā)生碰撞的幾率分別為：利用離散型隨機(jī)變量的抽樣方法，確定碰撞核種類：)()(,)()(,)()(111111mtmctcmtmbt

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

10、mbtmbeleleepeepeepeepeep利用離散型隨機(jī)變量的抽樣方法，確定反應(yīng)類型。在屏蔽問題中，中子與核反應(yīng)常只有彈性散射和吸收兩種類型，吸收截面為：這時(shí)，總截面為：發(fā)生彈性散射的幾率為：若 ，則為彈性散射；否則為吸收，發(fā)生吸收反應(yīng)意味著中子的歷史終止。)()(11mbtmbeleleepelp)()()(111mbcmbfmbaeee)()()(111mbambelmbteee(5) 確定碰撞后的能量與運(yùn)動(dòng)方向：如果中子被碰撞核吸收，則其輸運(yùn)歷史結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運(yùn)動(dòng)方向。中子能量 em 為：a是碰撞核的質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量；c

11、為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角。 可從質(zhì)心系中彈性散射角分布fc(c) 中抽樣產(chǎn)生。實(shí)驗(yàn)室系散射角l的余弦l為： 21)11(cos112aarrreecmmcccoscclaaa2112如果給出實(shí)驗(yàn)室系散射角余弦分布 fl(l)，可直接從 fl(l)中抽取l，此時(shí)能量em與l的關(guān)系式為：確定了實(shí)驗(yàn)室系散射角l后，再使用球面三角公式確定cosm ：其中為在0,2上均勻分布的方位角。cossinsincoscoscos11lmlmm222211) 1(llmmaaee至此，由sm-1完全可以確定sm。 因此，當(dāng)中子由源出發(fā)后，即s0確定后，重復(fù)步驟 (2)(5)，直到中子游動(dòng)歷史終

12、止。于是得到了一個(gè)中子的隨機(jī)游動(dòng)歷史 s0 ，s1 ，sm-1 ，sm，即也就是模擬了一個(gè)由源發(fā)出的中子的運(yùn)動(dòng)過程。mmmmmmeeeezzzzcos,cos,cos,cos,110110110以上模擬過程可分為兩大步：第一步確定粒子的初始狀態(tài)s0，第二步由狀態(tài)sm-1來確定狀態(tài)sm。這第二步又分為兩個(gè)過程：第一個(gè)過程是確定碰撞點(diǎn)位置zm ，稱為輸運(yùn)過程；第二個(gè)過程是確定碰撞后粒子的能量及運(yùn)動(dòng)方向，稱為碰撞過程。對(duì)于中子而言，碰撞過程是先確定散射角，進(jìn)而確定能量和運(yùn)動(dòng)方向；而對(duì)于光子，碰撞過程是先確定能量，再確定散射角以及運(yùn)動(dòng)方向。重復(fù)這兩個(gè)過程，直至粒子的歷史終止。這種模擬過程，是解任何類型

13、的粒子輸運(yùn)問題所共有的，它是蒙特卡羅方法解題的基本手段。在獲得中子的隨機(jī)游動(dòng)歷史后，我們要對(duì)所要計(jì)算的物理量進(jìn)行估計(jì)。對(duì)于屏蔽問題，我們要計(jì)算中子的穿透率?？疾烀總€(gè)中子的隨機(jī)游動(dòng)歷史，它可能穿透屏蔽（zma），可能被屏蔽發(fā)射回來（zm0），或者被吸收。設(shè)第 n 個(gè)中子對(duì)穿透的貢獻(xiàn)為n ，則如果我們共跟蹤了n 個(gè)中子，則穿透屏蔽的中子數(shù)為：3) 記錄結(jié)果，或者被吸收當(dāng)當(dāng)0, 0, 1mmnzaznnnn11則穿透屏蔽概率的近似值為：它是穿透率的一個(gè)無偏估計(jì)。我們稱這種直觀地模擬過程和估計(jì)方法為直接模擬方法。在置信水平 10.95 時(shí)， 的誤差為：其中 為n的均方差，由于n是一個(gè)服從二項(xiàng)分布的隨機(jī)

14、變量，所以或nnnnnnnp11)1(1)1 (npnnppn96. 1)1()1 ()1 ()1()1(22nnpppp為得到中子穿透屏蔽的能量、角分布，將能量、角度范圍分成若干個(gè)間隔：其中emax，emin分別表示能量的上、下限，對(duì)于穿透屏蔽的中子按其能量、方向分間隔記錄。設(shè)一穿透屏蔽的中子能量為em，其運(yùn)動(dòng)方向與z軸夾角為m，若能量em屬于第 i 個(gè)能量間隔ei，角度m屬于第 j 個(gè)角度間隔j，則分別在第 i 個(gè)能量計(jì)數(shù)器及第 j 個(gè)角度計(jì)數(shù)器中加 1。2010max01minjieeeee跟蹤 n 個(gè)中子后，則分別為穿透中子的能量分布和角分布。其中n1,i 和 n2,i 分別為第 i

15、個(gè)能量和第 j 個(gè)角度間隔的穿透中 子數(shù)。歸一后分別為：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 從模擬物理過程來說，直接模擬法是最簡單、也是最基本的方法。但是，在直接模擬法中，不管中子在屏蔽中經(jīng)過多少次碰撞，只要在介質(zhì)中被吸收，對(duì)穿透的貢獻(xiàn)就為零；因此在所跟蹤的粒子中絕大部分都對(duì)穿透沒有貢獻(xiàn)。而在許多屏蔽問題中，穿透率的數(shù)量級(jí)在10-6到10-8。進(jìn)一步，如果我們要求穿

16、透率達(dá)相對(duì)誤差小于1，即那么，n 要大到驚人的數(shù)量級(jí)1010到1012。顯然，這時(shí)用直接模擬法計(jì)算不是很有效。%1)1 (|pnpnpppnppp屏蔽物一般是由吸收強(qiáng)的介質(zhì)組成，因此在每次碰撞時(shí)，粒子很有可能被吸收而停止跟蹤?，F(xiàn)在改變模擬方法，在判斷碰撞類型 時(shí)，可以認(rèn)為粒子的 部分是彈性散射，而其余部分被吸收，即人為地把中子分成兩部分，一部分彈性散射，一部分吸收。彈性散射這部分繼續(xù)跟蹤；吸收部分則停止跟蹤。也就是說，我們利用中子權(quán)重的變化來反應(yīng)繼續(xù)彈性散射的部分。這就是簡單加權(quán)法的基本思想。1) 簡單加權(quán)法elp)()(11mbtmbeleleep顯然，在加權(quán)法中中子的權(quán)重w 已成為中子狀態(tài)

17、參數(shù)的組成部分。這時(shí)，中子歷史成為：對(duì)源中子，取w0=1。經(jīng)過碰撞中子權(quán)重的變化為：因子 稱為尚存因子。,cos,cos,cos,cos,110110110110mmmmmmmmwwwweeeezzzz)()(111mbtmbelmmeeww)()(11mbtmbelee這時(shí)，第 n 個(gè)中子對(duì)穿透的貢獻(xiàn)為：如果我們共跟蹤了n個(gè)中子，則穿透率p的無偏估計(jì)為：類似地，可以得到穿透中子的能量分布和角分布。只不過在對(duì)各計(jì)數(shù)器進(jìn)行的加 1 操作改為加wm。0, 0,mmmnzazw當(dāng)當(dāng)nnnnnp1)2(1簡單加權(quán)法的方差估計(jì)為：與直接模擬法相比，有注意到n1，有這表明簡單加權(quán)法的方差小于直接模擬法的方

18、差。這是因?yàn)榧訖?quán)法比直接模擬法減少了一次隨機(jī)抽樣。2) 簡單加權(quán)法的方差2)2(122)(1nnnnpnnnnnn1222)(122加權(quán)法的思想在蒙特卡羅方法中用途很廣泛。例如，對(duì)于具有中子增殖反應(yīng)，如裂變，(n,2n)，(n,3n) 反應(yīng)的中子輸運(yùn)問題，一個(gè)中子與核發(fā)生碰撞后，根據(jù)反應(yīng)的類型會(huì)產(chǎn)生不同數(shù)量的次級(jí)中子，每個(gè)次級(jí)中子又會(huì)產(chǎn)生新的次級(jí)中子，這樣鏈鎖反應(yīng) 下去，使得用直接模擬法模擬每一個(gè)中子是非常困難的。這種情況可以利用加權(quán)法來處理。3) 權(quán)重方法的其它應(yīng)用中子與核發(fā)生碰撞 后，產(chǎn)生的次級(jí)中子平均數(shù)為： 這里f 為裂變次級(jí)中子數(shù)。于是，碰撞后的權(quán)重為：而決定碰撞類型的幾率分別為：其中

19、加權(quán)法的思想，還可以應(yīng)用到連續(xù)分布情況和偏倚抽樣的問題tffnnnninel)3 ,()2,(32111mmwwtfftnntnntintel,3,2,)3 ,()2,(ffnnnninelt)3 ,()2,(3211 加權(quán)法雖然改進(jìn)了直接模擬法，但它同樣只關(guān)心中子是否穿透屏蔽這一信息，因此對(duì)每一個(gè)中子歷史的信息利用得很不充分。統(tǒng)計(jì)估計(jì)法能夠較多地利用中子的歷史信息，因而能得到更好的結(jié)果。一個(gè)中子，可能在介質(zhì)內(nèi)不發(fā)生碰撞而直接穿透屏蔽，也可能在介質(zhì)內(nèi)發(fā)生一次碰撞后再穿透屏蔽，或經(jīng)過二次碰撞穿透屏蔽，等等，這些事件是互不相容的，因此穿透概率p 可表示為：其中pm 是中子恰好經(jīng)過 m 次碰撞而穿透

20、屏蔽的概率。這表明，可以用求 pm (m=0,1, ) 的方法得到p。這樣，中子對(duì)穿透概率的貢獻(xiàn)就不只限于末次碰撞了。0mmpp01s0s1smmp1p0pmza0設(shè)中子的歷史為：根據(jù)該中子的歷史，我們可以估計(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 ) ，其

21、可能穿透的部分，正好是一個(gè)中子恰好經(jīng)過 m 次碰撞穿透的部分：這里的這種估計(jì)技巧，由于是對(duì)每次碰撞后的狀態(tài)，求其后未經(jīng)碰撞直接穿透的貢獻(xiàn)，因此該方法也稱為最后自由飛行估計(jì)。1, 1 , 000coscos)(expmmzaewpmmmmtmm其它，于是得到該中子對(duì)穿透的貢獻(xiàn)：如果我們共跟蹤了n個(gè)中子，則穿透率p的估計(jì)為：其方差估計(jì)為：10mmmppnnnnpnp1)3()(12)3(122)(1)(nnnpnpnp在直接模擬方法中，相對(duì)誤差為其中 為與置信水平 1相應(yīng)的量。如果構(gòu)造一個(gè)新的概率模型，使得該模型的穿透率p*與原模型的穿透率p之間存在關(guān)系：使用直接模擬方法 , 相對(duì)誤差為5. 指數(shù)變換法pnpnpppn)1 (pkp*pn 如果令*，即這意味著，達(dá)到同樣的相對(duì)誤差，跟蹤粒子的數(shù)目縮小 k 倍，從而減少 k 倍的計(jì)算量。指數(shù)變換法就是構(gòu)造一個(gè)新的概率模型的一個(gè)有效方法。knnpknpnpn*構(gòu)造如下偽過程：宏觀總截面為散射截面仍為el(e)。其中 emin、emax 分別為能量的下限和上限，為粒子的運(yùn)動(dòng)方向與 z 軸的夾角?？梢宰C明這個(gè)偽過程的穿透概率p* 與