JP2007192550A - 粒子輸送の計算方法、計算装置およびプログラム、核定数の計算方法、計算装置およびプログラム、ならびに、原子炉シミュレーション方法、原子炉シミュレータおよび原子炉シミュレーションプログラム - Google Patents

Abstract

【課題】任意形状・寸法の燃料集合体の核定数を短時間で評価できるようにする。
【解決手段】粒子輸送の計算方法において、低食い違い数列を発生する工程（Ｓ３）と、レイトレースの始点および方向を前記数列で指定する工程（Ｓ４）と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程（Ｓ５）と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程（Ｓ６）と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程（Ｓ７）と、を有する。ここで得られた粒子の分布から核定数を作成する（Ｓ８）。
【選択図】図１

Description

本発明は、粒子輸送の計算方法、計算装置およびプログラム、核定数の計算方法、計算装置およびプログラム、ならびに、原子炉シミュレーション方法、原子炉シミュレータおよび原子炉シミュレーションプログラムに関する。

原子炉の炉心設計に用いる原子炉シミュレーション装置は、一般的に、燃料棒や燃料集合体について格子配列等の規則的な体系を対象としている。この原子炉シミュレーション装置は、燃料集合体の核断面積の特性等を表す核定数を作成する二次元詳細計算部と、この核定数を用いて全炉心の核特性を評価する三次元シミュレータから構成されている。

二次元詳細計算部は、中性子輸送格子計算コードとも呼ばれる。二次元詳細計算部は、燃料集合体の長手方向（Ｚ方向）に垂直な、ある水平断面（Ｘ−Ｙ断面）に関して、燃料集合体の無限中性子増培係数、燃料棒毎の出力密度等（核特性）を、燃料集合体の燃焼度に応じて評価する。この評価に際しては、燃料棒の構造、配置、チャンネルボックス等の構造物も考慮する。ただし、燃料集合体がＸ−Ｙ方向にわたって、規則的に無限個配置されていると仮定し、Ｚ方向もそのＸ−Ｙ断面が無限に続いていると仮定している。

実際の炉心は、Ｘ−Ｙ方向もＺ方向も有限であるので、二次元詳細計算部で評価した核特性は、実際の炉心の核特性とは異なる。実際の炉心の核特性を評価するには、Ｚ方向がたとえば約１５ｃｍの燃料集合体毎に１つのメッシュとして扱う。すなわち、１５ｃｍ×１５ｃｍ×１５ｃｍの立方体を数値計算の１メッシュとして、この中では燃料や減速材等が均質に分布していると仮定する。実際の炉心を三次元体系で計算するには、多大な計算時間がかかるため、このような簡易化を行わざるをえない。この計算を行うプログラムを三次元シミュレータと呼ぶ。

このように二次元詳細計算と三次元シミュレータとに分割して評価するためには、多くの工夫が必要である。炉心の形状や燃料集合体の組成と形状を大きく変えると、これらの工夫が有効でなくなり、核特性評価手法の限界から、炉心や燃料集合体の設計が制約されることも多かった。

ところで、従来の沸騰水型原子炉（ＢＷＲ）の定型的燃料集合体は、たとえば、集合体の断面寸法が約１５ｃｍ×１５ｃｍで、燃料棒の配置は正方格子で燃料棒間隔も一定であった。二次元詳細計算部は、このような定型的な燃料集合体の核定数を求める。一方、将来炉の設計においては、燃料棒の配置は正方格子のみならず三角格子を採用したり、燃料集合体の寸法も大きくなる場合もある。このため、任意形状・寸法の燃料集合体の核定数評価が要求されている。任意形状・寸法の燃料集合体の核定数評価には非決定論的手法であるモンテカルロ法と決定論的手法であるCharacteristic法が用いられている。

群定数をベースにした場合には、モンテカルロ法とCharacteristic法は精度が同等である。モンテカルロ法では、精度が計算時間の平方根で決まるため、十分な精度を得るためには計算時間がかかるという欠点がある。サンプリングに準乱数を用いる準乱数モンテカルロ法は、計算精度が計算時間で決まるために、従来の確率的サンプリングを用いたモンテカルロ法よりも計算速度が速い長所がある。しかし、ランダムサンプリングにならないので、従来の中性子輸送モンテカルロ法の乱数を準乱数で置き換えるだけでは実現できない。Characteristic法は、衝突確率を求める際に燃料集合体を満遍なくフォローしていることの保証が難しい点が欠点である。

原子炉の核特性は、次式で表される中性子輸送ボルツマン方程式を解くことで求まる。

