蒙特卡羅方法在計算機上的實現(xiàn)

第五章蒙特卡羅方法在計算機上的實現(xiàn)

1.源分布抽樣過程

2.空間、能量和運動方向的隨機游動過程

3.記錄貢獻和分析結(jié)果過程

4.核截面數(shù)據(jù)的引用

5.蒙特卡羅程序結(jié)構(gòu)

作業(yè)

第五章蒙特卡羅方法在計算機上的實現(xiàn)

蒙特卡羅方法是隨著計算機的出現(xiàn)和發(fā)展

而逐步發(fā)展起來的。在計算機上能夠產(chǎn)生符合

要求的隨機數(shù)，實現(xiàn)對已知分布的抽樣，奠定

了蒙特卡羅方法在計算機上得以實現(xiàn)的基礎(chǔ)。

在計算機上使用蒙特卡羅方法解粒子輸運問題

大致包括三個過程：源分布抽樣過程，空間、

能量和運動方向的隨機游動過程以及記錄、分

析結(jié)果過程。

1.源分布抽樣過程

源分布抽樣的目的是產(chǎn)生粒子的初始狀態(tài)

5。=（知￡。，4）。下面我們介紹一些常見的特定

類型的源分布抽樣方法。

1)源粒子的位置常見分布的隨機抽樣

(1)圓內(nèi)均勻分布

設(shè)圓半徑為小，粒子在圓內(nèi)均勻分布時，從發(fā)射

點到中心的距離r的分布密度函數(shù)為：

2r

,2當0工廠工R0

，(吟=\人0

o其它

〃的抽樣方法為:

=凡-值虱。，2)

(2)圓環(huán)內(nèi)均勻分布

設(shè)圓環(huán)的內(nèi)半徑為火。,外半徑為與,則粒子在該圓

環(huán)內(nèi)均勻分布時，從發(fā)電寸點到中心的距離〃的分布密

度函數(shù)為：

當&W居

0其它

V

尸的抽樣方法為：

二（%一&）?max624）+凡曠=（R「Ro）七2+4

(3)球內(nèi)均勻分布

設(shè)球的半徑為七粒子在球內(nèi)均勻分布時，從發(fā)射

點到中心的距離r的分布密度函數(shù)為：

3r2

/0)={R3當OWE

0其它

一的抽樣方法為：

/=氏?值乂4，蜃，43)

在直角坐標系下，抽樣方法為：

Xo=Rpi，，為=尺,〃2，Zo=R.r/3

I

(4)球殼內(nèi)均勻分布

設(shè)球殼的內(nèi)半徑為火。,外半徑為在均勻分布時,

從發(fā)射點到中心的距離r的分布密度函數(shù)為：

當凡〈〃酒

f(叫=用一總

其它

〃的抽樣方法為:

<----------

R；+Ro&+火：

<----3-凡-火-］---

火;+Ro&+火；

x=42x=max(^2,^3)x=max?2,3,4)

r=(7?1-7?0)-x+7?0

在直角坐標系下，球殼內(nèi)點的坐標為：

x0=r-sin^cos^

%=/?sinOsino

z0-r-cosO

其中，〃由前面的抽樣方法確定，。、夕服從各向同性

分布，其抽樣方法為：

k

軸

