蒙特卡羅方法在粒子輸運(yùn)模擬中的應(yīng)用

1、蒙特卡羅方法在模擬粒子輸運(yùn)問(wèn)題上的應(yīng)用蒙特卡羅模擬方法方法一直是求解科學(xué)、工程和技術(shù)領(lǐng)域大量應(yīng)用問(wèn)題的常用數(shù)值方法，廣泛應(yīng)用于眾多科學(xué)領(lǐng)域。該方法被認(rèn)為是一種成熟的計(jì)算方法，然而，我們對(duì)大量問(wèn)題的了解還是不夠的。時(shí)至今日，盡管無(wú)法確定哪些問(wèn)題運(yùn)用蒙特卡羅模擬方法方法是有效的，但我們?nèi)〉玫囊粋€(gè)共識(shí)是，對(duì)于高維和復(fù)雜幾何區(qū)域的積分問(wèn)題，蒙特卡羅模擬方法求解的有效性是公認(rèn)的一個(gè)加速蒙特卡羅模擬方法收斂速度的常用算法是使用擬隨機(jī)數(shù)(QusaiRandom Numbers，簡(jiǎn)記為QRN)，它們?cè)跇颖究臻g是高度一致均勻分布的采用偽隨機(jī)數(shù)，次抽樣可將隨機(jī)誤差下降到O(N-1)，而Q蒙特卡羅模擬方法方法可將數(shù)

2、值積分“確定地”下降到O(N.1)。在其它應(yīng)用領(lǐng)域，QRN已用于加速蒙特卡羅模擬方法方法，并為許多計(jì)算科學(xué)家研究發(fā)展蒙特卡羅模擬方法自動(dòng)降方差技巧，研究損耗、微擾和幾何級(jí)數(shù)收斂方法成為今天蒙特卡羅模擬方法研究的新方向在核科學(xué)領(lǐng)域，國(guó)內(nèi)眾多單位和高校使用蒙特卡羅模擬方法程序，推動(dòng)了國(guó)內(nèi)蒙特卡羅模擬方法方法和軟件的發(fā)展目前我國(guó)迎來(lái)了核電發(fā)展的大好時(shí)光，發(fā)展具有自主知識(shí)產(chǎn)權(quán)的蒙特卡羅模擬方法粒子輸運(yùn)程序可謂是幾代人的愿望。模擬粒子輸運(yùn)問(wèn)題也是MC最具有代表性的例子之一。粒子的輸運(yùn)問(wèn)題帶有明顯的隨機(jī)性質(zhì)，粒子的輸運(yùn)過(guò)程是一個(gè)隨機(jī)過(guò)程。粒子的運(yùn)動(dòng)規(guī)律是根據(jù)大量粒子的運(yùn)動(dòng)狀況總結(jié)出來(lái)的，是一種統(tǒng)計(jì)規(guī)律。蒙

3、特卡羅模擬，實(shí)際上就是模擬相當(dāng)數(shù)量的粒子在介質(zhì)中運(yùn)動(dòng)的狀況，使粒子運(yùn)動(dòng)的統(tǒng)計(jì)規(guī)律得以重現(xiàn)。不過(guò)，這種模擬不是用實(shí)驗(yàn)方法，而是利用數(shù)值方法和技巧，即利用隨機(jī)數(shù)實(shí)現(xiàn)的。在查閱一些資料文獻(xiàn)之后，以蒙特卡羅方法解輻射屏蔽問(wèn)題為例說(shuō)明，蒙特卡羅方法在解決模擬粒子輸運(yùn)問(wèn)題上的應(yīng)用。輻射（光子和中子）屏蔽問(wèn)題是蒙特卡羅方法最早廣泛應(yīng)用的領(lǐng)域之一。下面將從物理直觀出發(fā)，說(shuō)明蒙特卡羅方法解決這類粒子輸運(yùn)問(wèn)題的基本方法和技巧。而這些方法和技巧對(duì)于諸如輻射傳播、多次散射和通量計(jì)算等一般粒子輸運(yùn)問(wèn)題都是適用的。在反應(yīng)堆工程和輻射的測(cè)量與應(yīng)用中，常常要用一些吸收材料做成屏蔽物擋住光子或中子。我們所關(guān)心的是經(jīng)過(guò)屏蔽后射線

