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

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

2、特卡洛方法的效率165. 蒙特卡洛方法應(yīng)用程序MCNP175.1 MCNP簡(jiǎn)述175.2 MCNP誤差的估計(jì)185.3 MCNP效率因素196. 結(jié)論19參考文獻(xiàn)201. 前言半個(gè)多世紀(jì)以來，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明，蒙特卡洛(Monte Carlo)方法作為一種獨(dú)立的方法被提出來，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。蒙特卡洛方法是一種計(jì)算方法，但與一般數(shù)值計(jì)算方法有很大區(qū)別。它是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。由于蒙特卡洛方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過程，解決一些數(shù)值方法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。蒙特卡洛模擬計(jì)算是解決中子在介質(zhì)中輸運(yùn)較為成熟

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

4、射輸運(yùn)問題的一種相當(dāng)成熟和有效的方法，而且它對(duì)于各種復(fù)雜問題，具有良好的通用性，實(shí)用性相當(dāng)廣泛，幾乎涉及核科學(xué)的各個(gè)領(lǐng)域。本文主要介紹蒙特卡洛的概念、原理和應(yīng)用及研究現(xiàn)狀。2. 蒙特卡洛方法概述蒙特卡洛方法又稱隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法。半個(gè)多世紀(jì)以來，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明 ，這種方法作為一種獨(dú)立的方法被提出來，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。蒙特卡洛方法是一種計(jì)算方法，但與一般數(shù)值計(jì)算方法有很大區(qū)別。它是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。由于蒙特卡洛方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過程，解決一些數(shù)值方法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。蒙特卡洛方法

5、的主要組成部分有：(1)概率密度函數(shù)(pdf) 必須給出描述一個(gè)物理系統(tǒng)的一組概率密度函數(shù)；(2)隨機(jī)數(shù)產(chǎn)生器能夠產(chǎn)生在區(qū)間0,1上均勻分布的隨機(jī)數(shù)；(3)抽樣規(guī)則如何從在區(qū)間0,1上均勻分布的隨機(jī)數(shù)出發(fā),隨機(jī)抽取服從給定的pdf的隨機(jī)變量；(4)模擬結(jié)果記錄記錄一些感興趣的量的模擬結(jié)果；(5)誤差估計(jì)必須確定統(tǒng)計(jì)誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；(6)減少方差的技術(shù)利用該技術(shù)可減少模擬過程中計(jì)算的次數(shù)；(7)并行和矢量化可以在先進(jìn)的并行計(jì)算機(jī)上運(yùn)行的有效算法2.1 蒙特卡洛方法的基本思想可以通俗地說，蒙特卡洛方法是用隨機(jī)試驗(yàn)的方法計(jì)算積分，即將所要計(jì)算的積分看作服從某種分布密度函

6、數(shù)f(r)的隨機(jī)變量(r)的數(shù)學(xué)期望通過某種試驗(yàn)，得到N個(gè)觀察值r1，r2，rN（用概率語言來說，從分布密度函數(shù)f(r)中抽取N個(gè)子樣r1，r2，rN，），將相應(yīng)的N個(gè)隨機(jī)變量的值g(r1)，g(r2)，g(rN)的算術(shù)平均值，作為積分的估計(jì)值（近似值）。為了得到具有一定精確度的近似解，所需試驗(yàn)的次數(shù)是很多的，通過人工方法作大量的試驗(yàn)相當(dāng)困難，甚至是不可能的。因此，蒙特卡洛方法的基本思想雖然早已被人們提出，卻很少被使用。本世紀(jì)四十年代以來，由于電子計(jì)算機(jī)的出現(xiàn)，使得人們可以通過電子計(jì)算機(jī)來模擬隨機(jī)試驗(yàn)過程，把巨大數(shù)目的隨機(jī)試驗(yàn)交由計(jì)算機(jī)完成，使得蒙特卡洛方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中