x=r-smOcos(p=r-242

Q^+A2^+A2^

(2必用3

y0=r-sm0smp=r-J；+

Z-.COSf宿

(5)圓柱內(nèi)均勻分布

圓柱內(nèi)均勻分布是指粒子發(fā)射點均勻地分布在底

半徑為尺高為2〃的圓柱內(nèi)。若固定圓柱的中心為

原點，圓柱的軸向為z軸，則分布密度函數(shù)為：

1當，+/4火2

?,\z\<H

f(x,y,z)=<2兀*HR

0其它

抽樣方法為:

<

XO=R，7,%=△,%，ZO=H.〃3

（6）點源分布

點源分布是指粒子由一固定點宙/；/；）發(fā)射,

其分布密度函數(shù)為：

/（x/,z）=S（x—耳）.況歹——Z；）

其中，/?）為狄拉克d函數(shù)，源粒子的抽樣方法為：

***

X=%,歹二歹0，z—Zo

在球坐標系中，粒子發(fā)射點到球心的距離〃的分布

密度函數(shù)為：

/（尸）=5（—3）

*

其中，%為點源到球心的距離。源粒子的位置抽樣為：

*

(7)球外平行束源分布

球外平行束源分布是指粒子平行入射到半徑為R

的球面上，或球外點源距離球很遠，可以近似地看作

平行束源。設(shè)〃為粒子發(fā)射點到球心的距離，其分布

密度函數(shù)為：

f(r)=3(r-R)

〃的抽樣方法為：

r=R

在直角坐標系中，抽樣方法為:

<1—

!<

為=火21，ZO=R?〃2

2)源粒子的能量常見分布的隨機抽樣

(1)單能源分布

單能源分布是指粒子的發(fā)射能量為一固定值乙,

其分布密度函數(shù)為：

f(E)=b(E-E°)

源粒子的能量為：

E=E。

(2)裂變中子譜分布

裂變中子譜分布的一般形式為：

E/A

f(E)=C-e--sh4BE,Emin<E<Emax

其中Z,B,C,Em[n，耳陋均為與元素有關(guān)的量。

對于鈾-235,

4=0.965,5=2.29,C=0.453,Emin=0,Em=（x）o

采用近似修正抽樣，抽樣方法為:

>

1：

M

<

■

E=E'=-上一暄

1-A23

其中，m-0.8746,M^O.2678,2=0.5543。

AmE

H、(E)=C-shV^-e叩{一力￡}

1—AXA2

止匕外，裂變譜分布也有以數(shù)值曲線形式給出的,

此時，用數(shù)值曲線抽樣方法抽取￡。

1

(3)麥克斯韋(Maxwell)譜分布

麥克斯韋譜分布的一般形式為：

2/73/2_

/位)=上=四?小,E>0

該分布的抽樣方法為

3V—e,.ln5上

<

E-_±

2B

1

3)源粒子運動方向常見分布的隨機抽樣

(1)各向同性分布

各向同性分布密度函數(shù)為：

一1

/(◎)=/(〃)?/?=「

力(。)=4

2In

其中，〃=cos仇。為運動方向與軸的夾角，夕為

方位角。

在直角坐標系下，各方向余弦"，V,墳為:

u=sinOcos。

v=sin8sine

w=cos6

其抽樣方法為：

（盤+/2堀+/2稔）244

u=sin3coscp=----------~~---

盤+/瑞+/2汽

v—singsin0=2/初3

長+///；+A2^

默—A?相一A2〃；

W=cosg=

盤+A2^+A2^

(2)半面各向同性分布

不妨設(shè)在e0的半面方向上各向同性發(fā)射粒子，

則在前述各向同性分布的抽樣方法中，用12代替叱就

能得到所需分布的抽樣。對于其它方向的情況，可用

類似的方法處理。

(3)球外平行束源分布

令〃=cos仇。為粒子運動方向的徑向夾角，貝IJ〃分

布密度函數(shù)為：

/(〃)=—2〃，-1<//<0

〃的抽樣方法為:

"一的看點)

(4)球外各向同性點源分布

設(shè)球外點源S到球心的距離為。0。點源S到球的

最大張角為

“1°)

則球外各向同性點源分布的抽樣方法是：」

先抽樣確定夕，再轉(zhuǎn)換成凡方古北巫『左大HV

*在直角坐標系下，取

cos-=1-(1-cos。)?4OS為Z軸，抽樣方法為：

u=sin。'

小R2—以s/j，

cos8=-v=0

Rw=-cos。'

4）次級粒子的源分布

在有關(guān)次級粒子（如裂變中子，中子生成光子，

光子生成中子）的輸運過程中，次級粒子源分布的抽

