蒙特卡洛方法在中子輸運(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 屏蔽問(wèn)題模型104.2 直接模擬方法114.2.1 狀態(tài)參數(shù)與狀態(tài)序列114.2.2 模擬運(yùn)動(dòng)過(guò)程124.2.3 記錄結(jié)果154.

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

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

4、解輻射輸運(yùn)問(wèn)題的一種相當(dāng)成熟和有效的方法,而且它對(duì)于各種復(fù)雜問(wèn)題,具有良好的通用性,實(shí)用性相當(dāng)廣泛,幾乎涉及核科學(xué)的各個(gè)領(lǐng)域.本文主要介紹蒙特卡洛的概念、原理和應(yīng)用及研究現(xiàn)狀.2 .蒙特卡洛方法概述蒙特卡洛方法又稱隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法.半個(gè)多世紀(jì)以來(lái),由于科學(xué)技術(shù)的開(kāi)展和電子計(jì)算機(jī)的創(chuàng)造,這種方法作為一種獨(dú)立的方法被提出來(lái),并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用.蒙特卡洛方法是一種計(jì)算方法,但與一般數(shù)值計(jì)算方法有很大區(qū)別.它是以概率統(tǒng)計(jì)理論為根底的一種方法.由于蒙特卡洛方法能夠比擬逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程,解決一些數(shù)值方法難以解決的問(wèn)題,因而該方法的應(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ù)可減少模擬過(guò)程中計(jì)算的次數(shù)；(7)并行和矢量化一可以在先進(jìn)的并行計(jì)算機(jī)上運(yùn)行的有效算法2.1 蒙特卡洛方法的根本思想可以通俗地說(shuō),蒙特卡洛方法是用隨機(jī)試驗(yàn)的方法計(jì)算積分,即將所要計(jì)算的積分看作服

6、從某種分布密度函數(shù)f(r)的隨機(jī)變量g(r)的數(shù)學(xué)期望:二g=0g(r)f(r)dr(0.1)通過(guò)某種試驗(yàn),得到N個(gè)觀察值n,rN(用概率語(yǔ)言來(lái)說(shuō),從分布密度函數(shù)f(r)中抽取N個(gè)子樣門(mén),2,rN,),將相應(yīng)的N個(gè)隨機(jī)變量的值g(n),1Ng(,g(rN)的算術(shù)平均值3n=,£g(G,作為積分的估計(jì)值(近似值).Ny為了得到具有一定精確度的近似解,所需試驗(yàn)的次數(shù)是很多的,通過(guò)人工方法作大量的試驗(yàn)相當(dāng)困難,甚至是不可能的.因此,蒙特卡洛方法的根本思想雖然早已被人們提出,卻很少被使用.本世紀(jì)四十年代以來(lái),由于電子計(jì)算機(jī)的出現(xiàn),使得人們可以通過(guò)電子計(jì)算機(jī)來(lái)模擬隨機(jī)試驗(yàn)過(guò)程,把巨大數(shù)目的隨機(jī)

7、試驗(yàn)交由計(jì)算機(jī)完成,使得蒙特卡洛方法得以廣泛地應(yīng)用,在現(xiàn)代化的科學(xué)技術(shù)中發(fā)揮應(yīng)有的作用.2.2 蒙特卡洛方法的收斂性、誤差蒙特卡洛方法作為一種計(jì)算方法,其收斂性與誤差是普遍關(guān)心的一個(gè)重要問(wèn)題.2.2.1 蒙特卡洛方法的收斂性由前面介紹可知,蒙特卡洛方法是由隨機(jī)變量X的簡(jiǎn)單子樣X(jué)1,X2,-1.NXn的算術(shù)平均值XN=一乙Xi.作為所求解的近似值.由大數(shù)止律可知,如X1,NYX2,Xn獨(dú)立同分布,且具有有限期望值,那么p1limxn=e(x)=i.即隨機(jī)變量X的簡(jiǎn)單子樣的算術(shù)平均值Xn,當(dāng)子樣數(shù)N充分大時(shí),以概率1收斂于它的期望值E(X)02.2.2 蒙特卡洛方法的誤差蒙特卡洛方法的近似值與真值

8、的誤差問(wèn)題,概率論的中央極限定理給出了答案.該定理指出,如果隨機(jī)變量序列Xi,X2,Xn獨(dú)立同分布,且具有有限非零的方差即0#仃2=f(x-E(X)2f(x)dx<g.f(X)是X的分布密度函數(shù)c那么_1疝-_I11.X2/2(0.2)(0.3)E(X)iNmPTiXN-E(X)<x存Ledt當(dāng)N充分大時(shí),有如下的近似式P1Xn-E(X)產(chǎn)e,2/2dt=1-a42其中a稱為置信度,1a稱為置信水平.這說(shuō)明,不等式Xn近似地以概率1-a成立,且誤差收斂速度的階為O(N/2).通常,蒙特卡洛方法的誤差e定義為(0.4)上式中%與置信度a是一一對(duì)應(yīng)的,根據(jù)問(wèn)題的要求確定出置信水平后,查

9、標(biāo)準(zhǔn)正態(tài)分布表,就可以確定出%.常用的a與心的對(duì)應(yīng)關(guān)系為：a=0.5%=0.6745;a=0.05%=0.96;a=0.003%=3.蒙特卡洛方法的誤差為概率誤差,這與其他數(shù)值計(jì)算方法是有區(qū)別的.誤差中的均方差是未知的,必須使用其估計(jì)值來(lái)代替以求出均方差偽由式(0.4)可知當(dāng)給定置信度a后,誤差e由6和N決定.要減小6,或者是增大N,或者是減小方差"在b固定的情況下,要把精度提升一個(gè)數(shù)量級(jí),試驗(yàn)次數(shù)N需增加兩個(gè)數(shù)量級(jí).因此,單純?cè)龃驨不是一個(gè)有效的方法.另一方面,如能減小估計(jì)的均方差o,比方降低一半,那誤差就減小一半,這相當(dāng)于N增大四倍的效果.因此降低方差的各種技巧,引起了人們的普遍

10、注意.1N1N7前XF?2(0.5)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)過(guò)程.蒙特卡洛方法可以局部代替物理實(shí)驗(yàn),甚至可以得到物理實(shí)驗(yàn)難以得到的結(jié)果.用蒙特卡洛方法解決實(shí)際問(wèn)題,可以直接從實(shí)際問(wèn)題本身出發(fā),而不從方程或數(shù)學(xué)表達(dá)式出發(fā).它有直觀、形象的特點(diǎn).(2)受幾何條件限制小.在計(jì)算s維空間中的任一區(qū)域Ds上的積分g=19(為雙27x()xx2dxs時(shí),無(wú)論區(qū)域Ds的形狀多么特殊,只要能給出Ds描述Ds的幾何特征的條件,就可以從Ds中均勻產(chǎn)生

11、N個(gè)點(diǎn)(姆,x2,婢),得n到積分的近似值gN=-DsSg(x,x2i),x,),其中Ds為區(qū)域Ds的體積.這是數(shù)Ny值方法難以作到的.(3)收斂速度與問(wèn)題的維數(shù)無(wú)關(guān).由誤差定義可知,在給定置信水平情況下,蒙特卡洛方法的收斂速度為O(N/2),與問(wèn)題本身的維數(shù)無(wú)關(guān).維數(shù)的變化,只引起抽樣時(shí)間及估計(jì)量計(jì)算時(shí)間的變化,不影響誤差.也就是說(shuō),使用蒙特卡洛方法時(shí),抽取的子樣總數(shù)N與維數(shù)s無(wú)關(guān).維數(shù)的增加,除了增加相應(yīng)的計(jì)算量外,不影響問(wèn)題的誤差.這一特點(diǎn),決定了蒙特卡洛方法對(duì)多維問(wèn)題的適應(yīng)性.(4)具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的水平.對(duì)于那些需要計(jì)算多個(gè)方案的問(wèn)題,使用蒙特卡洛方法有時(shí)不需要像常規(guī)

12、方法那樣逐個(gè)計(jì)算,而可以同時(shí)計(jì)一種情況,其他厚度的穿透概率在計(jì)算最厚一種情況時(shí)稍加處理便可同時(shí)得到.蔽層為均勻介質(zhì)的平板幾何,要計(jì)算假設(shè)干種厚度的穿透概率時(shí),例如,對(duì)于屏只需計(jì)算最厚的(5)誤差容易確定.對(duì)于一般計(jì)算方法,要給出計(jì)算結(jié)果與真值的誤差并不是一件容易的事情,而蒙特卡洛方法那么不然.根據(jù)蒙特卡洛方法的誤差公式,可以在計(jì)算所求量的同時(shí)計(jì)算出誤差.對(duì)干很復(fù)雜的蒙特卡洛方法計(jì)算問(wèn)題,也是容易確定的.(6)程序結(jié)構(gòu)簡(jiǎn)單,易于實(shí)現(xiàn).在計(jì)算機(jī)上進(jìn)行蒙特卡洛方法計(jì)算時(shí),程序結(jié)構(gòu)簡(jiǎn)單,分塊性強(qiáng),易于實(shí)現(xiàn).蒙特卡洛的主要缺點(diǎn)有：(1)收斂速度慢.如前所述,蒙特卡洛方法的收斂速度為O(N-/2),一般不

13、容易得到精確度較高的近似結(jié)果.對(duì)于維數(shù)少(三維以下)的問(wèn)題,不如其他方法好.(2)誤差具有概率性.由于蒙特卡洛方法的誤差是在一定置信水平下估計(jì)的,所以它的誤差具有概率性,而不是一般意義下的誤差.(3)在粒子輸運(yùn)問(wèn)題中,計(jì)算結(jié)果與系統(tǒng)大小有關(guān).經(jīng)驗(yàn)說(shuō)明,只有當(dāng)系統(tǒng)的大小與粒子的平均自由程可以相比擬時(shí)(一般在十個(gè)平均自由程左右),蒙特卡洛方法計(jì)算的結(jié)果較為滿意.但對(duì)于大系統(tǒng)或小概率事件的計(jì)算問(wèn)題,計(jì)算結(jié)果往往比真值偏低.而對(duì)于大系統(tǒng),數(shù)值方法那么是適用的.因此,在使用蒙特卡洛方法時(shí),可以考慮把蒙特卡洛方法與解析(或數(shù)值)方法相結(jié)合,取長(zhǎng)補(bǔ)短.2.4 蒙特卡洛方法的主要應(yīng)用范圍蒙特卡洛方法所特有的優(yōu)

14、點(diǎn),使得它的應(yīng)用范圍越來(lái)越廣.它的主要應(yīng)用范圍包括：粒子輸運(yùn)問(wèn)題,統(tǒng)計(jì)物理,典型數(shù)學(xué)問(wèn)題,真空技術(shù),激光技術(shù)以及醫(yī)學(xué),生物,探礦等方面.隨著科學(xué)技術(shù)的開(kāi)展,其應(yīng)用范圍將更加廣泛.蒙特卡洛方法在粒子輸運(yùn)問(wèn)題中的應(yīng)用范圍主要包括：實(shí)驗(yàn)核物理,反響堆物理,高能物理等方面.蒙特卡洛方法在實(shí)驗(yàn)核物理中的應(yīng)用范圍主要包括：通量及反響率,中子探測(cè)效率,光子探測(cè)效率,光子能量沉積譜及響應(yīng)函數(shù),氣體正比計(jì)數(shù)管反沖質(zhì)子譜,屢次散射與通量衰減修正等方面.3.隨機(jī)數(shù)隨機(jī)數(shù)是蒙特卡洛方法的主要組成局部之一.隨機(jī)數(shù)是指一個(gè)數(shù)列,其中的每一個(gè)體稱為隨機(jī)數(shù),其值與數(shù)列中的其它數(shù)無(wú)關(guān).在一個(gè)均勻分布的隨機(jī)數(shù)中,每一個(gè)體出現(xiàn)的概

15、率是均等的.物理中的很多過(guò)程需要隨機(jī)數(shù)確定,比方出射粒子的能量、方向等屬性,粒子與介質(zhì)的相互作用等等.所模擬的物理過(guò)程要求隨機(jī)數(shù)應(yīng)具有以下特性：1,隨機(jī)數(shù)序列應(yīng)是獨(dú)立的、互不相關(guān)的(uncorrelated):即序列中的任一子序列應(yīng)與其它的子序列無(wú)關(guān)；2 .長(zhǎng)的周期(longperiod):實(shí)際應(yīng)用中,隨機(jī)數(shù)都是用數(shù)學(xué)方法計(jì)算出來(lái)的,這些算法具有周期性,即當(dāng)序列到達(dá)一定長(zhǎng)度后會(huì)重復(fù)；3 .均勻分布的隨機(jī)數(shù)應(yīng)滿足均勻性(Uniformity):隨機(jī)數(shù)序列應(yīng)是均勻的、無(wú)偏的,即：如果兩個(gè)子區(qū)間的-面積I相等,那么落于這兩個(gè)子區(qū)間內(nèi)的隨機(jī)數(shù)的個(gè)數(shù)應(yīng)相等.4 .有效性(Efficiency):模擬結(jié)果

16、可靠,隨機(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可以使用物理方法產(chǎn)生隨機(jī)數(shù),用來(lái)作為隨機(jī)數(shù)發(fā)生器的物理源主要有兩種：一種是根據(jù)放射性物質(zhì)的放射性,另一種是利用計(jì)算機(jī)的固有噪聲.但在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)最實(shí)用、最常見(jiàn)的方法

17、是數(shù)學(xué)方法,即用如下遞推公式：nk=T(n,n1,nk),n=1,2,(3.1)產(chǎn)生隨機(jī)數(shù)序列.對(duì)于給定的初始值.&,去,泉,確定&+k,n=1,2.經(jīng)常使用的是k=1的情況.a)用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)有兩個(gè)特點(diǎn),即：遞推公式和初始值8,玄,反確定后,整個(gè)隨機(jī)數(shù)序列便被唯一確定.不滿足隨機(jī)數(shù)相互獨(dú)立的要求.b)由于隨機(jī)數(shù)序列是由遞推公式確定的,而在計(jì)算機(jī)上所能表示的0,1上的數(shù)又是有限的,因此,這種方法產(chǎn)生的隨機(jī)數(shù)序列就不可能不出現(xiàn)無(wú)限重復(fù).一旦出現(xiàn)這樣的n"n(n,門(mén),)使得彳*=由i=1,2,k成立隨機(jī)數(shù)序列便出現(xiàn)了周期性的循環(huán)現(xiàn)象.對(duì)于k=1的情況,只要有一個(gè)隨