ここで、ｒは中性子が位置する空間座標を示すベクトル、Ｅは中性子の運動エネルギー、ｗは中性子の運動方向を示すベクトルである。φ（r,E,Ω）は、r,E,Ωにおける中性子束であり、ボルツマン方程式の解である。式（１）の左辺第１項「Ω・∇」は方向ベクトルΩと微分演算子∇の内積を表すもので、展開すると次式になる。

Σ_T（r,E）は、中性子の巨視的全断面積である。Σ_S（r,E'→E,Ω'→Ω）は、エネルギーがE'で、運動方向がΩ'の中性子が散乱して、エネルギーがE、運動方向がΩに変わる確率を表す散乱断面積である。S_f（r,E,Ω）は、核分裂中性子源である。

式（１）の２番目の式として表されている核分裂中性子源S_f（r,E,Ω）は、臨界固有値ｋ、核分裂で発生する中性子のエネルギースペクトルχ（E）、核分裂で発生する中性子の数ν（E'）および巨視的核分裂断面積Σf（r,E）から構成されるが、中性子束φ（r,E,Ω）を含んでいるため、式（１）の１番目の式は、方程式の右辺が０となる非斉次線形微分方程式であることがわかる。

非斉次微分方程式は、基本的に非斉次連立方程式と同じようなものである。たとえば、次式で表される形になる連立方程式を非斉次方程式という。

この式の右辺を左辺に移項すると、次式のように、連立方程式の右辺が０となる方程式となる。

非斉次方程式が０以外の解を持つには、次式を満足すればよい。

この解λを固有値という。式（１）の場合、１番目の式の核分裂中性子源Ｓ_f（r,E,Ω）に、２番目の式を代入すると、式（５）においてλを１／ｋとみなしたことになる。

非斉次微分方程式には固有値が伴うが、原子炉の臨界問題ではこの固有値の逆数ｋを固有値と呼ぶことが多く、その物理的意味は中性子の増倍率である。固有値問題は、初めに核分裂源と固有値を仮定し、式（１）の１番目の式の右辺が０でないと仮定した斉次線形微分方程式を反復して解く解法がある。ここでは、核分裂中性子源が与えられているとする。

中性子輸送ボルツマン方程式には多くの解法があるが、基本的には次式で表されるグリーン関数の問題に帰着できる。

グリーン関数Ｇ（r,E,Ω｜r₀,E₀,Ω₀）は、空間r₀にエネルギーE₀、運動方向Ω₀の中性子が１個発生したときの中性子の分布で、空間の１点に単位入力が与えられたときの体系全体の応答である。グリーン関数Ｇがわかれば、中性子束φ（r,E,Ω）は、次式のように、グリーン関数Ｇと核分裂中性子源Ｓ_fの積分で求めることができる。

グリーン関数Ｇは、単位入力に対する応答である。実際の問題では、各点（r₀,E₀,Ω₀）に核分裂中性子源Ｓ_f（r₀,E₀,Ω₀）の入力があり、式（７）は単位入力の応答に実際の入力を積算したものになる。

グリーン関数を顕わに用いた方法にモンテカルロ法がある。モンテカルロ法の解法は、概ね次の通りである。

（Ａ）空間r₀に、エネルギーE₀、運動方向Ω₀の中性子を１個発生させる。

（Ｂ）空間r₀とエネルギーE₀から巨視的断面積Σ_T（r₀,E₀）が決まり、０から１の一様乱数と巨視的断面積を用いて、散乱までの飛行距離Ｌを求める。

（Ｃ）散乱を起こす空間r₁＝r₀＋LΩ₀での散乱断面積Σ_S（r₁,E₀→E₁,Ω₀→Ω₁）と一様乱数から散乱後のエネルギーE₁および運動方向Ω₁をサンプリングする。

（Ｄ）（Ｂ）から（Ｃ）の操作を終了条件まで繰り返す。

（Ｅ）（Ａ）から（Ｄ）について複数回繰り返す。

（Ｄ）における終了条件は、体系から中性子が漏れるとか、吸収されて中性子が消滅する等の直感的に判り易いものものあるが、他にも多くの高度な手法が提供されている。

（Ａ）から（Ｄ）の手続き中に、着目している空間r₀、エネルギーE₀、運動方向Ω₀の微小空間を通過した中性子の飛跡の長さ、もしくはそこでの散乱数を数え上げていき、（Ａ）から（Ｄ）を繰り返した回数で平均したものが、グリーン関数Ｇ（r,E,Ω｜r₀,E₀,Ω₀）である。