7、發(fā)揮應(yīng)有的作用。2.2 蒙特卡洛方法的收斂性、誤差蒙特卡洛方法作為一種計(jì)算方法，其收斂性與誤差是普遍關(guān)心的一個(gè)重要問題。2.2.1 蒙特卡洛方法的收斂性由前面介紹可知，蒙特卡洛方法是由隨機(jī)變量X的簡(jiǎn)單子樣X1，X2，XN的算術(shù)平均值.作為所求解的近似值。由大數(shù)定律可知，如X1，X2，XN獨(dú)立同分布，且具有有限期望值，則。即隨機(jī)變量X的簡(jiǎn)單子樣的算術(shù)平均值，當(dāng)子樣數(shù)N充分大時(shí)，以概率1收斂于它的期望值E(X)。2.2.2 蒙特卡洛方法的誤差蒙特卡洛方法的近似值與真值的誤差問題，概率論的中心極限定理給出了答案。該定理指出，如果隨機(jī)變量序列X1，X2，XN獨(dú)立同分布，且具有有限非零的方差2，即。f(

8、X)是X的分布密度函數(shù)。則當(dāng)N充分大時(shí)，有如下的近似式其中稱為置信度，1稱為置信水平。這表明，不等式近似地以概率1成立，且誤差收斂速度的階為。通常，蒙特卡洛方法的誤差定義為上式中與置信度是一一對(duì)應(yīng)的，根據(jù)問題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出。常用的與的對(duì)應(yīng)關(guān)系為：=0.5，=0.6745；=0.05，=0.96；=0.003，=3. 蒙特卡洛方法的誤差為概率誤差，這與其他數(shù)值計(jì)算方法是有區(qū)別的。誤差中的均方差是未知的，必須使用其估計(jì)值來代替以求出均方差。由式可知當(dāng)給定置信度后，誤差由和N決定。要減小，或者是增大N，或者是減小方差2。在固定的情況下，要把精度提高一個(gè)數(shù)量級(jí)，

9、試驗(yàn)次數(shù)N需增加兩個(gè)數(shù)量級(jí)。因此，單純?cè)龃驨不是一個(gè)有效的辦法。另一方面，如能減小估計(jì)的均方差，比如降低一半，那誤差就減小一半，這相當(dāng)于N增大四倍的效果。因此降低方差的各種技巧，引起了人們的普遍注意。2.3 蒙特卡洛方法的特點(diǎn)作為一種統(tǒng)計(jì)試驗(yàn)方法，蒙特卡洛方法因其優(yōu)點(diǎn)在諸多領(lǐng)域內(nèi)有著廣泛，但同時(shí)存在一些缺點(diǎn)。蒙特卡洛的主要優(yōu)點(diǎn)有：(1)能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過程。蒙特卡洛方法可以部分代替物理實(shí)驗(yàn)，甚至可以得到物理實(shí)驗(yàn)難以得到的結(jié)果。用蒙特卡洛方法解決實(shí)際問題，可以直接從實(shí)際問題本身出發(fā)，而不從方程或數(shù)學(xué)表達(dá)式出發(fā)。它有直觀、形象的特點(diǎn)。(2)受幾何條件限制小。在計(jì)

10、算s維空間中的任一區(qū)域Ds上的積分時(shí)，無論區(qū)域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產(chǎn)生N個(gè)點(diǎn)，得到積分的近似值，其中Ds為區(qū)域Ds的體積。這是數(shù)值方法難以作到的。(3)收斂速度與問題的維數(shù)無關(guān)。由誤差定義可知，在給定置信水平情況下，蒙特卡洛方法的收斂速度為，與問題本身的維數(shù)無關(guān)。維數(shù)的變化，只引起抽樣時(shí)間及估計(jì)量計(jì)算時(shí)間的變化，不影響誤差。也就是說，使用蒙特卡洛方法時(shí)，抽取的子樣總數(shù)N與維數(shù)s無關(guān)。維數(shù)的增加，除了增加相應(yīng)的計(jì)算量外，不影響問題的誤差。這一特點(diǎn)，決定了蒙特卡洛方法對(duì)多維問題的適應(yīng)性。(4)具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力。對(duì)于那些需要

