第五章蒙特卡羅方法在計(jì)算機(jī)上實(shí)現(xiàn)

1、l源分布抽樣過(guò)程源分布抽樣過(guò)程l空間、能量和運(yùn)動(dòng)方向的隨機(jī)游動(dòng)過(guò)程空間、能量和運(yùn)動(dòng)方向的隨機(jī)游動(dòng)過(guò)程l記錄貢獻(xiàn)和分析結(jié)果過(guò)程記錄貢獻(xiàn)和分析結(jié)果過(guò)程l核截面數(shù)據(jù)的引用核截面數(shù)據(jù)的引用l蒙特卡羅程序結(jié)構(gòu)蒙特卡羅程序結(jié)構(gòu)作作 業(yè)業(yè) 蒙特卡羅方法是隨著計(jì)算機(jī)的出現(xiàn)和發(fā)展而逐步發(fā)展起來(lái)的。在計(jì)算機(jī)上能夠產(chǎn)生符合要求的隨機(jī)數(shù)，實(shí)現(xiàn)對(duì)已知分布的抽樣，奠定了蒙特卡羅方法在計(jì)算機(jī)上得以實(shí)現(xiàn)的基礎(chǔ)。在計(jì)算機(jī)上使用蒙特卡羅方法解粒子輸運(yùn)問(wèn)題大致包括三個(gè)過(guò)程：源分布抽樣過(guò)程，空間、能量和運(yùn)動(dòng)方向的隨機(jī)游動(dòng)過(guò)程以及記錄、分析結(jié)果過(guò)程 。 源分布抽樣的目的是產(chǎn)生粒子的初始狀態(tài) 。下面我們介紹一些常見的特定類型的源分布抽樣