18、機(jī)數(shù)重復(fù),其后面的隨機(jī)數(shù)全部重復(fù),這與隨機(jī)數(shù)的要求是不相符的.由于這兩個(gè)問(wèn)題的存在,常稱用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)為偽隨機(jī)數(shù).3.1線性乘同余方法線性乘同余方法是由Lehmer在1951年提出來(lái)的,是一種最常用的產(chǎn)生偽隨機(jī)數(shù)的方法.乘同余方法中采用的遞推公式為In1=(aInc)modm(3.2)其中：I0為初始值,a為乘法器,c為增值,m為模數(shù),mod為取模運(yùn)算,(aIn+c)除以m后的余數(shù).&c、m皆為整數(shù).實(shí)型隨機(jī)序列：Inrnn0,1).In<m(3.3)float(m)rn0,1),In<m-1(3.4)float(m-1)上式中,獨(dú)立性和均勻性取決于參數(shù)a和c的選擇

19、.m應(yīng)盡可能地大,由于序列的周期不可能大于m.通常將m取為計(jì)算機(jī)所能表示的最大的整型量,在32位機(jī)上,319m=2=210o1961年,M.Greenbergerffi明,用線性乘同余方法產(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ù)序列a,宜,&全體作為子樣時(shí)的土勻性問(wèn)題.其中n為偽隨機(jī)數(shù)序列的最大容量.對(duì)于任意的0Wx0,l令Nn(x)表示偽隨機(jī)數(shù)序列a,W,&中適合不等式&&