11、計(jì)算多個(gè)方案的問題，使用蒙特卡洛方法有時(shí)不需要像常規(guī)方法那樣逐個(gè)計(jì)算，而可以同時(shí)計(jì)算所有的方案，其全部計(jì)算量幾乎與計(jì)算一個(gè)方案的計(jì)算量相當(dāng)。例如，對(duì)于屏蔽層為均勻介質(zhì)的平板幾何，要計(jì)算若干種厚度的穿透概率時(shí)，只需計(jì)算最厚的一種情況，其他厚度的穿透概率在計(jì)算最厚一種情況時(shí)稍加處理便可同時(shí)得到。(5)誤差容易確定。對(duì)于一般計(jì)算方法，要給出計(jì)算結(jié)果與真值的誤差并不是一件容易的事情，而蒙特卡洛方法則不然。根據(jù)蒙特卡洛方法的誤差公式，可以在計(jì)算所求量的同時(shí)計(jì)算出誤差。對(duì)干很復(fù)雜的蒙特卡洛方法計(jì)算問題，也是容易確定的。(6)程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)。在計(jì)算機(jī)上進(jìn)行蒙特卡洛方法計(jì)算時(shí)，程序結(jié)構(gòu)簡(jiǎn)單，分塊性強(qiáng)，

12、易于實(shí)現(xiàn)。蒙特卡洛的主要缺點(diǎn)有：(1)收斂速度慢。如前所述，蒙特卡洛方法的收斂速度為，一般不容易得到精確度較高的近似結(jié)果。對(duì)于維數(shù)少（三維以下）的問題，不如其他方法好。(2)誤差具有概率性。由于蒙特卡洛方法的誤差是在一定置信水平下估計(jì)的，所以它的誤差具有概率性，而不是一般意義下的誤差。(3)在粒子輸運(yùn)問題中，計(jì)算結(jié)果與系統(tǒng)大小有關(guān)。經(jīng)驗(yàn)表明，只有當(dāng)系統(tǒng)的大小與粒子的平均自由程可以相比較時(shí)（一般在十個(gè)平均自由程左右），蒙特卡洛方法計(jì)算的結(jié)果較為滿意。但對(duì)于大系統(tǒng)或小概率事件的計(jì)算問題，計(jì)算結(jié)果往往比真值偏低。而對(duì)于大系統(tǒng)，數(shù)值方法則是適用的。因此，在使用蒙特卡洛方法時(shí)，可以考慮把蒙特卡洛方法與解

13、析（或數(shù)值）方法相結(jié)合，取長(zhǎng)補(bǔ)短。2.4 蒙特卡洛方法的主要應(yīng)用范圍蒙特卡洛方法所特有的優(yōu)點(diǎn)，使得它的應(yīng)用范圍越來越廣。它的主要應(yīng)用范圍包括：粒子輸運(yùn)問題，統(tǒng)計(jì)物理，典型數(shù)學(xué)問題，真空技術(shù)，激光技術(shù)以及醫(yī)學(xué)，生物，探礦等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。蒙特卡洛方法在粒子輸運(yùn)問題中的應(yīng)用范圍主要包括：實(shí)驗(yàn)核物理，反應(yīng)堆物理，高能物理等方面。蒙特卡洛方法在實(shí)驗(yàn)核物理中的應(yīng)用范圍主要包括：通量及反應(yīng)率，中子探測(cè)效率，光子探測(cè)效率，光子能量沉積譜及響應(yīng)函數(shù)，氣體正比計(jì)數(shù)管反沖質(zhì)子譜，多次散射與通量衰減修正等方面。 3. 隨機(jī)數(shù)隨機(jī)數(shù)是蒙特卡洛方法的主要組成部分之一。隨機(jī)數(shù)是指一個(gè)數(shù)列，

