炉物理の研究 第74号 (2022年3月)＜第 15 回 炉物理部会賞受賞 記念寄稿＞
時間依存の非均質中性子輸送計算に関する計算手法の高度化
株式会社原子力エンジニアリング
解析サービス本部 解析技術グループ
辻田 浩介
1. はじめに
このたびは「時間依存非均質中性子輸送計算に関する計算手法の高度化」に関して、炉物理部会賞奨励賞をいただきまして、誠にありがとうございます。本研究成果は、名古屋大学における修士研究・博士研究の中で得られたものであり、本研究を進めるにあたりご指導いただきました山本章夫教授、遠藤知弘准教授のお力によるところが大きいです。
本稿では、本研究の概要についてご紹介致します。詳しい式展開は参考文献に記載されておりますのでそちらをご参照いただけましたら幸いです。
2. 研究概要
本 研 究 で は Characteristics 法 を 用 い た 時 間 依 存 中 性 子 輸 送 計 算 の 高 度 化 と し て 、Characteristics 法を用いた動特性計算が当時抱えていた以下の課題に関する検討を実施しました。
(1) 時間依存輸送方程式の角度中性子束時間微分項の非等方性を厳密に取り扱った場合、計算に膨大なメモリ量を要すること。
(2) Characteristics 法自体の計算量が非常に膨大であり、大規模な体系や詳細時間ステップの計算には膨大な時間を要すること。
1 つ目の課題は時間依存輸送方程式における角度中性子束時間微分項の取り扱いに関するものです。時間依存の中性子輸送方程式を以下に示しますが、従来法では本時間微分項に対して差分近似を適用するのが一般的です。
⃗⃗⃗, 𝑡)1 𝜕𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗ ∙ ∇𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) − Σ𝑡,𝑔 (𝑟⃗, 𝑡)𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) + 𝑄𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡)= −Ω𝑣𝑔𝜕𝑡(1)

しかし、Characteristics 法を用いた動特性計算で本時間微分項に対して時間差分を取ると、無数の中性子の飛行線について 1 つ前の時間ステップにおける角度中性子束の値を保持する必要がでてきてしまい、極めて膨大なメモリ量を要します。そこで、先行研究[1]では、下式のように角度中性子束時間微分項に対して等方近似を加えることで、メモリ量の大幅な低減を図っています。
（以下、等方近似法と記載します。
）
1 1 𝜕𝜙𝑔 (𝑟⃗, 𝑡)⃗⃗⃗ ∙ ∇𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) − Σ𝑡,𝑔 (𝑟⃗, 𝑡)𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) + 𝑄𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡)
= −Ω
4𝜋 𝑣𝑔
𝜕𝑡

(2)

炉物理の研究 第74号 (2022年3月)

