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

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

1.源分布抽樣過(guò)程

2.空間、能量和運(yùn)動(dòng)方向的隨機(jī)游動(dòng)過(guò)程

3.記錄貢獻(xiàn)和分析結(jié)果過(guò)程

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

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

作業(yè)

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

蒙特卡羅方法是隨著計(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ò)程。

1.源分布抽樣過(guò)程

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

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

類(lèi)型的源分布抽樣方法。

1)源粒子的位置常見(jiàn)分布的隨機(jī)抽樣

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

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

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

2r

,2當(dāng)0工廠(chǎng)工R0

，(吟=\人0

o其它

〃的抽樣方法為:

=凡-值虱。，2)

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

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

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

度函數(shù)為：

當(dāng)&W居

0其它

V

尸的抽樣方法為：

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

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

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

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

3r2

/0)={R3當(dāng)OWE

0其它

一的抽樣方法為：

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

在直角坐標(biāo)系下，抽樣方法為：

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

I

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

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

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

當(dāng)凡〈〃酒

f(叫=用一總

其它

〃的抽樣方法為:

<----------

R；+Ro&+火：

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

火;+Ro&+火；

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

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

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

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ā)射點(diǎn)均勻地分布在底

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

原點(diǎn)，圓柱的軸向?yàn)閦軸，則分布密度函數(shù)為：

1當(dāng)，+/4火2

?,\z\<H

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

0其它

抽樣方法為:

<

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

（6）點(diǎn)源分布

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

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

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

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

***

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

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

密度函數(shù)為：

/（尸）=5（—3）

*

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

*

(7)球外平行束源分布

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

的球面上，或球外點(diǎn)源距離球很遠(yuǎn)，可以近似地看作

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

密度函數(shù)為：

f(r)=3(r-R)

〃的抽樣方法為：

r=R

在直角坐標(biāo)系中，抽樣方法為:

<1—

!<

為=火21，ZO=R?〃2

2)源粒子的能量常見(jiàn)分布的隨機(jī)抽樣

(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)的量。

對(duì)于鈾-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ù)值曲線(xiàn)形式給出的,

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

1

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

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

2/73/2_

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

該分布的抽樣方法為

3V—e,.ln5上

<

E-_±

2B

1

3)源粒子運(yùn)動(dòng)方向常見(jiàn)分布的隨機(jī)抽樣

(1)各向同性分布

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

一1

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

力(。)=4

2In

其中，〃=cos仇。為運(yùn)動(dòng)方向與軸的夾角，夕為

方位角。

在直角坐標(biāo)系下，各方向余弦"，V,墳為:

u=sinOcos。

v=sin8sine

w=cos6

其抽樣方法為：

（盤(pán)+/2堀+/2稔）244

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

盤(pán)+/瑞+/2汽

v—singsin0=2/初3

長(zhǎng)+///；+A2^

默—A?相一A2〃；

W=cosg=

盤(pán)+A2^+A2^

(2)半面各向同性分布

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

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

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

類(lèi)似的方法處理。

(3)球外平行束源分布

令〃=cos仇。為粒子運(yùn)動(dòng)方向的徑向夾角，貝IJ〃分

布密度函數(shù)為：

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

〃的抽樣方法為:

"一的看點(diǎn))

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

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

最大張角為

“1°)

則球外各向同性點(diǎn)源分布的抽樣方法是：」

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

*在直角坐標(biāo)系下，取

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

u=sin。'

小R2—以s/j，

cos8=-v=0

Rw=-cos。'

4）次級(jí)粒子的源分布

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

光子生成中子）的輸運(yùn)過(guò)程中，次級(jí)粒子源分布的抽

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

（1）直接生成法

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

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

成的次級(jí)粒子進(jìn)行跟蹤。這種方法比較簡(jiǎn)單、直觀。

(2)離散分布法

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

方向分別記錄，得到次級(jí)粒子的空間、能量、運(yùn)動(dòng)方

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

焚巧，得到源分布的抽樣，對(duì)抽樣的源粒子進(jìn)行跟蹤、

記錄。

當(dāng)一個(gè)問(wèn)題需要用兩個(gè)以上的蒙卡程序處理時(shí)，

可采用這種方法。

2.空間、能量和運(yùn)動(dòng)方向的隨機(jī)游動(dòng)過(guò)程

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

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

層動(dòng)萬(wàn)向。冽+1。