14、其中的每一個(gè)體稱為隨機(jī)數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)。在一個(gè)均勻分布的隨機(jī)數(shù)中，每一個(gè)體出現(xiàn)的概率是均等的。物理中的很多過程需要隨機(jī)數(shù)確定，比如出射粒子的能量、方向等屬性，粒子與介質(zhì)的相互作用等等。所模擬的物理過程要求隨機(jī)數(shù)應(yīng)具有下列特性：1. 隨機(jī)數(shù)序列應(yīng)是獨(dú)立的、互不相關(guān)的(uncorrelated)：即序列中的任一子序列應(yīng)與其它的子序列無關(guān)；2. 長(zhǎng)的周期(long period)：實(shí)際應(yīng)用中，隨機(jī)數(shù)都是用數(shù)學(xué)方法計(jì)算出來的，這些算法具有周期性，即當(dāng)序列達(dá)到一定長(zhǎng)度后會(huì)重復(fù)；3. 均勻分布的隨機(jī)數(shù)應(yīng)滿足均勻性(Uniformity)：隨機(jī)數(shù)序列應(yīng)是均勻的、無偏的，即：如果兩個(gè)子區(qū)間的“面積

15、”相等，則落于這兩個(gè)子區(qū)間內(nèi)的隨機(jī)數(shù)的個(gè)數(shù)應(yīng)相等。4. 有效性（Efficiency):模擬結(jié)果可靠，隨機(jī)數(shù)的產(chǎn)生必須快速、有效，最好能夠進(jìn)行并行計(jì)算。為了產(chǎn)生隨機(jī)數(shù)，可以使用隨機(jī)數(shù)表。隨機(jī)數(shù)表是由0，1，9十個(gè)數(shù)字組成，每個(gè)數(shù)字以0.1的等概率出現(xiàn)，數(shù)字之間相互獨(dú)立。這些數(shù)字序列叫作隨機(jī)數(shù)字序列。如果要得到n位有效數(shù)字的隨機(jī)數(shù)，只需將表中每n個(gè)相鄰的隨機(jī)數(shù)字合并在一起，且在最高位的前邊加上小數(shù)點(diǎn)即可。例如，某隨機(jī)數(shù)表的第一行數(shù)字為7634258910，要想得到三位有效數(shù)字的隨機(jī)數(shù)依次為0.763，0.425，0.891?？梢允褂梦锢矸椒óa(chǎn)生隨機(jī)數(shù)，用來作為隨機(jī)數(shù)發(fā)生器的物理源主要有兩種：一種

16、是根據(jù)放射性物質(zhì)的放射性，另一種是利用計(jì)算機(jī)的固有噪聲。但在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)最實(shí)用、最常見的方法是數(shù)學(xué)方法，即用如下遞推公式： (3.1)產(chǎn)生隨機(jī)數(shù)序列。對(duì)于給定的初始值。1,2，k，確定n+k，=1，2。經(jīng)常使用的是k=1的情況。a) 用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)有兩個(gè)特點(diǎn)，即：遞推公式和初始值1,2，k確定后，整個(gè)隨機(jī)數(shù)序列便被唯一確定。不滿足隨機(jī)數(shù)相互獨(dú)立的要求。b) 由于隨機(jī)數(shù)序列是由遞推公式確定的，而在計(jì)算機(jī)上所能表示的0,1上的數(shù)又是有限的，因此，這種方法產(chǎn)生的隨機(jī)數(shù)序列就不可能不出現(xiàn)無限重復(fù)。一旦出現(xiàn)這樣的n，n(n< n)，使得成立隨機(jī)數(shù)序列便出現(xiàn)了周期性的循環(huán)現(xiàn)象。對(duì)于k=

17、1的情況，只要有一個(gè)隨機(jī)數(shù)重復(fù)，其后面的隨機(jī)數(shù)全部重復(fù)，這與隨機(jī)數(shù)的要求是不相符的。由于這兩個(gè)問題的存在，常稱用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)為偽隨機(jī)數(shù)。3.1 線性乘同余方法線性乘同余方法是由Lehmer在1951年提出來的，是一種最常用的產(chǎn)生偽隨機(jī)數(shù)的方法。乘同余方法中采用的遞推公式為 (3.2)其中：I0為初始值，a為乘法器，c為增值，m為模數(shù)，mod為取模運(yùn)算，除以m后的余數(shù)。a、c、m皆為整數(shù)。實(shí)型隨機(jī)序列： (3.3) (3.4)上式中，獨(dú)立性和均勻性取決于參數(shù)a和c的選擇。m應(yīng)盡可能地大，因?yàn)樾蛄械闹芷诓豢赡艽笥趍。通常將m取為計(jì)算機(jī)所能表示的最大的整型量，在32位機(jī)上，。1961年，M.

