Sn 法に用いる角度求積セットの数値積分精度 ∗千葉豪
平成 29 年 1 月 26 日

はじめに

1

Sn 法に用いる角度求積セットが果たすべき要件の一つに、中性子束の角度モーメント（用語として正確？）を可能な限り正確に計算することが挙げられる。中性子束の角度モーメントの計算が必要となる局面は、大きく以下の２つに分けられるであろう。
• 角度中性子束に依存する物理量、例えば反応率（0 次モーメント）や、平面からの中性子漏洩量（1 次モーメント）を評価したい場合。
• 角度微分散乱断面積が角度について関数形式で与えられている場合（散乱中性子源の評価のために必要となる）。
なお、Sn 法に用いる角度求積セットに関しては、Alcouﬀe のテキスト [1] が詳しく、かつ分かりやすいのでお薦めである。

一次元の場合2一次元体系における角度中性子束を ψ(µ) とし、これが以下のように L 次までの Legendre 多項式で記述できるものとする。
L∑
2l + 1
ψ(µ) =ϕl Pl (µ)(1)
2l=0
ここで Pl (µ) は l 次の Legendre 多項式を示す。展開係数 ϕl は以下のように定義される。
∫ϕl =ψ(µ)Pl (µ)dµ(2)Sn 法では、ϕl は以下のように求積で求められる。
ϕ̄l =J∑wj ψ(µj )Pl (µj )j=1

ここで、(µj , ωj ) は J 点の求積セットにおける j 番目の離散角度と重みを示す。
次に、ϕ̄l の精度について考えよう。

∗/Document/Study/SnQuadrature/1(3)式 (2) と式 (1) から以下の式が得られる。

ϕl =

L { ′
∑
2l + 1l′ =02}

∫ϕl ′Pl′ (µ)Pl (µ)dµ =L { ′∑
2l + 1l′ =022ϕl ′ ·δll′2l + 1}

(4)

また、式 (3) と式 (1) から以下の式が得られる。

ϕ̄l =L  ′∑2l + 1l′ =0

2ϕl ′J∑wj Pl′ (µj )Pl (µj )j=1(5)

従って、ϕ̄l の計算精度は、Legendre 多項式の直交性の求積による再現精度に置き換わることが分かる。
なお、Legendre 多項式の積はより高次の Legendre 多項式で展開することができ、例えば Pl と′Pl の積は (l + l′ ) 次以下の Legendre 多項式で展開できる。従って、Legendre 多項式の直交性の求積による再現精度は、以下の関係の求積による再現精度にさらに置き換えられることになる。
∫Pl (µ)dµ = 2δl0(6)一般的に一次元体系を対象とした Sn 法では、求積セットとして Gauss-Legendre セットが利用される。J 点の Gauss-Legendre セットを用いることで、(2J − 1) 次までの多項式に対する積分を厳密に求められることが分かっている 1 。Legendre 多項式 Pl (µ) の次数は l であることから、Gauss-Legendre 求積セットから計算される展開係数 ϕ̄l は、2(l + L) − 1 ≤ J が成立する限り、ϕlを厳密に再現することが言える。
以上について数値的に確認するため、実際に Gauss-Legendre 求積セットを用いて、Legendre 多項式の直交性に関する計算精度を評価した。求積セットを用いて計算したものは以下である。求積が正確に出来ていればゼロとなるべきものである。
√√2l + 12J2l′ + 1 ∑ωj Pl (µj )Pl′ (µj ) − δll′
2 j=1(7)角度分点が異なるいくつかの Gauss-Legendre 求積セットを用いて計算した結果を Fig. 1 に示す。
縦軸と横軸の数値がそれぞれ l、l′ に対応する。

1

これについては別稿「Gauss-Legendre 求積公式の導出」を参照のこと。

2(a) S4(b) S6990’err_s4’

100’err_s6’

10887710-2

610-2

610-4

510-4

54

410-6

310-632210-8

10010-100

1234

56

7810-81910-100123(c) S8456789(d) S10

990’err_s8’

100’err_s10’

1088

7710-2

610-2

610-4

510-4

54

410-6

310-6

32210-81010-810-10-10100

123456781090123456789Fig. 1: Reproduction error of orthgonality of the Legendre polynomials by the Gauss-Legendrequadratures with diﬀerent order3多次元の場合

3多次元の場合、一般的に、以下のように中性子束は角度について球面調和関数により展開される。

⃗ =
ψ(Ω)

lL
∑
2l + 1 ∑
l=04π

⃗ lm
Rlm (Ω)ϕ

(8)

m=−l

展開係数 ϕlm は以下で定義される。
∫ϕlm =4π⃗ lm (Ω)ψ(
⃗⃗
dΩRΩ)

(9)