（Ａ）における発生する空間r₀に、エネルギーE₀および運動方向Ω₀を核分裂中性子源Ｓ_f（r₀,E₀,Ω₀）に従って決めると、次式で表される重畳積分そのものになる。

何故なら、中性子の発生点を空間r₀、エネルギーE₀、運動方向Ω₀に限定した場合に、（Ｃ）のサンプリングで得られるものはグリーン関数Ｇ（r,E,Ω｜r₀,E₀,Ω₀）そのものであるからである。中性子の発生点を核分裂中性子源Ｓ_f（r₀,E₀,Ω₀）に従って変えると、Ｓ_f（r₀,E₀,Ω₀）が大きいところでは多数の中性子をサンプリングし、小さいところでは少数の中性子しかサンプリングしない。無限に近い個数の中性子をＳ_f（r₀,E₀,Ω₀）に従って発生すると、（Ｃ）のサンプリングで得られたものはグリーン関数をＳ_f（r₀,E₀,Ω₀）倍して足し合わせたものに近づき、式（８）と同じ操作になる。

モンテカルロ法の長所としては、任意の形状、組成の問題について厳密解に近い答えを出すこと、扱う次元が増えても計算時間の増加は少ないこと、共鳴反応のような中性子エネルギーで大きく変化する断面積でも正確に扱うことができること、統計的サンプリングに基づいているので、真値との乖離を標準偏差で見積もることができること、が挙げられる。

一方、モンテカルロ法の短所としては、統計的サンプリングを用いているので、多数の繰り返し回数が必要となり、基本的に計算時間がかかることが挙げられる。真値との乖離を示す標準偏差は、繰り返し回数の平方根に逆比例する。たとえば、１％の標準偏差を０．１％にまで改善するためには、１００倍の繰り返し回数が必要となり、計算時間も１００倍となる。この法則はクラメル・ラオの不等式として知られているもので、普遍的なものである。統計的サンプリングを基にするかぎり、この欠点から逃れることはできない。

また、統計的サンプリングに基づくモンテカルロ法は、中性子エネルギーを連続的に扱わなければならない問題や、組成・形状が複雑な問題には適しているが、より単純化された問題に限定すると過剰仕様となることが多い。

組成・形状が任意といっても燃料集合体の組成・形状は限定されてくることが多く、中性子エネルギーも連続的に考慮する必要がない場合もある。たとえば、中性子断面積をエネルギーに関して短冊状の階段関数で近似して表す群定数と呼ばれる定数が、核特性評価にはよく用いられている。

断面積を群定数で近似した場合にも、モンテカルロ法は有効であり、連続エネルギーで扱うよりも計算速度は改善されるが、標準偏差が繰り返し回数の平方根に逆比例することは避けられない。

特許文献１ないし特許文献３には、モンテカルロ法の計算時間を短縮できる方法として、相関サンプリングモンテカルロ法が開示されているが、精度が繰り返し回数（標本数）の平方根に逆比例することは変わらない。

サンプリングに準乱数を用いる準乱数モンテカルロ法は、計算精度が計算時間の平方根に反比例するのではなく、計算時間に反比例するため、従来の確率的サンプリングを用いたモンテカルロ法よりも計算時間が短い長所がある。この特性は、１９９０年代ごろから特に金融の分野で注目されてきた。金融の分野では数千、数万の資産を組み合わせた多次元ポートフォリオ問題の観点から、モンテカルロ法と準モンテカルロ法の適用が検討されてきた（たとえば特許文献４参照）。

準乱数モンテカルロ法を中性子輸送モンテカルロ法に適用する場合には、従来の中性子輸送モンテカルロ法の乱数を準乱数で置き換えるだけでは実現できない。なぜなら、準乱数はランダムに発生するものではないので、確率モデルに適用するにはサンプリングに準乱数を用いても確率事象をランダムなものではなく決定論的に表現するので、サンプリングの過程は、確率的事象にならないからである。この問題は、金融工学の世界でも指摘されていて、たとえばランダムウォークを準乱数でシミュレートしたときには、特に多次元問題で評価した結果が真値に収束しないことや、誤差が計算時間に反比例していかないことが、典型的な問題点として挙げられる。

次に、Characterisitc法について説明する。Characterisitic法も多くのバリエーションがあるが、ここではレイトレーシング手法に基づくGTRANコードについて説明する。

グリーン関数で、全く散乱を経ずに直進する中性子の成分のみを考える。グリーン関数を記述したボルツマン方程式（式（６））の左辺第三項は、散乱した成分を表す。散乱後、エネルギーも運動方向も変化しない確率は０なので、散乱項は０になる。したがって、このときのボルツマン方程式は、次式で表される。

