原子炉炉心をモデル化する方法及び対応するコンピュータプログラム製品
【課題】本発明は、より高い収束精度、よりよい計算エラー強さ及びより高い計算効率を示す原子炉モデル化方法を提供する。
【解決手段】原子炉炉心をモデル化するための方法であって、コンピュータで実行される計算のための格子（１２）のノードを構成するために炉心を立方体（１０）に分割するステップと、中性子束を、少なくとも１つの固有系の反復解手順を用いることによって計算するステップとを含み、固有系の反復の成分は、計算されるべきそれぞれの立方体（１０）についての、中性子束、中性子アウトカレント或るいは中性子インカレントのいずれかに対応する。制御パラメータは、中性子固有値μを摂動界面カレント方程式によって影響を与え且つ中性子固有値μを１に向って推進させるように変えられる。
【解決手段】原子炉炉心をモデル化するための方法であって、コンピュータで実行される計算のための格子（１２）のノードを構成するために炉心を立方体（１０）に分割するステップと、中性子束を、少なくとも１つの固有系の反復解手順を用いることによって計算するステップとを含み、固有系の反復の成分は、計算されるべきそれぞれの立方体（１０）についての、中性子束、中性子アウトカレント或るいは中性子インカレントのいずれかに対応する。制御パラメータは、中性子固有値μを摂動界面カレント方程式によって影響を与え且つ中性子固有値μを１に向って推進させるように変えられる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、原子炉の炉心をモデル化する方法、特に炉心内の中性子束を計算する方法に関する。
【背景技術】
【０００２】
そのようなモデル化方法の結果は、原子炉を建設し且つ始動する前に安全解析レポートを準備するために用いることができる。
【０００３】
これらの結果は、また、既存の原子炉に及び特に原子炉に装填された核燃料を管理するために有用である。特に、これらの結果は、いったい炉心設計がどのように発展すべきか又、炉心内の燃料集合体の位置、特に炉心内に導入されるべき新しい燃料集合体の位置をどのように決定すべきかを評価するために用いることが出来る。
【０００４】
そのようなモデル化方法は、コンピュータ計算によって実行される。この目的のために、炉心は、立方体に分割され、各立方体はデジタル計算を実行するための格子のノードを構成する。
【０００５】
通常、そのようなデジタル計算中に解かれるべき定常状態拡散方程式は、
【数１】
（１）
【数２】
（２）
になる。
【０００６】
ここに、λは第１中性子固有値であり、mは、ノードインデックスとも呼ばれる、立方体インデックスであり、Ｇは中性子エネルギー群の数であり、g、g’は中性子エネルギー群インデックスであり、uは立方体のデカルト軸インデックスであり、Σagm は、立方体mとエネルギー群gの巨視的吸収断面積を表し、Σfg'mは、立方体mとエネルギー群g’の巨視的分裂断面積を表し、Σgg'mは立方体mとエネルギー群gの巨視的減速断面積を表し、Φgmは中性子束を表し、ΣRgm，Φgmは対応する核反応（吸収、分裂）の反応率を表し、νは分裂毎に生ずる中性子の数であり、Ｘg は分裂から現れ、中性子エネルギーgでの分裂から出てくる中性子の関数であり、aumはデカルト軸uに沿う立方体mの幅であり、中性子アウトカレントjgul-mと jgur+m、中性子束Φgm及び中性子インカレントjgul+m と jgur-mの間の関係は、次式によって定義される。
【数３】
（３）
【０００７】
i=1, 2, 3で係数 cigumは、立方体mの特徴であり、且つノードの次元と巨視的断面積Σmに基づいている。
【０００８】
図１は、u=x,y及びzについて中性子インカレント jgul+mとjgur-m；u=x,y及びzについて中性子アウトカレントjgul-mと jgur+m；及び中性子束Φgmを示す立方体mの２次元での概略図である。インデックスlは、それぞれのデカルト軸線x,yについて立方体mの各左界面を、又rはそれぞれの各右界面を示す。インデックス+は、それぞれのデカルト軸線x,yについて左から右の向きを、又インデックス-は、それぞれに右から左への向きを表す。
【０００９】
定常状態拡散方程式（１）は、また、ノード展開法方程式に関してＮＥＭ方程式と名づけられる。
【００１０】
従来技術の方法では、計算努力のほとんどは、定常状態拡散方程式（１）に対応する大きな固有系の反復解にささげる部分に集中されている。
【００１１】
これらの計算努力をより少なくし、従って、固有系の解を加速させるために、粗メッシュ再平衡（Coarse Mesh Rebalancing :CMR）手順が用いられる。これらの手順では、与えられた反復のための中性子束及びカレントに、続く計算的に費用のかかる反復を追及する前に、補正因子を乗じる。乗法的補正は、正確な値λexactに近い第１中性子固有値λを有する固有スペクトルの非基本波長部分の存在を抑制するのに役立つ。
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、このやり方で実現される加速効果は、全炉心拡散レベルまで多レベル階層内のもっとも高い粗いメッシュレベルの数値近接に依存する。従って、そのようなＣＭＲ手順は、大変ゆっくりな収束又は収束のよどみすらもたらし、かくして、計算努力を増大させる。
【００１３】
本発明の目標は、より高い収束精度、よりよい計算エラー強さ及びより高い計算効率を提供する原子炉モデル化方法を提供することによって上述の問題を解決することであり、そのため、関係する中性子束計算を、短い計算時間内に且つとてもよい収束精度で得ることができる。
【課題を解決するための手段】
【００１４】
この目的のために、本発明は、請求項１に記載の方法に関する。
【００１５】
特定の実施形態では、本方法は、単独で又は組合せて取られる請求項２乃至１０の特徴を含むのがよい。
【００１６】
本発明は、また、請求項１１に記載のコンピュータープログラム製品に関する。
【００１７】
本発明は、単に例示によって与えられ且つ添付図面を参照して以下の詳細な説明を読むときによりよく理解される。
【図面の簡単な説明】
【００１８】
【図１】モデル化された原子炉炉心内のマクロ断面積、中性子インカレント、中性子束、中性子アウトカレントの間の関係の古典的表示である。
【図２】本発明の１つの側面による、原子炉炉心の分割及び立方体の界面と推進因子の関連を示す概略図である。
【図３】図２の立方体を出る中性子アウトカレントの等方性及び異方性部分の分解概略図である。
【図４】本発明の他の側面による反復解法処理からの多レベルＶサイクルの概略図である。
【図５】図４のボックスＶの詳細な図である。
【図６】本発明の他の側面による原子炉炉心の分割を示す概略図である。
【図７】従来技術方法及び本発明による方法の異なる側面毎に１組の固有系反復解処理の１組の収束曲線である。
【発明を実施するための形態】
【００１９】
以下の説明では、加圧水型原子炉（ＰＷＲ）の場合が考えられるが、本発明を他のタイプの原子炉に適用することも心に留めておくべきである。
【００２０】
本発明によるコンピュータ実行モデル化方法の第１ステップでは、原子炉の炉心が、従来技術の状態におけるように（図２に示す）立方体１０に分割される。各立方体１０は、数値計算がコンピュータにより実行される格子或いはネットワーク１２のノードに対応する。
【００２１】
表示を簡単にするために、実際には、３次元であることを心に留めておくべきである。
【００２２】
（図２の中心の）立方体インデックスmを有する立方体１０の隣りは、それぞれの立方体インデックス m’1(m)、m’2(m)、m’3(m)及びm’4(m)をつけた立方体１０である。以下の詳細な説明では、立方体１０を、立方体のそれぞれの立方体インデックスによって直接示す。
【００２３】
本発明によるコンピュータ実行モデル化方法の第２ステップでは、炉心内の中性子束Φgm は、定常状態拡散方程式（１）に対応する固有系の解によって計算される。この目的のために、反復解手順が用いられる。
【００２４】
除去演算子
、
非散乱演算子
、生成演算子
、インカレント演算子
が方程式（１）から、
【数４】
（４）
によって、定義される。
【００２５】
等方性アウトカレント世代演算子
及び単方向性カレント貫流演算子
は、
【数５】
（５）
によって定義される。
【００２６】
結合演算子
は、隣接する立方体１０についてアウトカレントをインカレントに結合し、
【数６】
（６）
によって定義される。
【００２７】
この結合演算子
は、例えば、立方体mについて、u方向左向きアウトカレントが、デカルト軸uの方向における左隣りのu方向左向きインカレントと等しく、又同様の等式が他の方向及び向きについて確かめられている事実を用いる。例えば、図２に示す隣りm’1(m)、m’2(m)、m’3(m)、m’4(m)が、それぞれ、方程式（6）のデカルト軸yの隣り(m+1)、デカルト軸xの隣り(m-1)、デカルト軸xの隣り(m+1)、デカルト軸yの隣り(m-1)であるとき、それぞれの隣りm’1(m)、m’2(m)、m’3(m)、m’4(m)から立方体mに至る中性子アウトカレントは、立方体mに向かう中性子インカレントである。
【００２８】
上記の演算子を用いると、方程式（１）及び（２）は
【数７】
（７）
と書かれ、
又は解かれるべき固有系である、
【数８】
と書かれる。
【００２９】
方程式（８）では、Φは中性子束縦ベクトルであり、各要素は、それぞれの立方体mに対し且つそれぞれのエネルギー群gに対する中性子束Φgmである（図１）。かくして、中性子束縦ベクトルΦの次元は、（GxM, 1）に等しく、ここに、Gはエネルギー群の数、Mは立方体１０の数である。 j(out) は、中性子アウトカレント縦ベクトルであり、各要素は、それぞれの立方体m、エネルギー群g及びuがx, y或いはzに等しいデカルト軸uについての、それぞれの中性子アウトカレント jgul-m, jgur+m,である。 j(in)は、中性子インカレント縦ベクトルであり、各要素は、それぞれの立方体m、エネルギー群g、デカルト軸uについての、それぞれの中性子インカレント jgul+m, jgur-mである。かくして、中性子アウトカレント縦ベクトル j(out)及び中性子インカレント縦ベクトルj(in)の次元は（6xGxM, 1）に等しい。以下の説明では、中性子アウトカレント縦ベクトル j(out) は、また、 joutと記載される。
【００３０】
かくして、式（８）によって定義される固有系の反復 (Φ, j(out), j(in))の成分は、立方体m毎に及びエネルギー群g毎に計算されるべき中性子束Φgm、中性子アウトカレント jgul-m, jgur+m、及び中性子インカレント jgul+m, jgur-mに対応し、これは、中性子束縦ベクトルΦ、中性子アウトカレント縦ベクトルj(out)及び中性子インカレント縦ベクトル j(in)のそれぞれの要素である。
【００３１】
反復解処理は、固有系を予備の固有系に調節するサブステップを含み、予備の固有系の反復の成分(j(out), j(in))は、中性子アウトカレント縦ベクトルj(out)及び中性子インカレント縦ベクトル j(in)のそれぞれの要素である、それぞれの立方体mに至る中性子インカレント jgul+m, jgur-m、或いはそれぞれの立方体mから来る中性子アウトカレント jgul-m, jgur+mの何れかにのみ対応する。
【００３２】
予備の固有系への固有系の調整は
【数９】
と書かれる方程式（８）の第１部分を変形することから始まる。
【００３３】
シンボリック演算子
を定義すると、中性子束縦ベクトルΦは、中性子インカレント縦ベクトル j(in) の関数で
【数１０】
（１０）
と表される。
【００３４】
この式を式７のアウトカレント平衡部分に、すなわち方程式８の第２部分に代入することによって、カレントのみの関係が得られる。
【数１１】
（１１）
【００３５】
次の演算子記号を定義した後に
【数１２】
（１２）
及び
【数１３】
（１３）
【００３６】
カレントのみの関係は、スペア固有系である、
【数１４】
（１４）
と書かれる。
【００３７】
中性子アウトカレント縦ベクトルj(out)と中性子インカレント縦ベクトルj(in)の間の関係
（ここに
は結合演算子である。）を用いると、スペア固有系が、スペア固有系の反復の成分(j(out)) が、中性子アウトカレント縦ベクトル j(out)の要素である、それぞれの立方体mから来る中性子アウトカレントjgul-m, jgur+m にのみ対応する。このスペア固有系は、
【数１５】
に対応する。
【００３８】
演算子
は、最初はちょうど１に等しくなく、 j(out)が正確な解 j(out)exactに収束するときのみ１に等しくなるので、もし同時に第１中性子固有値λが正確な値λexactに収束するならば、方程式１５は
【数１６】
（１６）
になる。
【００３９】
もし第１中性子固有値λが正確な値λexactに収束するならば、第２中性子固有値μは、１に向かって収束する。
【００４０】
数値最適化のために、方程式１６は、演算子Π-1を用いて前からの乗法により変形され、以下の方程式が得られる。
【数１７】
（１７）
【００４１】
項は、デカルト方向uに依存しない等方性の項である。
【００４２】
立方体インデックスm、デカルト軸u、u’、デカルト軸u、u’に沿うそれぞれの向きs、s’の方程式１７からそれぞれの項の式が、
〔数１８〕〔数１９〕〔数２０〕
（１８）（１９）（２０）
によって与えられる。
【００４３】
G、G’は、エネルギー群の組である。
と
は、それぞれの因子
と
を有する単一エネルギー演算子である。
は、スペクトル演算子
の因子である。
【００４４】
各組み合わせそれぞれ {u, s}、 {u’, s’}は、u、u’がx、y或いはzに等しく、s、s’がl或いはrに等しい立方体mの界面を定義することに気づくべきである。例えば、{x, l} は、デカルト軸xに沿う立方体mの左界面を定義する。
【００４５】
個々の項は
【数２１】
（２１）
【数２２】
【数２３】
（２３）
によって与えられる。
【００４６】
次に、方程式１７は、例えば、普通の反復解手順をガウス-ザイデル処理として用いることによって解かれる。
【００４７】
かくして、解は中性子アウトカレント縦ベクトルj(out)の要素 j(out)mgusについて計算される。これは、式１４による中性子インカレント縦ベクトル j(in)の解を決定でき、且つ最後に式１０による中性子束縦ベクトルΦの解を決定できる。
【００４８】
本モデル化方法によって得られた、計算された中性子束縦ベクトルΦは、例えば、核燃料を管理するために既存の原子炉炉心を制御するために用いられ、或いは新しい炉心を建設するために用いることができる。
【００４９】
モデル化方法は並列プロセッサで実行されてもよいし或いは単一プロセッサで実行されてもよい。上述のモデル化方法は、よりよい収束精度、よりよいコンピュータエラー強さ及びよりよい計算効率をもたらすことが証明されている。これは、その反復の成分のみが中性子アウトカレント j(out)に対応し、中性子束に依存しない、予備固有系の解に関連する。
【００５０】
さらに、本発明によるモデル化方法の収束精度は、1E-12まで改善され、これに対して、古典的なモデル化方法で得られる収束精度は1E-6に制限される。
【００５１】
他の実施形態では、解かれるべきスペア固有系の反復の成分は、中性子インカレント j(in)にのみ対応する。
【００５２】
計算エラー強さ及び計算効率をさらに改善するために、本発明の第２側面によれば、固有系は、先ず、いくつかの中性子エネルギー群の選択のための、固有系に対応する制限固有系で調整される。この選択は、また、いくつかのエネルギー群へのスペクトル制限と呼ばれる。次いで、制限固有系は、本発明の第１側面により解かれてもよい。
最後に、制限固有系の解は、固有系を解くために用いられる。
【００５３】
選択されたエネルギー群の数ＮＧＣは、エネルギー群の全数ngよりも小さい。ＮＧＣは、より少数の粗エネルギー群に併合される細かいエネルギー群の集合である
粗くされたスペクトル帯の数であってもよい。
【００５４】
ngは、例えば、８に等しく、ＮＧＣは、例えば、４、３、２又は１に等しくてもよい。
【００５５】
この第２の側面によると、それぞれの項が方程式１８乃至２０によって与えられて、推進因子
の方程式１７は、
【数２４】
（２４）
と表される。
ここに、G、G’は粗い中性子エネルギー群インデックスである。
【００５６】
等方性項
はデカルト軸u及び向きsによって特定される出て行くカレント方向と無関係であり、非等方性項とも呼ばれる貫流項
は、同じデカルト軸uに沿う同じ粗中性子エネルギー群Ｇ内で定義される。
【００５７】
立方体mを出るアウトカレントの等方性項及び非等方性項の分解を、図３に示す。
【００５８】
等方性項は全ての方向に同じであり、非等方性項は、立方体mの界面の間で変化する。
【００５９】
最初に、等方性項が容易に計算されるアウトカレントの分解とともに、いくつかのエネルギー群へのスペクトル制限について固有系を解くことは、よりよい計算エラー強さ及びよりよい計算効率をもたらす。これは、本発明の第２の側面によるスペクトル制限から生じる固有系の次元の減少に関連する。
【００６０】
エラー強さ及び計算効率をさらに改善するために、本発明の第３の側面によれば、複数のエネルギー群のための反復解処理は、図４に示されるように、頂レベル２１、第１中間レベル２２、４つの第２中間レベル２４、２６、２８、３０及び３つの底レベル３２、３４、３６を含む、多レベルＶサイクル２０の形態であるのがよい。
【００６１】
Ｖサイクルの頂レベル２１は、複数の中性子エネルギー群について、定常状態拡散方程式１、すなわちＮＥＭ方程式、に対応する固有系のための反復計算と、本発明の第２の側面によるいくつかのエネルギー群へのスペクトル制限のために制限固有系への固有系の調整を含む。頂レベル２１から生じる制限固有系は、次いで、頂レベル２１のすぐ下の第１中間レベル２２に送り込まれる。
【００６２】
第１中間レベル２２は、制限固有系を本発明の第１の側面によるスペア固有系に調整することを含み、スペア固有系の反復(j(out), j(in))の成分は、立方体mに入る中性子インカレントjgul+m, jgur-m或いは立方体mから去る、j(out)mgusとも記述される中性子アウトカレント jgul-m, jgur+m、のいずれかにのみ対応する。方程式１８乃至２０によって与えられる演算子
,
及び
の因子は、第１中間レベル２２のために陽関数表示で計算される。
【００６３】
各第２中間レベル２４、２６、２８、３０は、中性子エネルギー群の前の選択についての前のスペア固有系の、中性子エネルギー群の後の選択についての後のスペア固有系への調整を含み、後の選択での中性子エネルギー群の数は、前の選択での中性子エネルギー群の数より少ない。多レベルＶサイクル２０の下方向きで、第１中間レベル２２に続く第２中間レベル２４の前のスペア固有系は、第１中間レベル２２から生じるスペア固有系である。第１中間レベル２２に続かない各第２中間レベル２６、２８、３０の前のスペア固有系は、先行する第２中間レベル２４、２６、２８から生じる後のスペア固有系に対応する。言い換えれば、中性子エネルギー群の数は、第２中間レベルから次の第２中間レベルに下向きに減少する。最後の第２中間レベル３０では、中性子エネルギー群の数は、１に等しい。第２中間レベル２４、２６、２８、３０の中性子エネルギー群の数のこの制限は、後述するスペクトル圧縮代数によって、行われ、本発明の第２の側面によらなくてもよいことに気づくべきである。
【００６４】
底レベル３２、３４、３６は、Ｖサイクルの下方向きでは、最後の第２中間レベル３０で決定された単エネルギー群についての最後のスペア固有系の、従来技術手順の状態による解を含む。底レベル３４、３６は、立方体が先行する底レベル３２、３４の立方体よりも大きいそれぞれの粗い格子の、先行する底レベル３２、３４から生じる固有系の空間制限に対応する。
【００６５】
Ｖサイクルの底レベル３６では、単一エネルギー群についての固有系の解がコンピュータで計算される。
【００６６】
次いで、この解は、上方向きのＶサイクルの上レベルに再び導入される、そのため、それぞれのスペア固有系のために解が計算される。
【００６７】
上方向きのＶサイクルの頂レベル２１では、複数のエネルギー群についてＮＥＭ方程式の解が計算される。
【００６８】
図４の例示実施形態では、頂レベル２１でのＮＥＭ方程式について中性子エネルギー群の数ngは８に等しい。中間レベル２２、２４、２６、２８、３０について、それぞれの選択での中性子エネルギー群の数は、４、３、２及び１にそれぞれに等しい。図４に示すＶサイクルの３つの底レベル３２、３４、３６は、ＣＭＲ手順により１つのエネルギー群のための最後のスペア固有系の解を含む。
【００６９】
それぞれのスペア固有系の解に対応する４つの中間レベル２２、２４、２６、２８は、また、それぞれの参照記号SR4、SR3、SR2及びSR1（図５）によって示す。
【００７０】
ＳＲ４レベルの単一エネルギー演算子
の因子を陽関数表示的に計算する。引き続いて、ＳＲ３レベルの単一エネルギー演算子
の因子を、図５に示すように、レベルＳＲ４の単一エネルギー演算子の因子から導き出す、ここに各ボックス３７は、中性子エネルギー群を
【数２５】
（２５）
により表す。
【００７１】
式２５を、レベル２２と２４の間の矢印３８によって図５に示す。
【００７２】
引き続いて、ＳＲ２レベルの単一エネルギー演算子
の因子を、
【数２６】
（２６）
によりＳＲ３レベルの単一エネルギー演算子の因子から導かれる。
【００７３】
式２６を、レベル２４と２６の間の矢印４０によって図５に示す。最後に、ＳＲ１レベルの単一エネルギー演算子
の因子を、
【数２７】
（２７）
によりＳＲ２レベルの単一エネルギー演算子の因子から導き出す。
【００７４】
式２７を、レベル２６と２８の間の矢印４２によって図５に示す。
【００７５】
方程式２５乃至２７は、図５に示すスペクトル格子分割戦略に対応する単一エネルギー演算子
のスペクトル圧縮代数を定義する。異なる式にあっては単エネルギー演算子
の他のスペクトル凝縮代数を、他のスペクトル格子分割戦略を用いて、すなわちボックス３７の異なる配置を用いて定義してもよいことに気づくべきである。
【００７６】
単一エネルギー演算子
の因子を、下方向きにＳＲ４レベルからＳＲ１レベルまで単一エネルギー演算子
の因子について同様の方法で導き出す。ＳＲ４レベルのスペクトル演算子
の因子を、明確に計算する。引き続いて、ＳＲ３レベルのスペクトル演算子
の因子を、
【数２８】
（２８）
によってＳＲ４レベルのスペクトル演算子
の因子から導き出す。
【００７７】
引き続いて、ＳＲ２レベルのスペクトル演算子
の因子を、
【数２９】
（２９）
によってＳＲ３レベルのスペクトル演算子の因子から導き出す。
【００７８】
最後に、ＳＲ１レベルのスペクトル演算子
の因子を、
【数３０】
（３０）
によってＳＲ２レベルのスペクトル演算子の因子から導き出す。
【００７９】
方程式２８乃至３０は、単一エネルギー演算子
及び
の、図示しないスペクトル格子分割戦略に対応し、図５に示すスペクトル格子分割戦略と同様のスペクトル演算子
のスペクトル凝縮代数を定義する。スペクトル演算子
の他のスペクトル凝縮代数は、他のスペクトル格子分割戦略を用いて定義してもよいことを気づくべきである。
【００８０】
かくして、上方向きでのＶサイクルの頂レベル２１で中性子アウトカレント縦ベクトル j(out) の要素 j(out)mgus の解を計算する。この計算の解により、式１４による中性子インカレント縦ベクトル j(in) の解を決定することができ、最後に式１０による中性子束縦ベクトルΦの解を決定することができる。
【００８１】
次いで、原子炉炉心を、計算された中性子束縦ベクトルΦに基づいて建設し或いは運転する。
【００８２】
それぞれのＳＲ３、ＳＲ２及びＳＲ１レベルでの演算子
,
及び
について、複雑な算術演算を伴わない、スペクトル圧縮代数を用いることによって、演算子
,
及び
をとても安価に計算でき、かくして、モデル化方法のよりよい計算エラー強さ及びよりよい計算効率に至る。
【００８３】
本発明の第４の側面により、エラー強さ及び計算効率をさらに改善するために、立方体１０を第１カテゴリーと第２カテゴリーに分割する。
【００８４】
以下の説明では、図６に示すように、第１カテゴリーの立方体１０Ｒは、赤立方体と呼ばれ、第２カテゴリーの立方体１０Ｂは黒立方体と呼ばれるが、特別の限定的な意味を用語”黒”及び”赤”と関連させるべきではない。各赤立方体１０Ｒは、すぐ隣の立方体として黒立方体１０Ｂのみを有する。かくして、赤立方体１０Ｒのほとんどは、６つのすぐ隣りの黒立方体１０Ｂを有する。これは、立方体が考慮される立方体と共通の面を共有している、”すぐ”隣りによって理解されるべきである。
【００８５】
その結果、図６により例示されるように、格子１２は２次元表現の市松模様の暗い領域と明るい領域に関して視覚的に類似している。
【００８６】
次いで、立方体には例えば赤立方体１０Ｒによって始め、黒立方体１０Ｂによって終わる番号が付けられる。
【００８７】
以下の説明では、２つのカテゴリーでの立方体のこのような分割及びカテゴリーの次々の番号付けは、赤−黒順序付けと称される。
【００８８】
当該技術の辞書式順序付けの状態と比較して立方体の赤−黒順序付けの利点は、赤立方体１０Ｒについて、そのすぐ隣りの立方体が全て黒であって、逆の場合も同じである。
【００８９】
コンピュータが炉心内の中性子束を計算するために解かねばならない式２４は行列ベクトル形式で式数３１と書くことができ、
【数３１】
（３１）
【００９０】
典型的には０．９乃至０．９５の値からなり、そのため、積 ｙμがシフトを形成する、可変のパラメータｙを用いると、式３１の反復は次のシフト逆数陰関数形式で書かされる。
【数３２】
（３２）
【００９１】
ここに、s(n-1)は
【数３３】
（３３）
によって与えられる源である。
【００９２】
スペクトル演算子
及び単一エネルギー演算子
は希薄であり、赤立方体１０Ｒをすぐ隣りの黒立方体１０Ｂにのみ結合し、又その逆に結合する。かくして、赤黒順序付けは、赤の部分と黒の部分との間に次の便利な関係を可能にする。
【数３４】
（３４）
【００９３】
方程式３４が、解の赤の部分のみについて次の等式に変換される。
【数３５】
【００９４】
計算の観点から、赤黒順序付けの使用は、反復解手順中、もし、反復の最初の半分で、反復 d の、赤成分 dred が更新されるならば、次いで、反復のもう一つの部分中、全ての黒成分 dblackは、それらの隣りの赤成分dred に基づいて更新されることを意味する。換言すれば、各黒立方体１０Ｂの値は、全てのすぐ隣りの赤立方体１０Ｒの値に基づいて計算される。
【００９５】
赤黒順序付けを用いることによって、本発明の第１、第２又は第３の側面で記載された演算子
,
及び
の計算を改善する。かくして、赤黒順序付けは、計算効率を改善するために本発明の第１、第２或いは第３の側面を補完に用いられてもよい。
【００９６】
そのような赤黒順序付けは、本発明の第５側面を構成する特定の解手順と一緒に用いられるとき特に有用であることがわかる。
【００９７】
この手順は、特に高いエラー強さ安定化双共役勾配安定化（Bi-CGStab）手順である。この手順の適当な前置きの記載は、次の引用文献に見られる。
Y. Saad, "スパース線形システム反復法 (Iterative Methods for Sparse Linear Systems)", 第２版, 工業・応用数学会 (SIAM) (2003;)、
H.A. van der Vorst, "Bi-CGSTAB: 非対称の線形システムの解のためのBi-CGの高速且つ円滑な収束バリアント (a Fast and Smoothly Converging Variant of Bi-CG for the solution of nonsymmetric linear systems)", 日本応用数理学会 SIAM J. Sci. Stat. Comput. 13(2), pp.631-644 (1992)。
【００９８】
双共役勾配安定化（Bi-CGStab）手順は、与えられたΘとsにより、dの関数のための最小化原理から導かれ、そのために定常を方程式３５によって与えられる線形系をちょうど満たす特定の最適なdを中心とした小さな変分δdに関して適用し、且つ汎関数が最小値を仮定する。
【００９９】
かくして、式３５の効率的な解き方について、より直接的な考察によってではなく汎関数の最小化によって動かされる解手順を定義することが可能である。Bi-CGStab手順は、そのような最小原理から続き、反復での連続変化は、（Galerkin加重積分剰余のようである）関数の各変化が先の変化の全てに関して直交するように配列される。
【０１００】
本発明の第５側面による特殊なBi-CGStab手順が、解ベクトルd、解誤差 r ( r = s - Θ d)及び補助ベクトル r*, s及び p、で及び初期解推定値d0で下に与えられる。
【数３６】
（３６）
【０１０１】
上記Bi-CGStab数列は、普通Ｎ＜３のＮステップ後に打ち切られる。
【０１０２】
Bi-CGStab手順のために、予備調整が式３２のシフト逆陰関数形式に基づいており、ただし、典型的に0.9<ｙ<0.95である。
【０１０３】
シフトｙμのためのこの選択で、解処理中、ｙμが、系の最も小さい可能な固有値のｙ倍まで収束するので、演算子sは、非特異のままにするように保証される。
【０１０４】
方程式３２を解くことは、高等予備調整Krylov手順を展開するための完全なシナリオである、中性子源を有するわずかに臨界未満系の中性子束分布を決定することに数値的に等しい。
【０１０５】
赤黒順序付けを用いると、方程式３２は、上で説明したように、方程式３５に変換される。方程式３５によって与えられる予備調整系は、本発明の第６の側面を構成する。
【０１０６】
いったん方程式３５が、Bi-CGStab手順を介して解かれると、解黒部分が方程式３４の黒部分から計算される。
【０１０７】
格子上に投影されるように、隣りの（赤対黒）ノードの間の直接の中性子相互作用を意味する、この特定の予備調整法の使用は、解かれるべき系に予め含まれ、
【数３７】
（３７）
は、
【数３８】
（３８）
により与えられるゆえに、方程式３５の単位演算子をより支配的に形成する。
【０１０８】
この予備調整法は、空間結合、それ故に式の解を決定するノードの間の物理的な相互作用に関して、重要な情報、又は、その情報の大部分を、低計算コストで、うまく予め含める。
【０１０９】
他の実施形態では、反復解手順は、一般化した最小剰余法、また短縮GMRES、或いは平面ヤコビ法のようなザイデル処理、或いはクリロフ部分空間での反復手順である。
【０１１０】
本発明の第７側面によれば、改善された固有値計算のための変分の原理を次の説明に記載する。
【０１１１】
本発明の第１、第２或いは第３の側面について上述されたスペア固有系に対応する式（１６）は、第２中性子固有値μ及び中性子アウトカレント縦ベクトルjoutによって、次の形に書かされる。
【数３９】
（３９）
【０１１２】
ここに、c1は、立方体の特性により演算子
のための他の記号である。
【０１１３】
何が続くかにおいて、我々は、第１中性子固有値λがスカラー制御パラメータの役割を満たすと仮定する。もし全体の解が収束し、かくして安定化するならば、第２中性子固有値μが、１に近づく特性を用いると、変分原理又は摂動アプローチは、第１中性子固有値λを変化させて、第２中性子固有値μを１に向かわせ且つスペア固有系解の収束を加速することを含む。
【０１１４】
第１中性子固有値λのこの変分は、スペア固有系の各反復前に、すなわちスペア固有系の各第２、第３、第４或いは第５の反復前に、なされるのがよい。第１中性子固有値λの変分は、λがλ+δλにされるように偏差δλのλへの適用であり、これは摂動を起こした界面カレント方程式、
【数４０】
（４０）
ただし、演算子偏微分
【数４１】
（４１）
により第２中性子固有値μに影響を与える。
【０１１５】
摂動
が反復手順中ますます小さくなると期待されるので、摂動は方程式（４０）では都合よく無視され、任意の付加した重み関数
での積分は次式を得る。
【数４２】
（４２）
【０１１６】
第１中性子固有値λの変分が、第２中性子固有値μが１に向かって推進され、すなわちμ+δμが１に向かって推進されるようなものであるとき、近似
が置かれる。
この近似値及び方程式（４２）で、次の方程式が得られる。
【数４３】
（４３）
【０１１７】
等式
を用いると、δλが
【数４４】
（４４）
で与えられる。
【０１１８】
ここに、rout は、第２中性子固有値μのために置かれた目標値が１に等しいとき、界面アウトカレント平衡方程式、すなわち方程式（３９）の剰余を表す。剰余rout は、いかなる場合の各ＮＥＭ反復ステップ毎に計算される必要がある。
【０１１９】
この側面でスカラー制御パラメータと考慮される、第１中性子固有値λは、次いで、次式のように更新される。
【数４５】
（４５）
【０１２０】
ここに、
及び
は、それぞれ、第１中性子固有値λの先の値及び更新された値である。
【０１２１】
部分的な演算子微分係数
を決定し且つ式７を使うために、例えば中性子束Φは、
【数４６】
（４６）
と書かれる。
上式から：
【数４７】
が続く。
【０１２２】
次いで、近似導関数
は、
のような非散乱演算子
に対する行列の上対角部分を無視することによって、すなわち上方散乱、
【数４８】
（４８）
を無視することによって計算される。
【０１２３】
多くの計算の場合に、これは近似でもないが、もし上方散乱がモデル化されなければ、数学的に正確であることに気づくべきである。
【０１２４】
変分原理の好結果の適用にとって十分である
【数４９】
（４９）
を仮定する。
【０１２５】
を計算するために、それを先ず、
【数５０】
（５０）
と書き、
上式から：
【数５１】
（５１）
が続く。
【０１２６】
次いで、
【数５２】
（５２）
で演算子記号
を定義し、
を陰関数的に計算するために、テーラー式を
に適用して次式が得られる。
【数５３】
（５３）
【０１２７】
テーラー展開が、
【数５４】
（５４）
と並べ換える。
【０１２８】
方程式（５４）への連鎖律の適用により、導関数
が、
【数５５】
（５５）
によって与えられ、
これは
【数５６】
（５６）
を得る。
【０１２９】
導関数
を決定するのに効率的な反復スキームを可能にする最後の式は、
で方程式（５６）の前もっての乗法によって得られ、これは
【数５７】
（５７）
を与える。
【０１３０】
と書くことによって、
の計算が
【数５８】
（５８）
から得られる。
源の項sは、
【数５９】
（５９）
によって定義される。
【０１３１】
２つのエネルギー群の場合には、方程式（５８）は直接的に及び解析的に解かれる。
【０１３２】
２つ以上のエネルギー群の場合には、方程式（５８）は、
【数６０】
（６０）
の適用により、正確さが早期反復段階で要求されることなく、反復的に解かれる。
ここに、
及び
は、それぞれ、近似導関数
の以前の及び更新された値である。
【０１３３】
早期のとは、反復ステップの全必要数の最初の半分までの約３分の１を意味する。
【０１３４】
この変分アプローチは、有用である
の近似値のみを必要とするので、
の計算からの反復を、少し早く打ち切れることができることに気づくべきである。
【０１３５】
かくして、本発明のこの第７側面による変分原理は、方程式（４５）の分子
が収束のときますます小さくなり、分母
は適当な値にむしろ急速に収束するので、大変効果的である。この第７の側面では、分子は、それがゼロに向かって収束するので、推進原理を定義する。
【０１３６】
方程式（４５）による第１中性子固有値λを変えることは、第２中性子固有値μを１に向かって推進し、且つスペア固有系解の収束を加速させる。第１中性子固有値λのこの変分は、スペア固有系の各反復の前に、すなわちスペア固有系の各第２、第３、第４又は第５反復の前に行われてもよい。
【０１３７】
かくして、第１中性子固有値λの変分原理は、計算効率を改善するために本発明の第１、第２或いは第３の側面の補完に用いられてもよい。
【０１３８】
第７の側面の全ての上述した結果を、方法における計算効率を失うことなしに、格子１２の黒部分に減少させることができ、そのために、本発明の第４の側面で上述した赤黒順序付けがこの第７の側面の補完に用いられる。
【０１３９】
本発明の第８の側面により、固有値計算のための上述の変分原理は、典型的に、方程式（４）によって与えられる除去演算子
を変形させ、
と記される制御パラメータの計算のために一般化でき、変形除去演算子は
と記される。
【０１４０】
全ての位置での変形の例では、制御パラメータ
とは、硼素濃度である。選択された位置での変形の例では、制御パラメータ
は、一群の選択された制御棒の棒挿入深さを制御するためのパラメータである。
【０１４１】
と仮定する。
【０１４２】
すると、方程式（１６）は、第２中性子固有値μを用いて、次の形に書かれる：
【数６１】
（６１）
シンボリックスペクトル演算子
は
【数６２】
（６２）
によって定義される。
【０１４３】
もし制御パラメータ
との補正値を含む全体の解が収束し、かくして安定するならば、第２中性子固有値μが、１つに近づく性質を用いると、変分原理、又は摂動アプローチは、第２中性子固有値μを１に向かって推進させるためにこと制御パラメータ
を変化させることを含む。
【０１４４】
制御パラメータ
の変分は、
が
に推進させるように変分
の
への適用であり、これは、摂動界面カレント方程式
【数６３】
（６３）
部分演算子導関数：
【数６４】
（６４）
により第２中性子固有値μに影響を与える。
【０１４５】
摂動
が反復手順中だんだん小さくなると期待されるので、これらは方程式（６３）で都合よく無視され、任意の随伴重み関数
での積分は次の式を得る。
【数６５】
（６５）
【０１４６】
制御パラメータ
の変分で、第２中性子固有値μが１に向かって推進され、すなわち、そのμ+δμは１に向かって推進されるので、近似
が置かれる。この近似及び式（６５）を用いて、次の方程式が得られる。
【数６６】
（６６）
【０１４７】
等式
を用いると、
は、
【数６７】
（６７）
によって与えられる。
ここに、routは、界面アウトカレント平衡方程式、すなわち、方程式（３９）の残差を表し、第２中性子固有値μのために課された目標値は１に等しい。
【０１４８】
制御パラメータ
は次のように更新される：
【数６８】
（６８）
ここで
と
は、それぞれ、制御パラメータ
の以前の値及び更新された値である。
【０１４９】
変分原理の制御パラメータ
へのうまい適用にとって十分である、
が仮定される。
【０１５０】
を計算するために、テーラー公式の適用による上述の調整が、また、同じ方法で用いられる。
【０１５１】
本発明のこれら第８の側面は、本発明の上述の第７の側面と同じ利点を示し、計算効率を改善するために、本発明の上述の第７側面のためのものと同様の方法で、本発明の第１、第２或いは第３の側面の補完に用いられてもよい。
【０１５２】
図７は、追求したＮＥＭ反復ステップの数を横座標に、高速中性子エネルギー群の中性子束の最大の解の誤差である、高速中性子束誤差を縦座標にした、異なる固有系反復解手順の１組の収束曲線を表す。全てのエネルギー群がスペクトル的に結合されるので、他のエネルギー群の中性子束の最大の解の誤差は、高速中性子束誤差と全く同様である。かくして、異なる固有系反復解手順の間の比較がこれらの収束曲線によって可能である。
【０１５３】
曲線５０は、単一精度を用いた古典的な粗メッシュ再平衡（ＣＭＲ）手順用の収束曲線であり、曲線５１は、二重精度を用いた古典的なＣＭＲ手順の収束曲線である。単一及び二重の精度は、コンピュータプログラム製品における浮動小数点数を表すために用いられる古典的な精度である。
【０１５４】
曲線５２は、単一の精度を用いた単一エネルギー群のための、本発明の第１及び第７の側面による反復手順の収束曲線であり、曲線５３は、二重精度を用いた、本発明の同じ側面による反復手順の収束曲線である。
【０１５５】
曲線５４は、単一精度手順を用いた４つの粗エネルギー群のための、本発明の第１及び第２の側面による反復手順用の収束曲線であり、曲線５５は、二重精度を用いた、本発明の同じ側面による反復手順用の収束曲線である。
【０１５６】
曲線５６は、単一精度を用いた４つの粗エネルギー群のための、本発明の第１、第２及び第７の側面による反復手順用の収束曲線であり、曲線５７は、二重精度を用いた、本発明の同じ側面による反復手順用の収束曲線である。
【０１５７】
二重精度を用いた４つの粗エネルギー群のための、第１、第２及び第７の側面による手順（曲線５７）は、曲線５７について最大である収束曲線の勾配によって示されているように、５０回の追求ＮＥＭ反復ステップについての値が実質的に１０-13に等しい最も小さい誤差ともっとも良い計算効率の両方を有する手順である。
【０１５８】
二重精度を用いた単一エネルギー群のための、第１及び第７の側面による手順（曲線５３）は、１３０回の追求値がＮＥＭ反復ステップについて実質的に１０-9に等しい小さい誤差を出す。さらに、この手順（曲線５３）の計算効率は、曲線５５の勾配が曲線５３の勾配よりも大きいので、二重精度を用いた４つの粗エネルギー群のための、第１及び第２の側面による手順の計算効率（曲線５５）よりも価値が高い。
【０１５９】
二重精度を用いた４つの粗エネルギー群のための、第１及び第２の側面による手順は、４０回の追求ＮＥＭ反復ステップについて実質的に１０-7に等しい誤差を出す。
【０１６０】
かくして、二重精度は、高速中性子束誤差の値に大きな影響を及ぼすが、モデル化方法の計算効率に影響を与えることがない。
【０１６１】
本発明の第７側面は、また、計算効率に影響を与えるよりも、高速中性子束誤差の値に大きな影響を及ぼす。
【０１６２】
４つの粗エネルギー群のための、第１及び第２の側面は、最良の計算効率を提供し、図７は、100と10-5の間に含まれる誤差値に関する最も大きい勾配が、第１及び第２の側面に対応する曲線５４乃至５７の全てに関して得られることを示す。
【０１６３】
前述に開示されるように、本発明の第１乃至第８の側面は、よりよい計算効率を提供するエラー強さモデル化方法を達成するのを助け、そのため、関連する炉心シミュレーションを短い計算時間内に得ることができる。
【０１６４】
第１側面は、それ自体で、この目的を達成する上で役に立ち、かくして第２側面から第８側面まで別々に用いることができることを心に留めておくべきである。また、第２、第３、第７及び第８の側面は、第１の側面で或いは第４乃至第６側面のいずれかで必ずしも実行されない。
【技術分野】
【０００１】
本発明は、原子炉の炉心をモデル化する方法、特に炉心内の中性子束を計算する方法に関する。
【背景技術】
【０００２】
そのようなモデル化方法の結果は、原子炉を建設し且つ始動する前に安全解析レポートを準備するために用いることができる。
【０００３】
これらの結果は、また、既存の原子炉に及び特に原子炉に装填された核燃料を管理するために有用である。特に、これらの結果は、いったい炉心設計がどのように発展すべきか又、炉心内の燃料集合体の位置、特に炉心内に導入されるべき新しい燃料集合体の位置をどのように決定すべきかを評価するために用いることが出来る。
【０００４】
そのようなモデル化方法は、コンピュータ計算によって実行される。この目的のために、炉心は、立方体に分割され、各立方体はデジタル計算を実行するための格子のノードを構成する。
【０００５】
通常、そのようなデジタル計算中に解かれるべき定常状態拡散方程式は、
【数１】
（１）
【数２】
（２）
になる。
【０００６】
ここに、λは第１中性子固有値であり、mは、ノードインデックスとも呼ばれる、立方体インデックスであり、Ｇは中性子エネルギー群の数であり、g、g’は中性子エネルギー群インデックスであり、uは立方体のデカルト軸インデックスであり、Σagm は、立方体mとエネルギー群gの巨視的吸収断面積を表し、Σfg'mは、立方体mとエネルギー群g’の巨視的分裂断面積を表し、Σgg'mは立方体mとエネルギー群gの巨視的減速断面積を表し、Φgmは中性子束を表し、ΣRgm，Φgmは対応する核反応（吸収、分裂）の反応率を表し、νは分裂毎に生ずる中性子の数であり、Ｘg は分裂から現れ、中性子エネルギーgでの分裂から出てくる中性子の関数であり、aumはデカルト軸uに沿う立方体mの幅であり、中性子アウトカレントjgul-mと jgur+m、中性子束Φgm及び中性子インカレントjgul+m と jgur-mの間の関係は、次式によって定義される。
【数３】
（３）
【０００７】
i=1, 2, 3で係数 cigumは、立方体mの特徴であり、且つノードの次元と巨視的断面積Σmに基づいている。
【０００８】
図１は、u=x,y及びzについて中性子インカレント jgul+mとjgur-m；u=x,y及びzについて中性子アウトカレントjgul-mと jgur+m；及び中性子束Φgmを示す立方体mの２次元での概略図である。インデックスlは、それぞれのデカルト軸線x,yについて立方体mの各左界面を、又rはそれぞれの各右界面を示す。インデックス+は、それぞれのデカルト軸線x,yについて左から右の向きを、又インデックス-は、それぞれに右から左への向きを表す。
【０００９】
定常状態拡散方程式（１）は、また、ノード展開法方程式に関してＮＥＭ方程式と名づけられる。
【００１０】
従来技術の方法では、計算努力のほとんどは、定常状態拡散方程式（１）に対応する大きな固有系の反復解にささげる部分に集中されている。
【００１１】
これらの計算努力をより少なくし、従って、固有系の解を加速させるために、粗メッシュ再平衡（Coarse Mesh Rebalancing :CMR）手順が用いられる。これらの手順では、与えられた反復のための中性子束及びカレントに、続く計算的に費用のかかる反復を追及する前に、補正因子を乗じる。乗法的補正は、正確な値λexactに近い第１中性子固有値λを有する固有スペクトルの非基本波長部分の存在を抑制するのに役立つ。
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、このやり方で実現される加速効果は、全炉心拡散レベルまで多レベル階層内のもっとも高い粗いメッシュレベルの数値近接に依存する。従って、そのようなＣＭＲ手順は、大変ゆっくりな収束又は収束のよどみすらもたらし、かくして、計算努力を増大させる。
【００１３】
本発明の目標は、より高い収束精度、よりよい計算エラー強さ及びより高い計算効率を提供する原子炉モデル化方法を提供することによって上述の問題を解決することであり、そのため、関係する中性子束計算を、短い計算時間内に且つとてもよい収束精度で得ることができる。
【課題を解決するための手段】
【００１４】
この目的のために、本発明は、請求項１に記載の方法に関する。
【００１５】
特定の実施形態では、本方法は、単独で又は組合せて取られる請求項２乃至１０の特徴を含むのがよい。
【００１６】
本発明は、また、請求項１１に記載のコンピュータープログラム製品に関する。
【００１７】
本発明は、単に例示によって与えられ且つ添付図面を参照して以下の詳細な説明を読むときによりよく理解される。
【図面の簡単な説明】
【００１８】
【図１】モデル化された原子炉炉心内のマクロ断面積、中性子インカレント、中性子束、中性子アウトカレントの間の関係の古典的表示である。
【図２】本発明の１つの側面による、原子炉炉心の分割及び立方体の界面と推進因子の関連を示す概略図である。
【図３】図２の立方体を出る中性子アウトカレントの等方性及び異方性部分の分解概略図である。
【図４】本発明の他の側面による反復解法処理からの多レベルＶサイクルの概略図である。
【図５】図４のボックスＶの詳細な図である。
【図６】本発明の他の側面による原子炉炉心の分割を示す概略図である。
【図７】従来技術方法及び本発明による方法の異なる側面毎に１組の固有系反復解処理の１組の収束曲線である。
【発明を実施するための形態】
【００１９】
以下の説明では、加圧水型原子炉（ＰＷＲ）の場合が考えられるが、本発明を他のタイプの原子炉に適用することも心に留めておくべきである。
【００２０】
本発明によるコンピュータ実行モデル化方法の第１ステップでは、原子炉の炉心が、従来技術の状態におけるように（図２に示す）立方体１０に分割される。各立方体１０は、数値計算がコンピュータにより実行される格子或いはネットワーク１２のノードに対応する。
【００２１】
表示を簡単にするために、実際には、３次元であることを心に留めておくべきである。
【００２２】
（図２の中心の）立方体インデックスmを有する立方体１０の隣りは、それぞれの立方体インデックス m’1(m)、m’2(m)、m’3(m)及びm’4(m)をつけた立方体１０である。以下の詳細な説明では、立方体１０を、立方体のそれぞれの立方体インデックスによって直接示す。
【００２３】
本発明によるコンピュータ実行モデル化方法の第２ステップでは、炉心内の中性子束Φgm は、定常状態拡散方程式（１）に対応する固有系の解によって計算される。この目的のために、反復解手順が用いられる。
【００２４】
除去演算子
、
非散乱演算子
、生成演算子
、インカレント演算子
が方程式（１）から、
【数４】
（４）
によって、定義される。
【００２５】
等方性アウトカレント世代演算子
及び単方向性カレント貫流演算子
は、
【数５】
（５）
によって定義される。
【００２６】
結合演算子
は、隣接する立方体１０についてアウトカレントをインカレントに結合し、
【数６】
（６）
によって定義される。
【００２７】
この結合演算子
は、例えば、立方体mについて、u方向左向きアウトカレントが、デカルト軸uの方向における左隣りのu方向左向きインカレントと等しく、又同様の等式が他の方向及び向きについて確かめられている事実を用いる。例えば、図２に示す隣りm’1(m)、m’2(m)、m’3(m)、m’4(m)が、それぞれ、方程式（6）のデカルト軸yの隣り(m+1)、デカルト軸xの隣り(m-1)、デカルト軸xの隣り(m+1)、デカルト軸yの隣り(m-1)であるとき、それぞれの隣りm’1(m)、m’2(m)、m’3(m)、m’4(m)から立方体mに至る中性子アウトカレントは、立方体mに向かう中性子インカレントである。
【００２８】
上記の演算子を用いると、方程式（１）及び（２）は
【数７】
（７）
と書かれ、
又は解かれるべき固有系である、
【数８】
と書かれる。
【００２９】
方程式（８）では、Φは中性子束縦ベクトルであり、各要素は、それぞれの立方体mに対し且つそれぞれのエネルギー群gに対する中性子束Φgmである（図１）。かくして、中性子束縦ベクトルΦの次元は、（GxM, 1）に等しく、ここに、Gはエネルギー群の数、Mは立方体１０の数である。 j(out) は、中性子アウトカレント縦ベクトルであり、各要素は、それぞれの立方体m、エネルギー群g及びuがx, y或いはzに等しいデカルト軸uについての、それぞれの中性子アウトカレント jgul-m, jgur+m,である。 j(in)は、中性子インカレント縦ベクトルであり、各要素は、それぞれの立方体m、エネルギー群g、デカルト軸uについての、それぞれの中性子インカレント jgul+m, jgur-mである。かくして、中性子アウトカレント縦ベクトル j(out)及び中性子インカレント縦ベクトルj(in)の次元は（6xGxM, 1）に等しい。以下の説明では、中性子アウトカレント縦ベクトル j(out) は、また、 joutと記載される。
【００３０】
かくして、式（８）によって定義される固有系の反復 (Φ, j(out), j(in))の成分は、立方体m毎に及びエネルギー群g毎に計算されるべき中性子束Φgm、中性子アウトカレント jgul-m, jgur+m、及び中性子インカレント jgul+m, jgur-mに対応し、これは、中性子束縦ベクトルΦ、中性子アウトカレント縦ベクトルj(out)及び中性子インカレント縦ベクトル j(in)のそれぞれの要素である。
【００３１】
反復解処理は、固有系を予備の固有系に調節するサブステップを含み、予備の固有系の反復の成分(j(out), j(in))は、中性子アウトカレント縦ベクトルj(out)及び中性子インカレント縦ベクトル j(in)のそれぞれの要素である、それぞれの立方体mに至る中性子インカレント jgul+m, jgur-m、或いはそれぞれの立方体mから来る中性子アウトカレント jgul-m, jgur+mの何れかにのみ対応する。
【００３２】
予備の固有系への固有系の調整は
【数９】
と書かれる方程式（８）の第１部分を変形することから始まる。
【００３３】
シンボリック演算子
を定義すると、中性子束縦ベクトルΦは、中性子インカレント縦ベクトル j(in) の関数で
【数１０】
（１０）
と表される。
【００３４】
この式を式７のアウトカレント平衡部分に、すなわち方程式８の第２部分に代入することによって、カレントのみの関係が得られる。
【数１１】
（１１）
【００３５】
次の演算子記号を定義した後に
【数１２】
（１２）
及び
【数１３】
（１３）
【００３６】
カレントのみの関係は、スペア固有系である、
【数１４】
（１４）
と書かれる。
【００３７】
中性子アウトカレント縦ベクトルj(out)と中性子インカレント縦ベクトルj(in)の間の関係
（ここに
は結合演算子である。）を用いると、スペア固有系が、スペア固有系の反復の成分(j(out)) が、中性子アウトカレント縦ベクトル j(out)の要素である、それぞれの立方体mから来る中性子アウトカレントjgul-m, jgur+m にのみ対応する。このスペア固有系は、
【数１５】
に対応する。
【００３８】
演算子
は、最初はちょうど１に等しくなく、 j(out)が正確な解 j(out)exactに収束するときのみ１に等しくなるので、もし同時に第１中性子固有値λが正確な値λexactに収束するならば、方程式１５は
【数１６】
（１６）
になる。
【００３９】
もし第１中性子固有値λが正確な値λexactに収束するならば、第２中性子固有値μは、１に向かって収束する。
【００４０】
数値最適化のために、方程式１６は、演算子Π-1を用いて前からの乗法により変形され、以下の方程式が得られる。
【数１７】
（１７）
【００４１】
項は、デカルト方向uに依存しない等方性の項である。
【００４２】
立方体インデックスm、デカルト軸u、u’、デカルト軸u、u’に沿うそれぞれの向きs、s’の方程式１７からそれぞれの項の式が、
〔数１８〕〔数１９〕〔数２０〕
（１８）（１９）（２０）
によって与えられる。
【００４３】
G、G’は、エネルギー群の組である。
と
は、それぞれの因子
と
を有する単一エネルギー演算子である。
は、スペクトル演算子
の因子である。
【００４４】
各組み合わせそれぞれ {u, s}、 {u’, s’}は、u、u’がx、y或いはzに等しく、s、s’がl或いはrに等しい立方体mの界面を定義することに気づくべきである。例えば、{x, l} は、デカルト軸xに沿う立方体mの左界面を定義する。
【００４５】
個々の項は
【数２１】
（２１）
【数２２】
【数２３】
（２３）
によって与えられる。
【００４６】
次に、方程式１７は、例えば、普通の反復解手順をガウス-ザイデル処理として用いることによって解かれる。
【００４７】
かくして、解は中性子アウトカレント縦ベクトルj(out)の要素 j(out)mgusについて計算される。これは、式１４による中性子インカレント縦ベクトル j(in)の解を決定でき、且つ最後に式１０による中性子束縦ベクトルΦの解を決定できる。
【００４８】
本モデル化方法によって得られた、計算された中性子束縦ベクトルΦは、例えば、核燃料を管理するために既存の原子炉炉心を制御するために用いられ、或いは新しい炉心を建設するために用いることができる。
【００４９】
モデル化方法は並列プロセッサで実行されてもよいし或いは単一プロセッサで実行されてもよい。上述のモデル化方法は、よりよい収束精度、よりよいコンピュータエラー強さ及びよりよい計算効率をもたらすことが証明されている。これは、その反復の成分のみが中性子アウトカレント j(out)に対応し、中性子束に依存しない、予備固有系の解に関連する。
【００５０】
さらに、本発明によるモデル化方法の収束精度は、1E-12まで改善され、これに対して、古典的なモデル化方法で得られる収束精度は1E-6に制限される。
【００５１】
他の実施形態では、解かれるべきスペア固有系の反復の成分は、中性子インカレント j(in)にのみ対応する。
【００５２】
計算エラー強さ及び計算効率をさらに改善するために、本発明の第２側面によれば、固有系は、先ず、いくつかの中性子エネルギー群の選択のための、固有系に対応する制限固有系で調整される。この選択は、また、いくつかのエネルギー群へのスペクトル制限と呼ばれる。次いで、制限固有系は、本発明の第１側面により解かれてもよい。
最後に、制限固有系の解は、固有系を解くために用いられる。
【００５３】
選択されたエネルギー群の数ＮＧＣは、エネルギー群の全数ngよりも小さい。ＮＧＣは、より少数の粗エネルギー群に併合される細かいエネルギー群の集合である
粗くされたスペクトル帯の数であってもよい。
【００５４】
ngは、例えば、８に等しく、ＮＧＣは、例えば、４、３、２又は１に等しくてもよい。
【００５５】
この第２の側面によると、それぞれの項が方程式１８乃至２０によって与えられて、推進因子
の方程式１７は、
【数２４】
（２４）
と表される。
ここに、G、G’は粗い中性子エネルギー群インデックスである。
【００５６】
等方性項
はデカルト軸u及び向きsによって特定される出て行くカレント方向と無関係であり、非等方性項とも呼ばれる貫流項
は、同じデカルト軸uに沿う同じ粗中性子エネルギー群Ｇ内で定義される。
【００５７】
立方体mを出るアウトカレントの等方性項及び非等方性項の分解を、図３に示す。
【００５８】
等方性項は全ての方向に同じであり、非等方性項は、立方体mの界面の間で変化する。
【００５９】
最初に、等方性項が容易に計算されるアウトカレントの分解とともに、いくつかのエネルギー群へのスペクトル制限について固有系を解くことは、よりよい計算エラー強さ及びよりよい計算効率をもたらす。これは、本発明の第２の側面によるスペクトル制限から生じる固有系の次元の減少に関連する。
【００６０】
エラー強さ及び計算効率をさらに改善するために、本発明の第３の側面によれば、複数のエネルギー群のための反復解処理は、図４に示されるように、頂レベル２１、第１中間レベル２２、４つの第２中間レベル２４、２６、２８、３０及び３つの底レベル３２、３４、３６を含む、多レベルＶサイクル２０の形態であるのがよい。
【００６１】
Ｖサイクルの頂レベル２１は、複数の中性子エネルギー群について、定常状態拡散方程式１、すなわちＮＥＭ方程式、に対応する固有系のための反復計算と、本発明の第２の側面によるいくつかのエネルギー群へのスペクトル制限のために制限固有系への固有系の調整を含む。頂レベル２１から生じる制限固有系は、次いで、頂レベル２１のすぐ下の第１中間レベル２２に送り込まれる。
【００６２】
第１中間レベル２２は、制限固有系を本発明の第１の側面によるスペア固有系に調整することを含み、スペア固有系の反復(j(out), j(in))の成分は、立方体mに入る中性子インカレントjgul+m, jgur-m或いは立方体mから去る、j(out)mgusとも記述される中性子アウトカレント jgul-m, jgur+m、のいずれかにのみ対応する。方程式１８乃至２０によって与えられる演算子
,
及び
の因子は、第１中間レベル２２のために陽関数表示で計算される。
【００６３】
各第２中間レベル２４、２６、２８、３０は、中性子エネルギー群の前の選択についての前のスペア固有系の、中性子エネルギー群の後の選択についての後のスペア固有系への調整を含み、後の選択での中性子エネルギー群の数は、前の選択での中性子エネルギー群の数より少ない。多レベルＶサイクル２０の下方向きで、第１中間レベル２２に続く第２中間レベル２４の前のスペア固有系は、第１中間レベル２２から生じるスペア固有系である。第１中間レベル２２に続かない各第２中間レベル２６、２８、３０の前のスペア固有系は、先行する第２中間レベル２４、２６、２８から生じる後のスペア固有系に対応する。言い換えれば、中性子エネルギー群の数は、第２中間レベルから次の第２中間レベルに下向きに減少する。最後の第２中間レベル３０では、中性子エネルギー群の数は、１に等しい。第２中間レベル２４、２６、２８、３０の中性子エネルギー群の数のこの制限は、後述するスペクトル圧縮代数によって、行われ、本発明の第２の側面によらなくてもよいことに気づくべきである。
【００６４】
底レベル３２、３４、３６は、Ｖサイクルの下方向きでは、最後の第２中間レベル３０で決定された単エネルギー群についての最後のスペア固有系の、従来技術手順の状態による解を含む。底レベル３４、３６は、立方体が先行する底レベル３２、３４の立方体よりも大きいそれぞれの粗い格子の、先行する底レベル３２、３４から生じる固有系の空間制限に対応する。
【００６５】
Ｖサイクルの底レベル３６では、単一エネルギー群についての固有系の解がコンピュータで計算される。
【００６６】
次いで、この解は、上方向きのＶサイクルの上レベルに再び導入される、そのため、それぞれのスペア固有系のために解が計算される。
【００６７】
上方向きのＶサイクルの頂レベル２１では、複数のエネルギー群についてＮＥＭ方程式の解が計算される。
【００６８】
図４の例示実施形態では、頂レベル２１でのＮＥＭ方程式について中性子エネルギー群の数ngは８に等しい。中間レベル２２、２４、２６、２８、３０について、それぞれの選択での中性子エネルギー群の数は、４、３、２及び１にそれぞれに等しい。図４に示すＶサイクルの３つの底レベル３２、３４、３６は、ＣＭＲ手順により１つのエネルギー群のための最後のスペア固有系の解を含む。
【００６９】
それぞれのスペア固有系の解に対応する４つの中間レベル２２、２４、２６、２８は、また、それぞれの参照記号SR4、SR3、SR2及びSR1（図５）によって示す。
【００７０】
ＳＲ４レベルの単一エネルギー演算子
の因子を陽関数表示的に計算する。引き続いて、ＳＲ３レベルの単一エネルギー演算子
の因子を、図５に示すように、レベルＳＲ４の単一エネルギー演算子の因子から導き出す、ここに各ボックス３７は、中性子エネルギー群を
【数２５】
（２５）
により表す。
【００７１】
式２５を、レベル２２と２４の間の矢印３８によって図５に示す。
【００７２】
引き続いて、ＳＲ２レベルの単一エネルギー演算子
の因子を、
【数２６】
（２６）
によりＳＲ３レベルの単一エネルギー演算子の因子から導かれる。
【００７３】
式２６を、レベル２４と２６の間の矢印４０によって図５に示す。最後に、ＳＲ１レベルの単一エネルギー演算子
の因子を、
【数２７】
（２７）
によりＳＲ２レベルの単一エネルギー演算子の因子から導き出す。
【００７４】
式２７を、レベル２６と２８の間の矢印４２によって図５に示す。
【００７５】
方程式２５乃至２７は、図５に示すスペクトル格子分割戦略に対応する単一エネルギー演算子
のスペクトル圧縮代数を定義する。異なる式にあっては単エネルギー演算子
の他のスペクトル凝縮代数を、他のスペクトル格子分割戦略を用いて、すなわちボックス３７の異なる配置を用いて定義してもよいことに気づくべきである。
【００７６】
単一エネルギー演算子
の因子を、下方向きにＳＲ４レベルからＳＲ１レベルまで単一エネルギー演算子
の因子について同様の方法で導き出す。ＳＲ４レベルのスペクトル演算子
の因子を、明確に計算する。引き続いて、ＳＲ３レベルのスペクトル演算子
の因子を、
【数２８】
（２８）
によってＳＲ４レベルのスペクトル演算子
の因子から導き出す。
【００７７】
引き続いて、ＳＲ２レベルのスペクトル演算子
の因子を、
【数２９】
（２９）
によってＳＲ３レベルのスペクトル演算子の因子から導き出す。
【００７８】
最後に、ＳＲ１レベルのスペクトル演算子
の因子を、
【数３０】
（３０）
によってＳＲ２レベルのスペクトル演算子の因子から導き出す。
【００７９】
方程式２８乃至３０は、単一エネルギー演算子
及び
の、図示しないスペクトル格子分割戦略に対応し、図５に示すスペクトル格子分割戦略と同様のスペクトル演算子
のスペクトル凝縮代数を定義する。スペクトル演算子
の他のスペクトル凝縮代数は、他のスペクトル格子分割戦略を用いて定義してもよいことを気づくべきである。
【００８０】
かくして、上方向きでのＶサイクルの頂レベル２１で中性子アウトカレント縦ベクトル j(out) の要素 j(out)mgus の解を計算する。この計算の解により、式１４による中性子インカレント縦ベクトル j(in) の解を決定することができ、最後に式１０による中性子束縦ベクトルΦの解を決定することができる。
【００８１】
次いで、原子炉炉心を、計算された中性子束縦ベクトルΦに基づいて建設し或いは運転する。
【００８２】
それぞれのＳＲ３、ＳＲ２及びＳＲ１レベルでの演算子
,
及び
について、複雑な算術演算を伴わない、スペクトル圧縮代数を用いることによって、演算子
,
及び
をとても安価に計算でき、かくして、モデル化方法のよりよい計算エラー強さ及びよりよい計算効率に至る。
【００８３】
本発明の第４の側面により、エラー強さ及び計算効率をさらに改善するために、立方体１０を第１カテゴリーと第２カテゴリーに分割する。
【００８４】
以下の説明では、図６に示すように、第１カテゴリーの立方体１０Ｒは、赤立方体と呼ばれ、第２カテゴリーの立方体１０Ｂは黒立方体と呼ばれるが、特別の限定的な意味を用語”黒”及び”赤”と関連させるべきではない。各赤立方体１０Ｒは、すぐ隣の立方体として黒立方体１０Ｂのみを有する。かくして、赤立方体１０Ｒのほとんどは、６つのすぐ隣りの黒立方体１０Ｂを有する。これは、立方体が考慮される立方体と共通の面を共有している、”すぐ”隣りによって理解されるべきである。
【００８５】
その結果、図６により例示されるように、格子１２は２次元表現の市松模様の暗い領域と明るい領域に関して視覚的に類似している。
【００８６】
次いで、立方体には例えば赤立方体１０Ｒによって始め、黒立方体１０Ｂによって終わる番号が付けられる。
【００８７】
以下の説明では、２つのカテゴリーでの立方体のこのような分割及びカテゴリーの次々の番号付けは、赤−黒順序付けと称される。
【００８８】
当該技術の辞書式順序付けの状態と比較して立方体の赤−黒順序付けの利点は、赤立方体１０Ｒについて、そのすぐ隣りの立方体が全て黒であって、逆の場合も同じである。
【００８９】
コンピュータが炉心内の中性子束を計算するために解かねばならない式２４は行列ベクトル形式で式数３１と書くことができ、
【数３１】
（３１）
【００９０】
典型的には０．９乃至０．９５の値からなり、そのため、積 ｙμがシフトを形成する、可変のパラメータｙを用いると、式３１の反復は次のシフト逆数陰関数形式で書かされる。
【数３２】
（３２）
【００９１】
ここに、s(n-1)は
【数３３】
（３３）
によって与えられる源である。
【００９２】
スペクトル演算子
及び単一エネルギー演算子
は希薄であり、赤立方体１０Ｒをすぐ隣りの黒立方体１０Ｂにのみ結合し、又その逆に結合する。かくして、赤黒順序付けは、赤の部分と黒の部分との間に次の便利な関係を可能にする。
【数３４】
（３４）
【００９３】
方程式３４が、解の赤の部分のみについて次の等式に変換される。
【数３５】
【００９４】
計算の観点から、赤黒順序付けの使用は、反復解手順中、もし、反復の最初の半分で、反復 d の、赤成分 dred が更新されるならば、次いで、反復のもう一つの部分中、全ての黒成分 dblackは、それらの隣りの赤成分dred に基づいて更新されることを意味する。換言すれば、各黒立方体１０Ｂの値は、全てのすぐ隣りの赤立方体１０Ｒの値に基づいて計算される。
【００９５】
赤黒順序付けを用いることによって、本発明の第１、第２又は第３の側面で記載された演算子
,
及び
の計算を改善する。かくして、赤黒順序付けは、計算効率を改善するために本発明の第１、第２或いは第３の側面を補完に用いられてもよい。
【００９６】
そのような赤黒順序付けは、本発明の第５側面を構成する特定の解手順と一緒に用いられるとき特に有用であることがわかる。
【００９７】
この手順は、特に高いエラー強さ安定化双共役勾配安定化（Bi-CGStab）手順である。この手順の適当な前置きの記載は、次の引用文献に見られる。
Y. Saad, "スパース線形システム反復法 (Iterative Methods for Sparse Linear Systems)", 第２版, 工業・応用数学会 (SIAM) (2003;)、
H.A. van der Vorst, "Bi-CGSTAB: 非対称の線形システムの解のためのBi-CGの高速且つ円滑な収束バリアント (a Fast and Smoothly Converging Variant of Bi-CG for the solution of nonsymmetric linear systems)", 日本応用数理学会 SIAM J. Sci. Stat. Comput. 13(2), pp.631-644 (1992)。
【００９８】
双共役勾配安定化（Bi-CGStab）手順は、与えられたΘとsにより、dの関数のための最小化原理から導かれ、そのために定常を方程式３５によって与えられる線形系をちょうど満たす特定の最適なdを中心とした小さな変分δdに関して適用し、且つ汎関数が最小値を仮定する。
【００９９】
かくして、式３５の効率的な解き方について、より直接的な考察によってではなく汎関数の最小化によって動かされる解手順を定義することが可能である。Bi-CGStab手順は、そのような最小原理から続き、反復での連続変化は、（Galerkin加重積分剰余のようである）関数の各変化が先の変化の全てに関して直交するように配列される。
【０１００】
本発明の第５側面による特殊なBi-CGStab手順が、解ベクトルd、解誤差 r ( r = s - Θ d)及び補助ベクトル r*, s及び p、で及び初期解推定値d0で下に与えられる。
【数３６】
（３６）
【０１０１】
上記Bi-CGStab数列は、普通Ｎ＜３のＮステップ後に打ち切られる。
【０１０２】
Bi-CGStab手順のために、予備調整が式３２のシフト逆陰関数形式に基づいており、ただし、典型的に0.9<ｙ<0.95である。
【０１０３】
シフトｙμのためのこの選択で、解処理中、ｙμが、系の最も小さい可能な固有値のｙ倍まで収束するので、演算子sは、非特異のままにするように保証される。
【０１０４】
方程式３２を解くことは、高等予備調整Krylov手順を展開するための完全なシナリオである、中性子源を有するわずかに臨界未満系の中性子束分布を決定することに数値的に等しい。
【０１０５】
赤黒順序付けを用いると、方程式３２は、上で説明したように、方程式３５に変換される。方程式３５によって与えられる予備調整系は、本発明の第６の側面を構成する。
【０１０６】
いったん方程式３５が、Bi-CGStab手順を介して解かれると、解黒部分が方程式３４の黒部分から計算される。
【０１０７】
格子上に投影されるように、隣りの（赤対黒）ノードの間の直接の中性子相互作用を意味する、この特定の予備調整法の使用は、解かれるべき系に予め含まれ、
【数３７】
（３７）
は、
【数３８】
（３８）
により与えられるゆえに、方程式３５の単位演算子をより支配的に形成する。
【０１０８】
この予備調整法は、空間結合、それ故に式の解を決定するノードの間の物理的な相互作用に関して、重要な情報、又は、その情報の大部分を、低計算コストで、うまく予め含める。
【０１０９】
他の実施形態では、反復解手順は、一般化した最小剰余法、また短縮GMRES、或いは平面ヤコビ法のようなザイデル処理、或いはクリロフ部分空間での反復手順である。
【０１１０】
本発明の第７側面によれば、改善された固有値計算のための変分の原理を次の説明に記載する。
【０１１１】
本発明の第１、第２或いは第３の側面について上述されたスペア固有系に対応する式（１６）は、第２中性子固有値μ及び中性子アウトカレント縦ベクトルjoutによって、次の形に書かされる。
【数３９】
（３９）
【０１１２】
ここに、c1は、立方体の特性により演算子
のための他の記号である。
【０１１３】
何が続くかにおいて、我々は、第１中性子固有値λがスカラー制御パラメータの役割を満たすと仮定する。もし全体の解が収束し、かくして安定化するならば、第２中性子固有値μが、１に近づく特性を用いると、変分原理又は摂動アプローチは、第１中性子固有値λを変化させて、第２中性子固有値μを１に向かわせ且つスペア固有系解の収束を加速することを含む。
【０１１４】
第１中性子固有値λのこの変分は、スペア固有系の各反復前に、すなわちスペア固有系の各第２、第３、第４或いは第５の反復前に、なされるのがよい。第１中性子固有値λの変分は、λがλ+δλにされるように偏差δλのλへの適用であり、これは摂動を起こした界面カレント方程式、
【数４０】
（４０）
ただし、演算子偏微分
【数４１】
（４１）
により第２中性子固有値μに影響を与える。
【０１１５】
摂動
が反復手順中ますます小さくなると期待されるので、摂動は方程式（４０）では都合よく無視され、任意の付加した重み関数
での積分は次式を得る。
【数４２】
（４２）
【０１１６】
第１中性子固有値λの変分が、第２中性子固有値μが１に向かって推進され、すなわちμ+δμが１に向かって推進されるようなものであるとき、近似
が置かれる。
この近似値及び方程式（４２）で、次の方程式が得られる。
【数４３】
（４３）
【０１１７】
等式
を用いると、δλが
【数４４】
（４４）
で与えられる。
【０１１８】
ここに、rout は、第２中性子固有値μのために置かれた目標値が１に等しいとき、界面アウトカレント平衡方程式、すなわち方程式（３９）の剰余を表す。剰余rout は、いかなる場合の各ＮＥＭ反復ステップ毎に計算される必要がある。
【０１１９】
この側面でスカラー制御パラメータと考慮される、第１中性子固有値λは、次いで、次式のように更新される。
【数４５】
（４５）
【０１２０】
ここに、
及び
は、それぞれ、第１中性子固有値λの先の値及び更新された値である。
【０１２１】
部分的な演算子微分係数
を決定し且つ式７を使うために、例えば中性子束Φは、
【数４６】
（４６）
と書かれる。
上式から：
【数４７】
が続く。
【０１２２】
次いで、近似導関数
は、
のような非散乱演算子
に対する行列の上対角部分を無視することによって、すなわち上方散乱、
【数４８】
（４８）
を無視することによって計算される。
【０１２３】
多くの計算の場合に、これは近似でもないが、もし上方散乱がモデル化されなければ、数学的に正確であることに気づくべきである。
【０１２４】
変分原理の好結果の適用にとって十分である
【数４９】
（４９）
を仮定する。
【０１２５】
を計算するために、それを先ず、
【数５０】
（５０）
と書き、
上式から：
【数５１】
（５１）
が続く。
【０１２６】
次いで、
【数５２】
（５２）
で演算子記号
を定義し、
を陰関数的に計算するために、テーラー式を
に適用して次式が得られる。
【数５３】
（５３）
【０１２７】
テーラー展開が、
【数５４】
（５４）
と並べ換える。
【０１２８】
方程式（５４）への連鎖律の適用により、導関数
が、
【数５５】
（５５）
によって与えられ、
これは
【数５６】
（５６）
を得る。
【０１２９】
導関数
を決定するのに効率的な反復スキームを可能にする最後の式は、
で方程式（５６）の前もっての乗法によって得られ、これは
【数５７】
（５７）
を与える。
【０１３０】
と書くことによって、
の計算が
【数５８】
（５８）
から得られる。
源の項sは、
【数５９】
（５９）
によって定義される。
【０１３１】
２つのエネルギー群の場合には、方程式（５８）は直接的に及び解析的に解かれる。
【０１３２】
２つ以上のエネルギー群の場合には、方程式（５８）は、
【数６０】
（６０）
の適用により、正確さが早期反復段階で要求されることなく、反復的に解かれる。
ここに、
及び
は、それぞれ、近似導関数
の以前の及び更新された値である。
【０１３３】
早期のとは、反復ステップの全必要数の最初の半分までの約３分の１を意味する。
【０１３４】
この変分アプローチは、有用である
の近似値のみを必要とするので、
の計算からの反復を、少し早く打ち切れることができることに気づくべきである。
【０１３５】
かくして、本発明のこの第７側面による変分原理は、方程式（４５）の分子
が収束のときますます小さくなり、分母
は適当な値にむしろ急速に収束するので、大変効果的である。この第７の側面では、分子は、それがゼロに向かって収束するので、推進原理を定義する。
【０１３６】
方程式（４５）による第１中性子固有値λを変えることは、第２中性子固有値μを１に向かって推進し、且つスペア固有系解の収束を加速させる。第１中性子固有値λのこの変分は、スペア固有系の各反復の前に、すなわちスペア固有系の各第２、第３、第４又は第５反復の前に行われてもよい。
【０１３７】
かくして、第１中性子固有値λの変分原理は、計算効率を改善するために本発明の第１、第２或いは第３の側面の補完に用いられてもよい。
【０１３８】
第７の側面の全ての上述した結果を、方法における計算効率を失うことなしに、格子１２の黒部分に減少させることができ、そのために、本発明の第４の側面で上述した赤黒順序付けがこの第７の側面の補完に用いられる。
【０１３９】
本発明の第８の側面により、固有値計算のための上述の変分原理は、典型的に、方程式（４）によって与えられる除去演算子
を変形させ、
と記される制御パラメータの計算のために一般化でき、変形除去演算子は
と記される。
【０１４０】
全ての位置での変形の例では、制御パラメータ
とは、硼素濃度である。選択された位置での変形の例では、制御パラメータ
は、一群の選択された制御棒の棒挿入深さを制御するためのパラメータである。
【０１４１】
と仮定する。
【０１４２】
すると、方程式（１６）は、第２中性子固有値μを用いて、次の形に書かれる：
【数６１】
（６１）
シンボリックスペクトル演算子
は
【数６２】
（６２）
によって定義される。
【０１４３】
もし制御パラメータ
との補正値を含む全体の解が収束し、かくして安定するならば、第２中性子固有値μが、１つに近づく性質を用いると、変分原理、又は摂動アプローチは、第２中性子固有値μを１に向かって推進させるためにこと制御パラメータ
を変化させることを含む。
【０１４４】
制御パラメータ
の変分は、
が
に推進させるように変分
の
への適用であり、これは、摂動界面カレント方程式
【数６３】
（６３）
部分演算子導関数：
【数６４】
（６４）
により第２中性子固有値μに影響を与える。
【０１４５】
摂動
が反復手順中だんだん小さくなると期待されるので、これらは方程式（６３）で都合よく無視され、任意の随伴重み関数
での積分は次の式を得る。
【数６５】
（６５）
【０１４６】
制御パラメータ
の変分で、第２中性子固有値μが１に向かって推進され、すなわち、そのμ+δμは１に向かって推進されるので、近似
が置かれる。この近似及び式（６５）を用いて、次の方程式が得られる。
【数６６】
（６６）
【０１４７】
等式
を用いると、
は、
【数６７】
（６７）
によって与えられる。
ここに、routは、界面アウトカレント平衡方程式、すなわち、方程式（３９）の残差を表し、第２中性子固有値μのために課された目標値は１に等しい。
【０１４８】
制御パラメータ
は次のように更新される：
【数６８】
（６８）
ここで
と
は、それぞれ、制御パラメータ
の以前の値及び更新された値である。
【０１４９】
変分原理の制御パラメータ
へのうまい適用にとって十分である、
が仮定される。
【０１５０】
を計算するために、テーラー公式の適用による上述の調整が、また、同じ方法で用いられる。
【０１５１】
本発明のこれら第８の側面は、本発明の上述の第７の側面と同じ利点を示し、計算効率を改善するために、本発明の上述の第７側面のためのものと同様の方法で、本発明の第１、第２或いは第３の側面の補完に用いられてもよい。
【０１５２】
図７は、追求したＮＥＭ反復ステップの数を横座標に、高速中性子エネルギー群の中性子束の最大の解の誤差である、高速中性子束誤差を縦座標にした、異なる固有系反復解手順の１組の収束曲線を表す。全てのエネルギー群がスペクトル的に結合されるので、他のエネルギー群の中性子束の最大の解の誤差は、高速中性子束誤差と全く同様である。かくして、異なる固有系反復解手順の間の比較がこれらの収束曲線によって可能である。
【０１５３】
曲線５０は、単一精度を用いた古典的な粗メッシュ再平衡（ＣＭＲ）手順用の収束曲線であり、曲線５１は、二重精度を用いた古典的なＣＭＲ手順の収束曲線である。単一及び二重の精度は、コンピュータプログラム製品における浮動小数点数を表すために用いられる古典的な精度である。
【０１５４】
曲線５２は、単一の精度を用いた単一エネルギー群のための、本発明の第１及び第７の側面による反復手順の収束曲線であり、曲線５３は、二重精度を用いた、本発明の同じ側面による反復手順の収束曲線である。
【０１５５】
曲線５４は、単一精度手順を用いた４つの粗エネルギー群のための、本発明の第１及び第２の側面による反復手順用の収束曲線であり、曲線５５は、二重精度を用いた、本発明の同じ側面による反復手順用の収束曲線である。
【０１５６】
曲線５６は、単一精度を用いた４つの粗エネルギー群のための、本発明の第１、第２及び第７の側面による反復手順用の収束曲線であり、曲線５７は、二重精度を用いた、本発明の同じ側面による反復手順用の収束曲線である。
【０１５７】
二重精度を用いた４つの粗エネルギー群のための、第１、第２及び第７の側面による手順（曲線５７）は、曲線５７について最大である収束曲線の勾配によって示されているように、５０回の追求ＮＥＭ反復ステップについての値が実質的に１０-13に等しい最も小さい誤差ともっとも良い計算効率の両方を有する手順である。
【０１５８】
二重精度を用いた単一エネルギー群のための、第１及び第７の側面による手順（曲線５３）は、１３０回の追求値がＮＥＭ反復ステップについて実質的に１０-9に等しい小さい誤差を出す。さらに、この手順（曲線５３）の計算効率は、曲線５５の勾配が曲線５３の勾配よりも大きいので、二重精度を用いた４つの粗エネルギー群のための、第１及び第２の側面による手順の計算効率（曲線５５）よりも価値が高い。
【０１５９】
二重精度を用いた４つの粗エネルギー群のための、第１及び第２の側面による手順は、４０回の追求ＮＥＭ反復ステップについて実質的に１０-7に等しい誤差を出す。
【０１６０】
かくして、二重精度は、高速中性子束誤差の値に大きな影響を及ぼすが、モデル化方法の計算効率に影響を与えることがない。
【０１６１】
本発明の第７側面は、また、計算効率に影響を与えるよりも、高速中性子束誤差の値に大きな影響を及ぼす。
【０１６２】
４つの粗エネルギー群のための、第１及び第２の側面は、最良の計算効率を提供し、図７は、100と10-5の間に含まれる誤差値に関する最も大きい勾配が、第１及び第２の側面に対応する曲線５４乃至５７の全てに関して得られることを示す。
【０１６３】
前述に開示されるように、本発明の第１乃至第８の側面は、よりよい計算効率を提供するエラー強さモデル化方法を達成するのを助け、そのため、関連する炉心シミュレーションを短い計算時間内に得ることができる。
【０１６４】
第１側面は、それ自体で、この目的を達成する上で役に立ち、かくして第２側面から第８側面まで別々に用いることができることを心に留めておくべきである。また、第２、第３、第７及び第８の側面は、第１の側面で或いは第４乃至第６側面のいずれかで必ずしも実行されない。
【特許請求の範囲】
【請求項１】
原子炉炉心をモデル化するための方法であって、
コンピュータで実行される計算のための格子（１２）のノードを構成するために前記炉心を立方体（１０）に分割するステップと、
中性子束を、少なくとも１つの固有系の反復解手順を用いることによって計算するステップとを含み、固有系の反復の成分は、計算されるべきそれぞれの立方体（１０）についての、中性子束、中性子アウトカレント或るいは中性子インカレントのいずれかに対応し、
制御パラメータは、摂動界面カレント方程式によって中性子固有値μに影響を与え且つ中性子固有値μを１に向って推進させるように変えられる、ことを特徴とする前記方法。
【請求項２】
と記載される制御パラメータは、
によって与えられ、
ここに、
は、重み関数であり、Φは中性子束を示し、
rout は
によって与えられ、
は、結合演算子であり、 jout は中性子アウトカレントを示し、
は単方向性カレント貫流演算子であり、
c1は、立方体の特性に依存する演算子である、請求項１に記載の方法。
【請求項３】
と記載される制御パラメータは、硼素濃度である、請求項１又は２に記載の方法。
【請求項４】
と記載される制御パラメータは、選択された制御棒のグループについての棒挿入深さを制御するためのパラメータである、請求項１又は２に記載の方法。
【請求項５】
制御パラメータは、中性子固有値λであり、反復解手順中、中性子固有値λは中性子固有値μを１に向かって推進させるように変化し、中性子固有値λの変分は、λが、摂動界面カレント方程式によって中性子固有値μに影響を与えるλ+δλに推進されるように変分δλのλへの適用である、請求項１に記載の方法。
【請求項６】
前記中性子固有値λは、
による固有系解手順を反復する前に、予備調整され、
ここに、
は、重み関数であり、Φは中性子束を示し、
routは
によって与えられ、
は、結合演算子であり、 joutは中性子アウトカレントを示し、
は単方向性カレント貫流演算子であり、
c1 は、立方体の特性に依存する演算子である、請求項５に記載の方法。
【請求項７】
演算子偏微分
は、
の適用により反復的に解かれ、
ここで、
であり、
は、除去演算子を示し、
は、非散乱演算子
の下対角部分であり、
は生成演算子であり、
は、インカレント演算子である、請求項６に記載の方法。
【請求項８】
方法は、計算された中性子束に基づいて原子炉炉心を建設するステップを含む、請求項１乃至７の何れか１項に記載の方法。
【請求項９】
方法は、計算された中性子束に基づいて原子炉炉心を運転するステップを含む、請求項１乃至８の何れか１項に記載の方法。
【請求項１０】
方法は、コンピュータで実行される方法である、請求項１乃至９の何れか１項に記載の方法。
【請求項１１】
コンピュータプログラム製品はコンピュータ読み出し可能な媒体上に在り且つコンピュータで実行するためのコンピュータプログラム手段を含む、請求項１０に記載の方法。
【請求項１】
原子炉炉心をモデル化するための方法であって、
コンピュータで実行される計算のための格子（１２）のノードを構成するために前記炉心を立方体（１０）に分割するステップと、
中性子束を、少なくとも１つの固有系の反復解手順を用いることによって計算するステップとを含み、固有系の反復の成分は、計算されるべきそれぞれの立方体（１０）についての、中性子束、中性子アウトカレント或るいは中性子インカレントのいずれかに対応し、
制御パラメータは、摂動界面カレント方程式によって中性子固有値μに影響を与え且つ中性子固有値μを１に向って推進させるように変えられる、ことを特徴とする前記方法。
【請求項２】
と記載される制御パラメータは、
によって与えられ、
ここに、
は、重み関数であり、Φは中性子束を示し、
rout は
によって与えられ、
は、結合演算子であり、 jout は中性子アウトカレントを示し、
は単方向性カレント貫流演算子であり、
c1は、立方体の特性に依存する演算子である、請求項１に記載の方法。
【請求項３】
と記載される制御パラメータは、硼素濃度である、請求項１又は２に記載の方法。
【請求項４】
と記載される制御パラメータは、選択された制御棒のグループについての棒挿入深さを制御するためのパラメータである、請求項１又は２に記載の方法。
【請求項５】
制御パラメータは、中性子固有値λであり、反復解手順中、中性子固有値λは中性子固有値μを１に向かって推進させるように変化し、中性子固有値λの変分は、λが、摂動界面カレント方程式によって中性子固有値μに影響を与えるλ+δλに推進されるように変分δλのλへの適用である、請求項１に記載の方法。
【請求項６】
前記中性子固有値λは、
による固有系解手順を反復する前に、予備調整され、
ここに、
は、重み関数であり、Φは中性子束を示し、
routは
によって与えられ、
は、結合演算子であり、 joutは中性子アウトカレントを示し、
は単方向性カレント貫流演算子であり、
c1 は、立方体の特性に依存する演算子である、請求項５に記載の方法。
【請求項７】
演算子偏微分
は、
の適用により反復的に解かれ、
ここで、
であり、
は、除去演算子を示し、
は、非散乱演算子
の下対角部分であり、
は生成演算子であり、
は、インカレント演算子である、請求項６に記載の方法。
【請求項８】
方法は、計算された中性子束に基づいて原子炉炉心を建設するステップを含む、請求項１乃至７の何れか１項に記載の方法。
【請求項９】
方法は、計算された中性子束に基づいて原子炉炉心を運転するステップを含む、請求項１乃至８の何れか１項に記載の方法。
【請求項１０】
方法は、コンピュータで実行される方法である、請求項１乃至９の何れか１項に記載の方法。
【請求項１１】
コンピュータプログラム製品はコンピュータ読み出し可能な媒体上に在り且つコンピュータで実行するためのコンピュータプログラム手段を含む、請求項１０に記載の方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【公開番号】特開２０１１−６０２８４（Ｐ２０１１−６０２８４Ａ）
【公開日】平成２３年３月２４日（２０１１．３．２４）
【国際特許分類】
【外国語出願】
【出願番号】特願２０１０−１９７１０７（Ｐ２０１０−１９７１０７）
【出願日】平成２２年８月１７日（２０１０．８．１７）
【出願人】（５０９１０４５６２）アレヴァ エヌペ (28)
【Ｆターム（参考）】
【公開日】平成２３年３月２４日（２０１１．３．２４）
【国際特許分類】
【出願番号】特願２０１０−１９７１０７（Ｐ２０１０−１９７１０７）
【出願日】平成２２年８月１７日（２０１０．８．１７）
【出願人】（５０９１０４５６２）アレヴァ エヌペ (28)
【Ｆターム（参考）】
[ Back to top ]