20、lt;x,i=1,2,的個(gè)數(shù),那么(n)=sup|0/N-x|n(3.5)將偽隨機(jī)數(shù)序列己1,己2式從小至大重新排列,飛,壬一三,令<=0,"書(shū)=1,那么由n)的定義,容易證實(shí)I1產(chǎn)T一一mo-上|,|;-i|1>,很明顯,對(duì)于固定nn的n,n的值越小越好.它是描述偽隨機(jī)數(shù)序列均勻程度的根本量.對(duì)于任意1.1.隨機(jī)數(shù)序列,均有不等式6nX一成立.當(dāng)an=一成立時(shí),所對(duì)應(yīng)的偽隨機(jī)2n2n數(shù)序列為最正確分布.3.2.2偽隨機(jī)數(shù)的獨(dú)立性對(duì)于任意0«x,y«1,令Nnx,y表示8,2,&3,n,刈中適合不等式,<x,4.<y.的個(gè)數(shù),根據(jù)隨

21、機(jī)變量間相互獨(dú)立的定義和頻率近似概率的方法,Nn(x,y)"n)=supl-0<x,y<inNn(x)Nn(y)1nn(3.6)那么en標(biāo)志偽隨機(jī)數(shù)序列比已,晝的獨(dú)立程度,簡(jiǎn)稱為獨(dú)立偏度.對(duì)于固定的n,en的值越接近于零,偽隨機(jī)數(shù)序列的獨(dú)立性越好.4.蒙特卡洛方法在粒子輸運(yùn)上的應(yīng)用輻射光子和中子屏蔽問(wèn)題是蒙特卡洛方法最早廣泛應(yīng)用的領(lǐng)域之一.現(xiàn)主要從物理直觀出發(fā),說(shuō)明蒙特卡洛方法解決這類粒子輸運(yùn)問(wèn)題的根本方法和技巧.解決屏蔽問(wèn)題時(shí)可采取多種方法,如直接模擬方法、簡(jiǎn)單加權(quán)法、統(tǒng)計(jì)估計(jì)法、指數(shù)變換法等,這里只對(duì)直接模擬方法做介紹.4.1 屏蔽問(wèn)題模型在反響堆工程和輻射的測(cè)量與應(yīng)

