蒙特卡洛方法在中子輸運中的應(yīng)用

1、中子輸運理論與數(shù)值方法課程作業(yè)蒙特卡洛方法目錄1. 前言 32. 蒙特卡洛方法概述 32.1 蒙特卡洛方法的基本思想 42.2 蒙特卡洛方法的收斂性、誤差 42.2.1 蒙特卡洛方法的收斂性 42.2.2 蒙特卡洛方法的誤差 52.3 蒙特卡洛方法的特點 62.4 蒙特卡洛方法的主要應(yīng)用范圍 73. 隨機數(shù) 73.1 線性乘同余方法 93.2 偽隨機數(shù)序列的均勻性和獨立性 93.2.1 偽隨機數(shù)的均勻性 93.2.2 偽隨機數(shù)的獨立性 104. 蒙特卡洛方法在粒子輸運上的應(yīng)用 104.1 屏蔽問題模型 104.2 直接模擬方法 114.2.1 狀態(tài)參數(shù)與狀態(tài)序列 114.2.2 模擬運動過程

2、124.2.3 記錄結(jié)果 154.3 蒙特卡洛方法的效率 165. 蒙特卡洛方法應(yīng)用程序 MCNP 175.1 MCNP 簡述 175.2 MCNP誤差的估計 185.3 MCNP效 率因素 196. 結(jié)論 19參考文獻 201. 前言半個多世紀(jì)以來，由于科學(xué)技術(shù)的發(fā)展和電子計算機的發(fā)明，蒙特卡洛( Monte Carlo ) 方法作為一種獨立的方法被提出來，并首先在核武器的試驗與研 制中得到了應(yīng)用。 蒙特卡洛方法是一種計算方法， 但與一般數(shù)值計算方法有很大 區(qū)別。它是以概率統(tǒng)計理論為基礎(chǔ)的一種方法。 由于蒙特卡洛方法能夠比較逼真 地描述事物的特點及物理實驗過程， 解決一些數(shù)值方法難以解決的問

3、題， 因而該 方法的應(yīng)用領(lǐng)域日趨廣泛。 蒙特卡洛模擬計算是解決中子在介質(zhì)中輸運較為成熟、 有效的方法，對于原子能、 輻射防護、劑量學(xué)和輻射生物物理學(xué)等研究領(lǐng)域?qū)嶋H 問題的計算，都可以利用蒙特卡洛方法予以實現(xiàn)。粒子輸運過程可以用玻耳茲曼方程加以描述，然而，以此基礎(chǔ)上發(fā)展起來的 近似數(shù)值方法如擴散近似法、 離散坐標(biāo)方法在處理截面與能量相關(guān)以及散射各向 異性介質(zhì)、 復(fù)雜幾何條件問題時碰到了較大困難。 而蒙特卡洛方法在處理這類問 題時得心應(yīng)手，有很強的解題能力，并且近似較少，接近于真實情況。粒子輻射問題計算通常有輸運方程法、 蒙特卡洛法(MC法)、實驗測量法以及 經(jīng)驗法等幾種方法。 蒙特卡洛計算法又稱

4、隨機抽樣法或統(tǒng)計試驗法， 是基于計算 機模擬的思想， 抓住物理過程的數(shù)量和幾何特征， 進行數(shù)字模擬試驗， 該方法是 求解輻射輸運問題的一種相當(dāng)成熟和有效的方法， 而且它對于各種復(fù)雜問題， 具 有良好的通用性， 實用性相當(dāng)廣泛， 幾乎涉及核科學(xué)的各個領(lǐng)域。 本文主要介紹 蒙特卡洛的概念、原理和應(yīng)用及研究現(xiàn)狀。2. 蒙特卡洛方法概述 蒙特卡洛方法又稱隨機抽樣技巧或統(tǒng)計試驗方法。半個多世紀(jì)以來，由于科 學(xué)技術(shù)的發(fā)展和電子計算機的發(fā)明 ，這種方法作為一種獨立的方法被提出來， 并首先在核武器的試驗與研制中得到了應(yīng)用。 蒙特卡洛方法是一種計算方法， 但 與一般數(shù)值計算方法有很大區(qū)別。 它是以概率統(tǒng)計理論為