1）碰撞點(diǎn)位置的計(jì)算公式

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

的運(yùn)動(dòng)方向，則粒子第加+1次碰撞點(diǎn)的位置0+］為：

r

777+1=9n+L?Qm

即

X/x+l

ym+i=Zn+￡.v加

2加+1=2根+￡.%

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

的距離。

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

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

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

(1)直接抽樣方法

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

首先由自由程分布

/(夕)=「

中抽取P

夕=_1nq

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

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

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

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

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

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

z=0z=0

夕—￡第.％.(&)

則

i=0

z=0一」(&)

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

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

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

,max

邑人EG

z=0

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

(2)最大截面法

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

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

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

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

4二0

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

2f,max(￡〃)

I

1<〃.。團(tuán),Em)>

2-

L=L]

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

r

⑶限制抽樣法

當(dāng)介質(zhì)區(qū)域很小時(shí)，如使用直接抽樣法抽取輸運(yùn)

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

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

/(夕)=<匚"當(dāng)ovpv。"

0其它

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

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

夕=_1n

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

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

2）碰撞后能量5+1的隨機(jī)抽樣

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

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

能量紇人一般只與其碰撞前后運(yùn)動(dòng)方向的夾角（散射

鬲）有關(guān)。

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

（1）裂變中子譜

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

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

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

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

是：

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

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

E

F—____>21_

m+1(4+1)2

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

4c或4確定后，即可求出號(hào)+1。

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

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

系是：

F____________■

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

(A+

4+1

其中，?為第K個(gè)能級(jí)的閾能，底為第K個(gè)能級(jí)的激

友合肥里。

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

1

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

確定論后，再計(jì)算出￡

(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）碰撞后散射角的隨機(jī)抽樣

粒子碰撞后運(yùn)動(dòng)方向億什］的確定，一般與散射角

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

常見(jiàn)的散射角分布有如下幾種：

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

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

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

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

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

1+4〃

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

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

4________1+?J1—//Em?Nc

L+//紇）+24/l-

其中，。為第K個(gè)能級(jí)的閾能。

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

中子彈性散射角分布常以勒讓德多項(xiàng)式的展開(kāi)形

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

八上恪號(hào)力/(X)當(dāng)1昨1

0其它

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

該分布即為〃階勒讓德近似展開(kāi)。

勒讓德多項(xiàng)式由以下遞推公式確定：

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

鳥(niǎo)(%)=1

6(X)=X

考慮新的分布:

n

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

k=0

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

9n晉△711力/E)

縱二號(hào)-------------

￡甘有(/)]2

1=0乙

時(shí)，fa（x）依照勒讓德多項(xiàng)式展開(kāi)的前n項(xiàng)與/（X）的展

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

在實(shí)際問(wèn)題中，由于勒讓德多項(xiàng)式展開(kāi)項(xiàng)數(shù)不夠,

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

似分布：〃

，*（、）=￡咒-

其中:k=0

SI^Ikk=0

k=0

對(duì)于該近似分布，可用加抽樣方法進(jìn)行抽樣：

K—1K

、=打，當(dāng)

k=0k=0

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

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

告

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

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

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

/(X)=5(1

抽樣方法為:

4)碰撞后運(yùn)動(dòng)方向。團(tuán)+1的確定

實(shí)驗(yàn)室系散射角％確定后，依據(jù)不同的坐標(biāo)系的表

現(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］上均勻分布。

當(dāng)說(shuō)+噂一。時(shí)，不能使用上述公式，可用下

面的簡(jiǎn)單公式：

/用二be

v〃用二bd

叱用二叫

（2）確定@什］的球坐標(biāo)（“i，巴"1）

設(shè)的球坐標(biāo)分別為（/用其中，。為粒子運(yùn)動(dòng)

方向與z軸的夾角，夕為粒子運(yùn)動(dòng)方向在xy平面上投

影的方位角。則0什1的球坐標(biāo)（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）球形幾何的隨機(jī)游動(dòng)公式

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

而對(duì)球?qū)ΨQ(chēng)問(wèn)題，使用特殊形式更為方便。

（1）下次碰撞點(diǎn)的徑向位置〃+i的確定

兩次碰撞點(diǎn)間的距離上確定之后，下次碰撞點(diǎn)的

徑向位置G+i的計(jì)算公式為：

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

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

⑵粒子碰撞后瞬時(shí)運(yùn)動(dòng)方向的確定

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

