蒙特卡羅方法及應(yīng)用實(shí)驗(yàn)講義2016

1、蒙特卡羅方法及應(yīng)用實(shí)驗(yàn)講義東華理工大學(xué)核工系2016.815實(shí)驗(yàn)一 蒙特卡羅方法基本思想一、實(shí)驗(yàn)?zāi)康?、了解蒙特卡羅方法方法的基本思想；2、掌握蒙特卡羅方法計(jì)算面積、體積的方法；3、掌握由已知分布的隨機(jī)抽樣方法。二、實(shí)驗(yàn)原理Monte Carlo方法，又稱統(tǒng)計(jì)模擬方法或計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”進(jìn)行數(shù)值模擬的方法，一種采用統(tǒng)計(jì)抽樣理論近似求解物理或數(shù)學(xué)問(wèn)題的方法。如待求量可以表述成某些特征量的期望值、某些事件出現(xiàn)的概率或兩者的函數(shù)形式，那么可采用蒙特卡羅方法求解。在求解某些特征量的期望值或某些事件出現(xiàn)的概率時(shí)，必須構(gòu)建合符實(shí)際的數(shù)學(xué)模型。例如采用蒙特卡羅方法計(jì)算某函數(shù)所圍面積時(shí)，

2、構(gòu)建的數(shù)學(xué)模型是構(gòu)造一已知面積的可均勻抽樣區(qū)域，在該區(qū)域投點(diǎn)，由伯努利定理大數(shù)定理可知，進(jìn)入待求區(qū)域投點(diǎn)的頻率依概率1收斂于該事件出現(xiàn)的概率（面積之比）。由已知分布的隨機(jī)抽樣方法指的是由已知分布的總體中抽取簡(jiǎn)單子樣。具體方法很多，詳見(jiàn)教材第三章。三、實(shí)驗(yàn)內(nèi)容1、安裝所需計(jì)算工具（MATLAB、fortran、C+等）；2、學(xué)習(xí)使用rand(m,n)、unifrnd(a,b,m,n)函數(shù)3、求解下列問(wèn)題：3.0、蒲豐氏投針求圓周率。3.1、給定曲線y =2 x2 和曲線y3 = x2，曲線的交點(diǎn)為：P1( 1，1 )、P2( 1，1 )。曲線圍成平面有限區(qū)域，用蒙特卡羅方法計(jì)算區(qū)域面積；3.2、

3、計(jì)算所圍體積其中。4、對(duì)以下已知分布進(jìn)行隨機(jī)抽樣：4.1、；4.2、其中。四、實(shí)驗(yàn)報(bào)告編寫(xiě)1、給出各題的抽樣程序并解釋語(yǔ)句的含義；2、給出3.1和3.2抽樣結(jié)果誤差隨抽樣次數(shù)的關(guān)系圖，并解釋原因；表1 實(shí)驗(yàn)記錄表序號(hào)1234567試驗(yàn)次數(shù)1032×1045×1041.2×1051.5×1051.8×1062.0×107試驗(yàn)時(shí)間計(jì)算結(jié)果實(shí)驗(yàn)誤差3、給出4.1和4.2的抽樣框圖、試驗(yàn)累積頻率與理論累積頻率關(guān)系圖，并給出抽樣次數(shù)（>106）與抽樣時(shí)間。實(shí)驗(yàn)二 由已知分布的隨機(jī)抽樣方法一、實(shí)驗(yàn)?zāi)康?、掌握由已知分布的隨機(jī)抽樣方法。2、用

4、編程語(yǔ)言實(shí)現(xiàn)某具體隨機(jī)抽樣方法。二、實(shí)驗(yàn)原理由已知分布的隨機(jī)抽樣方法指的是由已知分布的總體中抽取簡(jiǎn)單子樣。具體方法很多，本實(shí)驗(yàn)綜合直接抽樣方法、挑選抽樣方法和替換抽樣方法，以散射方位角余弦分布的抽樣為例。實(shí)驗(yàn)原理詳見(jiàn)教材對(duì)應(yīng)章節(jié)。1.連續(xù)型分布的直接抽樣方法對(duì)于連續(xù)型分布，如果分布函數(shù)F(x) 的反函數(shù)F1(x)存在，則直接抽樣方法是：2.挑選抽樣方法為了實(shí)現(xiàn)從己知分布密度函數(shù)f(x)抽樣，選取與f(x)取值范圍相同的分布密度函數(shù)h(x)，如果則挑選抽樣方法為：3.替換法抽樣方法為了實(shí)現(xiàn)某個(gè)復(fù)雜的隨機(jī)變量 y 的抽樣，將其表示成若干個(gè)簡(jiǎn)單的隨機(jī)變量x1，x2，xn的函數(shù)得到 x1，x2，xn的