樣方法，主要可分為以下兩種：

（1）直接生成法

可將生成的次級粒子的位置、能量、方向、權(quán)重

等參數(shù)直接作為源分布的抽樣結(jié)果。也就是直接對生

成的次級粒子進行跟蹤。這種方法比較簡單、直觀。

(2)離散分布法

將生成的次級粒子的權(quán)重，按空間位置、能量、

方向分別記錄，得到次級粒子的空間、能量、運動方

向的離散的近似分布。再根據(jù)該分布，利用各種抽樣

焚巧，得到源分布的抽樣，對抽樣的源粒子進行跟蹤、

記錄。

當一個問題需要用兩個以上的蒙卡程序處理時，

可采用這種方法。

2.空間、能量和運動方向的隨機游動過程

粒子由狀態(tài)S加到狀態(tài)與+i時，需要

確定粒子的空間位置o+i，能量號+1和

層動萬向。冽+1。

1）碰撞點位置的計算公式

設(shè)o為粒子第加次碰撞點的位置，。機為碰撞后

的運動方向，則粒子第加+1次碰撞點的位置0+］為：

r

777+1=9n+L?Qm

即

X/x+l

ym+i=Zn+￡.v加

2加+1=2根+￡.%

其中3〃”昵，叱J為只”的方向余弦，L為兩次碰撞點間

的距離。

￡的分布密度函數(shù)為:

/⑷二邑(〃+1,￡,“)exp卜［Z(〃+14,Em￡>0

由/(￡)抽樣確定L的方法通常有三種：

(1)直接抽樣方法

確定上的直接抽樣方法是：

首先由自由程分布

/(夕)=「

中抽取P

夕=_1nq

再由下列關(guān)系式解出L。

夕=二“〃+/?4，&)"/

對于均勻介質(zhì)，有pInJ/---------------

4希1人l,田,/(&)2(&)

對于多層介質(zhì)，如果

￡的.?”夕＜￡然?％.(&)

z=0z=0

夕—￡第.％.(&)

則

i=0

z=0一」(&)

其中，羽，為粒子由r出發(fā)，沿a方向在順序經(jīng)過的

第，?個介質(zhì)區(qū)域內(nèi)走過的距離，EQ回)為第，個介質(zhì)

區(qū)域的宏觀總截面(1=1,2,…;1〃Q)。當

,max

邑人EG

z=0

時，意味著粒子穿出系統(tǒng)。

(2)最大截面法

對于多層介質(zhì)，或其他介質(zhì)密度與位置有關(guān)的問題，

在求A￡C=1,2,/〃皿)時，如果系統(tǒng)形狀復(fù)雜,

計算是非常煩雜的。在這種情況下，使用最大截面法

更方便。最大截面抽樣方法為：

4二0

A=L,——嶼一-----

2f,max(￡〃)

I

1<〃.。團,Em)>

2-

L=L]

其中\(zhòng)max(￡)=max￡M￡)

r

⑶限制抽樣法

當介質(zhì)區(qū)域很小時，如使用直接抽樣法抽取輸運

長度，粒子很容易穿出介質(zhì)，此時使用限制抽樣法確

定自由程個數(shù)P較好，P的分布密度函數(shù)為：

/(夕)=<匚"當ovpv。"

0其它

其中。加為粒子由％出發(fā)，沿與〃方向到達區(qū)域邊界的

自由程個數(shù)。P的抽樣方法是：

夕=_1n

然后用直接抽樣法中根據(jù)p計算￡的方法計算輸運長

度此時，粒子的權(quán)重需乘以糾偏因子(1-H%)。

2）碰撞后能量5+1的隨機抽樣

粒子在介質(zhì)中發(fā)生碰撞后，首先要確定與哪種原

子核發(fā)生何種反應(yīng)。粒子發(fā)生碰撞后（吸收除外）的

能量紇人一般只與其碰撞前后運動方向的夾角（散射

鬲）有關(guān)。