18、 Greenberger證明，用線性乘同余方法產(chǎn)生的隨機(jī)數(shù)序列具有周期m的條件是：(1)c和m互為質(zhì)數(shù)；(2)a-1是質(zhì)數(shù)p的倍數(shù)，其中p是a-1和m的公約數(shù)；(3)如果m是4的倍數(shù)，a-1也是4的倍數(shù)。3.2 偽隨機(jī)數(shù)序列的均勻性和獨(dú)立性3.2.1 偽隨機(jī)數(shù)的均勻性這里只考慮偽隨機(jī)數(shù)序列1,2，n全體作為子樣時(shí)的均勻性問題。其中n為偽隨機(jī)數(shù)序列的最大容量。對(duì)于任意的0x1，令Nn(x)表示偽隨機(jī)數(shù)序列1,2，n中適合不等式i< x,i=1,2,n的個(gè)數(shù)，則 (3.5)將偽隨機(jī)數(shù)序列1,2，n從小至大重新排列，令，則由(n)的定義，容易證明，很明顯，對(duì)于固定的，(n)的值越小越好。它是描

19、述偽隨機(jī)數(shù)序列均勻程度的基本量。對(duì)于任意隨機(jī)數(shù)序列，均有不等式成立。當(dāng)成立時(shí)，所對(duì)應(yīng)的偽隨機(jī)數(shù)序列為最佳分布。3.2.2 偽隨機(jī)數(shù)的獨(dú)立性對(duì)于任意，令表示(1,2),(2,3),(n,n+1)中適合不等式。的個(gè)數(shù)，根據(jù)隨機(jī)變量間相互獨(dú)立的定義和頻率近似概率的方法，令 (3.6)則(n)標(biāo)志偽隨機(jī)數(shù)序列1,2，n的獨(dú)立程度，簡(jiǎn)稱為獨(dú)立偏度。對(duì)于固定的n,(n)的值越接近于零，偽隨機(jī)數(shù)序列的獨(dú)立性越好。4. 蒙特卡洛方法在粒子輸運(yùn)上的應(yīng)用輻射（光子和中子）屏蔽問題是蒙特卡洛方法最早廣泛應(yīng)用的領(lǐng)域之一?，F(xiàn)主要從物理直觀出發(fā)，說明蒙特卡洛方法解決這類粒子輸運(yùn)問題的基本方法和技巧。解決屏蔽問題時(shí)可采取多

20、種方法，如直接模擬方法、簡(jiǎn)單加權(quán)法、統(tǒng)計(jì)估計(jì)法、指數(shù)變換法等，這里只對(duì)直接模擬方法做介紹。4.1 屏蔽問題模型在反應(yīng)堆工程和輻射的測(cè)量與應(yīng)用中，常常要用一些吸收材料做成屏蔽物擋住光子或中子。我們所關(guān)心的是經(jīng)過屏蔽后射線的強(qiáng)度及其能量分布，這就是屏蔽問題。當(dāng)屏蔽物的形狀復(fù)雜，散射各向異性，材料介質(zhì)不均勻 , 核反應(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ì)

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

22、間位置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ù)所求問題的類型和所用的方法來確定。對(duì)于無限平板幾何，取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)過 m 次碰撞后的狀態(tài)，其中rm：粒子在第 m 次碰撞點(diǎn)的位置Em：粒子第

23、 m 次碰撞后的能量m：粒子第 m 次碰撞后的運(yùn)動(dòng)方向tm：粒子到第 m 次碰撞時(shí)所經(jīng)歷的時(shí)間Wm：粒子第 m 次碰撞后的權(quán)重一個(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)序列描述，即S0，S1，SM-1，SM。或來描述。這里S0為粒子由源出發(fā)的狀態(tài)，稱為初態(tài)，SM 為粒子的終止?fàn)顟B(tài)。M稱為粒子運(yùn)動(dòng)的鏈長(zhǎng)。4.2.2 模擬運(yùn)動(dòng)過程這里以中子穿透均勻平板的模型來說明，這時(shí)狀態(tài)參數(shù)取S(z,E,cos)。模擬的步驟如下。(1)確定初始狀態(tài)S0：確定