この場合、運動方向は全く変わらないので、中性子の座標を表すベクトルrと中性子の発生点を表すベクトルr₀との差は、r - r₀ = L Ω₀となる。したがって、式（９）はΩ₀方向の飛行距離Lに関する１次元の微分方程式となる。この微分方程式は、次式で表される。

図６に衝突を起こさない中性子のグリーン関数Ｇ（L,E_g,Ω｜L₀,E_g,Ω₀）の変化を示す。ここで、巨視的全断面積によって領域を区分する。原子炉では、たとえば領域１が水、領域２が燃料被覆管、領域３が燃料ペレットに相当する。各領域内では、グリーン関数は次式で示されるように指数減衰の形となる。

次式で表される領域内で吸収された成分は、捕獲反応等で消滅する分と、散乱等によってエネルギーが変化したもの、方向が変化したものに分けられる。この描写は光が直進していくモデルと同じである。

散乱した中性子も考慮した連続的な分布関数としてのグリーン関数を用いる場合と、散乱しない中性子のグリーン関数のみを用いる場合との違いは、中性子束を求める時の積分範囲の違いである。原子炉の場合では、中性子の発生源が燃料ペレットのみなので、式（７）の積分範囲は燃料ペレットのみでよく、グリーン関数の始点もペレット領域のみとなる。

散乱しない中性子のグリーン関数のみを用いる場合には、散乱項を発生源と見なすので、燃料ペレット以外の領域も積分範囲に含まれ、グリーン関数の始点も全領域に及ぶ。この衝突を全く起こさない中性子に関するグリーン関数Ｇ（L,E_g,Ω｜L₀,E_g,Ω₀）で中性子束を表すと、式（７）の積分は次式で表され、右辺第２項に中性子束を含む漸化式となる。

実際に適用するには、これを離散化して燃料集合体内を有限の領域に分けて、そこの平均値で代表する。離散化には多くの手法が考えられる。

図７に、従来のcharacteristic法でのスキャンの例を示す。図７は４本の燃料棒が２×２の単位セルで配置したものであり、それぞれのセルについて水領域を４つ、被覆管、ペレットをそれぞれ１領域として、合計２４領域に分割されている。あるエネルギーE_gについて、角度を８成分（0°、45°、90°、135°、180°、225°、270°、315°）に分割したとすると、図７中の矢印の線で示すように各中性子の移動を１本の直線で表すことができる。すなわち、各中性子を１条の光線と見立てて、光で炉心内を走査するイメージとなる。走査することをレイトレースといいい、レイトレースをスキャンとも言う。これらのスキャンする線を増やすほど精度は上がる。グリーン関数の発生点の選択としては、これらの矢印の交点をとることも考えられる。

しかし、このように規則的に発生点を選ぶと、図７に示すように、粒状ガドリニヤ３１のような小さな領域は、全くスキャンできない。これを避けるには、ユーザが事前に発生点を設定してやればよいが、体系が複雑になると非常に煩雑になり、それのみに頼るのは現実的ではない。最適な発生点と方向を自動的に設定する手法が必要である。

次に、準モンテカルロ法について説明する。準モンテカルロ法の基本は、ランダムではなく決定論的数列を用いるところにある。積分領域Iにおける準モンテカルロ法による数値積分は、次式で表される。

準モンテカルロでは、積分範囲Iの元である点列x₁,…,x_Nをランダムではなく、一様に分布するように選ぶ。これらの点列の一様性は、数学的には食い違い性（discrepancy）で定義する。

このような数列は、正式には低食い違い数列（low−discrepancy sequences）または準乱数と呼ばれている。これらの数列は完全に決定論的なものでランダム成分は全くない。低食い違い数列の利点は、その一様性にある。

代表的な低食い違い数列には、Halton列、Sobol列、Faure列がある。これらの低食い違い数列は、van der Corputが１９３５年に提案した、次元s=1、基底b=2のvan der Corput列に由来している。１９６０年にHaltonは、このvan der Corput列を任意の次元に拡張した。これをHalton列という。１９６７年に、SobolがHaltonとは異なったアプローチで多次元の低食い違い数列を提案した。これをSobol列という。１９８２年にFaureがSobolのアイデアを拡張して、Sobol列を置き換える低食い違い数列を見いだした。これをFaure列という。Faure列は、Halton、Sobol列に比べて、多次元積分に関して極めて優れている。
特開平１１−２３７８７号公報 特開平１１−２９５４７２号公報 特開２００１−４７７９号公報 特開２００１−１１７９７４号公報

