小波尺度函数展开连续能量共振计算加速方法研究

核动力工程Nuclear Power Engineering Vo l.31.S2 Dec. 2 0 1 0
第31卷 增刊2
2 0 1 0年12月
文章编号：0258-0926(2010)S2-0016-05
小波尺度函数展开连续能量共振计算
加速方法研究
祖铁军，吴宏春，曹良志，杨伟炎，刘国明
（西安交通大学核科学与技术学院，西安，710049）
摘要：小波尺度函数展开连续能量共振计算方法可以处理复杂的共振计算问题。

该方法利用Daubechies 小波尺度函数对共振能量段中子注量率的能量变量进行离散，将中子注量率的求解转化为求解一系列的展开系数的方程组；由于要求解大量的展开系数，该方法的计算效率很低。

本研究采用离散方向概率方法求解展开系数。

离散方向概率方法结合了特征线方法和流耦合方法的优点，具有较高的计算精度和效率。

通过对压水堆栅元问题的数值检验，证明该方法能大大提高计算效率。

关键词：连续能量共振；Daubechies小波尺度函数；离散方向概率方法
中图分类号：TL323文献标识码：A
1 引言
传统的共振计算方法难以用于任意复杂几何多核素共振计算；本课题组提出的小波尺度函数展开连续能量共振计算方法[1]在快中子区和热中子区采用传统的多群方法，在共振区用Daubechies[2]小波尺度函数离散中子注量率的能量变量，将中子注量率的求解转化为展开系数的求解，重构中子注量率。

现有的WA VRESON[1]程序采用特征线（MOC）[3]方法求解小波尺度函数展开系数，在整个计算区域内进行特征线扫描，计算效率极低，限制了小波尺度函数展开连续能量共振计算方法在复杂问题中的应用。

本文采用离散方向概率方法[4]求解小波尺度函数展开系数。

首先将整个计算区域分解为若干个小的子区域，然后在子区域内进行特征线扫描，各子区域之间通过表面的角注量率连续进行耦合，在保持原有精度的同时提高计算效率。

2理论模型
2.1 小波尺度函数展开连续能量共振计算理论
模型
小波尺度函数由基函数通过平移和伸缩操作得到，具有良好的紧支性和正交性[1,2]，可以用来逼近剧烈变化的函数。

本文用小波尺度函数对共振区中子注量率能量变量进行离散：
()()()
∑
=
n
n i
n i
E
a
E,
,
,
,
,ϕ
φΩ
r
Ω
r（1）式中，()Ω
r,
,E
φ为中子注量率，cm-2·s-1；r为位置矢量，cm；Ω为方向矢量；E为能量，eV；
()Ω
r,
,n
i
a为小波尺度函数展开系数；()E n i,ϕ为Daubechie 小波尺度函数；小波尺度函数伸缩阶数2
2
2
,2,1−
+
=N
n i
"；N为Daubechie 小波阶数；i为展开阶数。