粒子碰撞后常見的能量分布有下面幾種情況。

（1）裂變中子譜

中子引起原子核裂變反應(yīng)時，裂變中子的能量服

從裂變譜分布。其抽樣方法可參考以前的介紹。

(2)中子彈性散射后能量的確定

中子彈性散射后，能量與質(zhì)心系散射角%的關(guān)系

是：

Em+i溶獷(1+2"+1)

能量與實驗室系散射角〃的關(guān)系是:

E

F—____>21_

m+1(4+1)2

其中，/為碰撞核的質(zhì)量，4=cos%,4=cos%。

4c或4確定后，即可求出號+1。

⑶中子非彈性散射后能量的確定

中子非彈性散射后,能量與質(zhì)心系散射角%的關(guān)

系是：

F____________■

Em+l-z/j_[、2K(l-以/&)+2431-以/&?%+U

(A+

4+1

其中，?為第K個能級的閾能，底為第K個能級的激

友合肥里。

如果確定了實驗室系散射角%,則根據(jù)下式

1

4c二AyjisIE"J/2(i-//￡〃J-]+筋+Mt

確定論后，再計算出￡

(4)光子康普頓(Compton)散射后能量的確定

光子發(fā)生康普頓散射后，其能量分布密度函數(shù)為：

///、1(a+1—%)111/I。

=----------+------r+F，1<X<1+26<

K(a)\a-xJxxx

其中，K(a)為歸一因子。

u(、12(i+l)[、141

K(a)=12—^(1+2^)+-+--------1

a2a2(1+2。)

x=a〔a',。和優(yōu)分別為光子散射前后的能量，以

加。理為單位，加°為電?青箏止質(zhì)量，C為先逮。

光子康普頓散射能量分布的抽樣方法為:

<

.A>

4。+29

1+2。

l+2a4x2=1+2。4

227("IP>

>7,a+1—X］、

+1

21aj4月

<

<

x=x

X=X]2

X的抽樣確定后，散射后的能量為:

，2二2Em

LF〃Z+1-a.加°。=—?mc

xQx

3）碰撞后散射角的隨機抽樣

粒子碰撞后運動方向億什］的確定，一般與散射角

有關(guān)。由已知分布抽樣確定散射角后，再確定。能+1。

常見的散射角分布有如下幾種：

（1）質(zhì)心系各向同性分布

散射角在質(zhì)心系服從各向同性分布時，其抽樣方

法為4c=〃=24-1。質(zhì)心系散射角生抽樣確定后，

需轉(zhuǎn)換成實驗室系散射角生。

在中子彈性散射情況下，轉(zhuǎn)換公式為:

1+4〃

其中A為碰撞核質(zhì)量，4。=cos%,4=cos%o

在中子非彈性散射情況下，轉(zhuǎn)換公式為：

4________1+?J1—//Em?Nc

L+//紇）+24/l-

其中，。為第K個能級的閾能。

(2)中子彈性散射勒讓德(Legendre)多項式分布

中子彈性散射角分布常以勒讓德多項式的展開形

式給定。散射角余弦x=cos。的分布密度函數(shù)為：

八上恪號力/(X)當1昨1

0其它

其中P/(x)為/階勒讓德多項式。

該分布即為〃階勒讓德近似展開。

勒讓德多項式由以下遞推公式確定：

(〃+1)只+1(x)—(2〃+1)吠(%)+叱T(X)=0

鳥(%)=1

6(X)=X

考慮新的分布:

n

fa(x)=￡R3(x-Xk)

k=0

當選取Xo，X],…x〃為P〃+i(x)=0的根，且

9n晉△711力/E)

縱二號-------------

￡甘有(/)]2

1=0乙

時，fa（x）依照勒讓德多項式展開的前n項與/（X）的展

開形式相同。因此，可以用力（X）作為/（X）的近似分布。

在實際問題中，由于勒讓德多項式展開項數(shù)不夠,

可能出現(xiàn)某個4為負值的現(xiàn)象。此時可以采用如下近

似分布：〃

，*（、）=￡咒-

其中:k=0

SI^Ikk=0

k=0

對于該近似分布，可用加抽樣方法進行抽樣：

K—1K

、=打，當

k=0k=0

此時，由于偏倚抽樣而引起的糾偏因子為WK，也就

是儻，粒子的權(quán)重要乘上限。

告

(3)光子康普頓散射角分布

光子的康普頓散射角與其散射前后的能量有關(guān)，它

的分布密度函數(shù)為：

/(X)=5(1

抽樣方法為:

4)碰撞后運動方向。團+1的確定

實驗室系散射角％確定后，依據(jù)不同的坐標系的表

現(xiàn)形式，有不同的確定方法。

(1)確定方向余弦%葉1，匕/1，

一bCHmUm-bdVm

Um+\+aUm

2

Ju2+v

Vmm

-bcwv+bdu

mnini

1m色+av

Vm+\m

%2+v：2

22

叱72+1%〃+以+叫

其中,

a=cos%,b=sin%-yll-a2,

c=cos/,d=sin/

方位角，在［0,2TT］上均勻分布。

當說+噂一。時，不能使用上述公式，可用下

面的簡單公式：

/用二be

v〃用二bd

叱用二叫

（2）確定@什］的球坐標（“i，巴"1）

設(shè)的球坐標分別為（/用其中，。為粒子運動

方向與z軸的夾角，夕為粒子運動方向在xy平面上投

影的方位角。則0什1的球坐標（g+1，9〃什1）分別由下式

確定：

C0S/+1=cos%cos%+sin%sin%cos/

sinasin/

sin(%_0〃J=

sin。

。/

COSL—cos，m”cosa1加1+1

cos@〃,+i—?！?？)=

sina,msin，m”+上l］

5）球形幾何的隨機游動公式

一般幾何的隨機游動公式可以應(yīng)用到球形幾何,

而對球?qū)ΨQ問題，使用特殊形式更為方便。

（1）下次碰撞點的徑向位置〃+i的確定

兩次碰撞點間的距離上確定之后，下次碰撞點的

徑向位置G+i的計算公式為：

G+i=G：+G+2L〃m.cos9m

設(shè)系統(tǒng)的外半徑為七如〃+1次，則粒子逃出系統(tǒng)。

⑵粒子碰撞后瞬時運動方向的確定

在球?qū)ΨQ系統(tǒng)中，粒子運動方向用其與徑向夾角