余弦來(lái)描述。使用球面三角公式，粒子碰撞后瞬時(shí)運(yùn)

動(dòng)方向與徑向夾痛余弦COS%+1的計(jì)算公式為：

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

sin

。m：vm

m+1

八,￡+r,“?cos，“

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

mr

m+1

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

%+1點(diǎn)進(jìn)入碰撞前瞬時(shí)運(yùn)動(dòng)方向與G+i徑向之間的夾

施。

6）點(diǎn)到給定邊界面的距離

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

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

距離問(wèn)題。在記錄對(duì)結(jié)果的貢獻(xiàn)時(shí)，也常使用類(lèi)似的

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

面的距離問(wèn)題，實(shí)際上是求直線(xiàn)（或半射線(xiàn)）與平面

或二次曲面的交點(diǎn)問(wèn)題。這是蒙特卡羅方法解粒子輸

運(yùn)的各種實(shí)際問(wèn)題時(shí)，所遇到的基本幾何問(wèn)題。

(1)點(diǎn)到平面的距離

點(diǎn)4=(%,%,Z。)沿方向4=(%,V。,%)的直線(xiàn)方程為:

r=「/?q

該直線(xiàn)到達(dá)方程為

ax+by+cz—d

的平面的距離為：

d_(QXo+byo+cZo)

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

au0+bv0+cw0

當(dāng)與平面平行時(shí)，即

auQ+bv°+cw0=0

直線(xiàn)與平面無(wú)交點(diǎn)。

如果/為負(fù)值，直線(xiàn)與平面也無(wú)交點(diǎn)。這時(shí)，粒子的

運(yùn)動(dòng)方向是背離平面的。

?

(2)點(diǎn)到球面的距離

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

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

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

將直線(xiàn)方程代入該球面方程，得到點(diǎn)分沿4方向到達(dá)

球面的距離I：

I=-5±VA

其中

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

A=B2—R；+R2

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

當(dāng)R04R時(shí)，即r0點(diǎn)在球內(nèi)，人三我/只有一個(gè)正

根：

I=—5+VA

當(dāng)火o〉R時(shí)，即乙點(diǎn)在球外，分以下三種情況：

a）若應(yīng)0,/無(wú)正實(shí)根，直線(xiàn)與球面無(wú)交點(diǎn)。

b）若3<0,A<0,/無(wú)實(shí)根，直線(xiàn)與球面無(wú)交點(diǎn)。

c）若3<0,A>0,/有兩個(gè)正實(shí)根，直線(xiàn)與球面有兩個(gè)交

點(diǎn)。

I=—5±VA

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

球面方程為r=Ro

當(dāng)〃0次時(shí)，即幾點(diǎn)在球內(nèi)，有一個(gè)交點(diǎn)：

2

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

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

當(dāng)r0>7?時(shí)，即乙點(diǎn)在球外，令

(7=-sin^0

當(dāng)cos%K)時(shí)，直線(xiàn)與球面無(wú)交點(diǎn)。

當(dāng)cos/<0時(shí)，若也R,則直線(xiàn)與球面無(wú)交點(diǎn)。

若d〈R,則有兩個(gè)交點(diǎn)：

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