将中子输运方程中的中子注量率用小波尺度函数展开，并利用小波尺度函数的正交性可得[1]：
()()
n
n
n i
n i
A
a
a,
,
,
,
,Ω
r
Ω
r
Ω+
∇⋅
()
∑
≠
−
+
=
n
l
n
l
n
i
l
l
A
a
F
S,
,
,Ω
r（2）
2
2
2
,2,1−
+
=N
l i
"（3）
()()()
l n
g
E
n
i
l i
n
l
A
E
E
E
E
Σ
A,
,
,
t
,
d
,=
=∫
∆
ϕ
ϕ
r（4）
()×
=∫
∆g
E
l i
l
E
S,ϕ
()
∫∫∞→
→
π4
'
,
'
,
,s
d'
d'
d'
,'
,
ΣE
E
E
E
E
Ω
Ω
r
Ω
Ω
r
φ（5）
收稿日期：2010-09-18；修回日期：2010-10-12
基金项目：国家自然科学基金项目（10875094）；863项目（2009AA050705）
祖铁军等：小波尺度函数展开连续能量共振计算加速方法研究 17
()()∫∆×=
g
E
l i l E E k F χϕ,π41
()()()E E E E ΣE d 'd 'd ',',','0π
40
f
ΩΩr r ∫∫∞φν （6）
式中，l 为小波尺度函数伸缩阶数；),(t E Σr 为宏观总截面，cm -1；ΩΩr →′→′,, ,s E E Σ为散射算子，cm -
1；
)(E χ为裂变中子能量份额；)(E ν为每次裂变产
生中子数；),(f E Σr 为宏观裂变截面，cm -
1；g
E ∆为第g 群能群宽度；n l A ,为总截面；l S 为散射源项；l
F 为裂变源项。

这些参数可以通过多群数据库和连续能量数据库数据计算得出，因此，在计算展开系数时为常数。

以上参数的计算可参见文献[1]。

方程（2）的二维积分形式表示为： ()
=xy xy n i Ωr a ,,
(
)
+⎟⎟⎠
⎞
⎜⎜⎝⎛′−×⎥⎥⎦⎤⎢⎢⎣⎡−+∫
∑≠S
n n
l n n l xy xy n i l l s A A Ωr a F S 0
,,,d sin exp ,θ
(
)
⎟
⎟⎠
⎞
⎜⎜⎝
⎛′′−
×′′∫
S
n n
xy xy
n i s A Ωr a 0,,d sin 1
exp ,θ （7） 式中，r xy 为位置矢量在x-y 平面上的投影，cm ；xy Ω为方向矢量Ω在x-y 平面上的投影；s 为xy 方向
中子径迹位置坐标在x-y 平面上的投影，cm ；S 为xy 方向中子入射位置到出射位置间距在x-y 平
面上的投影，cm ；n n
A ,′表示xy r ′处的总截面，n n A ,′′表示系统边界处的总截面，其中xy xy xy s Ωr r −=′，xy xy xy
s Ωr r −=′′；θ是方向xy Ω的极角。

2.2 离散方向概率方法的理论模型
离散方向概率方法（DDPM ）结合了穿透概
率方法和MOC 方法的优点。

该方法将计算系统分成若干子系统，再将子系统细分为许多节块。

各节块以及表面之间通过方向穿透概率和方向碰撞概率耦合。

该方法中主要有3个近似：①中子注量率的角度变量采用和离散纵标方法相同的方法进行离散；②子系统表面中子注量率在空间上均匀分布；③在节块内中子注量率在空间上均匀分布。

离散方向后二维中子积分输运方程可表示为如下形式：
()
()+⎟⎟⎠⎞⎜⎜⎝⎛′−×′=∫
s
n m
g m g xy n m g s ΣQ r 0
,t ,,d sin exp θφr r ()()⎥⎥⎦
⎤⎢⎢⎣⎡′′−×′′∫
s m
g n m n m g s Σ0,t ,,d sin 1
exp r r θφ （8） 式中，()
xy n m g r ,,φ为r xy 处方向),(n m 上第g 群中子
角注量；g Σ,t 为第g 群总截面，cm -
1；m 为方位
角ϕ 离散标号；n 为极角θ 离散标号；g Q 为第g 群源项，cm －
3·s －
1，源项的计算可参见文献[4]。

可以看出，方程（7）和方程（8）形式相同，因此只要将方程（8）中的参数替换为方程（7）中对应位置上的参数，即求得展开系数。

本文采用离散方向概率方法求解展开系数。

2.2.1 方向穿透概率 如图1所示，在子系统中，出射边界的注量率有2个来源：①来自入射边界
的中子注量率；②来自某方向m xy ,Ω特征线所经过的节块中的中子注量率。

图1 节块和边界的耦合关系
Fig. 1 Coupling of Meshes and Boundaries
in a Subsystem
沿特征线计算方程（8），可以得到子系统m xy ,Ω方向的出射中子注量率：
+=
∑=→LI
l n
m g l l l n m g n
m g l T
1
i in
,,,i o i ,,,out
,,,o φφ
∑=→I
i g
i l i n m g Q T
1
,,,,o （9）
式中，in ,,,i
n m g l φ为入射中子角注量率，cm -2·s -
1
；g i Q ,为节块i 内中子源项，cm －
3·s －
1；o i ,,,l l n m g T →为从
入射边l i 到出射边l o 的方向穿透概率；o ,,,l i n m g T →为从节块i 到出射边l o 的方向穿透概率。

方向穿透概率o i ,,,l l n m g T →的计算公式可以表示为：
核 动 力 工 程 V ol.31. S2. 2010
18 ∑∑
∈→⎥⎥⎦
⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛−=
i i
,,,,,t o
,)
(,o
i ,,,sin exp l l NI
n n
l ni n m ni g l m l b m l l n m g S Σw w T θ（10）
式中，)(,l b m w 为相邻特征线间宽度，cm ；o ,l m w 为出射边在m xy ,Ω方向总宽度，cm ；NI 为m xy ,Ω方向对出射边中子注量率有贡献的节块数；l ni n m S ,,,为特征线l 在节块ni 中的截距，cm 。

方向穿透概率o ,,,l i n m g T →的计算公式可以表示为：
=→o ,.,l i n m g T
∑
∈⎪⎩⎪⎨⎧×⎥⎥⎦
⎤
⎢⎢⎣
⎡⎟⎟⎠⎞⎜⎜⎝⎛−−i
l n l
i n m i g l m i g l b m S Σw Σw θsin exp 1,,,,,t o
,,,t )
(, ⎪⎭
⎪⎬⎫
⎟⎟⎠⎞⎜⎜⎝⎛−
∑
+=MI i mi n l mi n m mi g S Σ1,,,,,t sin exp θ （11） 式中，mi 是从节块i 开始对出射边中子注量率有贡献的节块数。

2.2.2 方向碰撞概率 在方程（9）中为了使用源迭代方法求解未知量，需要给出节块内中子角注量率的表达式。

在一个节块内满足中子守恒关系：
()
out ,,,in
,,,,,,,,,t sin n m g i n m g i i
n g i n m g i i g A Q Σφφθφ−+
= （12） 式中，i g Σ,,t 为节块i 内第g 群总截面，cm -
1；i
A 为节块面积，cm 2；in ,,,n m g i φ为m xy Ω,方向节块i 的入射中子注量率，cm -2·s -
1；out ,,,n m g i φ为m xy Ω,方向
节块i 的出射中子注量率，cm -2·s -
1。

节块i
的入射中子注量率in ,,,n m g i φ（即
i -1的出
射中子注量率）和出射中子注量率out ,,,n m g i φ可以采用和方程（9）相同的方法求得，不再赘述，代入
方程（12）后可以得到节块内中子角注量率的表达式：
∑∑
=→=→+
=
J
j g
j i j n m g l i n
m g l i l m g n m g i Q P
P 1
,,,,i
1
in
,,,i i ,,,,,φφ（13）
式中，i l n m g P →i ,,,为m xy ,Ω方向由入射边l i 到节块i 的碰撞概率；i j n m g P →,,,为m xy ,Ω方向由节块j 到节块i 的碰撞概率。

方向碰撞概率i l n m g P →i ,,,可以表示为：
×=
→i
i g n
i l n m g A ΣP ,,t i ,,,sin θ
×⎟⎟⎠
⎞⎜⎜⎝
⎛−
∑
∑
∈−=i
11,,,,,t )(,sin exp l l i nj n l
nj n m nj g l b m S Σw θ ⎥⎥⎦
⎤
⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛−−n l
i n m i g S Σθ
sin exp 1,,,,,t （14） 1≠j 时，方向碰撞概率i j n m g P →,,,可以表示
为：
×=
∑∈→j
l l b m i
j g i g n
i j n m g w
A ΣΣP )
(,,,t ,,t ,,,sin θ
×⎥⎥⎦⎤
⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛−−n l j n m j g S Σθsin exp 1 ,,,,,t ⎢⎢⎣⎡−⎟⎟⎠
⎞⎜
⎜⎝⎛−
∑
−+=11
,,,,,t sin exp i j mj n l mj n m mj g S Σθ ⎥⎥⎦
⎤
⎟⎟⎠⎞⎜⎜⎝
⎛−
∑
+=i j mj n l
mj n m mj g S Σ1,,,,,t sin exp θ （15） 在节块i 内的自碰撞概率为：
×=→i
i g i g n
i i n m g A ΣΣP ,,t ,,t ,,,sin θ
i g i l n l i n m i g l b m ΣS Σw ,,t ,,,,,t )(,1
1sin exp +⎥⎥⎦⎤⎢⎢⎣
⎡−⎟⎟⎠⎞⎜⎜⎝⎛−∑
∈θ （16） 另外，在这些概率之间存在互易关系和守恒关系，利用这些关系可以减少计算量。

3 数值结果
基于以上原理，开发了共振计算程序WRDDPM 。

为了验证基于DDPM 求解小波尺度函数展开系数的高效性和精确性，本文选择了PWR 栅元问题进行校验。

PWR 栅元的几何结构如图2所示。

在几何处理上，采用了和WA VRESON 相同的处理方式，将栅元分成7个节块，认为各节块内材料均匀。

0~π/2内选择4个均匀离散的幅角，对于极角，采用最优二极角。

特征线间距为 0.01 cm ；235U 的富集度3%；慢化剂为轻水；包壳为Zr 。

栅元中各核素的原子密度列于表1。

表2中列出了WRDDPM 和WA VRESON 计算得到的燃料区共振能群中子注量率，以及和
祖铁军等：小波尺度函数展开连续能量共振计算加速方法研究
19
图2 PWR栅元的几何描述cm
Fig. 2 Geometry Description of PWR Fuel Cell (cm) 表1 PWR栅元中个核素的原子核密度
Table 1 Atomic Density of Different Nuclides
in PWR Fuel Cell
原子核密度/1024 ·cm-3
材 料 235U 238U 16O 1H Zr 燃 料 0.0007014 0.0223920.0461868 ——
包 壳 ——— —0.04360
慢化剂 — — 0.03362 0.06723—
表2 燃料区共振能群平均中子注量率及相对误差①Table 2 Average Flux and Relative Error of Resonance Energy Group in Fuel Area
共振 能群 MCNP WA VRESON WRDDPM 误差1
/%
误差2
/%
15 0.281
2130.279 4710.279 4750.619 0.618
16 0.247
9190.249 6210.249 6290.687 0.690
17 0.241
5900.242 6630.242 6600.444 0.443
18 0.239
8260.240 7780.240 7740.397 0.395
19 0.235
7970.235 7570.235 7480.017 0.021
20 0.460
6780.45 9380.459 3660.282 0.285
21 0.444
6890.443 8780.443 8660.182 0.185
22 0.320
4340.319 9190.319 8980.161 0.167
23 0.207
6760.206 9540.206 9400.348 0.354
24 0.234
4520.233 3860.233 3650.455 0.463
25 0.218
8380.218 0210.218 0010.373 0.383
26 0.208
8960.20 7080.207 0710.869 0.873
27 0.327
4280.323 0720.323 051 1.330 1.337注：①误差1为WA VRESON计算中子注量率的相对误差；误差2为WRDDPM计算中子注量率的相对误差
MCNP计算结果的相对误差。

其中，用MCNP计算时，每代投入20000个粒子，共计算300代，其中前50代不参加结果统计。

从表中可以看出2种算法的结果非常接近，WA VDDPM与MCNP 的计算结果最大相对误差为1.337%。

表3中给出
表3 无限介质增殖系数（k inf）和计算时间
Table 3 k inf and Calculation Time
程序名称k inf计算时间/s
MCNP 1.38
496
1042 WAVRESON 1.383
085 6149 WRDDPM 1.382
977 51 误差1/%①0.135 —
误差2/%①0.143 —
注：①误差1为WA VRESON计算中子注量率的相对误差；
误差2为WRDDPM计算中子注量率的相对误差
了由2种方法求得的k inf，其中，WRDDPM的相
对误差为0.143%，而计算速度提高120倍。

由以上数值结果可知，WRDDPM和
WA VRESON的计算精度基本相同，但计算效率
明显提高。

由于DDPM和MOC相同，都是沿特
征线进行计算，因此可以得到相同的计算精度，
但是，DDPM与MOC的不同之处在于并非在整
个区域内进行扫描，而节块之间通过方向概率进
行耦合，从而大大提高计算效率。

4 结束语
采用离散方向概率方法求解共振区Daubechies小波尺度函数展开系数，并编制出相
应的计算程序WRDDPM，该程序可以达到现有Daubechies小波尺度函数展开连续能量共振计算
程序相同的计算精度，而计算效率明显提高。

从
而提高了Daubechies小波尺度函数展开连续能量
共振计算方法处理实际复杂工程问题的能力。

参考文献：
[1]杨伟炎，吴宏春，曹良志，等. 小波尺度函数展开在
连续能量中子输运方程计算中的应用[J]. 核动力工程,
2009, 30(5): 18-22.
[2]CHO Nam Zin, Park Chang Je. Wavelet Theory for
Solution of the Neutron Diffusion Equation [J]. Nuclear
Science and Engineering, 1996, 124: 417-430.
[3] CHEN Qi-chang, WU Hong-chun, CAO Liang-zhi. Auto
MOC--A 2D Neutron Transport Code for Arbitrary
Geometry Based on the Method of Characteristics and
Customization of AutoCAD[J]. Nuclear Engineering and
Design, 2008, 238: 2828-2833.
[4]刘国明. 非结构几何下中子积分输运方程的穿透概率
方法与离散方向概率方法研究[D]. 西安：西安交通大
学博士论文，2009.
核动力工程 V ol.31. S2. 2010 20
Study on Acceleration Method for Wavelets Scaling Function
Expansion Resonance Calculation
ZU Tie-jun, WU Hong-chun, CAO Liang-zhi, YANG Wei-yan, LIU Guo-ming
(School of Nuclear Science and Technology, Xi’an Jiaotong University, Xi’an, 710049, China)
Abstract: The continuous-energy resonance self-shielding calculation method based on the wavelets scaling function expansion method is a valuable potential method to solve the complex resonance problem. Daubechies’ wavelets scaling function expansion method is used to discretize the energy variable of neutron flux within the resonant energy range. In this method, the neutron transport equation is transformed to a set of coefficients equations of wavelets scaling function expansion. Calculation of the coefficients is so time-consuming that a powerful neutron transport calculation method is needed for better calculation efficiency. In this paper, the discrete direction probability method (DDPM) is employed as a tool for solving the wavelet scaling function expansion coefficients. The numerical results indicate that it is capable of computing the coefficients with high computational efficiency and precision.
Key words:Continuous-energy resonance self-shielding calculation, Daubechies’ wavelets scaling function, Discrete direction probability method
作者简介：
祖铁军（1985—），男，硕士研究生。

2009年毕业于西安交通大学核能系，获学士学位。

现从事核反应堆物理计算研究。

吴宏春（1964—），男，教授，博士生导师。

1994年毕业于西安交通大学核能系，获博士学位。

现从事核反应堆物理计算研究。

曹良志（1978—），男，副教授。

2005年毕业于于西安交通大学核能系，获博士学位。

现主要从事核反应堆物理计算研究。

（责任编辑：刘君）（上接第8页）
Subgroup Method for Resonance Self-Shielding Calculation and Research on Resonance Interference Effect
HUANG Shi-en1,2, WANG Kan1, YAO Dong2
(1. Department of Engineering Physics, Tsinghua University, Beijing, 100084, China;
2. Science and Technology on Reactor System Design Technology Laboratory, Chengdu, 610041, China)
Abstract: We combined the subgroup method and the method of characteristics for neutron resonance self-shielding calculation, and programmed a resonance calculation code SGMOC. This code is based on the WIMSD format multi-group data library. The numerical results of SGMOC are in good agreement with those of MCNP, and thus SGMOC has high calculation precision and geometry flexibility. Based on the SGMOC code, we studied two modification methods for the resonance interference effect calculation. The conditional probability method had a correction effect about 30~180pcm for k inf calculation of UO2 fuel cell. The method by virtue of NJOY code had a correction effect about 20~130pcm for k inf calculation of UO2 fuel cell.
Key words: Subgroup method, The method of characteristics, SGMOC, Resonance interference effect, Conditional probability
作者简介：
黄世恩（1984—），男，博士研究生。

2006年毕业于清华大学工程物理系，获学士学位。

现主要从事反应堆物理方向的科研与工程设计工作。

王侃（1965—），男，教授，博士生导师。

1993年毕业于西安交通大学能源与动力工程系反应堆工程与反应堆安全专业，获博士学位。

现主要从事核能科学与工程方向的科研与教学工作。

姚栋（1949—），男，研究员级高级工程师，博士生导师。

1976年毕业于北京大学物理系，现主要从事反应堆堆芯核设计及计算软件开发工作。

（责任编辑：刘 君）。