しかし、本近似の妥当性が定量的に評価された事例は（我々の知る限りでは）当時存在しなかったことから、本研究では等方近似法に起因する近似誤差の定量評価に着目しました。
ただし、角度中性子束時間微分項の非等方性を厳密に考慮した計算は、大規模体系では計算の実施そのものが困難であることから、まずは「実用的なメモリ量で、等方近似法よりも厳密に角度中性子束時間微分項の非等方性を取り扱う計算手法」の開発に着手しました。本
研究で開発した手法は以下の 2 つです。
・ 時空座標系を用いて n 次元体系の動特性計算を n+1 次元体系の定常計算問題に焼き直すことで、角度中性子束時間微分項に対し差分近似を適用することなく時間依存輸送方程式を解く方法 [2]・ 等方近似を適用しない場合に現れる前ステップの角度中性子束をメモリ上に保持せずに、1 ステップ前のバランス方程式から再計算する方法（再計算する 1 ステップ前のバランス方程式には等方近似法を使用）[3]
特に後者では、角度中性子束の時間微分項に1/𝑣𝑔 が掛かることを利用して、𝑣𝑔 Δ𝑡 > 1 を満たす動特性計算においては、実用的なメモリ量で等方近似誤差を大幅に低減した動特性計算が実現可能であること*、 またこれを用いて等方近似法起因の近似誤差が十分小さいことを確認しました[3]。
次に検討したのは Characteristics 法を用いた動特性計算の高速化です。本研究では、名古屋大学で開発された Multigrid Amplitude Function（MAF）法[4]の時間依存輸送計算への適用について検討しました。 MAF 法は、全中性子束を「空間に強く依存し時間に弱く依存する関数（形状関数）」と「空間に弱く依存し時間に強く依存する関数（振幅関数）」に分離し、両者の性質にあったメッシュ分割・時間ステップ分割の計算を組み合わせることで動特性計算の高速化を図る手法です。具体的には「詳細メッシュ・粗時間ステップの動特性計算」と「粗メッシュ・詳細時間ステップの動特性計算」を用いて時間・空間に関する精度を相互に補償することで、計算精度を劣化させること無く高速な計算を実現します。本研究ではCharacteristics 法を用いた動特性計算にこれを適用し、・計算コストの非常に高い Characteristics 法での計算は粗タイムステップ計算
・詳細タイムステップ計算は Characteristics 法よりも極めて高速な粗メッシュ拡散計算で行うことで、結果として従来法に比べ約 28.5 倍の高速化が達成しました[5]。なお、本成果については MAF 法そのものを開発された伴氏の功績によるところが大きいことを申し添えたいと思います。
さらに、Characteristics 法を用いた動特性計算で θ 法を適用した場合、前ステップの残余項（(4)式）も膨大なメモリ量を要する要因となるため、従来は前ステップの残余項を領域平均で均質化・等方化する等の近似が用いられることが多かった[1]のですが、MAF 法におけ*Δ𝑡はタイムステップ幅。軽水炉の熱群でも𝑣𝑔 はおよそ105 オーダーであり、Δ𝑡 = 1msec の非常に詳細な時間ステップであっても𝑣𝑔 Δ𝑡 > 102となるため、再計算される前ステップの角度中性子束に含まれる等方近似誤差が約 1/100 に低減される形になります。

炉物理の研究 第74号 (2022年3月)