5、基礎(chǔ)的一種方法。 由于 蒙特卡洛方法能夠比較逼真地描述事物的特點及物理實驗過程， 解決一些數(shù)值方 法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。蒙特卡洛方法的主要組成部分有：(1) 概率密度函數(shù)(pdf)必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù)；(2) 隨機數(shù)產(chǎn)生器一能夠產(chǎn)生在區(qū)間0,1上均勻分布的隨機數(shù)；(3) 抽樣規(guī)則一如何從在區(qū)間0,1上均勻分布的隨機數(shù)出發(fā)，隨機抽取服從給定的pdf的隨機變量；(4) 模擬結(jié)果記錄一記錄一些感興趣的量的模擬結(jié)果；(5) 誤差估計一必須確定統(tǒng)計誤差(或方差)隨模擬次數(shù)以及其它一些量的變化；(6) 減少方差的技術(shù)一利用該技術(shù)可減少模擬過程中計算的次數(shù)；(

6、7) 并行和矢量化一可以在先進的并行計算機上運行的有效算法2.1蒙特卡洛方法的基本思想可以通俗地說，蒙特卡洛方法是用隨機試驗的方法計算積分，即將所要計算 的積分看作服從某種分布密度函數(shù)f(r)的隨機變量g (r)的數(shù)學(xué)期望:g* ：.i0 g(r) f (r)dr(0.1)通過某種試驗，得到N個觀察值"，2,，rN (用概率語言來說，從分布密度 函數(shù)f(r)中抽取N個子樣1,2, ,“,)，將相應(yīng)的N個隨機變量的值g(r 1)，1 Ng(r 2)，g(r n)的算術(shù)平均值gN - ' g(ri)，作為積分的估計值(近似值)。N 7為了得到具有一定精確度的近似解，所需試驗的次數(shù)

7、是很多的，通過人工方 法作大量的試驗相當(dāng)困難，甚至是不可能的。因此，蒙特卡洛方法的基本思想雖 然早已被人們提出，卻很少被使用。本世紀(jì)四十年代以來，由于電子計算機的出 現(xiàn)，使得人們可以通過電子計算機來模擬隨機試驗過程，把巨大數(shù)目的隨機試驗交由計算機完成，使得蒙特卡洛方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中發(fā) 揮應(yīng)有的作用。2.2蒙特卡洛方法的收斂性、誤差蒙特卡洛方法作為一種計算方法，其收斂性與誤差是普遍關(guān)心的一個重要問 題。2.2.1蒙特卡洛方法的收斂性由前面介紹可知，蒙特卡洛方法是由隨機變量X的簡單子樣X1, X2,，Xn 1 N的算術(shù)平均值XN 、 xi .作為所求解的近似值。由大數(shù)定律可知

