第五章 蒙特卡罗方法在计算机上的实现
1,源分布抽样过程
2,空间、能量和运动方向的随机游动过程
3,记录贡献和分析结果过程
4,核截面数据的引用
5,蒙特卡罗程序结构
? 作 业
第五章 蒙特卡罗方法在计算机上的实现
蒙特卡罗方法是随着计算机的出现和发展
而逐步发展起来的。在计算机上能够产生符合
要求的随机数，实现对已知分布的抽样，奠定
了蒙特卡罗方法在计算机上得以实现的基础。
在计算机上使用蒙特卡罗方法解粒子输运问题
大致包括三个过程：源分布抽样过程，空间、
能量和运动方向的随机游动过程以及记录、分
析结果过程 。
1,源分布抽样过程
源分布抽样的目的是产生粒子的初始状态
。 下面我们介绍一些常见的特定
类型的源分布抽样方法 。
),,( 0000 ΩrS E?
1) 源粒子的位置常见分布的随机抽样
(1) 圆内均匀分布
设圆半径为 R0，粒子在圆内均匀分布时，从发射
点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
??
?
?
? ??
?
其它
当
0
02
)( 020
Rr
R
r
rf
),m ax ( 210 ???? Rr
(2) 圆环内均匀分布
设圆环的内半径为 R0，外半径为 R1，则粒子在该圆
环内均匀分布时，从发射点到中心的距离 r 的分布密
度函数为,
r 的抽样方法为,
??
?
?
? ??
??
其它
当
0
2
)( 102021
RrR
RR
r
rf
020103201
01
01
1
)(),m a x ()( RRRrRRRr
RR
RR
????????
?
?
?
???
?
≤ >
(3) 球内均匀分布
设球的半径为 R，粒子在球内均匀分布时，从发射
点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
在直角坐标系下，抽样方法为,
??
?
?
?
???
其它
当
0
0
3
)( 3
2
Rr
R
r
rf
),,m ax ( 321 ????? Rr
302010
2
3
2
2
2
1
,,
1
???
???
??????
???
RzRyRx
≤
>
(4) 球壳内均匀分布
设球壳的内半径为 R0，外半径为 R1，在均匀分布时，
从发射点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
??
?
?
?
??
??
其它
当
0
3
)( 103
0
3
1
2
RrR
RR
r
rf
001
432322
2
110
2
0
10
1
2
110
2
0
2
0
1
)(
),,m a x (),m a x (
3
3
RxRRr
xxx
RRRR
RR
RRRR
R
????
???
??
?
??
?
??????
?
?
≤
≤
>
>
在直角坐标系下，球壳内点的坐标为,
其中,r 由前面的抽样方法确定,θ,φ服从各向同性
分布，其抽样方法为,
?
??
??
c o s
s i ns i n
c o ss i n
0
0
0
??
??
??
rz
ry
rx
＞
≤
2
3
22
2
22
1
2
3
22
2
22
1
0
2
3
22
2
22
1
31
0
2
3
22
2
22
1
21
0
1
22
3
22
2
22
1
c o s
2
s i ns i n
2
c o ss i n
)(
???
???
?
???
??
??
???
??
??
????
AA
AA
rrz
AA
A
rry
AA
A
rrx
AA
??
??
????
??
????
??
????
???
(5) 圆柱内均匀分布
圆柱内均匀分布是指粒子发射点均匀地分布在底
半径为 R，高为 2H 的圆柱内。若固定圆柱的中心为
原点，圆柱的轴向为 z 轴，则分布密度函数为,
抽样方法为,
??
?
?
? ???
??
其它
当
0
||,
2
1
),,(
222
2 HzRyxHRzyxf ?
302010
2
2
2
1
,,
1
???
??
??????
??
HzRyRx
≤
>
(6) 点源分布
点源分布是指粒子由一固定点 发射，
其分布密度函数为,
其中,为狄拉克 δ函数，源粒子的抽样方法为,
在球坐标系中，粒子发射点到球心的距离 r 的分布
密度函数为,
其中,为点源到球心的距离。源粒子的位置抽样为,
)()()(),,( *0*0*0 zzyyxxzyxf ?????? ???
),,( *0*0*0 zyx
*0*0*0,,zzyyxx ???
*0rr?
*0r
)()( *0rrrf ?? ?
)(??
(7) 球外平行束源分布
球外平行束源分布是指粒子平行入射到半径为 R
的球面上，或球外点源距离球很远，可以近似地看作
平行束源。设 r 为粒子发射点到球心的距离,其分布
密度函数为,
r 的抽样方法为,
在直角坐标系中，抽样方法为,
)()( Rrrf ?? ?
Rr ?
2
0
2
0
2
0
2010
2
2
2
1
,
1
zyRx
RzRy
????
????
??
??
??
≤
>
2) 源粒子的能量常见分布的随机抽样
(1) 单能源分布
单能源分布是指粒子的发射能量为一固定值 E0,
其分布密度函数为,
源粒子的能量为,
)()( 0EEEf ?? ?
0EE ?
(2) 裂变中子谱分布
裂变中子谱分布的一般形式为,
其中 A,B,C,Emin,Emax 均为与元素有关的量。
对于铀 -235,
A=0.965,B=2.29,C=0.453,Emin=0,Emax=∞。
m axm i n,sh)( EEEBEeCEf AE ????? ?
采用近似修正抽样，抽样方法为,
其中,m≈0.8746,M1≈0.2678,λ≈0.5543。
此外，裂变谱分布也有以数值曲线形式给出的，
此时，用数值曲线抽样方法抽取 E 。
332
1
1
2
1
ln
1
)l n (
)(
?
?
??
?
?
A
A
EEAE
M
EH
m
?
???????
?
?
?
≤
≤
＞
＞
? ?EAmEBECAAEH ?? ??????? ???? e x psh1)( 21
(3) 麦克斯韦 (Maxwell) 谱分布
麦克斯韦谱分布的一般形式为,
该分布的抽样方法为
0,2)(
23
??? ?? EeEEf E???
＞
2
22
2
1
ln
2
3
ln
?
?
???
??
???
E
e
≤
3) 源粒子运动方向常见分布的随机抽样
(1) 各向同性分布
各向同性分布密度函数为,
其中,μ＝ cosθ,θ为运动方向与 z 轴的夹角,φ为
方位角。
?
??
?
??
2
1
)(,
2
1
)(
4
1
)()()(
21
21
??
???
ff
fff Ω
?
在直角坐标系下，各方向余弦 u,v,w 为,
其抽样方法为,
?
??
??
c o s
s ins in
c o ss in
?
?
?
w
v
u
＞
≤
2
3
22
2
22
1
2
3
22
2
22
1
2
3
22
2
22
1
31
2
3
22
2
22
1
21
1
22
3
22
2
22
1
c o s
2
s i ns i n
2
c o ss i n
)(
???
???
?
???
??
??
???
??
??
????
AA
AA
w
AA
A
v
AA
A
u
AA
??
??
??
??
??
??
??
???
(2) 半面各向同性分布
不妨设在 x≥0 的半面方向上各向同性发射粒子，
则在前述各向同性分布的抽样方法中，用 ξ2代替 η2就
能得到所需分布的抽样。对于其它方向的情况，可用
类似的方法处理。
(3) 球外平行束源分布
令 μ＝ cosθ,θ为粒子运动方向的径向夹角，则 μ分
布密度函数为,
μ的抽样方法为,
01,2)( ????? ???f
),m a x ( 21 ????u
(4) 球外各向同性点源分布
设球外点源 S 到球心的距离为 D0。点源 S 到球的
最大张角为 θ*,
则球外各向同性点源分布的抽样方法是,
先抽样确定,再转换成 θ。
0
22
0*c o s
D
RD ???
??
R
DR ?
?
???
??
??
?????
22
0
2
*
s in
c o s
)c o s1(1c o s
在直角坐标系下，取
OS 为 z 轴，抽样方法为,
?
?
???
?
??
c o s
0
s in
w
v
u
4) 次级粒子的源分布
在有关次级粒子（如裂变中子，中子生成光子，
光子生成中子）的输运过程中，次级粒子源分布的抽
样方法，主要可分为以下两种,
(1) 直接生成法
可将生成的次级粒子的位置、能量、方向、权重
等参数直接作为源分布的抽样结果。也就是直接对生
成的次级粒子进行跟踪。这种方法比较简单、直观。
(2) 离散分布法
将生成的次级粒子的权重，按空间位置、能量、
方向分别记录，得到次级粒子的空间、能量、运动方
向的离散的近似分布。再根据该分布，利用各种抽样
技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、
记录。
当一个问题需要用两个以上的蒙卡程序处理时，
可采用这种方法。
2,空间、能量和运动方向的随机游动过程
粒子由 状态 Sm到状态 Sm+1时, 需要
确定粒子的空间位置 rm+1,能量 Em+1和
运动方向 Ωm+1。
1) 碰撞点位置的计算公式
设 rm 为粒子第 m 次碰撞点的位置,Ωm 为碰撞后
的运动方向，则粒子第 m+1 次碰撞点的位置 rm+1 为,
即
其中 为 的方向余弦,L 为两次碰撞点间
的距离。
mmm
mmm
mmm
wLzz
vLyy
uLxx
???
???
???
?
?
?
1
1
1
mmm L Ωrr ???? 1
mmmm wvu Ω),,(
L 的分布密度函数为,
由 f (L) 抽样确定 L 的方法通常有三种,
(1) 直接抽样方法
确定 L 的直接抽样方法是,
首先由自由程分布
中抽取 ρ
再由下列关系式解出 L 。
? ? 0,),(e x p),()( 01 ??????? ?? LdlElELf L mmmtmmt Ωrr
? ???? L mmmt dlEl0 ),( Ωr?
?? ?? ef )(
?? ln??
对于均匀介质，有
对于多层介质，如果
则
其中,为粒子由 rm 出发，沿 Ωm 方向在顺序经过的
第 i 个介质区域内走过的距离,为第 i 个介质
区域的宏观总截面 ( i =1,2,…, Imax )。 当
时，意味着粒子穿出系统。
)(
ln
)( mtmt EEL ?????
??
)(
)(
)()(
,
1
0
,1
0
0
,
1
0
,
mIt
I
i
mitiI
i
i
I
i
miti
I
i
miti
E
EL
LL
ELEL
?
????
???
????????
?
?
??
?
?
?
?
?
?
?
?
?
iL?
)(,mit E?
?
?
????
m a x
0
,)(
I
i
miti EL?
(2) 最大截面法
对于多层介质，或其他介质密度与位置有关的问题，
在求 ( i =1,2,…, Imax ) 时，如果系统形状复杂，
计算是非常烦杂的。在这种情况下，使用最大截面法
更方便。最大截面抽样方法为,
其中
iL?
),(m a x)(m a x,EE tt rr ???
1
m a x,
1
2
m a x,
1
11
1
)(
),(
)(
ln
0
LL
E
EL
E
LL
L
mt
mmmt
mt
?
?
???
?
?
??
?
Ωr
?
?
≤
＞
(3) 限制抽样法
当介质区域很小时，如使用直接抽样法抽取输运
长度，粒子很容易穿出介质，此时使用限制抽样法确
定自由程个数 ρ较好,ρ的分布密度函数为,
其中 Dm 为粒子由 rm 出发，沿 Ωm 方向到达区域边界的
自由程个数。 ρ的抽样方法是,
然后用直接抽样法中根据 ρ计算 L 的方法计算输运长
度 L 。此时，粒子的权重需乘以纠偏因子 。
??
?
?
?
??
?? ?
?
其它
当
0
0
1)( mD
D
e
e
f m ??
?
])1(1l n [ ?? ????? ? mDe
)1( mDe??
2) 碰撞后能量 Em+1的随机抽样
粒子在介质中发生碰撞后，首先要确定与哪种原
子核发生何种反应。粒子发生碰撞后（吸收除外）的
能量 Em+1 一般只与其碰撞前后运动方向的夹角（散射
角）有关。
粒子碰撞后常见的能量分布有下面几种情况。
(1) 裂变中子谱
中子引起原子核裂变反应时，裂变中子的能量服
从裂变谱分布。其抽样方法可参考 以前的介绍 。
(2) 中子弹性散射后能量的确定
中子弹性散射后，能量与质心系散射角 θC的关系
是,
能量与实验室系散射角 θL的关系是,
其中,A 为碰撞核的质量,。
或 确定后，即可求出 Em+1。
? ?22221 1)1( LLmm AA EE ?? ??????
? ?12)1( 221 ????? Cmm AAA EE ?
LLCC ???? co s,co s ??
LC ??
(3) 中子非弹性散射后能量的确定
中子非弹性散射后，能量与质心系散射角 θC的关
系是,
其中,为第 K 个能级的阈能,为第 K 个能级的激
发态能量。
如果确定了实验室系散射角 θL，则根据下式
确定 后，再计算出 Em+1。
? ?? ?
KK
CmKmK
m
m
A
A
EAEA
A
E
E
??
???
?
?
?
?????
?
??
1
1121
)1(
2
21
K? K?
C?
? ?? ?1111 1 222 ??????? LLmKL
mK
C EAEA ??????
(4) 光子康普顿（ Compton）散射后能量的确定
光子发生康普顿散射后，其能量分布密度函数为,
其中,K(α) 为归一因子。
, 和 分别为光子散射前后的能量，以
m0c2 为单位,m0为电子静止质量,c 为光速。
?
?
?
?
? 211,1111
)(
1)(
32
2
???
?
?
?
?
?
?
?
?
????
?
??
?
?
?
??? x
xxxx
x
K
xf
22 )21(2
14
2
1)21ln ()1(21)(
????
??
???????
?
??
? ???K
?? ??x ???
光子康普顿散射能量分布的抽样方法为,
x 的抽样确定后, 散射后的能量为,
2
3
2
2
2
3
22
1
1
2
1
3
2
1
)1(
4
27
21
294
27
1
1
2
1
21
21
xx
x
x
x
xx
x
x
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
? ??
?
?
?
?
?
??
?
?
?
?
?
??
?
＞ ＞
＞ ≤
≤ ≤
x
Ecm
xcmE
m
m ???????
2
0
2
01
??
3) 碰撞后散射角的随机抽样
粒子碰撞后运动方向 Ωm+1的确定，一般与散射角
有关。由已知分布抽样确定散射角后，再确定 Ωm+1。
常见的散射角分布有如下几种,
(1) 质心系各向同性分布
散射角在质心系服从各向同性分布时，其抽样方
法为 。质心系散射角 θC抽样确定后,
需转换成实验室系散射角 θL。
12 ??? ??? C
在中子弹性散射情况下, 转换公式为,
其中 A 为碰撞核质量, 。
在中子非弹性散射情况下, 转换公式为,
其中, 为第 K 个能级的阈能 。
C
C
L AA
A
?
??
21
1
2 ??
??
? ? CmKmK
CmK
L
EAEA
EA
???
???
?????
????
1211
11
2
LLCC ???? co s,co s ??
K?
(2) 中子弹性散射勒让德 (Legendre) 多项式分布
中子弹性散射角分布常以勒让德多项式的展开形
式给定。散射角余弦 x＝ cosθ的分布密 度函数为,
其中 Pl(x) 为 l 阶勒让德多项式。
该分布即为 n 阶勒让德近似展开。
勒让德多项式由以下递推公式确定,
??
?
?
?
???
? ??
其它
当
0
1||)(
2
12
)( 0
xxPfl
xf
n
l
ll
xxP
xP
xnPxxPnxPn nnn
?
?
????? ??
)(
1)(
0)()()12()()1(
1
0
11
考虑新的分布,
当选取 x0,x1,… xn 为 Pn+1(x)＝ 0 的根，且
时,fa(x) 依照勒让德多项式展开的前 n 项与 f (x) 的展
开形式相同。因此，可以用 fa(x) 作为 f (x) 的近似分布。
?
?
?
?
?
??
?
n
l
kl
n
l
kll
k
xPl
xPfl
Φ
0
2
0
)]([
2
12
)(
2
12
?
?
??
n
k
kka xxΦxf
0
)()( ?
在实际问题中，由于勒让德多项式展开项数不够，
可能出现某个 为负值的现象。此时可以采用如下近
似分布,
其中,
对于该近似分布，可用加抽样方法进行抽样,
此时，由于偏倚抽样而引起的纠偏因子为 wK,也就
是说，粒子的权重要乘上 wK。
?
?
?
?
?
?
???
??
n
k
k
k
k
kn
k
k
k
k
n
k
kka
Φ|
Φ|
Φ
w
Φ|
Φ|
Ψ
xxΨxf
0
0
0
*
|
|
,
|
|
)()( ?
kΦ
??
?
?
?
???
K
k
k
K
k
kK ΨΨxx
0
1
0
,?当
(3) 光子康普顿散射角分布
光子的康普顿散射角与其散射前后的能量有关,它
的分布密度函数为,
抽样方法为,
)111()( LLf ????? ?????
??? 111 ????L
4) 碰撞后运动方向 Ωm+1的确定
实验室系散射角 θL确定后，依据不同的坐标系的表
现形式，有不同的确定方法。
(1) 确定方向余弦 um+1,vm+1,wm+1
mmmm
m
mm
mmm
m
m
mm
mmm
m
awvubcw
av
vu
bduvb c w
v
au
vu
bdvub c w
u
???
?
?
??
?
?
?
??
?
?
?
?
22
1
22
1
22
1
其中,
方位角 在 [0,2π] 上均匀分布。
当 时，不能使用上述公式，可用下
面的简单公式,
??
??
s i n,c o s
,1s i n,c o s 2
??
????
dc
aba LL
022 ?? mm vu
?
mm
m
m
aww
bdv
bcu
?
?
?
?
?
?
1
1
1
(2) 确定 Ωm+1的球坐标 (θm+1,φm+1)
设 Ωm的球坐标分别为 (θm,φm)，其中,θ为粒子运动
方向与 z 轴的夹角,φ为粒子运动方向在 x y 平面上投
影的方位角。则 Ωm+1的球坐标 (θm+1,φm+1) 分别由下式
确定,
1
1
1
1
1
1
s i ns i n
c o sc o sc o s
)c o s (
s i n
s i ns i n
)s i n (
c o ss i ns i nc o sc o sc o s
?
?
?
?
?
?
?
??
??
??
mm
mmL
mm
m
L
mm
LmLmm
??
???
??
?
??
??
??????
5) 球形几何的随机游动公式
一般几何的随机游动公式可以应用到球形几何，
而对球对称问题，使用特殊形式更为方便。
(1) 下次碰撞点的径向位置 rm+1的确定
两次碰撞点间的距离 L 确定之后，下次碰撞点的
径向位置 rm+1的计算公式为,
设系统的外半径为 R，如 rm+1≥R，则粒子逃出系统。
mmmm LL ?c o s2221 ?????? rrr
(2) 粒子碰撞后瞬时运动方向的确定
在球对称系统中，粒子运动方向用其与径向夹角
余弦来描述。使用球面三角公式，粒子碰撞后瞬时运
动方向与径向夹角余弦 cosθm+1的计算公式为,
其中,为在 [0,2π] 上均匀分布的方位角,为在
rm+1 点进入碰撞前瞬时运动方向与 rm+1 径向之间的夹
角。
1
1
1
c o s
c o s
s i ns i n
c o ss i ns i nc o sc o sc o s
?
?
?
??
??
???
????
m
mm
m
m
m
m
m
LmLmm
L
r
r
r
r
?
?
??
??????
?
m??
6) 点到给定边界面的距离
在抽样确定输运距离、判断粒子是否穿透系统时,
常遇到求由 rm 出发，沿 Ωm 方向到达某个区域表面的
距离问题。在记录对结果的贡献时，也常使用类似的
量。区域表面通常是平面或二次曲面。 求到达区域表
面的距离问题，实际上是求直线（或半射线）与平面
或二次曲面的交点问题。这是 蒙特卡罗方法解粒子输
运的各种实际问题时,所遇到的基本几何问题。
(1) 点到平面的距离
点 沿方向 的直线方程为,
该直线到达方程为
的平面的距离为,
当 与平面平行时，即
直线与平面无交点。
如果 l 为负值，直线与平面也无交点。这时，粒子的
运动方向是背离平面的。
),,( 0000 wvu?Ω
0000 ??? cwbvau
000
000 )(
cwbvau
czbyaxdl
??
????
dczbyax ???
00 Ωrr ??? l
),,( 0000 zyx?r
(2) 点到球面的距离
在三维直角坐标系中，设球心为 rc＝ (xc,yc,zc),
球半径为 R，则球面方程为,
将直线方程代入该球面方程，得到点 r0沿 Ω0方向到达
球面的距离 l,
其中
2222 )()()( Rzzyyxx ccc ??????
???? ?l
2
0
2
0
2
0
2
0
2
0
22
0
2
00000000
)()()(||
)()()()(
cccc
cccc
zzyyxxR
RR
wzzvyyuxx
????????
????
?????????
rr
Ωrr
?
?
当 R0≤R 时，即 r0 点在球内, Δ≥δ2,l 只有一个正
根,
当 R0＞ R 时，即 r0 点在球外，分以下三种情况,
a) 若 δ≥0,l 无正实根，直线与球面无交点。
b) 若 δ＜ 0,Δ＜ 0,l 无实根，直线与球面无交点。
c) 若 δ＜ 0,Δ≥0,l 有两个正实根，直线与球面有两个交
点。
???? ?l
???? ?l
在球坐标系中，不失一般性，设球心为 rc＝ 0，则
球面方程为 r＝ R。
当 r0≤R 时，即 r0 点在球内,有一个交点,
其中 θ0为 Ω0与 r0 的径向夹角。
当 r0＞ R 时，即 r0 点在球外,令
当 cosθ0≥0 时，直线与球面无交点。
当 cosθ0＜ 0 时，若 d≥R，则直线与球面无交点。
若 d＜ R，则有两个交点,
200200 )s in(c o s ?? ?????? rRrl
200200 )s in(c o s ?? ?????? rRrl
00 s in ??? rd
(3) 点到圆柱面的距离
设圆柱面的方程为,
其中 (xc,yc,0) 为圆柱的中心,R 为圆柱底半径。
点 r0沿 Ω0方向到达圆柱面的距离 l 为,
其中
222 )()( Ryyxx cc ????
2
01 w
l ? ???? ?
2
0
2
0
2
0
2
0
22
0
2
0000
)()(
)1)((
)()(
cc
cc
yyxxR
wRR
vyyuxx
????
?????
????
?
?
当 R0≤R 时,r0 点在圆柱内,如果,则 l 有
一个正根,
如果,即 Ω0平行于圆柱的对称轴，直线与
圆柱面无交点。
当 R0＞ R 时,r0 点在圆柱外，分以下三种情况,
a) 若 δ≥0,l 无正实根，直线与圆柱面无交点。
b) 若 δ＜ 0,Δ＜ 0,l 无实根，直线与圆柱面无交点。
c) 若 δ＜ 0,Δ≥0 且, l 有两个正实根，直线与圆柱
面有两个交点。
在 的情况下，直线与圆柱面不相交。
2
01 w
l ? ???? ?
2
01 w
l ? ???? ?
120 ?w
120 ?w
120 ?w
120 ?w
(4) 点到圆锥面的距离
设圆锥顶点在原点，以 z 轴为对称轴，则圆锥面
的方程为,
点 r0沿 Ω0方向到达圆锥面的距离 l 为,
其中
如果 Ω0与锥面某一母线平行，即,则
2222 zcyx ??
2
0
2 )1(1 wcl ??
???? ?
])1(1)[( 20220220202
00
2
0000
wczcyx
wzcvyux
???????
???
?
?
?2
2
0
22
0
2
0 zcyxl ????
1)1( 202 ?? wc
(5) 空腔处理
在粒子输运问题中，所考虑的系统常有空腔存在，
如中空的球壳,平板间的空隙等。粒子输运时，可有
两种处理空腔的方法,
a) 将空腔作为宏观总截面 Σt＝ 0 的区域,按通常的方法输
运。
b) 设 分别为由 rm 出发，沿 Ωm 方向到空腔区域的
近端和远端的交点，当粒子超过 时，以 为新的起
点，重新开始输运。
显然，这两种方法在统计上是等价的。
mm rr ???、
mr? mr?
7) 等效的边界条件
(1) 全反射边界
在反应堆活性区中，元件盒常常按正方形或六角
形排列。假定元件盒足够多，每个盒结构相同，那么
活性区中每个盒所占的栅胞的物理情况，可以代表整
个活性区中的状况。
进一步假定，元件盒是圆对称的，那么每个栅胞
中情况，可以用更小的单位（栅元）来反映。 比如对
六角形栅胞可取其 1/12 的 ΔOAB 来做代表； 正方形 栅
胞可用其 1/8 的 ΔO'A'B' 来做代表。这样一来问题就大
大简化 了。
现在的问题是怎样计算直角三角形栅元的物理量
（如通量）。用蒙特卡罗方法如何模拟中子在栅元内
的运动，反映出整个活性区对它的影响。
我们可把 OA',OB',A'B' 作为全反射边界来处理。
所谓全反射边界，就是当中子打到该边界上时，按镜
面反射的方式，从边界 上全部反射回来，中子的能量
与权重均不改变。
在这种边界上的反射条件，称之为全反射条件，
就是通常的镜面反射条件。
在全反射边界条件下，一条通过活性区若干个区
域的中子径迹，可以用栅元 ΔO'A'B' 中的一条折线轨
道来反映出来。
反过来，在直角三角形栅元 ΔO'A'B' 中任一条反
射成的折线轨道，都代表了中子在活性区内一条直线
轨道的作用。由于系统的对称性，在活性区内，凡是
与栅元内位置相当的地方，都有相同的物理情况，因
此栅元内各处的情况，当然代表了整个活性区的情况。
(2) 一般曲面全反射条件
对于一般曲面的全反射，设入射方向为 Ω，入射
点的内法线方向为 n,则反射方向 Ω' 为,
其中
设
则
? ?
nΩ
nnΩΩΩ
???
?????
?c o s2
2
nΩ ????c o s
),,( zyx nnn?n
)(c o s
c o s2
c o s2
c o s2
zyx
z
y
x
nwnvnu
nww
nvv
nuu
???????
????
????
????
?
?
?
?
(3) 平面全反射条件
设三角形栅元的横截面 ΔOAB 在 X-Y 平面上，
∠ OAB＝ θ。则边界 OA,OB,AB 上的反射都是平面
全反射。在任一与 X-Y 平面垂直且与 X 轴成 α角的平
面上，全反射条件为,
由此就可得到 OA,OB 和 AB 边上的全反射条件，对
于 OB 边,α=θ；对于 OA 边,α= 0；对于 AB 边，
α=π/ 2。
ww
vuv
vuu
??
?????
?????
??
??
2c o s2s in
2s in2c o s
(4) 反射层边界条件
对于具有大反射层的系统，如存放，运输和生产
裂变物质的仓 库、车厢和车间等，当中子从里面打到
四周墙上或反射层时，还要继续对它进行跟踪。这种
跟踪常常要花费很大的计算量，并且在结果中引起的
方差也比较大。如果在计算这种系统的不同方案中，
反 射层条件不变，那么这种大量重复的计算是很不经
济的。
中子射入反射层后，一部分被介质吸收，
只有一部分返回，由于中子的散射慢化，损失
一部分能量，因此反射回来的中子有一个能量
方向分布。显然，对这种反射层，不能应用全
反射条件。不过，我们仍然可以把它当做边界，
在边界上按反射层的物理作用来处理。
比如，如果反射层是一种平板几何，我们可以用
数值方法或蒙特卡罗方法，预先算好在各种不同入射
能量 E 下的反照率 β(E)，反射中子的能量分布
RE(E→ E' )。于是代替在反射层中眼踪中子，我们可
在反射层边界上作如下处理,
一旦中子打入反射 层，立即返回，反射后权重为
其中,E 为射入反射层中子的能量,W 为中子的权重。
反射后的能量 Eβ 由反射能谱 RE(E→ Eβ) 中抽样产生。
反射后的方向 Ωβ 由半平面各向同性分布或余弦分布中
抽样。反射后的中子位置为入射时的位置。
计算表明，对于大尺寸的反射层来说，这样的近
似，引 起的结果上的误差是可以忽略的，却能带来计
算量的大量节省。
)( EW ?? ??W
3,记录贡献与分析结果过程
在粒子输运问题中, 除了得到某些量的积
分结果外, 还需要得到这些量的方差, 协方差,
以及这些量的空间, 能量, 方向和时间的分布 。
这些量可以利用分类记录手续同时得到 。
1) 记录与结果
为了得到所求量的估计值，在粒子输运过程中需
进行记录，即求每个粒子对所求量的贡 献。
设模拟了 N 个粒子，所求量的估计值为,
其中 gi 为第 i 个粒子的总贡献。
?
?
?
N
i
iN gNg
1
1
记录的贡献由所求量决定。对于同一个所求量，
又随所用的蒙特卡罗技巧的不同而不同。 例如，所求
量是粒子穿透屏蔽概率，使用直接模拟法时，如粒子
穿透屏蔽，在叠加记录单元加,1” ( 初始值为零 )，否
则没有贡献。使用加权法时，如粒子穿透屏蔽，在叠
加记录单元加粒子的权重，否则没有贡献。使用统计
估计法时，粒子每发生一次碰撞 (包括零次碰撞 )，都
要记录贡献，等等。
2) 方差和协方差的估计
估计量 g 和 g' 的方差和协方差为,
它们可以用下式估计,
?
?
?
?
?
?
??
?
?
?
?
?
???
?
?
?
?
?
?
??
???
??
???
?
??
N
i
i
N
i
i
N
i
iigg
N
i
i
N
i
ig
g
N
g
N
gg
N
g
N
g
N
111
2
2
11
22
111
11
?
?
? ? ? ?
? ? ? ?? ?gEEgggE
EggE
gg
g
????
??
?
2
222
?
?
因此，要得到 和 的估计，只要对每一个历史记
录结果的 和 进行记录，并加以累加即可。
方差估计值 确定后，可得到误差
其中 为置信限，它随置信水平 而定。在通常
情况下，取 。
22 ggg ???
iii ggg ?2
2g?
N
g??? ??
??
96.1,95.01 ??? ???
??1
3) 位置、能量、方向、时间分布
在前面已经提到，用蒙特卡罗方法求某种
量的空间、能量、方向和时间分布，实质上是
得到这种分布的阶梯函数近似的估计值。而求
这种估计值是很方便的，只要将跟踪过程中所
得到的感兴趣量，按其状态的空间、能量、方
向、时间特征，分别记录其权 重，最后将这
些记录结果适当处理即可。
事先，将问题的空间、能量、方向（常按相对于
某个方向的夹角余弦）、时间范围，各分为如下不同
间隔,
再用一批存贮单元 {A} 记录相应间隔上阶梯函数
近似的累计值。;0;11;,,,
10
10
m a x01m i n
21
Tttt
EEEEE
VVV
L
K
J
I
?????
??????
?????
?
?
?
?
???
4,核截面数据的引用
用蒙特卡罗方法解粒子输运问题, 需要介质所包
含的各种原子核的核数据 。 以中子核数据为例, 需要
各种涉及到的核的微观总截面, 弹性散射, 非弹性散
射, n-2n 反应, 裂变, 俘获等截面；也需要这些反应
的相应能量, 角度分布, 次级粒子数, 以及其它关心
的粒子数及其能量, 方向分布 。 从输运方程中可以知
道, 有了这些数据, 问题就完全确定了；反映到蒙特
卡罗模拟中, 有了这些数据, 就能决定宏观总截面,
决定碰撞核的具体形式, 就能实现抽样和跟踪 。
在蒙特卡罗计算中, 引用的核数据有点截面, 分
段常数截面和群截面三种形式 。
1) 点截面形式
在跟踪粒子时，对粒子的每一种能量，先
从截面库中取出需要的核数据，再用插值（或
其它方式），求出相应能量的各种截面数据。
这种方法是直接的，也是比较精确的。不少通
用程序就是这样做的。这样做的条件是要有完
备的核截面数据库，计算机有大而快的存贮能
力。
2) 分段常数形式
将问题的能量范围（ Emin,Emax）分成许
多间隔，截面数据在每个间隔上看作与能量无
关，即截面取分段常数形式。这种引用截面数
据的好处是，数据量相对地少。问题是它要根
据问题的物理特点来划分能量间隔。而为了保
证误差较小，所取的间隔数一般是比较多的。
3) 群截面形式
解中子输运问题，常常采用多群近似方法。在多
群近似下，中子能量 E 用群指标 i 代替。为了实现蒙
特卡罗跟踪，只需要以下群截面,
显然，这种跟踪过程数据量小，程序简单。
群时的角分布；群散射到第点由第在
群的截面；群散射到第点由第在
群吸收截面；点的第在
群总截面；点的第在
jirr
jirr
irr
irr
:),|(
:)(
:)(
:)(
jif s
ji
s
i
a
i
t
?
?
?
?
?
?
5,蒙特卡罗程序结构
在电子计算机上, 蒙特卡罗方法解粒子输运问题
的程序, 一般都可分为：源抽样, 空间输运过程, 碰
撞过程, 记录过程和结果的处理与输出等部分 。
开始
数据预处理,各记录单元清零
取一个粒子历史
源分布抽样
输运过程
碰撞过程
历史终止否？
统计处理
做完给定历史数否？
结果的处理与输出
终止
记录过程
记录过程
记录过程
记录过程
至于粒子历史终止条件，根据问题的几何条件、
物理假定，处理方法，可归纳为以下几种,
1) 粒子从系统逃脱；
2) 粒子经碰撞被吸收；
3) 经俄国轮盘赌后，历史被终止；
4) 粒子能量低于给定能量（阈能） ;
5) 粒子位置越过某一界面；
6) 粒子飞行时间超过给定时间；
7) 粒子权重小于某个小量。
? 作 业
1)
1,源分布抽样过程
2,空间、能量和运动方向的随机游动过程
3,记录贡献和分析结果过程
4,核截面数据的引用
5,蒙特卡罗程序结构
? 作 业
第五章 蒙特卡罗方法在计算机上的实现
蒙特卡罗方法是随着计算机的出现和发展
而逐步发展起来的。在计算机上能够产生符合
要求的随机数，实现对已知分布的抽样，奠定
了蒙特卡罗方法在计算机上得以实现的基础。
在计算机上使用蒙特卡罗方法解粒子输运问题
大致包括三个过程：源分布抽样过程，空间、
能量和运动方向的随机游动过程以及记录、分
析结果过程 。
1,源分布抽样过程
源分布抽样的目的是产生粒子的初始状态
。 下面我们介绍一些常见的特定
类型的源分布抽样方法 。
),,( 0000 ΩrS E?
1) 源粒子的位置常见分布的随机抽样
(1) 圆内均匀分布
设圆半径为 R0，粒子在圆内均匀分布时，从发射
点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
??
?
?
? ??
?
其它
当
0
02
)( 020
Rr
R
r
rf
),m ax ( 210 ???? Rr
(2) 圆环内均匀分布
设圆环的内半径为 R0，外半径为 R1，则粒子在该圆
环内均匀分布时，从发射点到中心的距离 r 的分布密
度函数为,
r 的抽样方法为,
??
?
?
? ??
??
其它
当
0
2
)( 102021
RrR
RR
r
rf
020103201
01
01
1
)(),m a x ()( RRRrRRRr
RR
RR
????????
?
?
?
???
?
≤ >
(3) 球内均匀分布
设球的半径为 R，粒子在球内均匀分布时，从发射
点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
在直角坐标系下，抽样方法为,
??
?
?
?
???
其它
当
0
0
3
)( 3
2
Rr
R
r
rf
),,m ax ( 321 ????? Rr
302010
2
3
2
2
2
1
,,
1
???
???
??????
???
RzRyRx
≤
>
(4) 球壳内均匀分布
设球壳的内半径为 R0，外半径为 R1，在均匀分布时，
从发射点到中心的距离 r 的分布密度函数为,
r 的抽样方法为,
??
?
?
?
??
??
其它
当
0
3
)( 103
0
3
1
2
RrR
RR
r
rf
001
432322
2
110
2
0
10
1
2
110
2
0
2
0
1
)(
),,m a x (),m a x (
3
3
RxRRr
xxx
RRRR
RR
RRRR
R
????
???
??
?
??
?
??????
?
?
≤
≤
>
>
在直角坐标系下，球壳内点的坐标为,
其中,r 由前面的抽样方法确定,θ,φ服从各向同性
分布，其抽样方法为,
?
??
??
c o s
s i ns i n
c o ss i n
0
0
0
??
??
??
rz
ry
rx
＞
≤
2
3
22
2
22
1
2
3
22
2
22
1
0
2
3
22
2
22
1
31
0
2
3
22
2
22
1
21
0
1
22
3
22
2
22
1
c o s
2
s i ns i n
2
c o ss i n
)(
???
???
?
???
??
??
???
??
??
????
AA
AA
rrz
AA
A
rry
AA
A
rrx
AA
??
??
????
??
????
??
????
???
(5) 圆柱内均匀分布
圆柱内均匀分布是指粒子发射点均匀地分布在底
半径为 R，高为 2H 的圆柱内。若固定圆柱的中心为
原点，圆柱的轴向为 z 轴，则分布密度函数为,
抽样方法为,
??
?
?
? ???
??
其它
当
0
||,
2
1
),,(
222
2 HzRyxHRzyxf ?
302010
2
2
2
1
,,
1
???
??
??????
??
HzRyRx
≤
>
(6) 点源分布
点源分布是指粒子由一固定点 发射，
其分布密度函数为,
其中,为狄拉克 δ函数，源粒子的抽样方法为,
在球坐标系中，粒子发射点到球心的距离 r 的分布
密度函数为,
其中,为点源到球心的距离。源粒子的位置抽样为,
)()()(),,( *0*0*0 zzyyxxzyxf ?????? ???
),,( *0*0*0 zyx
*0*0*0,,zzyyxx ???
*0rr?
*0r
)()( *0rrrf ?? ?
)(??
(7) 球外平行束源分布
球外平行束源分布是指粒子平行入射到半径为 R
的球面上，或球外点源距离球很远，可以近似地看作
平行束源。设 r 为粒子发射点到球心的距离,其分布
密度函数为,
r 的抽样方法为,
在直角坐标系中，抽样方法为,
)()( Rrrf ?? ?
Rr ?
2
0
2
0
2
0
2010
2
2
2
1
,
1
zyRx
RzRy
????
????
??
??
??
≤
>
2) 源粒子的能量常见分布的随机抽样
(1) 单能源分布
单能源分布是指粒子的发射能量为一固定值 E0,
其分布密度函数为,
源粒子的能量为,
)()( 0EEEf ?? ?
0EE ?
(2) 裂变中子谱分布
裂变中子谱分布的一般形式为,
其中 A,B,C,Emin,Emax 均为与元素有关的量。
对于铀 -235,
A=0.965,B=2.29,C=0.453,Emin=0,Emax=∞。
m axm i n,sh)( EEEBEeCEf AE ????? ?
采用近似修正抽样，抽样方法为,
其中,m≈0.8746,M1≈0.2678,λ≈0.5543。
此外，裂变谱分布也有以数值曲线形式给出的，
此时，用数值曲线抽样方法抽取 E 。
332
1
1
2
1
ln
1
)l n (
)(
?
?
??
?
?
A
A
EEAE
M
EH
m
?
???????
?
?
?
≤
≤
＞
＞
? ?EAmEBECAAEH ?? ??????? ???? e x psh1)( 21
(3) 麦克斯韦 (Maxwell) 谱分布
麦克斯韦谱分布的一般形式为,
该分布的抽样方法为
0,2)(
23
??? ?? EeEEf E???
＞
2
22
2
1
ln
2
3
ln
?
?
???
??
???
E
e
≤
3) 源粒子运动方向常见分布的随机抽样
(1) 各向同性分布
各向同性分布密度函数为,
其中,μ＝ cosθ,θ为运动方向与 z 轴的夹角,φ为
方位角。
?
??
?
??
2
1
)(,
2
1
)(
4
1
)()()(
21
21
??
???
ff
fff Ω
?
在直角坐标系下，各方向余弦 u,v,w 为,
其抽样方法为,
?
??
??
c o s
s ins in
c o ss in
?
?
?
w
v
u
＞
≤
2
3
22
2
22
1
2
3
22
2
22
1
2
3
22
2
22
1
31
2
3
22
2
22
1
21
1
22
3
22
2
22
1
c o s
2
s i ns i n
2
c o ss i n
)(
???
???
?
???
??
??
???
??
??
????
AA
AA
w
AA
A
v
AA
A
u
AA
??
??
??
??
??
??
??
???
(2) 半面各向同性分布
不妨设在 x≥0 的半面方向上各向同性发射粒子，
则在前述各向同性分布的抽样方法中，用 ξ2代替 η2就
能得到所需分布的抽样。对于其它方向的情况，可用
类似的方法处理。
(3) 球外平行束源分布
令 μ＝ cosθ,θ为粒子运动方向的径向夹角，则 μ分
布密度函数为,
μ的抽样方法为,
01,2)( ????? ???f
),m a x ( 21 ????u
(4) 球外各向同性点源分布
设球外点源 S 到球心的距离为 D0。点源 S 到球的
最大张角为 θ*,
则球外各向同性点源分布的抽样方法是,
先抽样确定,再转换成 θ。
0
22
0*c o s
D
RD ???
??
R
DR ?
?
???
??
??
?????
22
0
2
*
s in
c o s
)c o s1(1c o s
在直角坐标系下，取
OS 为 z 轴，抽样方法为,
?
?
???
?
??
c o s
0
s in
w
v
u
4) 次级粒子的源分布
在有关次级粒子（如裂变中子，中子生成光子，
光子生成中子）的输运过程中，次级粒子源分布的抽
样方法，主要可分为以下两种,
(1) 直接生成法
可将生成的次级粒子的位置、能量、方向、权重
等参数直接作为源分布的抽样结果。也就是直接对生
成的次级粒子进行跟踪。这种方法比较简单、直观。
(2) 离散分布法
将生成的次级粒子的权重，按空间位置、能量、
方向分别记录，得到次级粒子的空间、能量、运动方
向的离散的近似分布。再根据该分布，利用各种抽样
技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、
记录。
当一个问题需要用两个以上的蒙卡程序处理时，
可采用这种方法。
2,空间、能量和运动方向的随机游动过程
粒子由 状态 Sm到状态 Sm+1时, 需要
确定粒子的空间位置 rm+1,能量 Em+1和
运动方向 Ωm+1。
1) 碰撞点位置的计算公式
设 rm 为粒子第 m 次碰撞点的位置,Ωm 为碰撞后
的运动方向，则粒子第 m+1 次碰撞点的位置 rm+1 为,
即
其中 为 的方向余弦,L 为两次碰撞点间
的距离。
mmm
mmm
mmm
wLzz
vLyy
uLxx
???
???
???
?
?
?
1
1
1
mmm L Ωrr ???? 1
mmmm wvu Ω),,(
L 的分布密度函数为,
由 f (L) 抽样确定 L 的方法通常有三种,
(1) 直接抽样方法
确定 L 的直接抽样方法是,
首先由自由程分布
中抽取 ρ
再由下列关系式解出 L 。
? ? 0,),(e x p),()( 01 ??????? ?? LdlElELf L mmmtmmt Ωrr
? ???? L mmmt dlEl0 ),( Ωr?
?? ?? ef )(
?? ln??
对于均匀介质，有
对于多层介质，如果
则
其中,为粒子由 rm 出发，沿 Ωm 方向在顺序经过的
第 i 个介质区域内走过的距离,为第 i 个介质
区域的宏观总截面 ( i =1,2,…, Imax )。 当
时，意味着粒子穿出系统。
)(
ln
)( mtmt EEL ?????
??
)(
)(
)()(
,
1
0
,1
0
0
,
1
0
,
mIt
I
i
mitiI
i
i
I
i
miti
I
i
miti
E
EL
LL
ELEL
?
????
???
????????
?
?
??
?
?
?
?
?
?
?
?
?
iL?
)(,mit E?
?
?
????
m a x
0
,)(
I
i
miti EL?
(2) 最大截面法
对于多层介质，或其他介质密度与位置有关的问题，
在求 ( i =1,2,…, Imax ) 时，如果系统形状复杂，
计算是非常烦杂的。在这种情况下，使用最大截面法
更方便。最大截面抽样方法为,
其中
iL?
),(m a x)(m a x,EE tt rr ???
1
m a x,
1
2
m a x,
1
11
1
)(
),(
)(
ln
0
LL
E
EL
E
LL
L
mt
mmmt
mt
?
?
???
?
?
??
?
Ωr
?
?
≤
＞
(3) 限制抽样法
当介质区域很小时，如使用直接抽样法抽取输运
长度，粒子很容易穿出介质，此时使用限制抽样法确
定自由程个数 ρ较好,ρ的分布密度函数为,
其中 Dm 为粒子由 rm 出发，沿 Ωm 方向到达区域边界的
自由程个数。 ρ的抽样方法是,
然后用直接抽样法中根据 ρ计算 L 的方法计算输运长
度 L 。此时，粒子的权重需乘以纠偏因子 。
??
?
?
?
??
?? ?
?
其它
当
0
0
1)( mD
D
e
e
f m ??
?
])1(1l n [ ?? ????? ? mDe
)1( mDe??
2) 碰撞后能量 Em+1的随机抽样
粒子在介质中发生碰撞后，首先要确定与哪种原
子核发生何种反应。粒子发生碰撞后（吸收除外）的
能量 Em+1 一般只与其碰撞前后运动方向的夹角（散射
角）有关。
粒子碰撞后常见的能量分布有下面几种情况。
(1) 裂变中子谱
中子引起原子核裂变反应时，裂变中子的能量服
从裂变谱分布。其抽样方法可参考 以前的介绍 。
(2) 中子弹性散射后能量的确定
中子弹性散射后，能量与质心系散射角 θC的关系
是,
能量与实验室系散射角 θL的关系是,
其中,A 为碰撞核的质量,。
或 确定后，即可求出 Em+1。
? ?22221 1)1( LLmm AA EE ?? ??????
? ?12)1( 221 ????? Cmm AAA EE ?
LLCC ???? co s,co s ??
LC ??
(3) 中子非弹性散射后能量的确定
中子非弹性散射后，能量与质心系散射角 θC的关
系是,
其中,为第 K 个能级的阈能,为第 K 个能级的激
发态能量。
如果确定了实验室系散射角 θL，则根据下式
确定 后，再计算出 Em+1。
? ?? ?
KK
CmKmK
m
m
A
A
EAEA
A
E
E
??
???
?
?
?
?????
?
??
1
1121
)1(
2
21
K? K?
C?
? ?? ?1111 1 222 ??????? LLmKL
mK
C EAEA ??????
(4) 光子康普顿（ Compton）散射后能量的确定
光子发生康普顿散射后，其能量分布密度函数为,
其中,K(α) 为归一因子。
, 和 分别为光子散射前后的能量，以
m0c2 为单位,m0为电子静止质量,c 为光速。
?
?
?
?
? 211,1111
)(
1)(
32
2
???
?
?
?
?
?
?
?
?
????
?
??
?
?
?
??? x
xxxx
x
K
xf
22 )21(2
14
2
1)21ln ()1(21)(
????
??
???????
?
??
? ???K
?? ??x ???
光子康普顿散射能量分布的抽样方法为,
x 的抽样确定后, 散射后的能量为,
2
3
2
2
2
3
22
1
1
2
1
3
2
1
)1(
4
27
21
294
27
1
1
2
1
21
21
xx
x
x
x
xx
x
x
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
? ??
?
?
?
?
?
??
?
?
?
?
?
??
?
＞ ＞
＞ ≤
≤ ≤
x
Ecm
xcmE
m
m ???????
2
0
2
01
??
3) 碰撞后散射角的随机抽样
粒子碰撞后运动方向 Ωm+1的确定，一般与散射角
有关。由已知分布抽样确定散射角后，再确定 Ωm+1。
常见的散射角分布有如下几种,
(1) 质心系各向同性分布
散射角在质心系服从各向同性分布时，其抽样方
法为 。质心系散射角 θC抽样确定后,
需转换成实验室系散射角 θL。
12 ??? ??? C
在中子弹性散射情况下, 转换公式为,
其中 A 为碰撞核质量, 。
在中子非弹性散射情况下, 转换公式为,
其中, 为第 K 个能级的阈能 。
C
C
L AA
A
?
??
21
1
2 ??
??
? ? CmKmK
CmK
L
EAEA
EA
???
???
?????
????
1211
11
2
LLCC ???? co s,co s ??
K?
(2) 中子弹性散射勒让德 (Legendre) 多项式分布
中子弹性散射角分布常以勒让德多项式的展开形
式给定。散射角余弦 x＝ cosθ的分布密 度函数为,
其中 Pl(x) 为 l 阶勒让德多项式。
该分布即为 n 阶勒让德近似展开。
勒让德多项式由以下递推公式确定,
??
?
?
?
???
? ??
其它
当
0
1||)(
2
12
)( 0
xxPfl
xf
n
l
ll
xxP
xP
xnPxxPnxPn nnn
?
?
????? ??
)(
1)(
0)()()12()()1(
1
0
11
考虑新的分布,
当选取 x0,x1,… xn 为 Pn+1(x)＝ 0 的根，且
时,fa(x) 依照勒让德多项式展开的前 n 项与 f (x) 的展
开形式相同。因此，可以用 fa(x) 作为 f (x) 的近似分布。
?
?
?
?
?
??
?
n
l
kl
n
l
kll
k
xPl
xPfl
Φ
0
2
0
)]([
2
12
)(
2
12
?
?
??
n
k
kka xxΦxf
0
)()( ?
在实际问题中，由于勒让德多项式展开项数不够，
可能出现某个 为负值的现象。此时可以采用如下近
似分布,
其中,
对于该近似分布，可用加抽样方法进行抽样,
此时，由于偏倚抽样而引起的纠偏因子为 wK,也就
是说，粒子的权重要乘上 wK。
?
?
?
?
?
?
???
??
n
k
k
k
k
kn
k
k
k
k
n
k
kka
Φ|
Φ|
Φ
w
Φ|
Φ|
Ψ
xxΨxf
0
0
0
*
|
|
,
|
|
)()( ?
kΦ
??
?
?
?
???
K
k
k
K
k
kK ΨΨxx
0
1
0
,?当
(3) 光子康普顿散射角分布
光子的康普顿散射角与其散射前后的能量有关,它
的分布密度函数为,
抽样方法为,
)111()( LLf ????? ?????
??? 111 ????L
4) 碰撞后运动方向 Ωm+1的确定
实验室系散射角 θL确定后，依据不同的坐标系的表
现形式，有不同的确定方法。
(1) 确定方向余弦 um+1,vm+1,wm+1
mmmm
m
mm
mmm
m
m
mm
mmm
m
awvubcw
av
vu
bduvb c w
v
au
vu
bdvub c w
u
???
?
?
??
?
?
?
??
?
?
?
?
22
1
22
1
22
1
其中,
方位角 在 [0,2π] 上均匀分布。
当 时，不能使用上述公式，可用下
面的简单公式,
??
??
s i n,c o s
,1s i n,c o s 2
??
????
dc
aba LL
022 ?? mm vu
?
mm
m
m
aww
bdv
bcu
?
?
?
?
?
?
1
1
1
(2) 确定 Ωm+1的球坐标 (θm+1,φm+1)
设 Ωm的球坐标分别为 (θm,φm)，其中,θ为粒子运动
方向与 z 轴的夹角,φ为粒子运动方向在 x y 平面上投
影的方位角。则 Ωm+1的球坐标 (θm+1,φm+1) 分别由下式
确定,
1
1
1
1
1
1
s i ns i n
c o sc o sc o s
)c o s (
s i n
s i ns i n
)s i n (
c o ss i ns i nc o sc o sc o s
?
?
?
?
?
?
?
??
??
??
mm
mmL
mm
m
L
mm
LmLmm
??
???
??
?
??
??
??????
5) 球形几何的随机游动公式
一般几何的随机游动公式可以应用到球形几何，
而对球对称问题，使用特殊形式更为方便。
(1) 下次碰撞点的径向位置 rm+1的确定
两次碰撞点间的距离 L 确定之后，下次碰撞点的
径向位置 rm+1的计算公式为,
设系统的外半径为 R，如 rm+1≥R，则粒子逃出系统。
mmmm LL ?c o s2221 ?????? rrr
(2) 粒子碰撞后瞬时运动方向的确定
在球对称系统中，粒子运动方向用其与径向夹角
余弦来描述。使用球面三角公式，粒子碰撞后瞬时运
动方向与径向夹角余弦 cosθm+1的计算公式为,
其中,为在 [0,2π] 上均匀分布的方位角,为在
rm+1 点进入碰撞前瞬时运动方向与 rm+1 径向之间的夹
角。
1
1
1
c o s
c o s
s i ns i n
c o ss i ns i nc o sc o sc o s
?
?
?
??
??
???
????
m
mm
m
m
m
m
m
LmLmm
L
r
r
r
r
?
?
??
??????
?
m??
6) 点到给定边界面的距离
在抽样确定输运距离、判断粒子是否穿透系统时,
常遇到求由 rm 出发，沿 Ωm 方向到达某个区域表面的
距离问题。在记录对结果的贡献时，也常使用类似的
量。区域表面通常是平面或二次曲面。 求到达区域表
面的距离问题，实际上是求直线（或半射线）与平面
或二次曲面的交点问题。这是 蒙特卡罗方法解粒子输
运的各种实际问题时,所遇到的基本几何问题。
(1) 点到平面的距离
点 沿方向 的直线方程为,
该直线到达方程为
的平面的距离为,
当 与平面平行时，即
直线与平面无交点。
如果 l 为负值，直线与平面也无交点。这时，粒子的
运动方向是背离平面的。
),,( 0000 wvu?Ω
0000 ??? cwbvau
000
000 )(
cwbvau
czbyaxdl
??
????
dczbyax ???
00 Ωrr ??? l
),,( 0000 zyx?r
(2) 点到球面的距离
在三维直角坐标系中，设球心为 rc＝ (xc,yc,zc),
球半径为 R，则球面方程为,
将直线方程代入该球面方程，得到点 r0沿 Ω0方向到达
球面的距离 l,
其中
2222 )()()( Rzzyyxx ccc ??????
???? ?l
2
0
2
0
2
0
2
0
2
0
22
0
2
00000000
)()()(||
)()()()(
cccc
cccc
zzyyxxR
RR
wzzvyyuxx
????????
????
?????????
rr
Ωrr
?
?
当 R0≤R 时，即 r0 点在球内, Δ≥δ2,l 只有一个正
根,
当 R0＞ R 时，即 r0 点在球外，分以下三种情况,
a) 若 δ≥0,l 无正实根，直线与球面无交点。
b) 若 δ＜ 0,Δ＜ 0,l 无实根，直线与球面无交点。
c) 若 δ＜ 0,Δ≥0,l 有两个正实根，直线与球面有两个交
点。
???? ?l
???? ?l
在球坐标系中，不失一般性，设球心为 rc＝ 0，则
球面方程为 r＝ R。
当 r0≤R 时，即 r0 点在球内,有一个交点,
其中 θ0为 Ω0与 r0 的径向夹角。
当 r0＞ R 时，即 r0 点在球外,令
当 cosθ0≥0 时，直线与球面无交点。
当 cosθ0＜ 0 时，若 d≥R，则直线与球面无交点。
若 d＜ R，则有两个交点,
200200 )s in(c o s ?? ?????? rRrl
200200 )s in(c o s ?? ?????? rRrl
00 s in ??? rd
(3) 点到圆柱面的距离
设圆柱面的方程为,
其中 (xc,yc,0) 为圆柱的中心,R 为圆柱底半径。
点 r0沿 Ω0方向到达圆柱面的距离 l 为,
其中
222 )()( Ryyxx cc ????
2
01 w
l ? ???? ?
2
0
2
0
2
0
2
0
22
0
2
0000
)()(
)1)((
)()(
cc
cc
yyxxR
wRR
vyyuxx
????
?????
????
?
?
当 R0≤R 时,r0 点在圆柱内,如果,则 l 有
一个正根,
如果,即 Ω0平行于圆柱的对称轴，直线与
圆柱面无交点。
当 R0＞ R 时,r0 点在圆柱外，分以下三种情况,
a) 若 δ≥0,l 无正实根，直线与圆柱面无交点。
b) 若 δ＜ 0,Δ＜ 0,l 无实根，直线与圆柱面无交点。
c) 若 δ＜ 0,Δ≥0 且, l 有两个正实根，直线与圆柱
面有两个交点。
在 的情况下，直线与圆柱面不相交。
2
01 w
l ? ???? ?
2
01 w
l ? ???? ?
120 ?w
120 ?w
120 ?w
120 ?w
(4) 点到圆锥面的距离
设圆锥顶点在原点，以 z 轴为对称轴，则圆锥面
的方程为,
点 r0沿 Ω0方向到达圆锥面的距离 l 为,
其中
如果 Ω0与锥面某一母线平行，即,则
2222 zcyx ??
2
0
2 )1(1 wcl ??
???? ?
])1(1)[( 20220220202
00
2
0000
wczcyx
wzcvyux
???????
???
?
?
?2
2
0
22
0
2
0 zcyxl ????
1)1( 202 ?? wc
(5) 空腔处理
在粒子输运问题中，所考虑的系统常有空腔存在，
如中空的球壳,平板间的空隙等。粒子输运时，可有
两种处理空腔的方法,
a) 将空腔作为宏观总截面 Σt＝ 0 的区域,按通常的方法输
运。
b) 设 分别为由 rm 出发，沿 Ωm 方向到空腔区域的
近端和远端的交点，当粒子超过 时，以 为新的起
点，重新开始输运。
显然，这两种方法在统计上是等价的。
mm rr ???、
mr? mr?
7) 等效的边界条件
(1) 全反射边界
在反应堆活性区中，元件盒常常按正方形或六角
形排列。假定元件盒足够多，每个盒结构相同，那么
活性区中每个盒所占的栅胞的物理情况，可以代表整
个活性区中的状况。
进一步假定，元件盒是圆对称的，那么每个栅胞
中情况，可以用更小的单位（栅元）来反映。 比如对
六角形栅胞可取其 1/12 的 ΔOAB 来做代表； 正方形 栅
胞可用其 1/8 的 ΔO'A'B' 来做代表。这样一来问题就大
大简化 了。
现在的问题是怎样计算直角三角形栅元的物理量
（如通量）。用蒙特卡罗方法如何模拟中子在栅元内
的运动，反映出整个活性区对它的影响。
我们可把 OA',OB',A'B' 作为全反射边界来处理。
所谓全反射边界，就是当中子打到该边界上时，按镜
面反射的方式，从边界 上全部反射回来，中子的能量
与权重均不改变。
在这种边界上的反射条件，称之为全反射条件，
就是通常的镜面反射条件。
在全反射边界条件下，一条通过活性区若干个区
域的中子径迹，可以用栅元 ΔO'A'B' 中的一条折线轨
道来反映出来。
反过来，在直角三角形栅元 ΔO'A'B' 中任一条反
射成的折线轨道，都代表了中子在活性区内一条直线
轨道的作用。由于系统的对称性，在活性区内，凡是
与栅元内位置相当的地方，都有相同的物理情况，因
此栅元内各处的情况，当然代表了整个活性区的情况。
(2) 一般曲面全反射条件
对于一般曲面的全反射，设入射方向为 Ω，入射
点的内法线方向为 n,则反射方向 Ω' 为,
其中
设
则
? ?
nΩ
nnΩΩΩ
???
?????
?c o s2
2
nΩ ????c o s
),,( zyx nnn?n
)(c o s
c o s2
c o s2
c o s2
zyx
z
y
x
nwnvnu
nww
nvv
nuu
???????
????
????
????
?
?
?
?
(3) 平面全反射条件
设三角形栅元的横截面 ΔOAB 在 X-Y 平面上，
∠ OAB＝ θ。则边界 OA,OB,AB 上的反射都是平面
全反射。在任一与 X-Y 平面垂直且与 X 轴成 α角的平
面上，全反射条件为,
由此就可得到 OA,OB 和 AB 边上的全反射条件，对
于 OB 边,α=θ；对于 OA 边,α= 0；对于 AB 边，
α=π/ 2。
ww
vuv
vuu
??
?????
?????
??
??
2c o s2s in
2s in2c o s
(4) 反射层边界条件
对于具有大反射层的系统，如存放，运输和生产
裂变物质的仓 库、车厢和车间等，当中子从里面打到
四周墙上或反射层时，还要继续对它进行跟踪。这种
跟踪常常要花费很大的计算量，并且在结果中引起的
方差也比较大。如果在计算这种系统的不同方案中，
反 射层条件不变，那么这种大量重复的计算是很不经
济的。
中子射入反射层后，一部分被介质吸收，
只有一部分返回，由于中子的散射慢化，损失
一部分能量，因此反射回来的中子有一个能量
方向分布。显然，对这种反射层，不能应用全
反射条件。不过，我们仍然可以把它当做边界，
在边界上按反射层的物理作用来处理。
比如，如果反射层是一种平板几何，我们可以用
数值方法或蒙特卡罗方法，预先算好在各种不同入射
能量 E 下的反照率 β(E)，反射中子的能量分布
RE(E→ E' )。于是代替在反射层中眼踪中子，我们可
在反射层边界上作如下处理,
一旦中子打入反射 层，立即返回，反射后权重为
其中,E 为射入反射层中子的能量,W 为中子的权重。
反射后的能量 Eβ 由反射能谱 RE(E→ Eβ) 中抽样产生。
反射后的方向 Ωβ 由半平面各向同性分布或余弦分布中
抽样。反射后的中子位置为入射时的位置。
计算表明，对于大尺寸的反射层来说，这样的近
似，引 起的结果上的误差是可以忽略的，却能带来计
算量的大量节省。
)( EW ?? ??W
3,记录贡献与分析结果过程
在粒子输运问题中, 除了得到某些量的积
分结果外, 还需要得到这些量的方差, 协方差,
以及这些量的空间, 能量, 方向和时间的分布 。
这些量可以利用分类记录手续同时得到 。
1) 记录与结果
为了得到所求量的估计值，在粒子输运过程中需
进行记录，即求每个粒子对所求量的贡 献。
设模拟了 N 个粒子，所求量的估计值为,
其中 gi 为第 i 个粒子的总贡献。
?
?
?
N
i
iN gNg
1
1
记录的贡献由所求量决定。对于同一个所求量，
又随所用的蒙特卡罗技巧的不同而不同。 例如，所求
量是粒子穿透屏蔽概率，使用直接模拟法时，如粒子
穿透屏蔽，在叠加记录单元加,1” ( 初始值为零 )，否
则没有贡献。使用加权法时，如粒子穿透屏蔽，在叠
加记录单元加粒子的权重，否则没有贡献。使用统计
估计法时，粒子每发生一次碰撞 (包括零次碰撞 )，都
要记录贡献，等等。
2) 方差和协方差的估计
估计量 g 和 g' 的方差和协方差为,
它们可以用下式估计,
?
?
?
?
?
?
??
?
?
?
?
?
???
?
?
?
?
?
?
??
???
??
???
?
??
N
i
i
N
i
i
N
i
iigg
N
i
i
N
i
ig
g
N
g
N
gg
N
g
N
g
N
111
2
2
11
22
111
11
?
?
? ? ? ?
? ? ? ?? ?gEEgggE
EggE
gg
g
????
??
?
2
222
?
?
因此，要得到 和 的估计，只要对每一个历史记
录结果的 和 进行记录，并加以累加即可。
方差估计值 确定后，可得到误差
其中 为置信限，它随置信水平 而定。在通常
情况下，取 。
22 ggg ???
iii ggg ?2
2g?
N
g??? ??
??
96.1,95.01 ??? ???
??1
3) 位置、能量、方向、时间分布
在前面已经提到，用蒙特卡罗方法求某种
量的空间、能量、方向和时间分布，实质上是
得到这种分布的阶梯函数近似的估计值。而求
这种估计值是很方便的，只要将跟踪过程中所
得到的感兴趣量，按其状态的空间、能量、方
向、时间特征，分别记录其权 重，最后将这
些记录结果适当处理即可。
事先，将问题的空间、能量、方向（常按相对于
某个方向的夹角余弦）、时间范围，各分为如下不同
间隔,
再用一批存贮单元 {A} 记录相应间隔上阶梯函数
近似的累计值。;0;11;,,,
10
10
m a x01m i n
21
Tttt
EEEEE
VVV
L
K
J
I
?????
??????
?????
?
?
?
?
???
4,核截面数据的引用
用蒙特卡罗方法解粒子输运问题, 需要介质所包
含的各种原子核的核数据 。 以中子核数据为例, 需要
各种涉及到的核的微观总截面, 弹性散射, 非弹性散
射, n-2n 反应, 裂变, 俘获等截面；也需要这些反应
的相应能量, 角度分布, 次级粒子数, 以及其它关心
的粒子数及其能量, 方向分布 。 从输运方程中可以知
道, 有了这些数据, 问题就完全确定了；反映到蒙特
卡罗模拟中, 有了这些数据, 就能决定宏观总截面,
决定碰撞核的具体形式, 就能实现抽样和跟踪 。
在蒙特卡罗计算中, 引用的核数据有点截面, 分
段常数截面和群截面三种形式 。
1) 点截面形式
在跟踪粒子时，对粒子的每一种能量，先
从截面库中取出需要的核数据，再用插值（或
其它方式），求出相应能量的各种截面数据。
这种方法是直接的，也是比较精确的。不少通
用程序就是这样做的。这样做的条件是要有完
备的核截面数据库，计算机有大而快的存贮能
力。
2) 分段常数形式
将问题的能量范围（ Emin,Emax）分成许
多间隔，截面数据在每个间隔上看作与能量无
关，即截面取分段常数形式。这种引用截面数
据的好处是，数据量相对地少。问题是它要根
据问题的物理特点来划分能量间隔。而为了保
证误差较小，所取的间隔数一般是比较多的。
3) 群截面形式
解中子输运问题，常常采用多群近似方法。在多
群近似下，中子能量 E 用群指标 i 代替。为了实现蒙
特卡罗跟踪，只需要以下群截面,
显然，这种跟踪过程数据量小，程序简单。
群时的角分布；群散射到第点由第在
群的截面；群散射到第点由第在
群吸收截面；点的第在
群总截面；点的第在
jirr
jirr
irr
irr
:),|(
:)(
:)(
:)(
jif s
ji
s
i
a
i
t
?
?
?
?
?
?
5,蒙特卡罗程序结构
在电子计算机上, 蒙特卡罗方法解粒子输运问题
的程序, 一般都可分为：源抽样, 空间输运过程, 碰
撞过程, 记录过程和结果的处理与输出等部分 。
开始
数据预处理,各记录单元清零
取一个粒子历史
源分布抽样
输运过程
碰撞过程
历史终止否？
统计处理
做完给定历史数否？
结果的处理与输出
终止
记录过程
记录过程
记录过程
记录过程
至于粒子历史终止条件，根据问题的几何条件、
物理假定，处理方法，可归纳为以下几种,
1) 粒子从系统逃脱；
2) 粒子经碰撞被吸收；
3) 经俄国轮盘赌后，历史被终止；
4) 粒子能量低于给定能量（阈能） ;
5) 粒子位置越过某一界面；
6) 粒子飞行时间超过给定时间；
7) 粒子权重小于某个小量。
? 作 业
1)