一方、求積により計算される展開係数 ϕ̄lm は以下のように書ける。

ϕ̄lm =J∑⃗ j )ψ(Ω⃗ j)
wj Rlm (Ω(10)j=1

一次元の場合と同様に考えると、ϕ̄lm の計算精度は、以下の球面調和関数の直交性の求積による再現精度に置き換わることが分かる。
∫4π⃗ lm (Ω)R⃗ l′ m′ (Ω)⃗ =dΩR

4πδll′ δmm′2l + 1(11)なお、球面調和関数の積はより高次の球面調和関数で展開することができ、例えば Rlm と Rl′ m′の積は (l + l′ ) 次以下の球面調和関数で展開されることが分かっている [2]。従って、球面調和関数の直交性の求積による再現精度は、以下の関係の求積による再現精度にさらに置き換えられることになる 2 。
∫⃗ lm (Ω)⃗ = 4πδl0 δm0dΩR(12)4π一次元の場合と同様に、求積セットによる球面調和関数の直交性に関する計算精度を評価した。
求積セットを用いて計算したものは以下である。求積が正確に出来ていればゼロとなるべきものである。
√√
J2l + 1 2l′ + 1 ∑⃗ j )Rl′ m′ (Ω⃗ j ) − δll′ δmm′
ωj Rlm (Ω(13)4π
4π j=12

例えば、文献 [3] の p.124 には以下の記述がある。
「the directions and weights of the angular quadrature should
be built so as to exactly integrate a maximum number of products of two spherical harmonics or, equivalently, tointegrate a maximum number of spherical harmonics」。

4

はじめに、Level-symmetric（LS）求積セットを用いて得た結果を Fig. 2 に示す。なお、X 軸（Y 軸）の数値は、式 (13) における l、m と Table 1 で示す対応をとる。

(a) LS4(b) LS83535100’err_ls4’100’err_ls8’

30

30
10-2

10-225

25
10-4

20
15

10-4

20
15

10-6

1010-6

10
10-8

10-8

5

510-10
51015

2025

30

10-1035510(c) LS121520253035(d) LS16

35

35

100

’err_ls12’

100’err_ls16’

30

30
10-2

10-2

25

25
10-4

20
15

10-4

20
15

10-6

1010-61010-810-85510-10
51015

20

25

30

10-10

35

51015

20

25

30

35

Fig. 2: Reproduction error of orthgonality of the spherical harmonics functions by the levelsymmetric quadratures with diﬀerent order

Table 1: Correspondence of index to the order of the spherical harmonics functions
Index
1234
5
6
7
8
9(l, m)(0,0)(1,-1)(1,0)(1,1)(2,-2)(2,-1)
(2,0)
(2,1)
(2,2)

Index
10
11
12
13
14
15
16
17
18

(l, m)(3,-3)(3,-2)(3,-1)
(3,0)
(3,1)
(3,2)
(3,3)
(4,-4)(4,-3)

Index
19
20
21
22
23
24
25
26
27

5(l, m)(4,-2)(4,-1)
(4,0)
(4,1)
(4,2)
(4,3)
(4,4)
(5,-5)(5,-4)

Index
28
2930
3132
33343536(l, m)(5,-3)(5,-2)(5,-1)
(5,0)
(5,1)
(5,2)
(5,3)
(5,4)
(5,5)

次に、Even-Odd（EO）求積セット [4] を用いて得た結果を Fig. 3 に示す。Level-Symmetric
セットと比較して直交性の再現性度が悪いが、これは Even-Odd 求積セットが、異なる媒質間で中性子束の角度依存性が滑らかとならない場合でも精度を保つ目的で、オクタントでのモーメント保存を条件として課していることに由来する。

(a) EO4(b) EO835

350’err_eo4’

100’err_eo8’

1030

30
10-2

10-22525
10-42015

10-42015

10-61010-610
10-810-8

5510-10510152025

30

10-10

35

510(c) EO121520253035(d) EO1635350’err_eo12’

100’err_eo16’

1030

30
10-210-2252510-4201510-4201510-61010-61010-810-85510-1051015202530

10-1035

5101520253035Fig. 3: Reproduction error of orthgonality of the spherical harmonics functions by the even-oddquadratures with diﬀerent order

6以上は、完全対称（Fully-symmetric）な求積セットによる計算結果であったが、完全対称性を満足しない（が、平面対称（Planar-symmetric）な）求積セットである Product quadrature（PQ）と Triangular quadrature（TQ）[5] についても同様の計算を行ったのでその結果を示す。PQ、TQ
におけるオクタント上での離散点の配置例を Fig. 4 に示す。

Fig. 4: Brief descriptions of product quadrature (left) and triangular quadrature (right)

また、各求積セットの離散点数について Table 2 にまとめる。

Table 2: Discrete points of various angular quadrature sets
Quadrature
LS4, EO4LS8, EO8
LS12, EO12
LS16, EO16
P2 T2 -PQ, DP2 T2 -PQ
P3 T3 -PQ, DP3 T3 -PQ
P5 T4 -PQ, DP5 T4 -PQ
P6 T6 -PQ, DP6 T6 -PQ
P2 T2 -TQ
P4 T4 -TQ
P6 T6 -TQ
P8 T8 -TQ

Number of discrete points per octant
310213649203631021367Gauss-Legendre と Chebyshev の PQ(PnTn-PQ) による計算結果を Fig. 5 に示す。直交性の再現性は Level symmetric セットより良好である。なお、PnTn-PQ は、「PQ による球面調和関数の求積」という観点からは最良のものとなる [5]。

(a) P2T2-PQ(b) P3T3-PQ

35

35

100’err_p2t2’

100’err_p3t3’

30

30
10-2

10-22525
10-42015

10-4201510-61010-61010-810-8

5510-1051015202530

10-1035510(c) P5T4-PQ

1520253035(d) P6T6-PQ3535100’err_p5t4’

100’err_p6t6’

3030
10-210-2252510-4201510-4201510-61010-61010-810-85510-1051015202530

10-1035

5101520253035

Fig. 5: Reproduction error of orthgonality of the spherical harmonics functions by the PnTn-PQ
with diﬀerent order

8次に、Gauss-double Legendre と Chebyshev の PQ（DPnTn-PQ）による計算結果を Fig. 6 に示す。PnTn-PQ と比較して再現精度が劣るが、これは DPn では極角方向について２分割して扱っていることに由来する（LS と EO の関係と類似と考えて良い）。

(a) DP2T2-PQ(b) DP3T3-PQ

35

35

0’err_dp2t2’

100’err_dp3t3’

1030

30
10-2

10-225

25
10-4

20
15

10-4

20
15

10-6

10

10-6

1010-810-85510-10
51015

2025

30

10-10

35

510(c) DP5T4-PQ

15

2025

30

35(d) DP6T6-PQ

35

35

100’err_dp5t4’

100’err_dp6t6’

30

30
10-2

10-2

25

25
10-4

20
15

10-4

20
15

10-6

10

10-6

10
10-8

10-8

5

510-10
51015

20

25

30

10-10

35

51015

20253035

Fig. 6: Reproduction error of orthgonality of the spherical harmonics functions by the DPnTnPQ quadratures with diﬀerent order

9さらに、Gauss-Legendre と Chebyshev の TQ（PnTn-TQ）による計算結果を Fig. 7 に示す。
この求積セットでは方位角方向の積分精度が悪化するため、PQ と比較して直交性の再現精度が極めて悪いことが分かる。

(a) P2T2-TQ(b) P4T4-TQ

35

35

0’err_tri_p2’

100’err_tri_p4’

1030

30
10-210-225

25
10-42015

10-42015

10-61010-61010-810-85510-10
51015

2025

30

10-1035510(c) P6T6-TQ

1520253035(d) P8T8-TQ

35

35

100’err_tri_p6’100’err_tri_p8’

3030
10-210-22525
10-4201510-4201510-61010-61010-810-85510-1051015202530

10-1035

5101520253035

Fig. 7: Reproduction error of orthgonality of the spherical harmonics functions by the PnTn-TQquadratures with diﬀerent order10なお、球面調和関数の求積という観点では、Lebedev 型 [6] のものが最良のものという記述が文献 [3] にある。以下、この文献の p.124 の記述を引用する。
「(Lebedev-type quadratures) oﬀers thebest performance ratio, i.e., the minimal number of directions for a given number of sphericalharmonics that are exactly integrated.」さらに、「the nodes of these quadratures are nearlyuniformly distributed」ともされており、Levedev 型は Ray-eﬀect に対しても有効な求積セットであると予想される。
最後に、Lebedev 型の求積セットに関する計算結果を Fig. 8 に示す。
「Precision-7」「Precision15」の全球面に対する離散点数はそれぞれ 26、86 である。

(a) Precision 7(b) Precision 15

35035-230

10-22510-42510-420-620

10-615-810

15

10-81010-101010-10

510-1210’err_leb_26’
30

1010-12105100’err_leb_86’-14

10-141051015202530

355101520253035Fig. 8: Reproduction error of orthgonality of the spherical harmonics functions by the rectangular DPnxn quadratures with diﬀerent order114

座標系を回転させたときの角度求積セットの数値積分精度⃗ の全球面の積分について、角度求積セットが l ≤ L 次まで厳密に再現でき球面調和関数 Rlm (Ω)
るとしよう。このとき、この求積セット、もしくは被積分関数である球面調和関数の座標軸をどの⃗ の求積は厳密に行われる。
ように回転させても、l ≤ L を満たす Rlm (Ω)
これについては、以下のように説明できよう。
L 次以下の球面調和関数で展開可能な球面上で定義される関数についてその座標系を回転させた場合、新しい座標系においてもやはり L 次以下の球面調和関数で展開可能と言える。従って、⃗ についてその座標系を回転させて定義される関数は、L 次以下の球面調和関数で展開可能
Rlm (Ω)⃗ の全球面の積分が角度求積セットで再現できるとすと言える。これより、l ≥ L 次までの Rlm (Ω)⃗ の座標系を回転させた関数の積分についても同様に再現できる
るならば、l ≤ L を満たす Rlm (Ω)
ことになる。
なお、「L 次以下の球面調和関数で展開可能な球面上で定義される関数についてその座標系を回転させた場合、新しい座標系においてもやはり L 次以下の球面調和関数で展開可能と言える」という点について、明確な説明が遠藤によって与えられている [2] ので、以下に引用する。
⃗ について、任意の角度で回転させた後に積分することを考える。そ『被積分関数である Rlm (Ω)⃗ に回転演算子を作用させて Ω⃗ ′ に回転させた後、Rlm (Ω⃗ ′ ) の数値積分を角度の操作は、立体角 Ω⃗ に回転演算子を作用させるという求積セットを使用して計算することになる。ここで、立体角 Ω⃗ = (Ωx , Ωy , Ωz ) のベクトルに回転行列を作用させることになるが、そのような変換をことは、Ω⃗ ′ の成分は Ωx 、Ωy 、Ωz の線形結合で表されることになる。L 次の球面調和関数行った後の角度 Ω⃗ は Ωx 、Ωy 、Ωz の L 次多項式で記述できることから、Ω⃗ ′ の各成分が Ωx 、Ωy 、Ωz の線形
Rlm (Ω)⃗ ′ ) は Ωx 、Ωy 、Ωz に関する L 次多項式で記述できる、つまり、L 次結合で表せるならば、Rlm (Ωまでの球面調和関数で展開できることになる。』

12

5

おわりに

本稿ではいくつかの角度求積セットについての計算結果を示したが、これら以外にも多くのセットが存在する。例えば、米国アルゴンヌ国立研究所の高速炉解析コードシステム SHARP[7] の離散座標法に基づく輸送ソルバー SN2ND には、本稿で採り上げたものに加えて Thurgood 型 [8]、COBE sky cube 型 [9]、Tegmark icosahedron 型 [10] といったものが実装されている。
本件についても、名古屋大学の遠藤知弘氏に様々なことをご教示いただいた。ここに深い謝意を表する。

参考文献[1] R.E. Alcouﬀe, R.D. O’Dell, “Transport calculations for nuclear reactors,” CRC Handbookof Nuclear Reactors Calculations, Vol.1, p.380, CRC Press, Boca Raton, Florida (1986).
[2] 遠藤知弘、私信.
[3] R. Sanchez, “Prospects in deterministic three-dimensional whole-core transport calculations,” Nucl. Eng. Technol., 44, p. 113-150 (2012).
[4] T. Endo, A. Yamamoto, “Development of new solid angle quadrature sets to satisfy evenand odd-moment conditions,” J. Nucl. Sci. Technol., 44, p. 1249-1258 (2007).
[5] R. Sanchez, J. Ragusa, “On the construction of Galerkin angular quadratures,” Nucl. Sci.
Eng., 169, p. 133-154 (2011).
[6] V.I. Lebedev, D.N. Laikov, “A quadrature formula for the sphere of the 131st algebraic
order of accuracy,” Doklade Mathematics, 59, p. 477 (1999).
[7] W.S. Yang, et al., “Neutronics modeling and simulation of SHARP for fast reactor analysis,”Nucl. Eng. Technol., 42, p. 520-544 (2010).
[8] C.P. Thurgood, et al., “The TN quadrature set for the discrete ordinates method,”
Transasction of the ASME, 117, p. 1068 (1995).
[9] F.K. Chan, E.M. O’Neill, “Feasibility study of a quadrilateralized spherical cube earth database,” EPRF Technocal Report 2-75, Computer Sciences Corporation (1975).
[10] M. Tegmark, “An icosahedron-based method for pixelizing the celestial sphere,” The Astrobnomical Journal, 470, L81 (1996).
[11] C. Rabiti, et al., “Spherical quadratures for the discrete ordinates method,” Trans. Am.
Nucl. Soc., 96, p. 475-476 (2007).
[12] C. Ahrens, C. Beylkin, “Rotationally invariant quadratures for the sphere,” Proc. R. Soc.
(London) Ser. A, 465, p. 3103-3125 (2009).

13