任意の体系・組成の燃料集合体の核定数を求めるには、モンテカルロ法が有効であるが、統計的サンプリングを用いているので、真値との乖離を示す標準偏差が繰り返し回数の平方根に逆比例する。したがって、要求精度が高くなるにつれて計算時間が急速に増加し、実用的でなくなる。

また、準乱数モンテカルロ法は、誤差が繰り返し回数（標本数）に反比例し、モンテカルロ法の計算時間の問題を解決できるが、通常のモンテカルロ法の乱数を準乱数に置き換えるだけでは実現できない。

また、Characteristic法は、計算時間はモンテカルロ法のような限界はないが、任意の体系・組成の燃料集合体では、中性子の発生点や運動方向を適切に選択することが難しい欠点がある。

そこで、本発明は、任意形状・寸法の体系における粒子の輸送計算を短時間で行うことができるようにすることを目的とする。さらに、その粒子として中性子を対象とした場合には、原子燃料の核定数を短時間で評価できるようにすることを目的とする。

上記目的を解決するため、本発明は、粒子輸送の計算方法において、低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、レイトレースの始点および方向を前記数列で指定する工程と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、を有することを特徴とする。

また、本発明は、原子燃料の核定数の計算方法において、低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、レイトレースの始点および方向を前記数列で指定する工程と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、前記粒子の分布から核定数を求める工程と、を有することを特徴とする。

また、本発明は、原子炉シミュレーション方法において、低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、レイトレースの始点および方向を前記数列で指定する工程と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、前記粒子の分布から核定数を求める工程と、前記核定数に基づいて、前記原子炉の三次元の特性を評価する工程と、を有することを特徴とする。

また、本発明は、粒子輸送の計算装置において、計算条件を入力する入力手段と、前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する演算手段と、前記粒子の分布を外部に出力する外部出力手段と、を有することを特徴とする。

また、本発明は、原子燃料の核定数の計算装置において、計算条件を入力する入力手段と、前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価し、前記粒子の分布から核定数を求める演算手段と、前記核定数を外部に出力する外部出力手段と、を有することを特徴とする。

また、本発明は、原子炉シミュレーション装置において、計算条件を入力する入力手段と、前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価し、前記粒子の分布から核定数を求める二次元計算部と、前記核定数および前記計算条件に基づいて、前記原子炉の三次元の特性を評価する三次元シミュレータと、を有することを特徴とする。

また、本発明は、粒子輸送の計算プログラムにおいて、コンピュータに、低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、レイトレースの始点および方向を前記数列で指定する手順と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、を実行させることを特徴とする。

また、本発明は、原子燃料集合体の核定数の計算プログラムにおいて、コンピュータに、低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、レイトレースの始点および方向を前記数列で指定する手順と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、前記粒子の分布から核定数を求める手順と、を実行させることを特徴とする。

また、本発明は、原子炉シミュレーションプログラムにおいて、コンピュータに、低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、レイトレースの始点および方向を前記数列で指定する手順と、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、前記粒子の分布から核定数を求める手順と、前記核定数に基づいて、前記原子炉の三次元の特性を評価する手順と、を実行させることを特徴とする。

本発明によれば、任意形状・寸法の体系における粒子の輸送計算を短時間で行なうことができる。さらに、その粒子として中性子を対象とした場合には、原子燃料の核定数を短時間で評価できる。

図２は、本発明の一実施形態の原子炉シミュレーション装置のブロック図である。

原子炉シミュレーション装置１は、燃料集合体の核断面積の特性等を表す核定数を作成する二次元詳細計算部２と、この核定数を用いて全炉心の核特性を評価する三次元シミュレータ３を有している。原子炉シミュレーション装置１は、たとえば１つの演算装置を有する１台のコンピュータであって、二次元詳細計算部２および三次元シミュレータ３は、それぞれの機能を実現するための計算プログラム（コード）によって動く１つの演算装置として実現することができる。

また、原子炉シミュレーション装置１は、入力条件等を入力するための入力手段５、計算結果を外部に出力するための外部出力部６および核定数を記憶しておく核定数記憶部４を有している。入力手段５は、たとえばキーボードやマウスである。表示手段６は、たとえば液晶ディスプレーやプリンタである。核定数記憶部４は、たとえばハードディスクである。