8、，如,人，,N iXn獨立同分布，且具有有限期望值，貝U一 Xm:P二 E(X)=1。即隨機變量X的簡單子樣的算術(shù)平均值Xn，當(dāng)子樣數(shù)N充分大時，以概率1收斂于它的期望值E(X)。222蒙特卡洛方法的誤差蒙特卡洛方法的近似值與真值的誤差問題， 概率論的中心極限定理給出了答 案。該定理指出，如果隨機變量序列 Xi，%,凡獨立同分布，且具有有限非 零的方差(T 2，即0乂2二(x-E(X)2 f (x)dx ： ： :。f(X)是X的分布密度函數(shù)。lNmiv廠対X e/2dt(0.2)當(dāng)N充分大時，有如下的近似式P(xn-E(X)c 倉卜話 t2/2 ,o e d1-：(0.3)n E(X) &l

9、t;其中a稱為置信度，1-a稱為置信水平。這表明，不等式X 近似地以概率1-a成立，且誤差收斂速度的階為 O(N/2)。通常，蒙特卡洛方法的誤差&定義為(0.4)上式中I與置信度a是一一對應(yīng)的，根據(jù)問題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出:，常用的a與I的對應(yīng)關(guān)系為：a =0.5，：.=0.6745 ; a =0.05， . =0.96 ; a =0.003， -.=3.蒙特卡洛方法的誤差為概率誤差,這與其他數(shù)值計算方法是有區(qū)別的。 誤差中的均方差c是未知的，必須使用 其估計值nnNiN' x；(- Xi)2(0.5)V N 7N 7來代替以求出均方差c。由式

10、(0.4)可知當(dāng)給定置信度a后，誤差&由c和N決定。要減小&，或者是增大N,或者是減小方差c 2。在c固定的情況下，要把精 度提高一個數(shù)量級，試驗次數(shù)N需增加兩個數(shù)量級。因此，單純增大N不是一個 有效的辦法。另一方面，如能減小估計的均方差c,比如降低一半，那誤差就減 小一半，這相當(dāng)于N增大四倍的效果。因此降低方差的各種技巧，引起了人們的 普遍注意。2.3蒙特卡洛方法的特點作為一種統(tǒng)計試驗方法，蒙特卡洛方法因其優(yōu)點在諸多領(lǐng)域內(nèi)有著廣泛，但 同時存在一些缺點。蒙特卡洛的主要優(yōu)點有：(1) 能夠比較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程。蒙特 卡洛方法可以部分代替物理實驗，

11、甚至可以得到物理實驗難以得到的結(jié)果。用蒙特卡洛方法解決實際問題，可以直接從實際問題本身出發(fā)，而不從方程或數(shù)學(xué)表 達式出發(fā)。它有直觀、形象的特點。 受幾何條件限制小。在 計算s維空間中的 任一區(qū)域DS上的積分 gg(xi ,x2 ,川兀dxdx2lldxs時，無論區(qū)域D的形狀多么特殊，只要能給出Ds描述D的幾何特征的條件，就可以從D中均勻產(chǎn)生N個點(xi(i),x2i)JH,xSi)，得到D N積分的近似值8n =邑遲g(xi(i),x2i)l,xSi)，其中D為區(qū)域DS的體積。這是數(shù)值 N y方法難以作到的。(3) 收斂速度與問題的維數(shù)無關(guān)。由誤差定義可知，在給定置信水平情況下，蒙特卡洛方法的

12、收斂速度為0(N，/2)，與問題本身的維數(shù)無關(guān)。維數(shù)的變化，只 引起抽樣時間及估計量計算時間的變化，不影響誤差。也就是說，使用蒙特卡洛方法時，抽取的子樣總數(shù)N與維數(shù)s無關(guān)。維數(shù)的增加，除了增加相應(yīng)的計算量 外，不影響問題的誤差。這一特點，決定了蒙特卡洛方法對多維問題的適應(yīng)性。(4) 具有同時計算多個方案與多個未知量的能力。對于那些需要計算多個方案的問題，使用蒙特卡洛方法有時不需要像常規(guī)方法那樣逐個計算，而可以同時 計算所有的方案，其全部計算量幾乎與計算一個方案的計算量相當(dāng)。 例如，對于 屏蔽層為均勻介質(zhì)的平板幾何，要計算若干種厚度的穿透概率時，只需計算最厚 的一種情況，其他厚度的穿透概率在計算

13、最厚一種情況時稍加處理便可同時得到。(5) 誤差容易確定。對于一般計算方法，要給出計算結(jié)果與真值的誤差并不 是一件容易的事情，而蒙特卡洛方法則不然。根據(jù)蒙特卡洛方法的誤差公式，可 以在計算所求量的同時計算出誤差。 對干很復(fù)雜的蒙特卡洛方法計算問題， 也是 容易確定的。(6) 程序結(jié)構(gòu)簡單，易于實現(xiàn)。在計算機上進行蒙特卡洛方法計算時，程序 結(jié)構(gòu)簡單，分塊性強，易于實現(xiàn)。蒙特卡洛的主要缺點有：(1) 收斂速度慢。如前所述，蒙特卡洛方法的收斂速度為0(N/2)，般不容易得到精確度較高的近似結(jié)果。對于維數(shù)少(三維以下)的問題，不如其他方 法好。(2) 誤差具有概率性。由于蒙特卡洛方法的誤差是在一定置信

14、水平下估計的， 所以它的誤差具有概率性，而不是一般意義下的誤差。(3) 在粒子輸運問題中，計算結(jié)果與系統(tǒng)大小有關(guān)。經(jīng)驗表明，只有當(dāng)系統(tǒng)的大小與粒子的平均自由程可以相比較時(一般在十個平均自由程左右)，蒙特卡洛方法計算的結(jié)果較為滿意。但對于大系統(tǒng)或小概率事件的計算問題，計算結(jié)果往往比真值偏低。而對于大系統(tǒng)，數(shù)值方法則是適用的。因此，在使用蒙特卡洛方法時，可以考慮把蒙特卡洛方法與解析(或數(shù)值)方法相結(jié)合，取長補短。2.4蒙特卡洛方法的主要應(yīng)用范圍蒙特卡洛方法所特有的優(yōu)點，使得它的應(yīng)用范圍越來越廣。它的主要應(yīng)用范 圍包括：粒子輸運問題，統(tǒng)計物理，典型數(shù)學(xué)問題，真空技術(shù)，激光技術(shù)以及醫(yī) 學(xué)，生物，探礦

15、等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。蒙特卡洛方法在粒子輸運問題中的應(yīng)用范圍主要包括：實驗核物理，反應(yīng)堆 物理，高能物理等方面。蒙特卡洛方法在實驗核物理中的應(yīng)用范圍主要包括：通量及反應(yīng)率，中子探測效率，光子探測效率，光子能量沉積譜及響應(yīng)函數(shù)，氣體 正比計數(shù)管反沖質(zhì)子譜，多次散射與通量衰減修正等方面。3. 隨機數(shù)隨機數(shù)是蒙特卡洛方法的主要組成部分之一。隨機數(shù)是指一個數(shù)列，其中的 每一個體稱為隨機數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)。在一個均勻分布的隨機數(shù)中， 每一個體出現(xiàn)的概率是均等的。物理中的很多過程需要隨機數(shù)確定，比如出射粒 子的能量、方向等屬性，粒子與介質(zhì)的相互作用等等。所模擬的物理過

16、程要求隨 機數(shù)應(yīng)具有下列特性：1. 隨機數(shù)序列應(yīng)是獨立的、互不相關(guān)的(un correlated):即序列中的任一 子序列應(yīng)與其它的子序列無關(guān)；2. 長的周期(Iong period):實際應(yīng)用中，隨機數(shù)都是用數(shù)學(xué)方法計算出來的，這些算法具有周期性，即當(dāng)序列達到一定長度后會重復(fù)；3. 均勻分布的隨機數(shù)應(yīng)滿足均勻性(Un iformity):隨機數(shù)序列應(yīng)是均勻的、 無偏的，即：如果兩個子區(qū)間的“面積”相等，則落于這兩個子區(qū)間內(nèi) 的隨機數(shù)的個數(shù)應(yīng)相等。4. 有效性(Efficiency):模擬結(jié)果可靠，隨機數(shù)的產(chǎn)生必須快速、有效，最好能夠進行并行計算。為了產(chǎn)生隨機數(shù)，可以使用隨機數(shù)表。隨機數(shù)表是由

17、0, 1,，9十個數(shù)字組成，每個數(shù)字以0.1的等概率出現(xiàn)，數(shù)字之間相互獨立。這些數(shù)字序列叫作隨 機數(shù)字序列。如果要得到n位有效數(shù)字的隨機數(shù)，只需將表中每n個相鄰的隨機數(shù) 字合并在一起，且在最高位的前邊加上小數(shù)點即可。例如，某隨機數(shù)表的第一行數(shù)字為7634258910，要想得到三位有效數(shù)字的隨機數(shù)依次為0.763，0.425，0.891?？梢允褂梦锢矸椒óa(chǎn)生隨機數(shù)，用來作為隨機數(shù)發(fā)生器的物理源主要有兩種： 一種是根據(jù)放射性物質(zhì)的放射性，另一種是利用計算機的固有噪聲。但在計算機 上產(chǎn)生隨機數(shù)最實用、最常見的方法是數(shù)學(xué)方法，即用如下遞推公式：nk 訂(n，川，n2)，門=1,2, “ |(3.1)產(chǎn)

18、生隨機數(shù)序列。對于給定的初始值。E 1, E 2，E k,確定E n+k，n =1, 2。經(jīng)常使用的是k=1的情況。a)用數(shù)學(xué)方法產(chǎn)生的隨機數(shù)有兩個特點，即：遞推公式和初始值E1, E 2,E k確定后，整個隨機數(shù)序列便被唯一確定。不滿足隨機數(shù)相互獨立的要 求。b)由于隨機數(shù)序列是由遞推公式確定的，而在計算機上所能表示的0,1上的數(shù)又是有限的，因此，這種方法產(chǎn)生的隨機數(shù)序列就不可能不出現(xiàn)無限 重復(fù)。一旦出現(xiàn)這樣的n，, n (n ' <n)，使得嚴(yán)叮i=1,2,|H，k成 立隨機數(shù)序列便出現(xiàn)了周期性的循環(huán)現(xiàn)象。對于k=1的情況，只要有一個隨機數(shù)重復(fù)，其后面的隨機數(shù)全部重復(fù)，這與隨

19、機數(shù)的要求是不相符的。由于這兩個問題的存在，常稱用數(shù)學(xué)方法產(chǎn)生的隨機數(shù)為偽隨機數(shù)。3.1線性乘同余方法線性乘同余方法是由Lehme在1951年提出來的，是一種最常用的產(chǎn)生偽隨機數(shù)的方法。乘同余方法中采用的遞推公式為In 1 二(aIn c)mod m(3.2)其中：10為初始值，a為乘法器，c為增值，仃為模數(shù)，mo(為取模運算，(aIn c)除以m后的余數(shù)。a、c、仃皆為整數(shù)。實型隨機序列：Infloat (m)0,1), In ： m(3.3)float (m -1)0,1),In 5-1(3.4)上式中，獨立性和均勻性取決于參數(shù)a和c的選擇。m應(yīng)盡可能地大，因為序列的周期不可能大于m通常將