余弦來描述。使用球面三角公式，粒子碰撞后瞬時運

動方向與徑向夾痛余弦COS%+1的計算公式為：

cos?！ㄈ~1=cos。：cos%+sine；sin%cos/

sin

。m：vm

m+1

八,￡+r,“?cos，“

cosO----------------

mr

m+1

其中，力為在［0,2用上均勻分布的方位角，現(xiàn)為在

%+1點進入碰撞前瞬時運動方向與G+i徑向之間的夾

施。

6）點到給定邊界面的距離

在抽樣確定輸運距離、判斷粒子是否穿透系統(tǒng)時，

常遇到求由G出發(fā)，沿方向到達某個區(qū)域表面的

距離問題。在記錄對結(jié)果的貢獻時，也常使用類似的

量。區(qū)域表面通常是平面或二次曲面。求到達區(qū)域表

面的距離問題，實際上是求直線（或半射線）與平面

或二次曲面的交點問題。這是蒙特卡羅方法解粒子輸

運的各種實際問題時，所遇到的基本幾何問題。

(1)點到平面的距離

點4=(%,%,Z。)沿方向4=(%,V。,%)的直線方程為:

r=「/?q

該直線到達方程為

ax+by+cz—d

的平面的距離為：

d_(QXo+byo+cZo)

/=----------------

au0+bv0+cw0

當與平面平行時，即

auQ+bv°+cw0=0

直線與平面無交點。

如果/為負值，直線與平面也無交點。這時，粒子的