二次元詳細計算部２は、中性子輸送格子計算コードとも呼ばれる。二次元詳細計算部２は、燃料仕様等、すなわち、ウランの濃縮度やガドリニアなどの可燃性毒物の濃度、燃料棒の構造、配置、チャンネルボックス等の構造が入力されて、計算を行う。この計算は、たとえば燃料集合体の長手方向（Ｚ方向）に垂直な、ある水平断面（Ｘ−Ｙ断面）に関して、燃料集合体の無限中性子増培係数、燃料棒毎の出力密度等（核特性）を、燃料集合体の燃焼度に応じて評価する。

二次元詳細計算部２で計算された結果のうち、三次元シミュレータ３で用いられるものは核定数として核定数記憶部４に保存される。

二次元詳細計算部２は、燃料集合体がＸ−Ｙ方向にわたって、規則的に無限個配置されていると仮定し、Ｚ方向にはそのＸ−Ｙ断面が無限に続いていると仮定している。しかし、実際の炉心は、Ｘ−Ｙ方向もＺ方向も有限であるので、二次元詳細計算部で評価した核特性は、実際の炉心の核特性とは異なる。そこで、二次元詳細計算部２で作成した核定数４を用い、さらに運転条件等を入力して、三次元シミュレータ３によって炉心内の三次元の特性を評価する。

炉心の核特性を評価するには、Ｚ方向がたとえば約１５ｃｍの燃料集合体毎に１つのメッシュとして扱う。すなわち、１５ｃｍ×１５ｃｍ×１５ｃｍの立方体を数値計算の１メッシュとして、この中では燃料や減速材等が均質に分布していると仮定して計算を行う。

図１は、本発明の一実施形態の二次元詳細計算部の計算の流れ図である。

まず、ユーザは、燃料集合体の形状と組成を入力する（工程Ｓ１）。すると、コンピュータは磁気ディスク４１に格納された群定数と組成・形状データを取り込む（工程Ｓ２）。

その後、準乱数を発生させ（工程Ｓ３）、その準乱数（低食い違い数列）を用いてレイトレースのスキャン範囲を自動的に設定する（工程Ｓ４）。その後に、レイトレースの処理を実施し（工程Ｓ５）、グリーン関数を作成する（工程Ｓ６）。そのグリーン関数を用いて、中性子輸送ボルツマン方程式のソルバーによってボルツマン方程式を解き（工程Ｓ７）、核定数を作成する（工程Ｓ８）。

本実施形態では、準乱数として、ｓ次元のFaure列を用いる。次に、このFaure列の生成方法を説明する。

まず、ｒを、２とｓのどちらか大きい方の整数以上で最小の素数とし、規定ｒに関する任意の整数ｎを以下の多項式で記述する。

ここで、ａ^ｉ _ｊ（ｎ）は０以上ｒ未満の整数である。ｎ番目のFaure数列における、第1次のFaure数は、次式で表される。

これ以降の、第ｋ次のFaure数は、第（ｋ−１）次のFaure数から以下の漸化式で表される係数ａ^ｋ _ｊ（ｎ）によって表される。

ここで、Ｃ^ｉ _ｊは次式で表される二項係数である。

係数ａ^ｋ _ｊ（ｎ）を用いて、この第ｋ次のFaure数は次のようになる。

このようにして、ｎ番目のS個のFaure数列を決めることができ、それをφ_ｎ＝（φ^１ _ｎ、φ^２ _ｎ，…，φ^Ｓ _ｎ）と記述する。

本実施形態では、３次のFaure列を用いてレイトレースの始点と方向を設定する。すなわち、φ_ｎ＝（φ^１ _ｎ、φ^２ _ｎ，φ^３ _ｎ）を用いて、中性子発生点（x_n,y_n）および運動方向θnを設定する。x_n,y_n,θ_nは次式で与える。

x_n = L_x φ¹ _n + x₀
y_n = L_y φ² _n + y₀
θ_n = N φ³ _n Δθ
ここで、（x₀,y₀）は原点座標、L_xおよびL_yはスキャン範囲のＸ方向およびＹ方向の長さである。Nは角度の分割数で、たとえばN＝８である。Δθは角度の分解能で、たとえばΔθ＝４５°である。

図３は、本実施形態の二次元詳細計算部で設定した中性子発生点および運動方向の例である。図中の黒丸は中性子発生点１１を示していて、矢印は中性子移動方向１２を示している。この例は、１辺が１０ｃｍの正方形の中に１０００個の中性子を発生させたものであり、図３はその１辺が３ｃｍの正方形を拡大したものである。このため、図３にプロットされている中性子は９０個である。

このようにして発生させた中性子のレイトレースを行う（工程Ｓ４）。図３において、領域は水２１、被覆管２１およびペレット２３に区分されている。中性子移動方向１２に沿って、中性子の軌跡がどのように各領域を通過するかを調べる。