2、方法。),(0000rSEl 圓內(nèi)均勻分布設(shè)圓半徑為R0，粒子在圓內(nèi)均勻分布時(shí)，從發(fā)射點(diǎn)到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)002)(020RrRrrf),max(210 Rrl 圓環(huán)內(nèi)均勻分布設(shè)圓環(huán)的內(nèi)半徑為R0，外半徑為R1，則粒子在該圓環(huán)內(nèi)均勻分布時(shí)，從發(fā)射點(diǎn)到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)02)(102021RrRRRrrf02010320101011)(),max()(RRRrRRRrRRRRl 球內(nèi)均勻分布設(shè)球的半徑為R，粒子在球內(nèi)均勻分布時(shí)，從發(fā)射點(diǎn)到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：在直角坐標(biāo)系下，抽樣

3、方法為：其它當(dāng)003)(32RrRrrf),max(321 Rr302010232221,1RzRyRxl 球殼內(nèi)均勻分布設(shè)球殼的內(nèi)半徑為R0，外半徑為R1，在均勻分布時(shí)，從發(fā)射點(diǎn)到中心的距離 r 的分布密度函數(shù)為： r 的抽樣方法為：其它當(dāng)03)(1030312RrRRRrrf001432322211020101211020201)(),max(),max(33RxRRrxxxRRRRRRRRRRR在直角坐標(biāo)系下，球殼內(nèi)點(diǎn)的坐標(biāo)為： 其中，r 由前面的抽樣方法確定，、服從各向同性分布，其抽樣方法為：cossinsincossin000rzryrx23222221232222210232222

4、21310232222212101223222221cos2sinsin2cossin)(AAAArrzAAArryAAArrxAAl 圓柱內(nèi)均勻分布圓柱內(nèi)均勻分布是指粒子發(fā)射點(diǎn)均勻地分布在底半徑為 R，高為 2H 的圓柱內(nèi)。若固定圓柱的中心為原點(diǎn)，圓柱的軸向?yàn)?z 軸，則分布密度函數(shù)為： 抽樣方法為：其它當(dāng)0| ,21),(2222HzRyxHRzyxf3020102221,1HzRyRxl 點(diǎn)源分布 點(diǎn)源分布是指粒子由一固定點(diǎn) 發(fā)射，其分布密度函數(shù)為：其中， 為狄拉克函數(shù)，源粒子的抽樣方法為：在球坐標(biāo)系中，粒子發(fā)射點(diǎn)到球心的距離 r 的分布密度函數(shù)為：其中， 為點(diǎn)源到球心的距離。源粒子的位

5、置抽樣為：)()()(),(*0*0*0zzyyxxzyxf),(*0*0*0zyx*0*0*0,zzyyxx*0rr *0r)()(*0rrrf)(l 球外平行束源分布球外平行束源分布是指粒子平行入射到半徑為 R 的球面上，或球外點(diǎn)源距離球很遠(yuǎn)，可以近似地看作平行束源。設(shè) r 為粒子發(fā)射點(diǎn)到球心的距離 , 其分布密度函數(shù)為： r 的抽樣方法為：在直角坐標(biāo)系中，抽樣方法為：)()(RrrfRr 20202020102221,1zyRxRzRyl 單能源分布單能源分布是指粒子的發(fā)射能量為一固定值 E0 ，其分布密度函數(shù)為 ：源粒子的能量為：)()(0EEEf0EE l 裂變中子譜分布裂變中子譜分

6、布的一般形式為： 其中A，B，C，Emin，Emax 均為與元素有關(guān)的量。對(duì)于鈾-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。maxmin,sh)(EEEBEeCEfAE采用近似修正抽樣，抽樣方法為：其中，m0.8746，M10.2678，0.5543。此外，裂變譜分布也有以數(shù)值曲線形式給出的，此時(shí)，用數(shù)值曲線抽樣方法抽取 E 。3321121ln1)ln()(AAEEAEMEHmEAmEBECAAEHexpsh1)(21l 麥克斯韋(Maxwell) 譜分布麥克斯韋譜分布的一般形式為：該分布的抽樣方法為0,2)(23EeEEfE22221ln23lnEel

7、 各向同性分布各向同性分布密度函數(shù)為： 其中，cos，為運(yùn)動(dòng)方向與 z 軸的夾角，為方位角。21)(,21)(41)()()(2121fffff在直角坐標(biāo)系下，各方向余弦 u，v，w 為：其抽樣方法為：cossinsincossinwvu2322222123222221232222213123222221211223222221cos2sinsin2cossin)(AAAAwAAAvAAAuAAl 半面各向同性分布不妨設(shè)在 x0 的半面方向上各向同性發(fā)射粒子，則在前述各向同性分布的抽樣方法中，用2代替2就能得到所需分布的抽樣。對(duì)于其它方向的情況，可用類似的方法處理。l 球外平行束源分布令cos

8、，為粒子運(yùn)動(dòng)方向的徑向夾角，則分布密度函數(shù)為：的抽樣方法為：01,2)(f),max(21ul 球外各向同性點(diǎn)源分布設(shè)球外點(diǎn)源 S 到球心的距離為D0。點(diǎn)源 S 到球的最大張角為*，則球外各向同性點(diǎn)源分布的抽樣方法是：先抽樣確定 ，再轉(zhuǎn)換成。0220*cosDRD RDR2202*sincos)cos1 (1cos在直角坐標(biāo)系下，取 OS 為 z 軸，抽樣方法為：cos0sinwvu 在有關(guān)次級(jí)粒子（如裂變中子，中子生成光子，光子生成中子）的輸運(yùn)過(guò)程中，次級(jí)粒子源分布的抽樣方法，主要可分為以下兩種：l 直接生成法 可將生成的次級(jí)粒子的位置、能量、方向、權(quán)重等參數(shù)直接作為源分布的抽樣結(jié)果。也就是

9、直接對(duì)生成的次級(jí)粒子進(jìn)行跟蹤。這種方法比較簡(jiǎn)單、直觀。l 離散分布法 將生成的次級(jí)粒子的權(quán)重，按空間位置、能量、方向分別記錄，得到次級(jí)粒子的空間、能量、運(yùn)動(dòng)方向的離散的近似分布。再根據(jù)該分布，利用各種抽樣技巧，得到源分布的抽樣，對(duì)抽樣的源粒子進(jìn)行跟蹤、記錄。 當(dāng)一個(gè)問(wèn)題需要用兩個(gè)以上的蒙卡程序處理時(shí)，可采用這種方法。 粒子由狀態(tài)Sm到狀態(tài)Sm+1時(shí)，需要確定粒子的空間位置 rm+1，能量 Em+1和運(yùn)動(dòng)方向m+1。 設(shè) rm 為粒子第 m 次碰撞點(diǎn)的位置，m 為碰撞后的運(yùn)動(dòng)方向，則粒子第 m+1 次碰撞點(diǎn)的位置 rm+1 為：即其中 為 的方向余弦，L 為兩次碰撞點(diǎn)間的距離。mmmmmmmmm

10、wLzzvLyyuLxx111mmmL rr1mmmmwvu),(L 的分布密度函數(shù)為：由 f (L) 抽樣確定 L 的方法通常有三種：l 直接抽樣方法確定 L 的直接抽樣方法是：首先由自由程分布中抽取(1)再由下列關(guān)系式解出 L 。0,),(exp),()(01LdlElELfLmmmtmmtrrLmmmtdlEl0),(r ef)(ln對(duì)于均勻介質(zhì)，有對(duì)于多層介質(zhì)，如果則其中， 為粒子由 rm 出發(fā)，沿m 方向在順序經(jīng)過(guò)的第 i 個(gè)介質(zhì)區(qū)域內(nèi)走過(guò)的距離， 為第 i 個(gè)介質(zhì)區(qū)域的宏觀總截面 ( i =1，2，Imax )。 當(dāng)時(shí)，意味著粒子穿出系統(tǒng)。)(ln)(mtmtEEL)()()()(

11、,10,100,10,mItIimitiIiiIimitiIimitiEELLLELELiL)(,mitEmax0,)(IimitiELl 最大截面法對(duì)于多層介質(zhì)，或其他介質(zhì)密度與位置有關(guān)的問(wèn)題，在求 ( i =1，2，Imax ) 時(shí)，如果系統(tǒng)形狀復(fù)雜，計(jì)算是非常煩雜的。在這種情況下，使用最大截面法更方便。最大截面抽樣方法為：其中iL),(max)(max,EEttrr1max,12max,1111)(),()(ln0LLEELELLLmtmmmtmtrl 限制抽樣法 當(dāng)介質(zhì)區(qū)域很小時(shí)，如使用直接抽樣法抽取輸運(yùn)長(zhǎng)度，粒子很容易穿出介質(zhì)，此時(shí)使用限制抽樣法確定自由程個(gè)數(shù)較好，的分布密度函數(shù)為：

12、其中 Dm 為粒子由rm 出發(fā)，沿m 方向到達(dá)區(qū)域邊界的自由程個(gè)數(shù)。的抽樣方法是：然后用直接抽樣法中根據(jù)計(jì)算 L 的方法計(jì)算輸運(yùn)長(zhǎng)度 L 。此時(shí)，粒子的權(quán)重需乘以糾偏因子 。其它當(dāng)001)(mDDeefm)1 (1lnmDe)1 (mDe 粒子在介質(zhì)中發(fā)生碰撞后，首先要確定與哪種原子核發(fā)生何種反應(yīng)。粒子發(fā)生碰撞后（吸收除外）的能量 Em+1 一般只與其碰撞前后運(yùn)動(dòng)方向的夾角（散射角）有關(guān)。 粒子碰撞后常見的能量分布有下面幾種情況。l 裂變中子譜(1) 中子引起原子核裂變反應(yīng)時(shí)，裂變中子的能量服從裂變譜分布。其抽樣方法可參考以前的介紹。l 中子彈性散射后能量的確定 中子彈性散射后，能量與質(zhì)心系散

13、射角C的關(guān)系是：能量與實(shí)驗(yàn)室系散射角L的關(guān)系是：其中，A 為碰撞核的質(zhì)量， ?；?確定后，即可求出 Em+1。222211) 1(LLmmAAEE12) 1(221CmmAAAEELLCCcos,cosLCl 中子非彈性散射后能量的確定 中子非彈性散射后，能量與質(zhì)心系散射角C的關(guān)系是：其中， 為第 K 個(gè)能級(jí)的閾能， 為第 K 個(gè)能級(jí)的激發(fā)態(tài)能量。 如果確定了實(shí)驗(yàn)室系散射角L，則根據(jù)下式確定 后，再計(jì)算出 Em+1。KKCmKmKmmAAEAEAAEE11121) 1(221KKC11111222LLmKLmKCEAEAl 光子康普頓（Compton）散射后能量的確定光子發(fā)生康普頓散射后，其能

14、量分布密度函數(shù)為：其中， K() 為歸一因子。 ， 和 分別為光子散射前后的能量，以 m0c2 為單位，m0為電子靜止質(zhì)量，c 為光速。211,1111)(1)(322xxxxxxKxf22)21 (21421)21ln() 1(21)(Kx光子康普頓散射能量分布的抽樣方法為：x 的抽樣確定后，散射后的能量為：232223221121321) 1(427212942711212121xxxxxxxxxxEcmxcmEmm 20201粒子碰撞后運(yùn)動(dòng)方向m+1的確定，一般與散射角有關(guān)。由已知分布抽樣確定散射角后，再確定m+1。常見的散射角分布有如下幾種：l 質(zhì)心系各向同性分布散射角在質(zhì)心系服從各向

15、同性分布時(shí)，其抽樣方法為 。質(zhì)心系散射角C抽樣確定后，(1)需轉(zhuǎn)換成實(shí)驗(yàn)室系散射角L。12 C在中子彈性散射情況下，轉(zhuǎn)換公式為：其中 A 為碰撞核質(zhì)量， 。在中子非彈性散射情況下，轉(zhuǎn)換公式為：其中， 為第 K 個(gè)能級(jí)的閾能。CCLAAA2112CmKmKCmKLEAEAEA1211112LLCCcos,cosKl 中子彈性散射勒讓德 (Legendre) 多項(xiàng)式分布 中子彈性散射角分布常以勒讓德多項(xiàng)式的展開形式給定。散射角余弦 xcos的分布密 度函數(shù)為：其中 Pl(x) 為 l 階勒讓德多項(xiàng)式。該分布即為 n 階勒讓德近似展開。勒讓德多項(xiàng)式由以下遞推公式確定：其它當(dāng)01|)(212)(0 x

16、xPflxfnlllxxPxPxnPxxPnxPnnnn)(1)(0)()() 12()() 1(1011考慮新的分布：當(dāng)選取 x0，x1， xn 為 Pn+1(x)0 的根，且時(shí)，fa(x) 依照勒讓德多項(xiàng)式展開的前 n 項(xiàng)與 f (x) 的展開形式相同。因此，可以用 fa(x) 作為 f (x) 的近似分布。nlklnlkllkxPlxPfl020)(212)(212nkkkaxxxf0)()(在實(shí)際問(wèn)題中，由于勒讓德多項(xiàng)式展開項(xiàng)數(shù)不夠，可能出現(xiàn)某個(gè) 為負(fù)值的現(xiàn)象。此時(shí)可以采用如下近似分布：其中：對(duì)于該近似分布，可用加抽樣方法進(jìn)行抽樣：此時(shí)，由于偏倚抽樣而引起的糾偏因子為 wK ，也就是說(shuō)

17、，粒子的權(quán)重要乘上wK。nkkkkknkkkknkkka|w|xxxf000*|,|)()(kKkkKkkKxx010,當(dāng)l 光子康普頓散射角分布光子的康普頓散射角與其散射前后的能量有關(guān) , 它的分布密度函數(shù)為：抽樣方法為：)111 ()(LLf111L實(shí)驗(yàn)室系散射角L確定后，依據(jù)不同的坐標(biāo)系的表現(xiàn)形式，有不同的確定方法。l 確定方向余弦 um+1，vm+1，wm+1mmmmmmmmmmmmmmmmmmawvubcwavvubduvbcwvauvubdvubcwu221221221其中，方位角 在 0, 2 上均勻分布。 當(dāng) 時(shí)，不能使用上述公式，可用下面的簡(jiǎn)單公式：sin,cos,1sin,

18、cos2dcabaLL022mmvummmmawwbdvbcu111l 確定m+1的球坐標(biāo) (m+1，m+1)設(shè)m的球坐標(biāo)分別為 (m,m)，其中，為粒子運(yùn)動(dòng)方向與 z 軸的夾角， 為粒子運(yùn)動(dòng)方向在 x y 平面上投影的方位角。則m+1的球坐標(biāo) (m+1,m+1) 分別由下式確定：111111sinsincoscoscos)cos(sinsinsin)sin(cossinsincoscoscosmmmmLmmmLmmLmLmm 一般幾何的隨機(jī)游動(dòng)公式可以應(yīng)用到球形幾何，而對(duì)球?qū)ΨQ問(wèn)題，使用特殊形式更為方便。l 下次碰撞點(diǎn)的徑向位置 rm+1的確定 兩次碰撞點(diǎn)間的距離 L 確定之后，下次碰撞點(diǎn)的