運動方向是背離平面的。

?

(2)點到球面的距離

在三維直角坐標系中，設(shè)球心為〃=(%,”/0)，

球半徑為七則球面方程為：

Ofy+Q—yj+(z-zj=R2

將直線方程代入該球面方程，得到點分沿4方向到達

球面的距離I：

I=-5±VA

其中

S=%-〃)?4=(%-^>0+(%-紇)V。+(Z°-Zc)%

A=B2—R；+R2

R；=匕一〃「=(%—%)2+(歹0-"P+Go—Z)

當R04R時，即r0點在球內(nèi)，人三我/只有一個正

根：

I=—5+VA

當火o〉R時，即乙點在球外，分以下三種情況：

a）若應(yīng)0,/無正實根，直線與球面無交點。

b）若3<0,A<0,/無實根，直線與球面無交點。

c）若3<0,A>0,/有兩個正實根，直線與球面有兩個交

點。

I=—5±VA

在球坐標系中，不失一般性，設(shè)球心為〃=0,則

球面方程為r=Ro

當〃0次時，即幾點在球內(nèi)，有一個交點：

2

I=-rQ-cos^0+JR2—(%-sin^0)

其中綜為Oo與勺的徑向夾角。

當r0>7?時，即乙點在球外，令

(7=-sin^0

當cos%K)時，直線與球面無交點。

當cos/<0時，若也R,則直線與球面無交點。

若d〈R,則有兩個交點：

/二一q?cos'±J/?2一伍.sin4)2

(3)點到圓柱面的距離

設(shè)圓柱面的方程為：

)2+(f丁"

其中為圓柱的中心，火為圓柱底半徑。

點勺沿4方向到達圓柱面的距離I為:

,一3±VA

/=----%-

1一%

其中

3二(%-凡)%+(%-

R；=(Xo~Xc)2+。0-以)2

當火0逆時，打點在圓柱內(nèi)，如果用則/有

一個正根：

7-3+VA

/=----5-

1一%

如果用=1,即4平行于圓柱的對稱軸，直線與

圓柱面無交點。

當凡〉火時，尸0點在圓柱外，分以下三種情況：

a）若后0,/無正實根，直線與圓柱面無交點。

b）若3<0,A<0,/無實根，直線與圓柱面無交點。

c）若d<0,AN0且晡聲，有兩個正實根，直線與圓柱

面有兩個交點。

7-5±VA

/=----%-

1一%

在晡=1的情況下，直線與圓柱面不相交。

(4)點到圓錐面的距離

設(shè)圓錐頂點在原點，以z軸為對稱軸，則圓錐面

的方程為：

x2+y2=c2z2

點打沿4方向到達圓錐面的距離1為:

一5±VA

1-(1+。2)若

其中

d=x^+y^-c2z^

△=3?一(x；+M-c2z1)[!-(1+c2)Wg]

如果4與錐面某一母線平行，即(1+/)晡=i,則

Xo+^o-^o

1=

(5)空腔處理

在粒子輸運問題中，所考慮的系統(tǒng)常有空腔存在，

如中空的球殼，平板間的空隙等。粒子輸運時，可有

兩種處理空腔的方法：

a)將空腔作為宏觀總截面％=0的區(qū)域，按通常的方法輸

運。

b)設(shè)4、寸分別為由％出發(fā)，沿。加方向到空腔區(qū)域的

近端和遠端的交點，當粒子超過匕時，以廠為新的起

點，重新開始輸運。

顯然，這兩種方法在統(tǒng)計上是等價的。

7）等效的邊界條件

（1）全反射邊界

在反應(yīng)堆活性區(qū)中，元件盒常常按正方形或六角

形排列。假定元件盒足夠多，每個盒結(jié)構(gòu)相同，那么

活性區(qū)中每個盒所占的柵胞的物理情況，可以代表整

進一步假定，元件盒是圓對稱的，那么每個柵胞

中情況，可以用更小的單位（柵元）來反映。比如對