5、抽樣后，即可確定 y 的抽樣，這種方法叫作替換法抽樣。三、實(shí)驗(yàn)內(nèi)容蒲豐氏問(wèn)題的算法如何產(chǎn)生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均勻分布的，其分布密度函數(shù)為：類似地，的分布密度函數(shù)為：因此，產(chǎn)生任意的（x,）的過(guò)程就變成了由f1(x)抽樣x及由f2()抽樣的過(guò)程了。由此得到： 其中1，2均為（0,1）上均勻分布的隨機(jī)變量。每次投針試驗(yàn)，實(shí)際上變成在計(jì)算機(jī)上從兩個(gè)均勻分布的隨機(jī)變量中抽樣得到（x,），然后定義描述針與平行線相交狀況的隨機(jī)變量s(x,)，為 如果投針次，則 是針與平行線相交概率的估計(jì)值。事實(shí)上， 于是有四、實(shí)驗(yàn)報(bào)告1、給出源程序程序并解釋語(yǔ)句的含義；2、作出抽樣框

6、圖、試驗(yàn)累積頻率與理論累積頻率關(guān)系圖，并給出抽樣次數(shù)（>106）與抽樣時(shí)間。實(shí)驗(yàn)三 MCNP方法在實(shí)驗(yàn)核物理中的應(yīng)用一、實(shí)驗(yàn)?zāi)康?、了解MCNP程序運(yùn)行流程；2、掌握MCNP輸入文件編寫(xiě)規(guī)范；3、理解模擬內(nèi)容、并能編寫(xiě)輸入文件、運(yùn)行，并獲得計(jì)算結(jié)果；二、實(shí)驗(yàn)原理MCNP是一種常見(jiàn)的粒子輸運(yùn)模擬軟件，軟件的安裝、運(yùn)行和輸入文件編寫(xiě)方法詳見(jiàn)相關(guān)參考資料。MCNP輸入文件編寫(xiě)完成后，先確認(rèn)輸入模型是否正確，在DOS環(huán)境下進(jìn)行，打開(kāi)運(yùn)行DOS環(huán)境，進(jìn)行以下操作：DOS命令操作命令含義Mcnp i=name.inp o=name.o打開(kāi)畫(huà)圖框PX vx輸出模型在x=vx面上的切面PY vy輸出模型

7、在y=vy面上的切面PZ vz輸出模型在z=vz面上的切面FACTOR m將輸出圖放大1/m倍Extent a b切面沿兩坐標(biāo)軸方向分別放大ORIGIN X Y Z定義畫(huà)圖中心位置（X,Y,Z）三、實(shí)驗(yàn)內(nèi)容1、學(xué)習(xí)MCNP程序常見(jiàn)各種運(yùn)行方法；2、編寫(xiě)以下問(wèn)題的輸入文件；2.1對(duì)課堂講解的實(shí)例，模擬溴化鑭探測(cè)器對(duì)點(diǎn)源的能譜，實(shí)驗(yàn)做一遍。2.2有一HPGe探測(cè)器，結(jié)構(gòu)如圖1所示。分別給出位于探測(cè)器軸心、距離探測(cè)器晶體中心25cm處的137Cs源、60Co、131I源對(duì)應(yīng)特征射線的探測(cè)效率（計(jì)算時(shí)相應(yīng)特征射線的源粒子至少為107個(gè)），并給出三者混合源（活度比為1:1:3）的能譜圖（源發(fā)射總粒子數(shù)大