各領域内ではグリーン関数は、式（１２）のように指数減衰の形となる。これを基に、それぞれの中性子の軌跡に沿ったグリーン関数を作成する（工程Ｓ６）。

このグリーン関数を用いて、式（１０）で表されるボルツマン方程式を解き、粒子の分布を評価する（工程７）。この解を整理することにより、ある燃料のあるＸ−Ｙ平面の核定数４ができる（工程Ｓ８）。

必要に応じて二次元詳細計算部２によって計算を繰り返し、必要な核定数４を準備する。その後、運転条件等の入力を与えることにより、三次元シミュレータ３は核定数４を用いて炉心三次元計算を行う。炉心三次元計算の結果は、表示手段６によって表示される。

従来のCharacteristic法では衝突確率を求める際に、中性子の運動方向や発生空間を指定しなければならない。その指定を誤ると、結果は真値と大きく乖離する。中性子の発生点を規則的に選択すると、微小領域などをスキャンできない可能性がある。また、擬似乱数を用いる場合で、任意形状の場合には、適切なスキャン範囲を自動的に設定できないこともある。そこでレイトレースのスキャン範囲が設定された後に、適切にスキャン範囲が設定されているか否かを判定し、適切にスキャン範囲が設定されていない場合には、ユーザが手入力で範囲を設定することが必要になる。

しかし、本実施形態では、これらの指定に準乱数（低食い違い数列）を用いることで、任意の場合について燃料集合体を満遍なくフォローしていることを保証できるため、スキャン範囲が適切か否かの判定や、手入力による範囲の設定は不要である。

また、Faure列を用いた数値積分には、以下のような利点がある。

サンプル数Nとすると、通常のモンテカルロ法ではN^-1/2に誤差が比例するのに対して、準モンテカルロ法の誤差は(logN)^(s-1)/2/N^-1に比例して変化する。通常のモンテカルロ法に用いる擬似乱数の周期は比較的短い場合が多いが、Faure列のような低食い違い数列には周期はない。また、高速な収束性とそれに伴う高精度を保証できる。

図４は、擬似乱数と準乱数によって多数の二次元の点をプロットした図であって、（ａ）は擬似乱数によるもので、（ｂ）は準乱数によるものである。これらは、擬似乱数および準乱数によって、２つの系列の数列x₁,x₂,…,x_nとy₁,y₂,…,y_nを発生し、（x₁,y₁）,（x₂,y₂）,…,（x_n,y_n）の点をプロットしたものである。

擬似乱数はランダムに面を埋めていくが一様性に劣るために、ランダムに隙間ができているのに対して、Faure準乱数では規則的に隙間ができ一様に面を埋めていくことがわかる。

ここで、擬似乱数を用いたモンテカルロ法と準乱数を用いた準モンテカルロ法による標本数と誤差の関係を調べるため、次式で表されるπの数値積分を行った。なお、モンテカルロ法によるπの計算は円の数値積分の例としてよく引用される。

ここで、ξ_ｉは０から１の間をとる実数である。モンテカルロ法ではξ_ｉとして０から１の一様乱数を用いる。準モンテカルロ法では低食い違い数列（準乱数）を用いる。

図５は、数値積分で求めたπの誤差を標本数Ｎの関数として表したものである。図５の１３はモンテカルロ法による数値積分の結果、１４は準モンテカルロ法による数値積分の結果を示したものである。なお、ここで、πの誤差とは、式（１６）の数値積分で求めた値とπの真値の差の絶対値である。

図５から、モンテカルロ法では誤差が標本数の平方根に反比例するが、準モンテカルロ法では標本数に反比例することがわかる。

低食い違い数列としては、Faure列の代わりに、Halton列、Sobol列、または、Niederreiter列を用いてもよい。また、Faure列、Halton列、Sobol列、または、Niederreiter列の変形を用いてもよい。

また、低食い違い数列の代わりに、擬似乱数を用いてもよい。なお、この場合には誤差は標本数の平方根に反比例するが、準乱数の規則性が問題となる場合に、これを解決することができる。

なお、本発明は、上述の実施形態に限定されず、様々な形態で実施することができる。たとえば、二次元詳細計算部２と三次元シミュレータ３を別々のコンピュータ上に実現してもよい。また、核定数記憶部４は、取り外し可能な記憶媒体であってもよい。二次元詳細計算部２の計算結果を外部出力できるようにしてもよい。表示手段６に、記憶装置を接続して、計算結果を保存できるようにしてもよい。