22、用中,常常要用一些吸收材料做成屏蔽物擋住光子或中子.我們所關(guān)心的是經(jīng)過(guò)屏蔽后射線的強(qiáng)度及其能量分布,這就是屏蔽問(wèn)題.當(dāng)屏蔽物的形狀復(fù)雜,散射各向異性,材料介質(zhì)不均勻,核反響截面與能量、位置有關(guān)時(shí),難以用數(shù)值方法求解,用蒙特卡洛方法能夠得到滿意的結(jié)果.粒子的輸運(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ī)律.蒙特卡洛模擬,實(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ù)來(lái)實(shí)現(xiàn)的為方便起見(jiàn),選用平板屏蔽模型,在厚度為a,長(zhǎng)、寬無(wú)限的

23、平板左側(cè)放置一個(gè)強(qiáng)度,具有能量、方向分布的輻射源S,見(jiàn)圖4.1.求粒子穿透屏蔽概率(穿透率)及其能量、方向分布.穿透率就是由源發(fā)出的平均一個(gè)粒子穿透屏蔽的數(shù)目.同時(shí),假定粒子在兩次碰撞之間按直線運(yùn)動(dòng),且粒子之間的相互作用可以忽略.屏蔽物圖4.1屏蔽問(wèn)題模型4.2 直接模擬方法直接模擬方法就是直接從物理問(wèn)題出發(fā),模擬粒子的真實(shí)物理過(guò)程.4.2.1 狀態(tài)參數(shù)與狀態(tài)序列粒子在介質(zhì)中的運(yùn)動(dòng)的狀態(tài),可用一組參數(shù)來(lái)描述,稱之為狀態(tài)參數(shù).它通常包括：粒子的空間位置r,能量E和運(yùn)動(dòng)方向Q,以$=(r,E,)表示.有時(shí)還需要其他的參數(shù),如粒子的時(shí)間t和附帶的權(quán)重W,這時(shí)狀態(tài)參數(shù)為S'=(r,E,Q,t,