8、于3×108個(gè)）。圖1 HPGe探測(cè)器結(jié)構(gòu)圖四、實(shí)驗(yàn)報(bào)告1、給出2.1和2.2的MCNP輸入文件并解釋每一行的含義；2、分別運(yùn)行實(shí)例，給出實(shí)驗(yàn)結(jié)果，并對(duì)結(jié)果進(jìn)行分析。實(shí)驗(yàn)四 MCNP模擬計(jì)算射線造成的劑量一、實(shí)驗(yàn)?zāi)康?、掌握應(yīng)用軟件MCNP、應(yīng)用范圍以及在輻射劑量計(jì)算和防護(hù)中的作用；2、進(jìn)一步掌握MCNP程序基本用法；3、利用MCNP解決一個(gè)簡(jiǎn)單的求解射線在空氣、組織等效材料(肌肉)中造成的劑量沉積的計(jì)算問(wèn)題，并進(jìn)行結(jié)果分析，得出結(jié)論；4、利用MCNP程序解決實(shí)際工作中碰到的實(shí)際問(wèn)題；二、實(shí)驗(yàn)內(nèi)容1、學(xué)習(xí)MCNP程序的基本組成、操作方法以及問(wèn)題描述文件的寫(xiě)法；2、利用MCNP程序計(jì)算

9、簡(jiǎn)單的射線源在空氣、肌肉模型中的劑量沉積分布，并對(duì)計(jì)算結(jié)果進(jìn)行分析并繪圖，得出結(jié)論，調(diào)整數(shù)據(jù)重新計(jì)算，并與理論計(jì)算結(jié)果進(jìn)行比較；三、內(nèi)容簡(jiǎn)介1、MCNP程序的計(jì)算流程如下圖1所示：圖1 MCNP程序的計(jì)算流程2、MCNP輸入文件inp的基本結(jié)構(gòu) 通過(guò)這個(gè)文件描述并建立一個(gè)蒙特卡羅計(jì)算問(wèn)題，對(duì)問(wèn)題的幾何結(jié)構(gòu)、材料、記數(shù)要求等給以描述，如果需要，便可直接運(yùn)行。該文件的格式如下：(optional可選) 信息塊 空行分隔符 標(biāo)題卡 柵元卡1 柵元卡2 。 柵元卡n 空行分隔符 曲面卡1 曲面卡2 。 曲面卡n 空行分隔符 數(shù)據(jù)卡1 數(shù)據(jù)卡2 。 數(shù)據(jù)卡n 空行分隔符 (optional) 其它選擇項(xiàng)

10、 (optional)其中柵元卡用來(lái)描述由不同的封閉曲面分割的立體空間區(qū)域，并用獨(dú)有的數(shù)字ID號(hào)加以標(biāo)示，同時(shí)在各個(gè)柵元卡中說(shuō)明包圍該區(qū)域的曲面類型(曲面卡)、填充該區(qū)域的材料類型(材料卡)以及對(duì)應(yīng)的材料密度等；曲面卡是用來(lái)描述不同類型曲面的，并用獨(dú)有的數(shù)字ID號(hào)加以標(biāo)示，最終曲面卡被應(yīng)用在柵元卡中，并利用交(與)、聯(lián)(或)、補(bǔ)(非)這些邏輯運(yùn)算符號(hào)聯(lián)合不同曲面組成所需要的復(fù)雜的柵元。在mcnp中支持的常見(jiàn)曲面類型見(jiàn)參考文獻(xiàn)3,4。數(shù)據(jù)卡類型很多，主要有粒子類型標(biāo)識(shí)卡mode、重要性卡imp、通用源卡sdef、粒子計(jì)數(shù)器卡Fn、材料描述卡Mn以及粒子截?cái)嗫?nps或ctme)等，數(shù)據(jù)卡的類型涉

11、及到了方方面面，類型很多，具體請(qǐng)見(jiàn)參考文獻(xiàn)3,4。圖3 例子的幾何示意圖下面利用一個(gè)簡(jiǎn)單的例子來(lái)配合說(shuō)明mcnp中的輸入卡(inp)的編寫(xiě)格式。3、一個(gè)簡(jiǎn)單的說(shuō)明例子為說(shuō)明如何填寫(xiě)INP文件，這里例舉一個(gè)簡(jiǎn)單問(wèn)題。如圖3所示，在一個(gè)邊長(zhǎng)10cm的石墨立方體3中有兩個(gè)半徑0.5cm的球形空間，球1中充滿氧氣，球2是鐵球。在球1中置一14MeV各向同性中子點(diǎn)源，計(jì)算球2外表面與能量相關(guān)的中子通量。建立的INP文件如下：SAMPLE PROBLEM INPUT DECK1 1 0.0014 -72 2 7.86 -83 3 1.60 1 2 3 4 5 6 7 84 0 -1:2:3:-4:5:-6