24、粒子的初始狀態(tài)，實(shí)際上就是要從中子源的空間位置、能量和方向分布中抽樣。設(shè)源分布為，則分別從各自的分布中抽樣確定初始狀態(tài)。(2)確定下一個(gè)碰撞點(diǎn)：已知狀態(tài)Sm-1，要確定狀態(tài)Sm，首先要確定下一個(gè)碰撞點(diǎn)的位置zm。在相鄰兩次碰撞之間，中子的輸運(yùn)長(zhǎng)度l服從如下分布： (4.1)對(duì)于平板模型，l服從分布： (4.2)其中，t為介質(zhì)的中子宏觀總截面。積分稱為粒子輸運(yùn)的自由程數(shù)。顯然，粒子輸運(yùn)的自由程數(shù)服從指數(shù)分布，因此從f(l)中抽樣確定l，就是要從積分方程中解出l。對(duì)于單一介質(zhì)，則下一個(gè)碰撞點(diǎn)的位置為 (4.3)如果zma，則中子穿透屏蔽，若zm0, 則中子被反射出屏蔽。這兩種情況，均視為中子歷史終

25、止。(3)確定被碰撞的原子核：通常介質(zhì)由幾種原子核組成，中子與核碰撞時(shí)，要確定與哪一種核碰撞。設(shè)介質(zhì)由A、B、C 三種原子核組成，其核密度分別為NA、NB、NC，則介質(zhì)的宏觀總截面為： (4.4)其中分別為核A、B、C的宏觀總截面。其定義如下：，分別表示(·)核的宏觀總截面、核密度和微觀總截面。由于中子截面表示中子與核碰撞可能性的大小，因此，很自然地，中子與A、B、C核發(fā)生碰撞的幾率分別為： (4.5)若，則中子與A核碰撞；若，則中子與B核碰撞；若，則與C核碰撞。(4)確定碰撞類型：確定了碰撞的核(比如B核)后，就要進(jìn)一步確定碰撞類型。中子與核的反應(yīng)類型有彈性散射、非彈性散射、(n,

26、2n)反應(yīng)，裂變和俘獲等，它們的微觀截面分別為，則有 (4.6)各種反應(yīng)發(fā)生的幾率分別為 (4.7)利用離散型隨機(jī)變量的抽樣方法，確定反應(yīng)類型。在屏蔽問題中，中子與核反應(yīng)常只有彈性散射和吸收兩種類型，吸收截面為：。這時(shí)，總截面為： (4.8)發(fā)生彈性散射的幾率為：若，則為彈性散射；否則為吸收，發(fā)生吸收反應(yīng)意味著中子的歷史終止。(5)確定碰撞后的能量與運(yùn)動(dòng)方向：如果中子被碰撞核吸收，則其輸運(yùn)歷史結(jié)束。如果發(fā)生彈性散射，需要確定散射后中子的能量和運(yùn)動(dòng)方向。中子能量Em為： (4.9)其中。A是碰撞核的質(zhì)量與中子質(zhì)量之比，一般就取元素的原子量；C為質(zhì)心系中中子散射前后方向間的夾角，即偏轉(zhuǎn)角?？蓮馁|(zhì)心

27、系中彈性散射角分布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 ：。各角度關(guān)系如圖4.2所示。圖4.2 角度關(guān)系示意圖至此，由Sm-1完全可以確定Sm。因此，當(dāng)中子由源出發(fā)后，即S0確定后，重復(fù)步驟 (2)(5)，直到中子游動(dòng)歷史終止。于是得到了一個(gè)中子的隨機(jī)游動(dòng)歷史S0，S1，SM-1，SM，即也就是模擬了一個(gè)由源發(fā)出的中子的運(yùn)動(dòng)過程。4.2.3 記錄結(jié)果在獲得中子的隨機(jī)游動(dòng)歷史后，我們要對(duì)所要計(jì)算的物理量進(jìn)行估計(jì)。對(duì)于屏蔽問題，