19、徑向位置 rm+1的計(jì)算公式為：設(shè)系統(tǒng)的外半徑為R，如 rm+1R，則粒子逃出系統(tǒng)。mmmmLLcos2221rrrl 粒子碰撞后瞬時(shí)運(yùn)動(dòng)方向的確定 在球?qū)ΨQ系統(tǒng)中，粒子運(yùn)動(dòng)方向用其與徑向夾角余弦來(lái)描述。使用球面三角公式，粒子碰撞后瞬時(shí)運(yùn)動(dòng)方向與徑向夾角余弦 cosm+1的計(jì)算公式為：其中， 為在 0, 2 上均勻分布的方位角， 為在 rm+1 點(diǎn)進(jìn)入碰撞前瞬時(shí)運(yùn)動(dòng)方向與 rm+1 徑向之間的夾角。111coscossinsincossinsincoscoscosmmmmmmmmLmLmmLrrrrm在抽樣確定輸運(yùn)距離、判斷粒子是否穿透系統(tǒng)時(shí)， 常遇到求由 rm 出發(fā)，沿m 方向到達(dá)某個(gè)區(qū)域表

20、面的距離問(wèn)題。在記錄對(duì)結(jié)果的貢獻(xiàn)時(shí)，也常使用類似的量。區(qū)域表面通常是平面或二次曲面。 求到達(dá)區(qū)域表面的距離問(wèn)題，實(shí)際上是求直線（或半射線）與平面或二次曲面的交點(diǎn)問(wèn)題。這是 蒙特卡羅方法解粒子輸運(yùn)的各種實(shí)際問(wèn)題時(shí) , 所遇到的基本幾何問(wèn)題。l 點(diǎn)到平面的距離點(diǎn) 沿方向 的直線方程為：該直線到達(dá)方程為的平面的距離為：當(dāng) 與平面平行時(shí)，即直線與平面無(wú)交點(diǎn)。如果 l 為負(fù)值，直線與平面也無(wú)交點(diǎn)。這時(shí)，粒子的運(yùn)動(dòng)方向是背離平面的。),(0000wvu0000cwbvau000000)(cwbvauczbyaxdldczbyax00rrl),(0000zyxrl 點(diǎn)到球面的距離在三維直角坐標(biāo)系中，設(shè)球心為

21、 rc(xc, yc, zc) ，球半徑為 R，則球面方程為：將直線方程代入該球面方程，得到點(diǎn) r0沿0方向到達(dá)球面的距離 l ：其中2222)()()(Rzzyyxxcccl2020202020220200000000)()()(|)()()()(cccccccczzyyxxRRRwzzvyyuxxrrrr當(dāng) R0R 時(shí)，即 r0 點(diǎn)在球內(nèi) ，2，l 只有一個(gè)正根：當(dāng) R0R 時(shí)，即 r0 點(diǎn)在球外，分以下三種情況：l 若0，l 無(wú)正實(shí)根，直線與球面無(wú)交點(diǎn)。l 若0，0，l 無(wú)實(shí)根，直線與球面無(wú)交點(diǎn)。a) 若0，0，l 有兩個(gè)正實(shí)根，直線與球面有兩個(gè)交點(diǎn)。ll在球坐標(biāo)系中，不失一般性，設(shè)球心