六角形柵胞可取其1/12的AOZB來做代表；正方形柵

胞可用其1/8的AOZ2來做代表。這樣一來問題就大

大簡化了。

現(xiàn)在的問題是怎樣計算直角三角形柵元的物理量

（如通量）。用蒙特卡羅方法如何模擬中子在柵元內(nèi)

的運動，反映出整個活性區(qū)對它的影響。

我們可把。T、0B'、/9作為全反射邊界來處理。

所謂全反射邊界，就是當中子打到該邊界上時，按鏡

面反射的方式，從邊界上全部反射回來，中子的能量

與權(quán)重均不改變。

在這種邊界上的反射條件，稱之為全反射條件，

就是通常的鏡面反射條件。

在全反射邊界條件下，一條通過活性區(qū)若干個區(qū)

域的中子徑跡，可以用柵元A。為聲中的一條折線軌

道來反映出來。

反過來，在直角三角形柵元A。力2中任一條反

射成的折線軌道，都代表了中子在活性區(qū)內(nèi)一條直線

軌道的作用。由于系統(tǒng)的對稱性，在活性區(qū)內(nèi)，凡是

與柵元內(nèi)位置相當?shù)牡胤?，都有相同的物理情況，因

此柵元內(nèi)各處的情況，當然代表了整個活性區(qū)的情況。

(2)一般曲面全反射條件

對于一般曲面的全反射，設(shè)入射方向為入射

點的內(nèi)法線方向為〃，則反射方向2為：

。-￡2-2（p-n）n

=q+2cos9?九

cosO=T2"

u="+2cos6?〃X

v=v+2cos9“

w=w+2cos。?〃二

cos。=一（"?〃、+v"+w?電）

(3)平面全反射條件

設(shè)三角形柵元的橫截面A0Z5在XT平面上，

ZOAB=0o則邊界。4、0B、48上的反射都是平面

全反射。在任一與x-y平面垂直且與x軸成Q角的平

面上，全反射條件為：

u二》?cos2a+v?sin2a

M=〃?sin2。-v?cos2a

wf=w

O

由此就可得到O/、OB和48邊上的全反射條件，對

于05邊，a=0；對于。4邊，a=0；對于Z5邊，

。=兀/2。

(4)反射層邊界條件

對于具有大反射層的系統(tǒng)，如存放，運輸和生產(chǎn)

裂變物質(zhì)的倉庫、車廂和車間等，當中子從里面打到

四周墻上或反射層時，還要繼續(xù)對它進行跟蹤。這種

跟蹤常常要花費很大的計算量，并且在結(jié)果中引起的

方差也比較大。如果在計算這種系統(tǒng)的不同方案中，

反射層條件不變，那么這種大量重復(fù)的計算是很不經(jīng)

濟的。

中子射入反射層后，一部分被介質(zhì)吸收，

只有一部分返回，由于中子的散射慢化，損失

一部分能量，因此反射回來的中子有一個能量

方向分布。顯然，對這種反射層，不能應(yīng)用全

反射條件。不過，我們?nèi)匀豢梢园阉斪鲞吔?

在邊界上按反射層的物理作用來處理。

比如，如果反射層是一種平板幾何，我們可以用

數(shù)值方法或蒙特卡羅方法，預(yù)先算好在各種不同入射

能量E下的反照窣打(￡),反J寸中子的能量分希