24、W)0狀態(tài)參數(shù)通常要根據(jù)所求問(wèn)題的類型和所用的方法來(lái)確定.對(duì)于無(wú)限平板幾何,取S=(z,E,cos,其)中z為粒子的位置坐標(biāo),a為粒子的運(yùn)動(dòng)方向與Z軸的夾角.對(duì)于球?qū)ΨQ幾何,取S=(r,E,cos,乳肘r表示粒子所在位置到球心的距離,效粒子的運(yùn)動(dòng)方向與其所在位置的徑向夾角.粒子第m次碰撞后的狀態(tài)參數(shù)為Sm=七,Em,Qm),或S:=七,Em,Qm,tm)Wn),它表示一個(gè)由源發(fā)出的粒子,在介質(zhì)中經(jīng)過(guò)m次碰撞后的狀態(tài),其中rm：粒子在第m次碰撞點(diǎn)的位置Em：粒子第m次碰撞后的能量Qm：粒子第m次碰撞后的運(yùn)動(dòng)方向tm：粒子到第m次碰撞時(shí)所經(jīng)歷的時(shí)間Wm：粒子第m次碰撞后的權(quán)重一個(gè)由源發(fā)出的粒子在介

25、質(zhì)中運(yùn)動(dòng),經(jīng)過(guò)假設(shè)干次碰撞后,直到其運(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,SM0或r0,n,.工齡"E0,E1,Em/,EmkQ0,Q,Qmj.,Qm/來(lái)描述.這里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)過(guò)程這里以中子穿透均勻平板的模型來(lái)說(shuō)明,這時(shí)狀態(tài)參數(shù)取S=(z,E,cos磯摸擬的步驟如下.(1)確定初始狀態(tài)S.：確定粒子的初始狀態(tài),實(shí)際上就是要從中子源的空間位置、能量和方向分布中抽

26、樣.設(shè)源分布為f(z0,E0,cosct0)=3其止伍)f3(cos«0),那么分別從各自的分布中抽樣確定初始狀態(tài).(2)確定下一個(gè)碰撞點(diǎn)：狀態(tài)Sm-1,要確定狀態(tài)Sm,首先要確定下一個(gè)碰撞點(diǎn)的位置Zm.在相鄰兩次碰撞之間,中子的輸運(yùn)長(zhǎng)度l服從如下分布：f(l)=Et(rm+lQm,Em/)expj0工t(rm_L+rQm,Em/)dl'(4.1)對(duì)于平板模型,l服從分布：f(l)t(Zmlxp故)/mE,m三t(Zm二lcos：m,EmJ)dl/(4.2)l其中,X為介質(zhì)的中子宏觀總截面.積分J0,(rm+l''Qm,Em)dl'稱為粒子輸運(yùn)的自由程