20、m取為計算機所能表示的最大的整型量，在 32位機上,319m =22 10。 1961年，M. Greenberger證明，用線性乘同余方法產(chǎn)生的隨機數(shù)序列具有周期m的條件是：(1)c和n互為質(zhì)數(shù)；(2)a-1是質(zhì)數(shù)p的倍數(shù)，其中p 是a-1和m的公約數(shù)；(3)如果n是4的倍數(shù)，a-1也是4的倍數(shù)。3.2偽隨機數(shù)序列的均勻性和獨立性3.2.1偽隨機數(shù)的均勻性這里只考慮偽隨機數(shù)序列E 1, E 2，E n全體作為子樣時的均勻性問題。其中 n為偽隨機數(shù)序列的最大容量。對于任意的 0w x< 1，令M(x)表示偽隨機數(shù)序列 E 1, E 2，E n中適合不等式E i< x,i=1,2,n

21、的個數(shù)，則(3.5)將偽隨機數(shù)序列E 1, E 2，E n從小至大重新排列1七2 一川一 n ,令0 = 0, n 1 -1，則由S (n)的定義，容易證明、(n) = max I i丄l，l i J 丄I，很II。蜜 In nJ明顯，對于固定的n, S (n)的值越小越好。它是描述偽隨機數(shù)序列均勻程度的一 11基本量。對于任意隨機數(shù)序列，均有不等式、：(n) _ 成立。當(dāng)；(n)成立時， 2n2n所對應(yīng)的偽隨機數(shù)序列為最佳分布。3.2.2偽隨機數(shù)的獨立性對于任意 0 一 X, y 一1，令 NnX，y)表示(E 1, E 2),( E 2, E 3),(E n, E n+1)中適合不等式i