本発明の一実施形態の二次元詳細計算部の計算の流れ図である。 本発明の一実施形態の原子炉シミュレーション装置のブロック図である。 本発明の一実施形態の二次元詳細計算部で設定した中性子発生点および運動方向の例である。 擬似乱数と準乱数によって多数の二次元の点をプロットした図であって、（ａ）は擬似乱数によるもので、（ｂ）は準乱数によるものである。 数値積分で求めたπの誤差を標本数Ｎの関数として表したグラフである。 衝突を起こさない中性子のグリーン関数の変化の例を示すグラフである。 従来のcharacteristic法でのスキャンの例である。

符号の説明

１…原子炉シミュレーション装置、２…二次元詳細計算部、３…三次元シミュレータ、４…核定数記憶部、５…入力手段、６…表示手段、１１…中性子発生点、１２…中性子移動方向、１３…モンテカルロ法による数値積分の結果、１４…は準モンテカルロ法による数値積分の結果、２１…水、２２…被覆管、２３…ペレット、３１…粒状ガドリニア、４１…磁気ディスク

Claims (10)

1. 粒子輸送の計算方法において、
   低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、
   レイトレースの始点および方向を前記数列で指定する工程と、
   前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、
   前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、
   前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、
   を有することを特徴とする粒子輸送の計算方法。
2. 前記低食い違い数列は、Faure列、変形されたFaure列、Sobol列、変形されたSobol列、Halton列、変形されたHalton列、Niederreiter列および変形されたNiederreiter列のいずれかであることを特徴とする請求項１記載の粒子輸送の計算方法。
3. 原子燃料の核定数の計算方法において、
   低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、
   レイトレースの始点および方向を前記数列で指定する工程と、
   前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、
   前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、
   前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、
   前記粒子の分布から核定数を求める工程と、
   を有することを特徴とする核定数の計算方法。
4. 原子炉シミュレーション方法において、
   低食い違い数列および擬似乱数のいずれかの数列を発生する工程と、
   レイトレースの始点および方向を前記数列で指定する工程と、
   前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する工程と、
   前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する工程と、
   前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する工程と、
   前記粒子の分布から核定数を求める工程と、
   前記核定数に基づいて、前記原子炉の三次元の特性を評価する工程と、
   を有することを特徴とする原子炉シミュレーション方法。
5. 粒子輸送の計算装置において、
   計算条件を入力する入力手段と、
   前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する演算手段と、
   前記粒子の分布を外部に出力する外部出力手段と、
   を有することを特徴とする粒子輸送の計算装置。
6. 原子燃料の核定数の計算装置において、
   計算条件を入力する入力手段と、
   前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価し、前記粒子の分布から核定数を求める演算手段と、
   前記核定数を外部に出力する外部出力手段と、
   を有することを特徴とする核定数の計算装置。
7. 原子炉シミュレーション装置において、
   計算条件を入力する入力手段と、
   前記計算条件に基づき、低食い違い数列および擬似乱数のいずれかの数列を発生し、レイトレースの始点および方向を前記数列で指定し、前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定し、前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成し、前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価し、前記粒子の分布から核定数を求める二次元計算部と、
   前記核定数および前記計算条件に基づいて、前記原子炉の三次元の特性を評価する三次元シミュレータと、
   を有することを特徴とする原子炉シミュレーション装置。
8. 粒子輸送の計算プログラムにおいて、コンピュータに、
   低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、
   レイトレースの始点および方向を前記数列で指定する手順と、
   前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、
   前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、
   前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、
   を実行させるための粒子輸送の計算プログラム。
9. 原子燃料集合体の核定数の計算プログラムにおいて、コンピュータに、
   低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、
   レイトレースの始点および方向を前記数列で指定する手順と、
   前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、
   前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、
   前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、
   前記粒子の分布から核定数を求める手順と、
   を実行させるための核定数の計算プログラム。
10. 原子炉シミュレーションプログラムにおいて、コンピュータに、
    低食い違い数列および擬似乱数のいずれかの数列を発生する手順と、
    レイトレースの始点および方向を前記数列で指定する手順と、
    前記始点および前記方向を用いてレイトレースし、粒子の軌跡を決定する手順と、
    前記軌跡に沿って、粒子の輸送方程式のグリーン関数を作成する手順と、
    前記グリーン関数を用いて粒子輸送ボルツマン方程式を解き、粒子の分布を評価する手順と、
    前記粒子の分布から核定数を求める手順と、
    前記核定数に基づいて、前記原子炉の三次元の特性を評価する手順と、
    を実行させるための原子炉シミュレーションプログラム。

Images