4、的強(qiáng)度及其能量分布，這就是屏蔽問(wèn)題。當(dāng)屏蔽物的形狀復(fù)雜，散射各向異性，材料介質(zhì)不均勻 , 核反應(yīng)截面與能量、位置有關(guān)時(shí)，難以用數(shù)值方法求解，用蒙特卡羅方法模擬卻能夠得到滿意的結(jié)果。為方便起見，選用平板屏蔽模型，在厚度為 a，長(zhǎng)、寬無(wú)限的平板左側(cè)放置一個(gè)強(qiáng)度已知，具有已知能量、方向分布的輻射源 S 。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源發(fā)出的平均一個(gè)粒子穿透屏蔽的數(shù)目。同時(shí)，假定粒子在兩次碰撞之間按直線運(yùn)動(dòng) , 且粒子之間的相互作用可以忽略。采用直接從物理問(wèn)題出發(fā)，模擬粒子的真實(shí)物理過(guò)程的方法。粒子在介質(zhì)中的運(yùn)動(dòng)的狀態(tài)，可用一組參數(shù)來(lái)描述，稱之為狀態(tài)參數(shù)。它通常包括：粒

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

6、態(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ù)。一個(gè)由源發(fā)出的粒子在介質(zhì)中運(yùn)動(dòng)，經(jīng)過(guò)若干次碰撞后，直到其運(yùn)動(dòng)歷史結(jié)束（如逃出系統(tǒng)或被吸收等）。假定粒子在兩次碰撞之間按直線運(yùn)動(dòng)，其運(yùn)動(dòng)方向與能量均不改變，則粒子在介質(zhì)中的運(yùn)動(dòng)過(guò)程可用以下碰撞點(diǎn)的狀態(tài)序列 描述： S0 ，S1 ，SM-1 ，SM或者更詳細(xì)些 , 用來(lái)描述。這里 S0 為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，SM 為粒子的終止?fàn)顟B(tài)。M 稱為粒子運(yùn)動(dòng)的鏈長(zhǎng)

7、。這樣的序列稱為粒子隨機(jī)運(yùn)動(dòng)的歷史，模擬一個(gè)粒子的運(yùn)動(dòng)過(guò)程，就變成確定狀態(tài)序列的問(wèn)題。接著，我們模擬運(yùn)動(dòng)過(guò)程。為簡(jiǎn)單起見，這里以中子穿透均勻平板的模型來(lái)說(shuō)明，這時(shí)狀態(tài)參數(shù) 取 S( z , E , cos)。模擬的步驟如下：(1) 確定初始狀態(tài) S0 ：確定粒子的初始狀態(tài)，實(shí)際上就是要從中子源的空間位置、能量和方向分布中抽樣。設(shè)源分布為則分別從各自的分布中抽樣確定初始狀態(tài)。對(duì)于平板情況，抽樣得到 z00。確定下一個(gè)碰撞點(diǎn) ：已知狀態(tài)Sm-1，要確定狀態(tài)Sm，首先要確定下一個(gè)碰撞點(diǎn)的位置 zm。在相鄰兩次碰撞之間，中子的輸運(yùn)長(zhǎng)度 l 服從如下分布：對(duì)于平板模型，l 服從分布：其中，t 為介質(zhì)的中

8、子宏觀總截面，積分 稱為粒子輸運(yùn)的自由程數(shù)，系統(tǒng)的大小通常就是用系統(tǒng)的自由程數(shù)表示的。顯然，粒子輸運(yùn)的自由程數(shù)服從指數(shù)分布，因此從 f ( l ) 中抽樣確定 l，就是要從積分方程 中解出 l。對(duì)于單一介質(zhì)則下一個(gè)碰撞點(diǎn)的位置如果 zma，則中子穿透屏蔽，若 zm0, 則中子被反射出屏蔽。這兩種情況，均視為中子歷史終止。確定被碰撞的原子核 ：通常介質(zhì)由幾種原子核組成，中子與核碰撞時(shí)，要確定與哪一種核碰撞。設(shè)介質(zhì)由A、B、C 三種原子核組成，其核密度分別為NA、NB、NC，則介質(zhì)的宏觀總截面為：其中 分別為核A、B、C 的宏觀總截面。其定義如下：分別表示(·)核的宏觀總截面、核密度和微