る時間離散化誤差は（十分なサイズの粗メッシュで振幅関数の空間離散化が行われていれば）ほぼ振幅関数の時間離散化手法によって決定されるため、「形状関数の離散化に完全陰解法」、「振幅関数の離散化に θ 法」を採用することによって、前ステップの残余項に対する近似を完全に排除しつつ、θ 法による良好な時間積分性能を得ることができることも確認できました。[5]1 1 𝜙𝑔 (𝑟⃗, 𝑡𝑛+1 ) − 𝜙𝑔 (𝑟⃗, 𝑡𝑛 )⃗⃗⃗, 𝑡𝑛 ) + 𝜃𝑅𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡𝑛+1 )= (1 − 𝜃)𝑅𝑔 (𝑟⃗, Ω4𝜋 𝑣𝑔Δ𝑡(3)⃗⃗⃗, 𝑡) = −Ω⃗⃗⃗ ∙ ∇𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) − Σ𝑡,𝑔 (𝑟⃗, 𝑡)𝜓𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡) + 𝑄𝑔 (𝑟⃗, Ω⃗⃗⃗, 𝑡)
𝑅𝑔 (𝑟⃗, Ω(4)しかし、かなりの高速化が実現できたとはいえ、Characteristics 法を用いる以上、計算量は依然としてかなりのものであり、実時間アプリケーションが実現可能なレベルの高速化には届きません。そこで博士研究では、「Characteristics 法の結果を再現する Reduced Order
Model（ROM）の開発」について検討を進めています。
その中で注目したのが Proper OrthogonalDecomposition（POD）です[6] 。POD では、解析対象の問題を表現することに特化した直交基底を特異値分解で事前に求めておき、中性子束分布の直交基底展開（(5)式）を考えます。
これにより、中性子のバランス方程式は「直交基底の展開係数に関する式（(6)式）」に変換されるため、大幅な自由度の削減が実現できます。
𝜙⃗⃗ = 𝐔𝜑⃗⃗ (𝜙⃗⃗ ∶ 中性子束, 𝐔 ∶ 直交基底, 𝜑⃗⃗ ∶ 展開係数)

(5)

11 𝐓
𝐅𝜙⃗⃗ → 𝐔 𝐓 𝐀𝐔𝜑⃗⃗ =𝐔 𝐅𝐔𝜑⃗⃗𝑘eff𝑘eff(6)𝐀𝜙⃗⃗ =

̂𝜑→ 𝐀⃗⃗ =1𝐅̂𝜑⃗⃗𝑘eff

ただし、POD を適用するためには上式のようにバランス方程式を行列表現する必要があるのですが、Characteristics 法のように無数の中性子パスを考慮する手法でそれらすべてを行列表現することはメモリ量の観点から非現実的です。そこで本研究では、中性子束を一定（もしくは線形）とみなす非構造メッシュ領域毎の中性子バランス方程式を Characteristics
法の計算結果から再構成することで領域単位でのバランス方程式を行列で表現し、POD を適用することを考えました[7]。 本 ROM の詳細な検証は現在進行中ですが、C5G7-TD ベンチマーク問題（2D）において、Characteristics 法+MAF 法でも約 4 時間(16 スレッドの並列計算)かかった計算と同等の計算精度を、本 ROM ではシングルスレッド計算で約 3 秒†と非常に高速に再現できることが分かっております。
3. おわりに
今では国内外の様々な機関にて Characteristics 法を用いた動特性計算コードが開発されていますが、本研究を始めた当初は今ほど数が多くありませんでした。そのため、修士研究・博士研究で検討に使用した解析コードは全て一から自作・検証したものになります。この炉
心計算コードを”自分で書く”機会は炉心解析に関する理論やコード開発技術を学ぶ上で非†直交基底の構築にかかる計算時間を除く。

炉物理の研究 第74号 (2022年3月)

常に良い機会でしたし、毎日少しずつコードが完成に近づいていく姿を見るのは非常に楽しいものでした。また、コード開発に際しては、本稿では紙面の関係上書ききれなかった数
多くの計算理論や高速化のテクニックを、先生方を始め、研究室の先輩方からも教えていただきました。先生方、先輩方からのご指導無しでは本成果は得られなかったことを、この場をお借りして深く御礼申し上げます。今後も炉物理分野の技術発展に少しでも貢献できるよう、日々精進して参りたいと思います。
参考文献
1.

Hursin M, Downar TJ, Thomas J. PWR control rod ejection analysis with the method ofcharacteristic code DeCART. Proc PHYSOR 2008; 2008 Sep 14–19; Interlaken, Switzerland.

2.

Tsujita K, Endo T, Yamamoto A, Kinetic Calculation Method in Space-Time Frame Using
Characteristics Line. Trans Am Nucl Soc. 2012;106,743–746.

3.

Tsujita K, Endo T, Yamamoto A, Higher order treatment on temporal derivative of angular flux
for time dependent MOC. Proc M&C 2013; 2013 May 5–9; Sum Valley, Idaho.

4.

Ban Y, Endo T, Yamamoto A. A unified approach for numerical calculation of space-dependent
kinetic equation. J Nucl Sci Technol. 2012;49(5):496–515.

5.

Tsujita K, Endo T, Yamamoto A. Application of the multigrid amplitude function method fortime-dependent MOC based on the linear source approximation. J Nucl Sci Technol.
2020;57(6):646–662.

6.

Tsujita K, Endo T, Yamamoto A. Fast reproduction of time-dependent diffusion calculations using
the reduced order model based on the proper orthogonal and singular value decompositions. JNucl Sci Technol. 2021;58(2):173–183.

7.

Tsujita K, Endo T, Yamamoto A. Fast reproduction of time-dependent MOC calculations using
the reduced order model based on the proper orthogonal and singular value decompositions. Trans
Am Nucl Soc. 2020;123:1349–1353.

炉物理の研究 第74号 (2022年3月)＜第15回 炉物理部会賞受賞 記念寄稿＞
研究題目：燃料デブリによる臨界事故時の影響評価に関する研究（研究実施時）東京工業大学 環境・社会理工学院（現）国立研究開発法人 日本原子力研究開発機構
安全研究センター 臨界安全研究グループ
福田 航大

はじめにこの度は炉物理部会賞奨励賞という栄えある賞を頂く機会に恵まれましたこと、大変光栄に感じております。この場を借りて炉物理部会賞関係者の皆様にお礼申し上げます。
受賞対象となった研究は、私が東京工業大学修士課程・博士課程に在学している間、先導原子力研究所小原研究室にて取り組んだ研究です。５年間に亘り丁寧に指導してくださった小原徹
教授をはじめ、様々なご助言をくださった西山潤助教授、以下に紹介する投稿論文で共著者として共に研究をさせて頂いた村本武司さんと Tuya Delgersaikhan に心よりお礼申し上げます。また学会発表の場などで他大学の先生方や研究者の方々から頂いた本研究に対する真摯なコメン
ト無くしては本研究を深めることができなかったと身に染みて感じております。皆様、ありがとうございました。

研究の概要
さて本研究は、福島第一原子力発電所廃止措置での燃料デブリ取出しの際の臨界事故のリスクを適切に評価することは作業時の安全確保のため必要である、という背景のもと、燃料デブリによる臨界事故時の影響評価上重要な因子を明らかにすることを目的としたものです。
燃料デブリ体系における動特性解析並びに臨界事故解析は、非常に複雑なものです。様々な性状を持つ燃料デブリが弱く結合し（さらに時間と共にその位置が複雑に変化する可能性もあり）一点炉近似モデルでは扱うことができません。動特性解析において重要な要素である反応度フィードバック効果は燃料温度上昇のみならず、放射線分解ガスをはじめとしたその他の効果をどのように扱うべきかという問題があります。また、いざ影響評価を行おうとすると水による遮蔽の効果と反応度の影響を同時に考慮する必要があります。
本研究は、このような複雑な条件下での解析を精度よく行う上で重要なパラメータを抽出しそれらが結果に与える影響を明らかにしたもので、臨界安全解析上重要な知見を得ることができたのではないかと考えています。

炉物理の研究 第74号 (2022年3月)

大変恐縮ですが、研究の詳細は以下の論文に纏められておりますので、興味を持って頂いた方には是非ご覧いただければと思います。

◼

Kodai Fukuda, Dergelsaikan Tuya, Jun Nishiyama and Toru Obara, “Radiation Dose Analysis
in Criticality Accident of Fuel Debris in Water,” Nuclear Science and Engineering. Vol.194, No.
3, pp.181-189 (2020); https://doi.org/10.1080/00295639.2019.1665459.

◼

Kodai Fukuda, Jun Nishiyama and Toru Obara, “Reactivity feedback effect on supercriticaltransient analysis of fuel debris,” Nuclear Science and Engineering. Vol.194, No. 7, pp.493-507(2020); https://doi.org/10.1080/00295639.2020.1743580.

◼

Kodai Fukuda, Jun Nishiyama and Toru Obara, “Supercritical Transient Analysis for RampReactivity Insertion Using Multi-region Integral Kinetics Code,” Nuclear Science andEngineering; Vol. 195, pp.453-463 (2021); https://doi.org/10.1080/00295639.2020.1847979.

◼

Kodai Fukuda, Jun Nishiyama and Toru Obara, “Impact of Change of Neutron Generation Timeon Released Energy in Falling Fuel Debris at Supercritical Condition,” Journal of Nuclear
Science and Technology; Vol. 59, Issue 2, pp.230-239 (2021);
https://doi.org/10.1080/00223131.2021.1963344.

最後に炉物理部会のホームページに掲載されている過去の部会賞受賞者一覧を拝見すると、炉物理界で輝かしいご活躍をされている方々のお名前が並んでいます。冒頭に申し上げたように、この研究は小原先生をはじめとした多数の方々のお力添えがあって進めることができたものです。ですので、正直なところ、私なんかがここに名前を載せていただいてもいいのだろうかという気持ちはあります。しかし、その賞の名に恥じぬように研究を頑張っていこうという気持ちで今後過ごしていきたいと考えています。至らない点ばかりですが、今後ともご指導・ご鞭撻のほどよろしくお願い致します。
福田 航大
2022 年 2 月