22、為 rc0，則球面方程為 rR。當(dāng) r0R 時(shí)，即 r0 點(diǎn)在球內(nèi) ，有一個(gè)交點(diǎn)：其中0為0與 r0 的徑向夾角。當(dāng) r0R 時(shí)，即 r0 點(diǎn)在球外 ，令當(dāng) cos00 時(shí)，直線與球面無(wú)交點(diǎn)。當(dāng) cos00 時(shí)，若 dR，則直線與球面無(wú)交點(diǎn)。若 dR，則有兩個(gè)交點(diǎn)：200200)sin(cosrRrl200200)sin(cosrRrl00sin rdl 點(diǎn)到圓柱面的距離設(shè)圓柱面的方程為：其中 (xc, yc, 0) 為圓柱的中心，R 為圓柱底半徑。點(diǎn) r0沿0方向到達(dá)圓柱面的距離 l 為：其中222)()(Ryyxxcc201wl2020202022020000)()()1)()()(ccc

23、cyyxxRwRRvyyuxx 當(dāng) R0R 時(shí)，r0 點(diǎn)在圓柱內(nèi) ，如果 ，則 l 有一個(gè)正根： 如果 ，即0平行于圓柱的對(duì)稱軸，直線與圓柱面無(wú)交點(diǎn)。 當(dāng) R0R 時(shí)，r0 點(diǎn)在圓柱外，分以下三種情況：l 若0，l 無(wú)正實(shí)根，直線與圓柱面無(wú)交點(diǎn)。l 若0，0，l 無(wú)實(shí)根，直線與圓柱面無(wú)交點(diǎn)。l 若0，0 且 ，l 有兩個(gè)正實(shí)根，直線與圓柱面有兩個(gè)交點(diǎn)。l在 的情況下，直線與圓柱面不相交。201wl201wl120w120w120w120wl 點(diǎn)到圓錐面的距離 設(shè)圓錐頂點(diǎn)在原點(diǎn)，以 z 軸為對(duì)稱軸，則圓錐面的方程為：點(diǎn) r0沿0方向到達(dá)圓錐面的距離 l 為：其中如果0與錐面某一母線平行，即 ，則

24、2222zcyx202)1 (1wcl)1 (1)(202202202020020000wczcyxwzcvyux22022020zcyxl1)1 (202wcl 空腔處理 在粒子輸運(yùn)問(wèn)題中，所考慮的系統(tǒng)常有空腔存在，如中空的球殼 , 平板間的空隙等。粒子輸運(yùn)時(shí)，可有兩種處理空腔的方法：l 將空腔作為宏觀總截面t0 的區(qū)域 , 按通常的方法輸運(yùn)。l 設(shè) 分別為由 rm 出發(fā)，沿m 方向到空腔區(qū)域的近端和遠(yuǎn)端的交點(diǎn)，當(dāng)粒子超過(guò) 時(shí)，以 為新的起點(diǎn)，重新開始輸運(yùn)。 顯然，這兩種方法在統(tǒng)計(jì)上是等價(jià)的。mmrr 、mrmr l 全反射邊界在反應(yīng)堆活性區(qū)中，元件盒常常按正方形或六角形排列。假定元件盒足夠

25、多，每個(gè)盒結(jié)構(gòu)相同，那么活性區(qū)中每個(gè)盒所占的柵胞的物理情況，可以代表整個(gè)活性區(qū)中的狀況。 進(jìn)一步假定，元件盒是圓對(duì)稱的，那么每個(gè)柵胞中情況，可以用更小的單位（柵元）來(lái)反映。比如對(duì)六角形柵胞可取其 1/12 的OAB 來(lái)做代表；正方形柵胞可用其 1/8 的OAB 來(lái)做代表。這樣一來(lái)問(wèn)題就大大簡(jiǎn)化了。 現(xiàn)在的問(wèn)題是怎樣計(jì)算直角三角形柵元的物理量（如通量）。用蒙特卡羅方法如何模擬中子在柵元內(nèi)的運(yùn)動(dòng)，反映出整個(gè)活性區(qū)對(duì)它的影響。我們可把OA、OB、AB 作為全反射邊界來(lái)處理。所謂全反射邊界，就是當(dāng)中子打到該邊界上時(shí)，按鏡面反射的方式，從邊界 上全部反射回來(lái)，中子的能量與權(quán)重均不改變。 在這種邊界上的反

26、射條件，稱之為全反射條件，就是通常的鏡面反射條件。 在全反射邊界條件下，一條通過(guò)活性區(qū)若干個(gè)區(qū)域的中子徑跡，可以用柵元OAB 中的一條折線軌 道來(lái)反映出來(lái)。 反過(guò)來(lái)，在直角三角形柵元OAB 中任一條反射成的折線軌道，都代表了中子在活性區(qū)內(nèi)一條直線軌道的作用。由于系統(tǒng)的對(duì)稱性，在活性區(qū)內(nèi)，凡是與柵元內(nèi)位置相當(dāng)?shù)牡胤?，都有相同的物理情況，因此柵元內(nèi)各處的情況，當(dāng)然代表了整個(gè)活性區(qū)的情況。l 一般曲面全反射條件 對(duì)于一般曲面的全反射，設(shè)入射方向?yàn)?，入射點(diǎn)的內(nèi)法線方向?yàn)?n ，則反射方向 為：其中設(shè)則nnncos22ncos),(zyxnnnn)(coscos2cos2cos2zyxzyxnwnvnu

27、nwwnvvnuul 平面全反射條件 設(shè)三角形柵元的橫截面OAB 在 X-Y 平面上，OAB。則邊界 OA、OB、AB 上的反射都是平面全反射。在任一與 X-Y 平面垂直且與 X 軸成角的平面上，全反射條件為：由此就可得到OA、OB 和 AB 邊上的全反射條件，對(duì)于 OB 邊，=；對(duì)于 OA 邊，= 0；對(duì)于 AB 邊，=/ 2。wwvuvvuu2cos2sin2sin2cosl 反射層邊界條件 對(duì)于具有大反射層的系統(tǒng)，如存放，運(yùn)輸和生產(chǎn)裂變物質(zhì)的倉(cāng) 庫(kù)、車廂和車間等，當(dāng)中子從里面打到四周墻上或反射層時(shí)，還要繼續(xù)對(duì)它進(jìn)行跟蹤。這種跟蹤常常要花費(fèi)很大的計(jì)算量，并且在結(jié)果中引起的方差也比較大。如果

28、在計(jì)算這種系統(tǒng)的不同方案中，反 射層條件不變，那么這種大量重復(fù)的計(jì)算是很不經(jīng)濟(jì)的。 中子射入反射層后，一部分被介質(zhì)吸收，只有一部分返回，由于中子的散射慢化，損失一部分能量，因此反射回來(lái)的中子有一個(gè)能量方向分布。顯然，對(duì)這種反射層，不能應(yīng)用全反射條件。不過(guò)，我們?nèi)匀豢梢园阉?dāng)做邊界，在邊界上按反射層的物理作用來(lái)處理。 比如，如果反射層是一種平板幾何，我們可以用數(shù)值方法或蒙特卡羅方法，預(yù)先算好在各種不同入射能量 E 下的反照率(E)，反射中子的能量分布 RE(EE )。于是代替在反射層中眼蹤中子，我們可在反射層邊界上作如下處理：一旦中子打入反射 層，立即返回，反射后權(quán)重為其中，E 為射入反射層中子