12、¬空行1 PZ -52 PZ 53 PY 54 PY -55 PX 56 PX -57 S 0 4 2.5 .58 S 0 4 4 .5¬空行MODE pIMP:p 1 1 1 0SDEF POS=0 4 2.5 ERG=14 F2:n 8E0 1E-5 1E-4 1E-3 .01 .1 1 2 3 4 5 6 7 8 9 10 11 12 13M1 8016 1M2 26000 1M3 6000 1NPS 100000¬空行本例中沒(méi)有信息塊，第一行是標(biāo)題卡，之后至空格前為柵元塊。柵元卡上依次填寫(xiě)柵元號(hào)、材料號(hào)、密度和構(gòu)成柵元界面的曲面號(hào)(帶正負(fù)號(hào))，這里定義了4

13、個(gè)柵元：柵元1由球面7圍成，里面填充材料1(16O2氣體)，密度是0.0014g/cm3；柵元2由球面8圍成，填充材料2(鐵)，密度7.86g/cm3；柵元3由平面1、2、3、4、5、6圍成，不包括球面7、8以內(nèi)的空間，填充材料3(石墨)，密度1.6g/cm3；柵元4是柵元3以外的空間，為真空。曲面卡上需要填寫(xiě)曲面號(hào)、曲面類型和曲面參數(shù)，本例中定義了8個(gè)曲面，前6個(gè)為與原點(diǎn)距離5cm垂直于各坐標(biāo)軸的平面，后兩個(gè)是半徑0.5cm的球面，球心分別在(0,-4,-2.5)和(0,4,4)。數(shù)據(jù)塊中指定了問(wèn)題類型、源、記數(shù)方式、材料和運(yùn)行粒子數(shù)，各卡數(shù)據(jù)項(xiàng)的意義如下：MODE卡問(wèn)題類型是中子輸運(yùn)IMP

14、卡4個(gè)柵元的中子重要性分別是 1 1 1 0SDEF卡位于(0,-4,-2.5)、能量14MeV的各向同性點(diǎn)源F2卡在曲面8上做中子通量記數(shù)E0卡對(duì)記數(shù)能量分區(qū)，114MeV之間間隔為1MeV，1MeV10-5MeV之間間隔為一個(gè)數(shù)量級(jí)M1卡材料1是16O核素M2卡材料2是Fe元素M3卡材料3是C元素NPS卡運(yùn)行源粒子數(shù)100000以上例子僅用于說(shuō)明INP文件格式，有關(guān)各輸入卡的詳細(xì)內(nèi)容，具體使用方法見(jiàn)參考文獻(xiàn)3,4。四、實(shí)驗(yàn)步驟1、編寫(xiě)對(duì)應(yīng)任務(wù)的輸入文件，如實(shí)例Dose.inp。2、運(yùn)行實(shí)例Dose，得到計(jì)算結(jié)果文件Dose.o，對(duì)Dose.o文件進(jìn)行分析，了解并熟悉計(jì)算過(guò)程，得到并分析計(jì)算

15、結(jié)果；3、理論計(jì)算空氣比釋動(dòng)能值，并與上述計(jì)算結(jié)果計(jì)算得到的空氣比釋動(dòng)能值進(jìn)行比較，分析差異，給出分析結(jié)果。4、在上述例子的基礎(chǔ)上，對(duì)計(jì)算空間進(jìn)行柵元細(xì)分，得到更為精細(xì)的計(jì)算結(jié)果，并按照前面順序?qū)τ?jì)算結(jié)果進(jìn)行劑量沉積分析。五、問(wèn)題求解1、問(wèn)題1個(gè)1Ci的Cs-137源， 計(jì)算距離其30cm處的空氣比釋動(dòng)能K、肌肉材料的吸收劑量D。2、設(shè)計(jì)要求：（1） 所設(shè)計(jì)的源為各向同性通用源：即用SDEF卡設(shè)置源；（2） 計(jì)算光子在肌肉和空氣中所沉積的能量；（3） 所描述的幾何有一定數(shù)量的柵元，幾何形狀可以按照自己的興趣選擇；3、光子作用的材料肌肉材料如下：1001. -0.102 6000. -0.143 7014. -0.034 8016. -0.710 11023. -0.001 15031. -0.002 16032. -0.003 17000. -0.001 19000. -0.004空氣材料如下 6000. -0.000124 7014. -0.7553 8016. -0.2318 18000. -0.01289894、輸入卡的編寫(xiě)5、運(yùn)行結(jié)果六、思考問(wèn)題1、什么是