RE(ETE')。于是代替在反射層中眼蹤中子，我們可

在反射層邊界上作如下處理：

一旦中子打入反射層，立即返回，反射后權(quán)重為

其中，￡為射入反射層中子的能量，沙為中子的權(quán)重。

反射后的能量與由反射能譜火￡(￡一為)中抽樣產(chǎn)生。

反射后的方向應(yīng)由半平面各向同性分布或余弦分布中

抽樣。反射后的中子位置為入射時的位置。

計算表明，對于大尺寸的反射層來說，這樣的近

似，引起的結(jié)果上的誤差是可以忽略的，卻能帶來計

算量的大量節(jié)省。

3.記錄貢獻與分析結(jié)果過程

在粒子輸運問題中，除了得到某些量的積

分結(jié)果外，還需要得到這些量的方差、協(xié)方差、

以及這些量的空間、能量、方向和時間的分布。

這些量可以利用分類記錄手續(xù)同時得到。

1）記錄與結(jié)果

為了得到所求量的估計值，在粒子輸運過程中需

進行記錄，即求每個粒子對所求量的貢獻。

設(shè)模擬了N個粒子，所求量的估計值為:

_1N

孰=后卒

其中&為第，個粒子的總貢獻。

記錄的貢獻由所求量決定。對于同一個所求量，

又隨所用的蒙特卡羅技巧的不同而不同。例如，所求

量是粒子穿透屏蔽概率，使用直接模擬法時，如粒子

穿透屏蔽，在疊加記錄單元加“1”（初始值為零），否

則沒有貢獻。使用加權(quán)法時，如粒子穿透屏蔽，在疊

加記錄單元加粒子的權(quán)重，否則沒有貢獻。使用統(tǒng)計

估計法時，粒子每發(fā)生一次碰撞（包括零次碰撞），都

要記錄貢獻，等等。

2）方差和協(xié)方差的估計

估計量g和@的方差和協(xié)方差為:

城=Mg）（￡g）2

咤=￡（gg'）—（￡g）（￡g'）

它們可以用下式估計:

N、2

1N1

Z=11NZ=1

1"(\NV]N

gig"2gi

N心INi=i人N1)

因此，要得至Ijb：和b；,,的估計，只要對每一個歷史記

錄結(jié)果的,和g,￡進行記錄，并加以累加即可。

方差估計值U確定后，可得到誤差

其中凡為置信限，它隨置信水平1-。而定。在通常

情況下，取1—。=0.95,%=1.96。

3）位置、能量、方向、時間分布

在前面已經(jīng)提到，用蒙特卡羅方法求某種

量的空間、能量、方向和時間分布，實質(zhì)上是

得到這種分布的階梯函數(shù)近似的估計值。而求

這種估計值是很方便的，只要將跟蹤過程中所

得到的感興趣量，按其狀態(tài)的空間、能量、方

向、時間特征，分別記錄其權(quán)重，最后將這

些記錄結(jié)果適當處理即可。

事先，將問題的空間、能量、方向（常按相對于

某個方向的夾角余弦）、時間范圍，各分為如下不同

間隔：

匕匕…，匕;

￡min=Ej<E、<E°=Max；

—l=4o<〃i<…=1；

0=Zo<tx<-?-<tL=T;

再用一批存貯單元｛A｝記錄相應(yīng)間隔上階梯函數(shù)

近似的累計值。

4.核截面數(shù)據(jù)的引用

用蒙特卡羅方法解粒子輸運問題，需要介質(zhì)所包

含的各種原子核的核數(shù)據(jù)。以中子核數(shù)據(jù)為例，需要

各種涉及到的核的微觀總截面、彈性散射、非彈性散

射、n-2n反應(yīng)、裂變、俘獲等截面；也需要這些反應(yīng)

的相應(yīng)能量、角度分布、次級粒子數(shù)，以及其它關(guān)心

的粒子數(shù)及其能量、方向分布。從輸運方程中可以知

道，有了這些數(shù)據(jù)，問題就完全確定了；反映到蒙特

卡羅模擬中，有了這些數(shù)據(jù)，就能決定宏觀總截面，

決定碰撞核的具體形式，就能實現(xiàn)抽樣和跟蹤。

在蒙特卡羅計算中，引用的核數(shù)據(jù)有點截面、分

段常數(shù)截面和群截面三種形式。

1）點截面形式

在跟蹤粒子時，對粒子的每一種能量，先

從截面庫中取出需要的核數(shù)據(jù)，再用插值（或

其它方式），求出相應(yīng)能