28、我們要計(jì)算中子的穿透率。考察每個(gè)中子的隨機(jī)游動(dòng)歷史，它可能穿透屏蔽（zMa），可能被屏蔽發(fā)射回來（zM0），或者被吸收。設(shè)第n個(gè)中子對(duì)穿透的貢獻(xiàn)為n，則如果我們共跟蹤了N 個(gè)中子，則穿透屏蔽的中子數(shù)為：。則穿透屏蔽概率的近似值為： (4.10)我們稱這種直觀地模擬過程和估計(jì)方法為直接模擬方法。在置信水平10.95時(shí)，的誤差為： (4.11)其中為n的均方差，由于n是一個(gè)服從二項(xiàng)分布的隨機(jī)變量，所以 (4.12)或 (4.13)為得到中子穿透屏蔽的能量、角分布，將能量、角度范圍分成若干個(gè)間隔：，。其中Emax，Emin分別表示能量的上、下限，對(duì)于穿透屏蔽的中子按其能量、方向分間隔記錄。設(shè)一穿透屏

29、蔽的中子能量為EM，其運(yùn)動(dòng)方向與Z軸夾角為M，若能量EM屬于第 i 個(gè)能量間隔Ei，角度M屬于第 j 個(gè)角度間隔j，則分別在第 i 個(gè)能量計(jì)數(shù)器及第 j 個(gè)角度計(jì)數(shù)器中加1。跟蹤N個(gè)中子后，則 (4.14) (4.15)分別為穿透中子的能量分布和角分布。其中N1,i和N2,i分別為第 i 個(gè)能量和第 j 個(gè)角度間隔的穿透中子數(shù)。歸一后分別為： (4.16) (4.17)4.3 蒙特卡洛方法的效率衡量蒙特卡洛技巧的好壞，除了看其方差大小外，還要看其所需費(fèi)用（計(jì)算時(shí)間）多少，即從該技巧的效率Ef（方差與費(fèi)用乘積的倒數(shù)）全面考慮： (4.18)其中2為方差，T為所需費(fèi)用。Ef大時(shí)，所用方法的效率高；

30、否則，效率低。在一般情況下，直接模擬方法、簡(jiǎn)單加權(quán)法、統(tǒng)計(jì)估計(jì)法、指數(shù)變換法等方法中有些方法雖然減小了方差，卻增加了費(fèi)用。例如，加權(quán)法、統(tǒng)計(jì)估計(jì)法雖然較直接模擬方法減小了方差，卻使每個(gè)粒子的運(yùn)動(dòng)鏈長(zhǎng)增加，或記錄貢獻(xiàn)的計(jì)算時(shí)間增加。因此，不能認(rèn)為方差小的方法一定好，要從方法的效率全面考慮。在有些情況下，直接模擬方法仍然是一個(gè)被廣泛使用的方法。5. 蒙特卡洛方法應(yīng)用程序MCNP5.1 MCNP簡(jiǎn)述MCNP(A General Monte Carlo Code for Neutron and Particle Transport)是一套通用的、三維空間中連續(xù)能量中子、光子和帶電粒子（離子）聯(lián)合輸運(yùn)過

31、程模擬程序，在軍事和工業(yè)領(lǐng)域有著廣泛應(yīng)用。是基于蒙特卡洛方法的用于計(jì)算三維復(fù)雜幾何結(jié)構(gòu)中的中子、光子、電子或者耦合中子、光子、電子輸運(yùn)問題的通用軟件包，也具有計(jì)算核臨界系統(tǒng)（包括次臨界和超臨界系統(tǒng)）本征值問題的能力。該軟件包通過FORTRAN語言編程實(shí)現(xiàn)。MCNP程序具有超強(qiáng)的幾何處理能力，幾何系統(tǒng)由幾何空間單元(cell)組成，而幾何空間單元的界面(surface)由平面、二次曲面及特殊的四次橢圓環(huán)曲面組成。幾何空間單元中的材料由包括同位素在內(nèi)的多種核素組成，使用精確的點(diǎn)截面參數(shù)，對(duì)特定的評(píng)價(jià)庫(ENDF/B-IV,V,V,VI庫或ENDL851庫)，考慮了該庫給出的所有中子反應(yīng)類型。在截面