(3)點(diǎn)到圓柱面的距離

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

)2+(f丁"

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

點(diǎn)勺沿4方向到達(dá)圓柱面的距離I為:

,一3±VA

/=----%-

1一%

其中

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

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

當(dāng)火0逆時(shí)，打點(diǎn)在圓柱內(nèi)，如果用則/有

一個(gè)正根：

7-3+VA

/=----5-

1一%

如果用=1,即4平行于圓柱的對(duì)稱(chēng)軸，直線(xiàn)與

圓柱面無(wú)交點(diǎn)。

當(dāng)凡〉火時(shí)，尸0點(diǎn)在圓柱外，分以下三種情況：

a）若后0,/無(wú)正實(shí)根，直線(xiàn)與圓柱面無(wú)交點(diǎn)。

b）若3<0,A<0,/無(wú)實(shí)根，直線(xiàn)與圓柱面無(wú)交點(diǎn)。

c）若d<0,AN0且晡聲，有兩個(gè)正實(shí)根，直線(xiàn)與圓柱

面有兩個(gè)交點(diǎn)。

7-5±VA

/=----%-

1一%

在晡=1的情況下，直線(xiàn)與圓柱面不相交。

(4)點(diǎn)到圓錐面的距離

設(shè)圓錐頂點(diǎn)在原點(diǎn)，以z軸為對(duì)稱(chēng)軸，則圓錐面

的方程為：

x2+y2=c2z2

點(diǎn)打沿4方向到達(dá)圓錐面的距離1為:

一5±VA

1-(1+。2)若

其中

d=x^+y^-c2z^

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

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

Xo+^o-^o

1=

(5)空腔處理

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

如中空的球殼，平板間的空隙等。粒子輸運(yùn)時(shí)，可有

兩種處理空腔的方法：

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

運(yùn)。

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

近端和遠(yuǎn)端的交點(diǎn)，當(dāng)粒子超過(guò)匕時(shí)，以廠(chǎng)為新的起

點(diǎn)，重新開(kāi)始輸運(yùn)。

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

7）等效的邊界條件

（1）全反射邊界

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

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

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

進(jìn)一步假定，元件盒是圓對(duì)稱(chēng)的，那么每個(gè)柵胞

中情況，可以用更小的單位（柵元）來(lái)反映。比如對(duì)

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

胞可用其1/8的AOZ2來(lái)做代表。這樣一來(lái)問(wèn)題就大

大簡(jiǎn)化了。

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

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

的運(yùn)動(dòng)，反映出整個(gè)活性區(qū)對(duì)它的影響。

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

所謂全反射邊界，就是當(dāng)中子打到該邊界上時(shí)，按鏡

面反射的方式，從邊界上全部反射回來(lái)，中子的能量

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

在這種邊界上的反射條件，稱(chēng)之為全反射條件，

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

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

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

道來(lái)反映出來(lái)。

反過(guò)來(lái)，在直角三角形柵元A。力2中任一條反

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

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

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

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

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

對(duì)于一般曲面的全反射，設(shè)入射方向?yàn)槿肷?/p>

點(diǎn)的內(nèi)法線(xiàn)方向?yàn)椤ǎ瑒t反射方向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邊上的全反射條件，對(duì)

于05邊，a=0；對(duì)于。4邊，a=0；對(duì)于Z5邊，

。=兀/2。

(4)反射層邊界條件

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

裂變物質(zhì)的倉(cāng)庫(kù)、車(chē)廂和車(chē)間等，當(dāng)中子從里面打到

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

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

方差也比較大。如果在計(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下的反照窣打(￡),反J寸中子的能量分希

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

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

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

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

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

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

抽樣。反射后的中子位置為入射時(shí)的位置。

計(jì)算表明，對(duì)于大尺寸的反射層來(lái)說(shuō)，這樣的近

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

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

3.記錄貢獻(xiàn)與分析結(jié)果過(guò)程

在粒子輸運(yùn)問(wèn)題中，除了得到某些量的積

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

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

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

1）記錄與結(jié)果

為了得到所求量的估計(jì)值，在粒子輸運(yùn)過(guò)程中需

進(jìn)行記錄，即求每個(gè)粒子對(duì)所求量的貢獻(xiàn)。

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

_1N

孰=后卒

其中&為第，個(gè)粒子的總貢獻(xiàn)。

記錄的貢獻(xiàn)由所求量決定。對(duì)于同一個(gè)所求量，

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

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

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

則沒(méi)有貢獻(xiàn)。使用加權(quán)法時(shí)，如粒子穿透屏蔽，在疊

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

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

要記錄貢獻(xiàn)，等等。

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

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

城=Mg）（￡g）2

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

它們可以用下式估計(jì):

N、2

1N1

Z=11NZ=1

1"(\NV]N

gig"2gi

N心INi=i人N1)

因此，要得至Ijb：和b；,,的估計(jì)，只要對(duì)每一個(gè)歷史記

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

方差估計(jì)值U確定后，可得到誤差

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

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

3）位置、能量、方向、時(shí)間分布

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

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

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

這種估計(jì)值是很方便的，只要將跟蹤過(guò)程中所

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

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

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

事先，將問(wèn)題的空間、能量、方向（常按相對(duì)于

某個(gè)方向的夾角余弦）、時(shí)間范圍，各分為如下不同

間隔：

匕匕…，匕;

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

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

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

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

近似的累計(jì)值。

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

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

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

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

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

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

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

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

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

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

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

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

1）點(diǎn)截面形式

在跟蹤粒子時(shí)，對(duì)粒子的每一種能量，先

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

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