29、的能量，W 為中子的權(quán)重。反射后的能量 E 由反射能譜 RE(EE) 中抽樣產(chǎn)生。反射后的方向 由半平面各向同性分布或余弦分布中抽樣。反射后的中子位置為入射時(shí)的位置。 計(jì)算表明，對(duì)于大尺寸的反射層來(lái)說(shuō)，這樣的近似，引 起的結(jié)果上的誤差是可以忽略的，卻能帶來(lái)計(jì)算量的大量節(jié)省。)(EWW 在粒子輸運(yùn)問(wèn)題中，除了得到某些量的積分結(jié)果外，還需要得到這些量的方差、協(xié)方差、以及這些量的空間、能量、方向和時(shí)間的分布。這些量可以利用分類記錄手續(xù)同時(shí)得到。 為了得到所求量的估計(jì)值，在粒子輸運(yùn)過(guò)程中需進(jìn)行記錄，即求每個(gè)粒子對(duì)所求量的貢 獻(xiàn)。 設(shè)模擬了 N 個(gè)粒子，所求量的估計(jì)值為：其中 gi 為第 i 個(gè)粒子的總

30、貢獻(xiàn)。NiiNgNg11 記錄的貢獻(xiàn)由所求量決定。對(duì)于同一個(gè)所求量，又隨所用的蒙特卡羅技巧的不同而不同。 例如，所求量是粒子穿透屏蔽概率，使用直接模擬法時(shí)，如粒子穿透屏蔽，在疊加記錄單元加“1” ( 初始值為零 )，否則沒有貢獻(xiàn)。使用加權(quán)法時(shí)，如粒子穿透屏蔽，在疊加記錄單元加粒子的權(quán)重，否則沒有貢獻(xiàn)。使用統(tǒng)計(jì)估計(jì)法時(shí)，粒子每發(fā)生一次碰撞 (包括零次碰撞)，都要記錄貢獻(xiàn)，等等。估計(jì)量 g 和 g 的方差和協(xié)方差為：它們可以用下式估計(jì)：NiiNiiNiiiggNiiNiiggNgNggNgNgN11122112211111 gEEgggEEggEggg2222因此，要得到 和 的估計(jì)，只要對(duì)每一個(gè)歷史記錄結(jié)果的 和 進(jìn)行記錄，并加以累加即可。方差估計(jì)值 確定后，可得到誤差其中 為置信限，它隨置信水平 而定。在通常情況下，取 。22gggiiiggg22gNg96. 1,95. 011 在前面已經(jīng)提到，用蒙特卡羅方法求某種量的空間、能量、方向和時(shí)間分布，實(shí)質(zhì)上是得到這種分布的階梯函數(shù)近似的估計(jì)值。而求這種估計(jì)值是很方便的，只要將跟蹤過(guò)程中所得到的感興趣量，按其狀態(tài)的空間、能量、方向、時(shí)間特征，分別記錄其權(quán) 重，最后將這些記錄結(jié)果適當(dāng)處理即可。 事先，將問(wèn)題的空間、能量、方向（常按相對(duì)于某個(gè)方向的夾角余弦）、時(shí)間