27、數(shù).顯然,粒子輸運(yùn)的自由程數(shù)服從指數(shù)分布,因此從f中抽樣確定l,l就是要從積分方程J01t(rm+l''Qm,Em,)dl'=Tn上中解出l.對(duì)于單一介質(zhì)lnl=-,那么下一個(gè)碰撞點(diǎn)的包置為三t(Eml)(4.3)lnZm=Zmjlcos-m4=Zm4"COS-m4-t(Em)如果Zm刃,那么中子穿透屏蔽,假設(shè)Zm&0,那么中子被反射出屏蔽.這兩種情況,均視為中子歷史終止.3確定被碰撞的原子核：通常介質(zhì)由幾種原子核組成,中子與核碰撞時(shí),要確定與哪一種核碰撞.設(shè)介質(zhì)由A、B、C三種原子核組成,其核密度分別為Na、Nb、Nc,那么介質(zhì)的宏觀總截面為：ItE

28、m=E：Em三：Em言Em,4.4其中工A,工B,工C分別為核A、B、C的宏觀總截面.其定義如下：IEm=N產(chǎn)tEm,工Em、»,、仃品分別表示核的宏觀總截面、核密度和微觀總截面.由于中子截面表示中子與核碰撞可能性的大小,因此,很自然地,中子與A、B、C核發(fā)生碰撞的幾率分別為：PA=力EmLPB=Em.,PC=D4.5三tEmJ'Em.XEm假設(shè)EPa,那么中子與A核碰撞；假設(shè)CEPa+R,那么中子與B核碰撞；假設(shè)之a(chǎn)Pa+Pb,那么與C核碰撞.4確定碰撞類型：確定了碰撞的核比方B核后,就要進(jìn)一步確定碰撞類型.中子與核的反響類型有彈性散射、非彈性散射、n,2n反響,裂變和俘獲

29、等,它們的微觀截面分別為端EmJEm,、.:,2nEmJ邸Em工和仃;Em,那么有BBBBBB5(Em4)=Jl(Em)+5n(Em4)+b(n,2n)(Em4)+bf(Em)+Dc(Em)(4.6)各種反響發(fā)生的幾率分別為匕PnRn,2n)PfPC=4(Em,)WtB(Em)=：iB(Em)；、B(Eml)(n,2n)(Em4)/-tB(Em4)=；葉昌)；邸昌1)(4.7)利用離散型隨機(jī)變量的抽樣方法,確定反響類型.在屏蔽問(wèn)題中,中子與核反響二4(Em)；±tB(Eml)常只有彈性散射和吸收兩種類型,吸收截面為：仃；Em=O：Em十b：Em.這時(shí),總截面為："Em=；：

30、Em二；Em4.8發(fā)生彈性散射的幾率為：Pel=°BEmJ假設(shè)WPei,那么為彈性散射；否那么為吸收,二tEm發(fā)生吸收反響意味著中子的歷史終止.5確定碰撞后的能量與運(yùn)動(dòng)方向：如果中子被碰撞核吸收,那么其輸運(yùn)歷史結(jié)束.如果發(fā)生彈性散射,需要確定散射后中子的能量和運(yùn)動(dòng)方向.中子能量Em為：Em=停但11-rcos%4.9其中=4二2.A是碰撞核的質(zhì)量與中子質(zhì)量之比,一般就取元素的原子量；色A1為質(zhì)心系中中子散射前前方向間的夾角,即偏轉(zhuǎn)角.%=cos%可從質(zhì)心系中彈性散射角分布fc©中抽樣產(chǎn)生.實(shí)驗(yàn)室系散射角4的余弦乩為：Li1+AC產(chǎn)L-'c二°1A2A%如果