22、:：: x, i 1 :：: y。的個數(shù)，根據(jù)隨機變量間相互獨立的定義和頻率近似概率的方法，令(n)-sup10 空,y 21Nn(x,y)nNn(x) Nn(y)|n n(3.6)則& (n)標(biāo)志偽隨機數(shù)序列E 1, E 2，E n的獨立程度，簡稱為獨立偏度。對于固 定的n, & (n)的值越接近于零，偽隨機數(shù)序列的獨立性越好。4. 蒙特卡洛方法在粒子輸運上的應(yīng)用輻射(光子和中子)屏蔽問題是蒙特卡洛方法最早廣泛應(yīng)用的領(lǐng)域之一。現(xiàn)主要從物理直觀出發(fā)，說明蒙特卡洛方法解決這類粒子輸運問題的基本方法和技 巧。解決屏蔽問題時可采取多種方法，如直接模擬方法、簡單加權(quán)法、統(tǒng)計估計 法、指

23、數(shù)變換法等，這里只對直接模擬方法做介紹。4.1屏蔽問題模型在反應(yīng)堆工程和輻射的測量與應(yīng)用中，常常要用一些吸收材料做成屏蔽物擋 住光子或中子。我們所關(guān)心的是經(jīng)過屏蔽后射線的強度及其能量分布，這就是屏蔽問題。當(dāng)屏蔽物的形狀復(fù)雜，散射各向異性，材料介質(zhì)不均勻，核反應(yīng)截面與能量、位置有關(guān)時，難以用數(shù)值方法求解，用蒙特卡洛方法能夠得到滿意的結(jié) 果。粒子的輸運問題帶有明顯的隨機性質(zhì)， 粒子的輸運過程是一個隨機過程。 粒 子的運動規(guī)律是根據(jù)大量粒子的運動狀況總結(jié)出來的，是一種統(tǒng)計規(guī)律。蒙特卡洛模擬，實際上就是模擬相當(dāng)數(shù)量的粒子在介質(zhì)中運動的狀況，使粒子運動的統(tǒng) 計規(guī)律得以重現(xiàn)。不過，這種模擬不是用實驗方法，

24、而是利用數(shù)值方法和技巧， 即利用隨機數(shù)來實現(xiàn)的。為方便起見，選用平板屏蔽模型，在厚度為 a,長、寬無限的平板左側(cè)放置 一個強度已知，具有已知能量、方向分布的輻射源 S,見圖4.1。求粒子穿透屏 蔽概率(穿透率)及其能量、方向分布。穿透率就是由源發(fā)出的平均一個粒子穿 透屏蔽的數(shù)目。同時，假定粒子在兩次碰撞之間按直線運動，且粒子之間的相互作用可以忽略。S7屏蔽物圖4.1屏蔽問題模型4.2直接模擬方法直接模擬方法就是直接從物理問題出發(fā)，模擬粒子的真實物理過程。4.2.1狀態(tài)參數(shù)與狀態(tài)序列粒子在介質(zhì)中的運動的狀態(tài)，可用一組參數(shù)來描述，稱之為狀態(tài)參數(shù)。它通 常包括：粒子的空間位置r，能量E和運動方向Q,

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