32、數(shù)據(jù)文件中收集了多種評(píng)價(jià)庫的數(shù)據(jù)。對(duì)熱中子還配備了相應(yīng)的截面數(shù)據(jù)，可按自由氣體模型或S模型處理。對(duì)光子考慮了相干和非相干散射，并處理了光電吸收后可能的熒光發(fā)射或電子對(duì)產(chǎn)生。MCNP3版(1983 年)和3A版(1985年)發(fā)行后，這一軟件就成為用蒙特卡洛方法模擬核過程最流行的通用程序，程序在計(jì)算輻射能量沉積和輻射計(jì)量等方面取得成功。88年出版的 MCNP3B程序具有重復(fù)構(gòu)造和結(jié)構(gòu)的能力，能夠解決特征譜線的問題，可以很好地模擬中子和光子的聯(lián)合輸運(yùn)問題，使用的主要核數(shù)據(jù)庫是ENDF/B-4。91年MCNP4版問世，這時(shí)程序可以模擬中子、光子、帶電粒子(離子)的聯(lián)合輸運(yùn)過程，可以模擬探測(cè)器的測(cè)量結(jié)果

33、。MCNP4版使用了更新的ENDF/B-6評(píng)價(jià)核數(shù)據(jù)庫，加入了脈沖中子源功能等。MCNP5版（2003年）提高了彩色描點(diǎn)能力（64種顏色），提高了處理中性粒子照相問題能力，為源增加了新選項(xiàng)，并對(duì)廣泛應(yīng)用的windows系統(tǒng)有了更好的支持。該程序是目前國際上在核技術(shù)領(lǐng)域中應(yīng)用最廣泛、效果較佳、具有通用性的蒙特卡洛模擬計(jì)算程序，許多核反應(yīng)蒙特卡洛專用程序都引用該程序的核心部分。MCNP程序涉及面如此之多，關(guān)鍵是通過讀入一個(gè)經(jīng)用戶創(chuàng)建的稱為INP的輸入文件來進(jìn)行計(jì)算。該文件必須遵循按照柵元卡(card)的格式進(jìn)行組織，指定描述空間問題的信息，具體地有：(1)空間幾何體的描述說明；(2)幾何體的使用材

34、料描述和交叉區(qū)域的選擇估計(jì)；(3)中子、光子以及電子這3種粒子源的位置和特性說明；(4)必要的回答卡和標(biāo)記卡的類型；(5)任何必需的冗余量消除技術(shù)，以提高計(jì)算效率。目前，MCNP以其靈活、通用的特點(diǎn)以及強(qiáng)大的功能被廣泛應(yīng)用于輻射防護(hù)與射線測(cè)定、輻射屏蔽設(shè)計(jì)優(yōu)化、反應(yīng)堆設(shè)計(jì)、（次）臨界裝置實(shí)驗(yàn)、醫(yī)學(xué)以及檢測(cè)器設(shè)計(jì)與分析等學(xué)科領(lǐng)域，并得到一致認(rèn)可。5.2 MCNP誤差的估計(jì)蒙特卡洛方法的結(jié)果方法的結(jié)果代表被抽樣的許多歷史過程貢獻(xiàn)的平均值，假定P(x)是選擇一個(gè)隨機(jī)步的幾率密度函數(shù)，x是這個(gè)隨機(jī)步產(chǎn)生的被估計(jì)的記錄值，其平均值記為： (5.1)E(x)近似期望值可以通過MC方法得到：，其中N是粒子數(shù)目，xi是從P(x)中第i個(gè)歷史的值。從加強(qiáng)大數(shù)定理：，x的方差是離散度的度量，定義為：，為標(biāo)準(zhǔn)差，MCNP方法可以估計(jì)這個(gè)