31、給出實(shí)驗(yàn)室系散射角余弦分布fLQ,可直接從fLQ中抽取乩,此時(shí)能量E2Em與比的關(guān)系式為：Em='2"+jA21+支.確定了實(shí)驗(yàn)室系散射角小后,再使用球面三角公式確定cosm:cos%=cosamucos9L+sin«msinLcos?.各角度關(guān)系如圖4.2所示圖4.2角度關(guān)系示意圖至此,由Sm-1完全可以確定Sm.因此,當(dāng)中子由源出發(fā)后,即So確定后,重復(fù)步于是得到了一個(gè)中子的隨機(jī)游動(dòng)歷史So,驟25,直到中子游動(dòng)歷史終止Si,Sm-i,Sm,即fZ0,Z1,ZMA,ZmE.,E1,em,Em<cosa0,cos%,coscM,coso(M)也就是模擬了一個(gè)

32、由源發(fā)出的中子的運(yùn)動(dòng)過(guò)程.4.2.3記錄結(jié)果在獲得中子的隨機(jī)游動(dòng)歷史后,我們要對(duì)所要計(jì)算的物理量進(jìn)行估計(jì).對(duì)于屏蔽問(wèn)題,我們要計(jì)算中子的穿透率.考察每個(gè)中子的隨機(jī)游動(dòng)歷史,它可能穿n個(gè)中子對(duì)穿透屏蔽ZM刃,可能被屏蔽發(fā)射回來(lái)ZM<0,或者被吸收.設(shè)第透的奉獻(xiàn)為小,那么n_11,當(dāng)Zm之a(chǎn)n=0,當(dāng)ZmM0如果我們共跟蹤了N個(gè)中子,那么穿透屏蔽的中子數(shù)為:NNi八nn1那么穿透屏蔽概率的近似值為:即二NiNNni(4.10)我們稱這種直觀地模擬過(guò)程和估計(jì)方法為直接模擬方法.在置信水平1-0.95時(shí),叩的誤差為:(4.11)(4.12)即-P<其中仃為中的均方差,由于汴是一個(gè)服從二項(xiàng)分

33、布的隨機(jī)變量,所以二2=P(1-P)(4.13)為得到中子穿透屏蔽的能量、角分布,將能量、角度范圍分成假設(shè)干個(gè)間隔:Emin=E|<<E1<E0=Emax,0=口.<口j=可2.其中Emax,Emin分別表示能量的上、下限,對(duì)于穿透屏蔽的中子按其能量、方向分間隔記錄.設(shè)一穿透屏蔽的中子能量為Em,其運(yùn)動(dòng)方向與Z軸夾角為cm,假設(shè)能量Em屬于第i個(gè)能量問(wèn)隔出i,角度0M屬于第j個(gè)角度間隔A.那么分別在第i個(gè)能量計(jì)數(shù)器及第j個(gè)角度計(jì)數(shù)器中加1.跟蹤N個(gè)中子后,那么Nii(4.14)Ji=1,2,IN石j=1,2,J分別為穿透中子的能量分布和角分布.其中N1,i和N2,i分別

34、為第i個(gè)能量和第j個(gè)角(4.15)度間隔的穿透中子數(shù).歸一后分別為:3*PN1,ig)PN1,i(4.16)p(1)*白PN2N2,jj=1,2,J(4.17)4.3蒙特卡洛方法的效率衡量蒙特卡洛技巧的好壞,除了看其方差大小外,還要看其所需費(fèi)用計(jì)算時(shí)間多少,即從該技巧的效率Ef方差與費(fèi)用乘積的倒數(shù)全面考慮：Ef4.18二T其中,為方差,T為所需費(fèi)用.Ef大時(shí),所用方法的效率高；否那么,效率低.在一般情況下,直接模擬方法、簡(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ì)

35、算時(shí)間增加.因此,不能認(rèn)為方差小的方法一定好,要從方法的效率全面考慮.在有些情況下,直接模擬方法仍然是一個(gè)被廣泛使用的方法.5.蒙特卡洛方法應(yīng)用程序一MCNP5.1MCNP簡(jiǎn)述MCNPAGeneralMonteCarloCodeforNeutronandParticleTransport是一套通用的、三維空間中連續(xù)能量中子、光子和帶電粒子離子聯(lián)合輸運(yùn)過(guò)程模擬程序,在軍事和工業(yè)領(lǐng)域有著廣泛應(yīng)用.是基于蒙特卡洛方法的用于計(jì)算三維復(fù)雜幾何結(jié)構(gòu)中的中子、光子、電子或者耦合中子、光子、電子輸運(yùn)問(wèn)題的通用軟件包,也具有計(jì)算核臨界系統(tǒng)包括次臨界和超臨界系統(tǒng)本征值問(wèn)題的水平.該軟件包通過(guò)FORTRAN語(yǔ)言編程

