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

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

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

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

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

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

6、圖、試驗累積頻率與理論累積頻率關(guān)系圖，并給出抽樣次數(shù)（>106）與抽樣時間。實驗三 MCNP方法在實驗核物理中的應(yīng)用一、實驗?zāi)康?、了解MCNP程序運行流程；2、掌握MCNP輸入文件編寫規(guī)范；3、理解模擬內(nèi)容、并能編寫輸入文件、運行，并獲得計算結(jié)果；二、實驗原理MCNP是一種常見的粒子輸運模擬軟件，軟件的安裝、運行和輸入文件編寫方法詳見相關(guān)參考資料。MCNP輸入文件編寫完成后，先確認(rèn)輸入模型是否正確，在DOS環(huán)境下進(jìn)行，打開運行DOS環(huán)境，進(jìn)行以下操作：DOS命令操作命令含義Mcnp i=name.inp o=name.o打開畫圖框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定義畫圖中心位置（X,Y,Z）三、實驗內(nèi)容1、學(xué)習(xí)MCNP程序常見各種運行方法；2、編寫以下問題的輸入文件；2.1對課堂講解的實例，模擬溴化鑭探測器對點源的能譜，實驗做一遍。2.2有一HPGe探測器，結(jié)構(gòu)如圖1所示。分別給出位于探測器軸心、距離探測器晶體中心25cm處的137Cs源、60Co、131I源對應(yīng)特征射線的探測效率（計算時相應(yīng)特征射線的源粒子至少為107個），并給出三者混合源（活度比為1:1:3）的能譜圖（源發(fā)射總粒子數(shù)大

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

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

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

11、及到了方方面面，類型很多，具體請見參考文獻(xiàn)3,4。圖3 例子的幾何示意圖下面利用一個簡單的例子來配合說明mcnp中的輸入卡(inp)的編寫格式。3、一個簡單的說明例子為說明如何填寫INP文件，這里例舉一個簡單問題。如圖3所示，在一個邊長10cm的石墨立方體3中有兩個半徑0.5cm的球形空間，球1中充滿氧氣，球2是鐵球。在球1中置一14MeV各向同性中子點源，計算球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¬空行本例中沒有信息塊，第一行是標(biāo)題卡，之后至空格前為柵元塊。柵元卡上依次填寫柵元號、材料號、密度和構(gòu)成柵元界面的曲面號(帶正負(fù)號)，這里定義了4

13、個柵元：柵元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以外的空間，為真空。曲面卡上需要填寫曲面號、曲面類型和曲面參數(shù)，本例中定義了8個曲面，前6個為與原點距離5cm垂直于各坐標(biāo)軸的平面，后兩個是半徑0.5cm的球面，球心分別在(0,-4,-2.5)和(0,4,4)。數(shù)據(jù)塊中指定了問題類型、源、記數(shù)方式、材料和運行粒子數(shù)，各卡數(shù)據(jù)項的意義如下：MODE卡問題類型是中子輸運IMP

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

15、結(jié)果；3、理論計算空氣比釋動能值，并與上述計算結(jié)果計算得到的空氣比釋動能值進(jìn)行比較，分析差異，給出分析結(jié)果。4、在上述例子的基礎(chǔ)上，對計算空間進(jìn)行柵元細(xì)分，得到更為精細(xì)的計算結(jié)果，并按照前面順序?qū)τ嬎憬Y(jié)果進(jìn)行劑量沉積分析。五、問題求解1、問題1個1Ci的Cs-137源， 計算距離其30cm處的空氣比釋動能K、肌肉材料的吸收劑量D。2、設(shè)計要求：（1） 所設(shè)計的源為各向同性通用源：即用SDEF卡設(shè)置源；（2） 計算光子在肌肉和空氣中所沉積的能量；（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、輸入卡的編寫5、運行結(jié)果六、思考問題1、什么是