9、觀總截面。由于中子截面表示中子與核碰撞可能性的大小，因此，很自然地，中子與A、B、C 核發(fā)生碰撞的幾率分別為：利用離散型隨機(jī)變量的抽樣方法，確定碰撞核種類：確定碰撞類型 ：確定了碰撞的核(比如B核)后，就要進(jìn)一步確定碰撞類型。中子與核的反應(yīng)類型有彈性散射、非彈性散射、(n,2n)反應(yīng)，裂變和俘獲等，它們的微觀截面分別為則有各種反應(yīng)發(fā)生的幾率分別為利用離散型隨機(jī)變量的抽樣方法，確定反應(yīng)類型。在屏蔽問(wèn)題中，中子與核反應(yīng)常只有彈性散射和吸收兩種類型，吸收截面為：這時(shí)，總截面為：發(fā)生彈性散射的幾率為：若 ，則為彈性散射；否則為吸收，發(fā)生吸收反應(yīng)意味著中子的歷史終止。如果中子被碰撞核吸收，則其輸運(yùn)歷史結(jié)

10、束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運(yùn)動(dòng)方向。中子能量 Em 為：A是碰撞核的質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量；C 為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角。 可從質(zhì)心系中彈性散射角分布fC(C) 中抽樣產(chǎn)生。實(shí)驗(yàn)室系散射角L的余弦L為：如果給出實(shí)驗(yàn)室系散射角余弦分布 fL(L)，可直接從 fL(L)中抽取L，此時(shí)能量Em與L的關(guān)系式為：確定了實(shí)驗(yàn)室系散射角L后，再使用球面三角公式確定cosm ：其中為在0,2上均勻分布的方位角。至此，由Sm-1完全可以確定Sm。 因此，當(dāng)中子由源出發(fā)后，即S0確定后，重復(fù)步驟 (2)(5)，直到中子游動(dòng)歷史終止。于是得到了一個(gè)中子的隨

11、機(jī)游動(dòng)歷史 S0 ，S1 ，SM-1 ，SM，即也就是模擬了一個(gè)由源發(fā)出的中子的運(yùn)動(dòng)過(guò)程。以上模擬過(guò)程可分為兩大步：第一步確定粒子的初始狀態(tài)S0，第二步由狀態(tài)Sm-1來(lái)確定狀態(tài)Sm。這第二步又分為兩個(gè)過(guò)程：第一個(gè)過(guò)程是確定碰撞點(diǎn)位置zm ，稱為輸運(yùn)過(guò)程；第二個(gè)過(guò)程是確定碰撞后粒子的能量及運(yùn)動(dòng)方向，稱為碰撞過(guò)程。對(duì)于中子而言，碰撞過(guò)程是先確定散射角，進(jìn)而確定能量和運(yùn)動(dòng)方向；而對(duì)于光子，碰撞過(guò)程是先確定能量，再確定散射角以及運(yùn)動(dòng)方向。重復(fù)這兩個(gè)過(guò)程，直至粒子的歷史終止。這種模擬過(guò)程，是解任何類型的粒子輸運(yùn)問(wèn)題所共有的，它是蒙特卡羅方法解題的基本手段。參考文獻(xiàn)：謝仲生，鄧力中子輸運(yùn)理論數(shù)值計(jì)算方法M西北工業(yè)大學(xué)出版社，2004：3640上官丹驊，許海燕，鄧力基于Monte Carlo方法的中子標(biāo)識(shí)輸運(yùn)計(jì)算J清華大學(xué)學(xué)報(bào)，2007，47(S1)：10621064鄧力，張文勇MCNP程序在MPI下的并行化及完善J數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，2003，24(3)：161166竹生東，鄧力，李樹，等堆外核儀表系統(tǒng)(RPN)的預(yù)設(shè)效驗(yàn)系數(shù)理論計(jì)算J核動(dòng)力工程，2004，25(2)：152155鄧力，謝仲生碳氧比能譜測(cè)井的蒙特卡羅模擬J地