36、實(shí)現(xiàn).MCNP程序具有超強(qiáng)的幾何處理水平,幾何系統(tǒng)由幾何空間單元cell組成,而幾何空間單元的界面surface由平面、二次曲面及特殊的四次橢圓環(huán)曲面組成.幾何空間單元中的材料由包括同位素在內(nèi)的多種核素組成,使用精確的點(diǎn)截面參數(shù),對(duì)特定的評(píng)價(jià)庫(kù)ENDF/B-IV,V,V,VI庫(kù)或ENDL851庫(kù),考慮了該庫(kù)給出的所有中子反響類型.在截面數(shù)據(jù)文件中收集了多種評(píng)價(jià)庫(kù)的數(shù)據(jù).對(duì)熱中子還配備了相應(yīng)的截面數(shù)據(jù),可按自由氣體模型或S模型處理.對(duì)光子考慮了相干和非相干散射,并處理了光電吸收后可能的熒光發(fā)射或電子對(duì)產(chǎn)生.MCNP3版1983年和3A版1985年發(fā)行后,這一軟件就成為用蒙特卡洛方法模擬核過(guò)程最流

37、行的通用程序,程序在計(jì)算輻射能量沉積和輻射計(jì)量等方面取得成功.88年出版的MCNP3B程序具有重復(fù)構(gòu)造和結(jié)構(gòu)的水平,能夠解決特征丫譜線的問(wèn)題,可以很好地模擬中子和光子的聯(lián)合輸運(yùn)問(wèn)題,使用的主要核數(shù)據(jù)庫(kù)是ENDF/B-4.91年MCNP4版問(wèn)世,這時(shí)程序可以模擬中子、光子、帶電粒子離子的聯(lián)合輸運(yùn)過(guò)程,可以模擬探測(cè)器的測(cè)量結(jié)果.MCNP4版使用了更新的ENDF/B-6評(píng)價(jià)核數(shù)據(jù)庫(kù),參加了脈沖中子源功能等.MCNP5版200孫提升了彩色描點(diǎn)水平64種顏色,提升了處理中性粒子照相問(wèn)題水平,為源增加了新選項(xiàng),并對(duì)廣泛應(yīng)用的windows系統(tǒng)有了更好的支持.該程序是目前國(guó)際上在核技術(shù)領(lǐng)域中應(yīng)用最廣泛、效果

38、較佳、具有通用性的蒙特卡洛模擬計(jì)算程序,許多核反響蒙特卡洛專用程序都引用該程序的核心局部.MCNP程序涉及面如此之多,關(guān)鍵是通過(guò)讀入一個(gè)經(jīng)用戶創(chuàng)立的稱為INP的輸入文件來(lái)進(jìn)行計(jì)算.該文件必須遵循根據(jù)柵元卡card的格式進(jìn)行組織,指定描述空間問(wèn)題的信息,具體地有：1空間幾何體的描述說(shuō)明；2幾何體的使用材料描述和交叉區(qū)域的選擇估計(jì)；3中子、光子以及電子這3種粒子源的位置和特性說(shuō)明；4必要的答復(fù)卡和標(biāo)記卡的類型；5任何必需的冗余量消除技術(shù),以提升計(jì)算效率.目前,MCNP以其靈活、通用的特點(diǎn)以及強(qiáng)大的功能被廣泛應(yīng)用于輻射防護(hù)與射線測(cè)定、輻射屏蔽設(shè)計(jì)優(yōu)化、反響堆設(shè)計(jì)、(次)臨界裝置實(shí)驗(yàn)、醫(yī)學(xué)以及檢測(cè)器設(shè)計(jì)與分析等學(xué)科領(lǐng)域,并得到一致認(rèn)可.5.2MCNP誤差的估計(jì)蒙特卡洛方法的結(jié)果方法的結(jié)果代表被抽樣的許多歷史過(guò)程奉獻(xiàn)的平均值,假定P(x)是選擇一個(gè)隨機(jī)步的幾率密度函數(shù),x是這個(gè)隨機(jī)步產(chǎn)生的被估計(jì)的記錄值,其平均值記為：E(x)=xP(x)dx(5.1),一、一、,一1NE(x)近似期望值可以通