26、經(jīng)過m次碰撞后的狀態(tài)，其中rm：粒子在第 m次碰撞點的位置Em：粒子第m次碰撞后的能量Q m粒子第m次碰撞后的運動方向tm：粒子到第 m次碰撞時所經(jīng)歷的時間wm：粒子第m次碰撞后的權(quán)重一個由源發(fā)出的粒子在介質(zhì)中運動，經(jīng)過若干次碰撞后，直到其運動歷史結(jié) 束(如逃出系統(tǒng)或被吸收等)。假定粒子在兩次碰撞之間按直線運動，其運動方 向與能量均不改變，則粒子在介質(zhì)中的運動過程可用以下碰撞點的狀態(tài)序列描述,即 SO, S1,，SM-1, SM或'ro ,斤，川，rMrME。，E1，HI, Eg, EmI Qo,甜，III,©M _1，Qm 來描述。這里So為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，S

27、m為粒子的終止?fàn)顟B(tài)。M稱 為粒子運動的鏈長。422模擬運動過程這里以中子穿透均勻平板的模型來說明，這時狀態(tài)參數(shù)取S= (z,E,cos a )。模擬的步驟如下。(1) 確定初始狀態(tài)So：確定粒子的初始狀態(tài)，實際上就是要從中子源的空間位置、能量和方向分布 中抽樣。設(shè)源分布為f(Zo,Eo,COS： o) = fi(Zo)f2(Eo)f3(COS： o)，則分別從各自的分 布中抽樣確定初始狀態(tài)。(2) 確定下一個碰撞點：已知狀態(tài)Sni，要確定狀態(tài)Sm首先要確定下一個碰撞點的位置Zmo在相鄰兩次 碰撞之間，中子的輸運長度I服從如下分布：f(l)Ft(rmI Qm4,Em4)expon(瞎T 編亠Em

28、j)dl(4.1)對于平板模型，l服從分布：f(l)=打(Zm4 l COS： m-Em4)expj.o 1t(Zm4 l COS ： m J, Em/)dl 】 件2) 其中，工t為介質(zhì)的中子宏觀總截面。積分 J；£(rm+r 'Qm_l, Em_l)dl *稱為粒子輸 運的自由程數(shù)。顯然，粒子輸運的自由程數(shù)服從指數(shù)分布，因此從f(l)中抽樣確l定l，就是要從積分方程,oit(rmd l Qmd,Emd)d-ln中解出l。對于單一介質(zhì)丨 也，則下一個碰撞點的位置為Et（EmJzm = zm J '丨 COS'm J = zm JIn匕偏二）cos<mJ

29、(4.3)如果Zm> a,則中子穿透屏蔽，若Zm< 0,則中子被反射出屏蔽。這兩種情況，均 視為中子歷史終止。（3）確定被碰撞的原子核：通常介質(zhì)由幾種原子核組成，中子與核碰撞時，要確定與哪一種核碰撞。設(shè) 介質(zhì)由A、B、C三種原子核組成，其核密度分別為 N、NB、NC,則介質(zhì)的宏觀總 截面為：1t(Em4 亍(Em J 寸(Em)叮(Em)(4.4)其中空叮，寸分別為核A、B、C的宏觀總截面。其定義如下：即（Em）=N（）V>（EmJ，U）（Em）、N（）.、二（Em）分別表示（）核的宏觀總截面、核密度和微觀總截面。由于中子截面表示中子與核碰撞可能性的大小， 因此，很自然地，中

30、子與A、B、C核發(fā)生碰撞的幾率分別為：Pa =工/(Em)PRlt(Em4/B亍(Em)lt(Em4)卩。=?。‥m） C 3t（EmJ(4.5)若飛Pa ,則中子與A核碰撞；若Pa Pb，則中子與B核碰撞；若Pa Pb，則與C核碰撞。（4）確定碰撞類型：確定了碰撞的核（比如B核）后，就要進一步確定碰撞類型。中子與核的反應(yīng)類型有彈性散射、非彈性散射、（n,2n）反應(yīng)，裂變和俘獲等，它們的微觀截面分 別為匚：（Em、二iBgJ、強,2n）（EmJ、匚和匚，則有屮屋皿“二誌乍皿"V（Emd）飛爲(wèi)n）（EmJ V（EmJ V（EmJ（4.6）各種反應(yīng)發(fā)生的幾率分別為Pnpn,2n)PfFC

31、利用離散型隨機變量的抽樣方法,el （Em）.；t （ Em 1）in （Em_!）'-t （ Em 1） 眾2n）（EmQYtB（Em）；-f （Em）.；t （ Em J ）-'c （Em/）：；-t （Em）確定反應(yīng)類型。在屏蔽問題中，中子與核反應(yīng)(4.7)常只有彈性散射和吸收兩種類型,吸收截面為：二（Em）*；©）；叮耳）。這時，總截面為:發(fā)生彈性散射的幾率為：巳若乞Pel，則為彈性散射；否則為吸收,6 （Em）(4.8)5B(EmJJB(EmJ y：(Em)發(fā)生吸收反應(yīng)意味著中子的歷史終止。（5）確定碰撞后的能量與運動方向：如果中子被碰撞核吸收，則其輸運歷史

32、結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運動方向。中子能量 Em為：Em 二 E21 r 1-r cosc（4.9）其中=（AW）2。A是碰撞核的質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量； A+10 C為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角。匕二cosc可從質(zhì)心系中彈性散射角分布fc（卩C）中抽樣產(chǎn)生。實驗室系散射角0 L的余弦卩L為:二 1+A 巴L _ 1 A2 2Ac如果給出實驗室系散射角余弦分布f L（卩l(xiāng)），可直接從f l（卩l(xiāng)）中抽取卩L，此時能量E與卩L的關(guān)系式為：Em二占京 ' A2 -1 叮2。確定了實驗室系散射 角0 L后，再使用球面三角公式確定co

33、s a m :cos m 二 cos m4 cosl sin ： m4 sinL cos 。各角度關(guān)系如圖4.2所示圖4.2角度關(guān)系示意圖至此，由缶完全可以確定S。因此，當(dāng)中子由源出發(fā)后，即S。確定后，重復(fù)步驟 ，直到中子游動歷史終止。于是得到了一個中子的隨機游動歷史 $,Smi， Sm,即卩'Z0,Z1,III,ZM 1,ZME。，E1,III,Em丄Emicos%cos%,III,cos m 1,coso(m也就是模擬了一個由源發(fā)出的中子的運動過程。4.2.3記錄結(jié)果在獲得中子的隨機游動歷史后，我們要對所要計算的物理量進行估計。對于 屏蔽問題，我們要計算中子的穿透率?？疾烀總€中子的

34、隨機游動歷史，它可能穿 透屏蔽(Zm> a)，可能被屏蔽發(fā)射回來(Zm< 0)，或者被吸收。設(shè)第n個中子對穿 透的貢獻為n n,貝U0,當(dāng)Zm -a當(dāng)Zm空0N如果我們共跟蹤了 N個中子，則穿透屏蔽的中子數(shù)為：n。則穿透屏蔽n =1概率的近似值為：P二叫n(4.10)N N n ¥我們稱這種直觀地模擬過程和估計方法為直接模擬方法。在置信水平1-a=0.95時，F(xiàn)?(1)的誤差為：工_1.96旳N N(4.11)其中匚為n n的均方差，由于n n是一個服從二項分布的隨機變量，所以(4.12)宀P（仁P）P用1）（仁即）(4.13)為得到中子穿透屏蔽的能量、角分布，將能量、角

35、度范圍分成若干個間隔:Emin= E| ： 111 : Ei ；： Eq Emax ，。0 ：1 ：丨 11 ": 2。其中 Emax， Enin分別表示能量的上、下限，對于穿透屏蔽的中子按其能量、方向分間隔記錄。設(shè)一穿透屏 蔽的中子能量為曰，其運動方向與Z軸夾角為a m,若能量E屬于第i個能量間隔 E，角度a m屬于第j個角度間隔 a j，則分別在第i個能量計數(shù)器及第j 個角度計數(shù)器中加1。跟蹤N個中子后，則i =1,2川 |,1N2,j(4.14)Nj =1,2，Hl,J(4.15)分別為穿透中子的能量分布和角分布。其中N,i和分別為第i個能量和第個角度間隔的穿透中子數(shù)。歸一后分

36、別為:? 1)* N“N1,iPN1,i =ir = Nnii(1)*= 1,2，川，1(4.16)N2,jj =1,2,|,J(4.17)4.3蒙特卡洛方法的效率衡量蒙特卡洛技巧的好壞，除了看其方差大小外，還要看其所需費用（計算時間）多少，即從該技巧的效率 E （方差與費用乘積的倒數(shù)）全面考慮：1Ef二飛（4.佝CT T其中/為方差，T為所需費用。Ef大時，所用方法的效率高；否則，效率低在一般情況下，直接模擬方法、簡單加權(quán)法、統(tǒng)計估計法、指數(shù)變換法等方法中 有些方法雖然減小了方差，卻增加了費用。例如，加權(quán)法、統(tǒng)計估計法雖然較直 接模擬方法減小了方差， 卻使每個粒子的運動鏈長增加， 或記錄貢獻

37、的計算時間 增加。因此，不能認為方差小的方法一定好，要從方法的效率全面考慮。在有些 情況下，直接模擬方法仍然是一個被廣泛使用的方法。5. 蒙特卡洛方法應(yīng)用程序 MCNP5.1 MCNP 簡述MCNP（A General Monte Carlo Code for Neutron and Particle Transport） 是一套通用的、三維空間中連續(xù)能量中子、光子和帶電粒子（離子）聯(lián)合輸運過 程模擬程序， 在軍事和工業(yè)領(lǐng)域有著廣泛應(yīng)用。 是基于蒙特卡洛方法的用于計算 三維復(fù)雜幾何結(jié)構(gòu)中的中子、光子、電子或者耦合中子、光子、電子輸運問題的 通用軟件包， 也具有計算核臨界系統(tǒng) （包括次臨界和超臨

38、界系統(tǒng)） 本征值問題的 能力。該軟件包通過FORTRAN言編程實現(xiàn)。MCN程序具有超強的幾何處理能力，幾何系統(tǒng)由幾何空間單元（cell）組成，而幾何空間單元的界面 （surface） 由平面、二次曲面及特殊的四次橢圓環(huán)曲面組 成。幾何空間單元中的材料由包括同位素在內(nèi)的多種核素組成， 使用精確的點截 面參數(shù)，對特定的評價庫（ENDF/B-IV,V,V,VI庫或ENDL85庫），考慮了該庫給出 的所有中子反應(yīng)類型。 在截面數(shù)據(jù)文件中收集了多種評價庫的數(shù)據(jù)。 對熱中子還 配備了相應(yīng)的截面數(shù)據(jù)，可按自由氣體模型或 S模型處理。對光子考慮了相干和 非相干散射，并處理了光電吸收后可能的熒光發(fā)射或電子對產(chǎn)生

39、。MCNP版（1983年）和3A版（1985年）發(fā)行后，這一軟件就成為用蒙特卡洛方法 模擬核過程最流行的通用程序， 程序在計算輻射能量沉積和輻射計量等方面取得 成功。88年出版的MCNP3程序具有重復(fù)構(gòu)造和結(jié)構(gòu)的能力，能夠解決特征丫譜 線的問題， 可以很好地模擬中子和光子的聯(lián)合輸運問題， 使用的主要核數(shù)據(jù)庫是 ENDF/B-4 91年MCNP版問世，這時程序可以模擬中子、光子、帶電粒子（離子）的聯(lián)合輸運過程，可以模擬探測器的測量結(jié)果。MCNP版使用了更新的ENDF/B-6評價核數(shù)據(jù)庫，加入了脈沖中子源功能等。MCNR版（2003年）提高了彩色描點能力（ 64種顏色），提高了處理中性粒子照相問題

40、能力，為源增加了新選項，并 對廣泛應(yīng)用的windows系統(tǒng)有了更好的支持。該程序是目前國際上在核技術(shù)領(lǐng)域 中應(yīng)用最廣泛、 效果較佳、 具有通用性的蒙特卡洛模擬計算程序， 許多核反應(yīng)蒙 特卡洛專用程序都引用該程序的核心部分。MCN程序涉及面如此之多，關(guān)鍵是通過讀入一個經(jīng)用戶創(chuàng)建的稱為INP的輸入文件來進行計算。該文件必須遵循按照柵元卡(card)的格式進行組織，指定描 述空間問題的信息，具體地有：(1)空間幾何體的描述說明；(2)幾何體的使用材 料描述和交叉區(qū)域的選擇估計；(3)中子、光子以及電子這3種粒子源的位置和特 性說明； 必要的回答卡和標(biāo)記卡的類型；(5)任何必需的冗余量消除技術(shù)，以 提高計算效率。目前，MCN以其靈活、通用的特點以及強大的功能被廣泛應(yīng)用于 輻射防護與射線測定、輻射屏蔽設(shè)計優(yōu)化、反應(yīng)堆設(shè)計、(次)臨界裝置實驗、醫(yī)學(xué)以及檢測器設(shè)計與分析等學(xué)科領(lǐng)域，并得到一致認可。5.2 MCNP誤差的估計蒙特卡洛方法的結(jié)果方法的結(jié)果代表被抽樣的許多歷史過程貢獻的平均值， 假定P(x)是選擇一個隨機步的幾率密度函數(shù)，x是這個隨機步產(chǎn)生的被估計的記 錄值，其平均值記為：E(x) = xP(x)dx(5.1)單NE(x)近似期望值可以通過MCT法得到：x二丄，其中N是粒子數(shù)目，Xi是從N yP(x)中第i個歷史的值從加強大數(shù)定理：P n