(57)【要約】【課題】サブチャンネル解析コード等のマクロコードを用いてサブクール沸騰領域に関する解析を適切に行えるサブクール沸騰特性予測装置を提供すること。【解決手段】マクロコード処理装置１２からサブクール沸騰の解析対象のプロセス量を解析対象情報入力手段１３を介して入力し、初期設定手段１４において解析領域を設定すると共に、その解析領域に二相流解析のための粒子を配置する。解析領域は解析領域出力手段１５からマクロコード処理装置１２に出力される。二相流解析手段１６は初期粒子配置を入力すると共に、マクロコード処理装置１２から解析領域の壁面温度を入力して二相流解析を行う。そして、変換手段１７により二相流解析結果をマクロコード処理装置１２で処理可能なサブクール沸騰特性量に変換して出力する。【選択図】  図１
2005050042.tif 000002
【特許請求の範囲】
【請求項１】
気体および液体を粒子で模擬して、粒子の位置および速度を陽的計算ステップと陰的計算ステップとを用いて計算することにより気液二相流の解析を行う二相流解析方法において、
前記陽的計算ステップは、粒子の周辺の粒子についてカーネル関数を足し合わせて計算する粒子数密度を計算する粒子数密度の計算ステップと、
液体粒子に働く表面張力を計算する表面張力項の計算ステップ、
粒子の運動量拡散効果を計算する運動量拡散項の計算ステップと、
前記表面張力と前記拡散効果と粒子に働く重力により粒子の位置および速度を計算する位置速度の陽的仮更新ステップと、
粒子間の熱の移動および粒子の相変化を計算する熱および相変化計算ステップとを有し、
前記陰的計算ステップは、液体粒子を気体粒子と仮定して圧力方程式を解いて気体粒子の圧力を求めるステップと、
この気体粒子の圧力を求めるステップの後に液体と気体の境界の圧力条件として気体粒子の圧力を与えて圧力方程式を解いて液体粒子の圧力を求めるステップと、
粒子間の圧力の勾配による効果により粒子の速度の変化量および粒子の位置の変化量を計算する圧力勾配項の計算ステップと、
前記粒子の速度の変化量と前記粒子の位置の変化量により粒子の位置および速度を計算する位置速度の陰的仮更新ステップとを有し、
前記表面張力項の計算ステップは、周囲にある液体粒子の割合を用いて表面張力の補正を行い、
前記圧力勾配項の計算ステップは、気体粒子と液体粒子の間に働く圧力勾配項の補正を行うことを特徴とする二相流解析方法。
【請求項２】
前記陽的計算ステップは、液体の蒸発により発生する気体を模擬する気体粒子に仮想的な距離を与えて計算することを特徴とする請求項１に記載の二相流解析方法。
【請求項３】
気液二相流の解析を行うためにコンピュータを、気体および液体を模擬した粒子の位置および速度を計算する陽的計算手段および陰的計算手段として機能させる二相流解析プログラムにおいて、
前記陽的計算手段は、粒子の周辺の粒子についてカーネル関数を足し合わせて計算する粒子数密度を計算する粒子数密度の計算手段と、
液体粒子に働く表面張力を計算する表面張力項の計算手段、
粒子の運動量拡散効果を計算する運動量拡散項の計算手段と、
前記表面張力と前記拡散効果と粒子に働く重力により粒子の位置および速度を計算する位置速度の陽的仮更新手段と、
粒子間の熱の移動および粒子の相変化を計算する熱および相変化計算手段とを有し、
前記陰的計算手段は、液体粒子を気体粒子と仮定して圧力方程式を解いて気体粒子の圧力を求める手段と、
この気体粒子の圧力を求めるステップの後に液体と気体の境界の圧力条件として気体粒子の圧力を与えて圧力方程式を解いて液体粒子の圧力を求める手段と、
粒子間の圧力の勾配による効果により粒子の速度の変化量および粒子の位置の変化量を計算する圧力勾配項の計算手段と、
前記粒子の速度の変化量と前記粒子の位置の変化量により粒子の位置および速度を計算する位置速度の陰的仮更新手段プとを有し、
前記表面張力項の計算手段は、周囲にある液体粒子の割合を用いて表面張力の補正を行い、
前記圧力勾配項の計算手段は、気体粒子と液体粒子の間に働く圧力勾配項の補正を行うことを特徴とする二相流解析プログラム。
【請求項４】
前記陽的計算手段は、液体の蒸発により発生する気体を模擬する気体粒子に仮想的な距離を与えて計算することを特徴とする請求項３に記載の二相流解析プログラム。
【請求項５】
解析対象をメッシュに分割してサブクール沸騰の解析を行うためのマクロコード処理装置から前記サブクール沸騰の解析対象のプロセス量を入力する解析対象情報入力手段と、
前記解析対象情報入力手段から入力された解析対象のプロセス量に基づいて解析領域を設定すると共にその解析領域に二相流解析のための粒子を配置する初期設定手段と、
前記初期設定手段で設定した解析領域を前記マクロコード処理装置に出力する解析領域出力手段と、
前記解析領域出力手段を介して前記初期設定手段で配置された粒子配置を入力すると共に前記マクロコード処理装置から解析領域の壁面温度を入力し二相流解析を行う二相流解析手段と、
前記二相流解析手段の二相流解析結果を前記マクロコード処理装置で処理可能なサブクール沸騰特性量に変換する変換手段とを備えたことを特徴とするサブクール沸騰特性量予測装置。
【請求項６】
前記変換手段で変換されたサブクール沸騰特性量に基づいてサブクール沸騰特性量に関する相関式を演算する相関式作成手段を備えたことを特徴とする請求項５記載のサブクール沸騰特性量予測装置。
【請求項７】
前記二相流解析手段は、前記初期設定手段で設定された初期粒子配列を入力すると共に前記マクロコード処理装置から壁面温度を入力して表面張力を計算し、粒子速度・位置座標の陽的な仮の更新を行って相変化計算を行い、次に仮想付加次元計算を行い、圧力方程式を解法して気液間相互作用の高精度化アルゴリズム処理を行い、圧力項による粒子の速度・位置座標の陰的更新を行い、圧力補正が収束するまで繰り返し演算を行い、圧力補正が収束したときは二相流解析結果を出力することを特徴とする請求項５記載のサブクール沸騰特性量予測装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、流路断面に蒸気と過冷液とが共存する状態におけるサブクール沸騰の特性を計測する二相流解析方法、二相流解析プログラム、およびサブクール沸騰特性量予測装置に関する。
【０００２】
【従来の技術】
サブクール沸騰領域は、流路断面内に蒸気と過冷液とが共存するため熱力学的非平衡状態にあり、管長に沿うボイド率分布は、伝熱と流れの絡んだ複雑な過程で決まる。例えば、水冷却型原子炉においてはボイド率が系の核的反応度に影響を有するので、この領域のボイド率を予測することが必要となる。
【０００３】
サブチャンネル解析コード等のマクロコードを用いてサブクール沸騰領域に関する解析を実施する場合、ペクレ数、スタントン数と、ボイド率、サブクール度等のサブクール沸騰特性量に関する相関式が必要となる。そこで、必要となる相関式は実験結果から作成するようにしている。また、解析によりサブクール沸騰特性量に関する相関式を得るには、流体解析対象の空間に計算格子を設定し流体解析を行うことになる。
【０００４】
【発明が解決しようとする課題】
ところが、実験結果から得られた相関式は、限られた体系形状や流体の熱的条件に対してのみ適用可能であり、体系形状や流体の熱的条件が異なる場合には、その都度実験を行い相関式を求めなければならない。
【０００５】
また、解析により相関式を求めるには、相変化を伴う気液が流体解析の対象である場合には、設定した計算格子の一つ一つの内部において、液体が気体の流れの中でどのような形状をとるかを考慮しなければならない。従って、実験から決定される流動様式や気液の間の抵抗等に関する多くの構成方程式が必要となる。
そのため、それら構成方程式が得られていない条件においてはサブクール沸騰特性量間の関連付けを行うことは困難である。
【０００６】
このように、実験結果からサブクール沸騰特性量に関する相関式を作成することは、複雑な体系に対しては困難であり、その結果、適用範囲が限定された相関式しか得られていない。従って、流路の形状が異なるとその都度実験を行うことが必要となる。特に、原子炉燃料集合体のスペーサ部に見られるような複雑構造に対しては、スペーサの形状を変更する度に実験が必要となるが、これは時間的にもコストの面からも実現は困難である。
【０００７】
本発明の目的は、サブチャンネル解析コード等のマクロコードを用いてサブクール沸騰領域に関する解析を適切に行うことができる二相流解析方法、二相流解析プログラム、およびサブクール沸騰特性量予測装置を提供することである。
【０００８】
【課題を解決するための手段】
本発明は、解析対象をメッシュに分割してサブクール沸騰の解析を行うためのマクロコード処理装置からサブクール沸騰の解析対象のプロセス量を解析対象情報入力手段を介して入力し、初期設定手段において解析領域を設定すると共に、その解析領域に二相流解析のための粒子を配置し、解析領域出力手段から初期設定手段で設定した解析領域をマクロコード処理装置に出力する。一方、初期設定手段で配置された粒子配置を解析領域出力手段を介して二相流解析手段に入力すると共に、マクロコード処理装置から解析領域の壁面温度を二相流解析手段に入力し、二相流解析手段で気液二相流の解析を行う。そして、変換手段により二相流解析手段の二相流解析結果をマクロコード処理装置で処理可能なサブクール沸騰特性量に変換して出力する。また、必要に応じて、変換手段で変換されたサブクール沸騰特性量に基づいてサブクール沸騰特性量に関する相関式を演算する相関式作成手段を備える。
【０００９】
【発明の実施の形態】
以下、本発明の実施の形態を説明する。図１は本発明の実施の形態に係わるサブクール沸騰特性量予測装置１１のブロック構成図である。以下の説明では解析対象として、沸騰水型原子炉における燃料集合体の場合について説明する。
【００１０】
マクロコード処理装置１２は、解析対象である燃料集合体をメッシュに分割してサブクール沸騰の解析を行うものであり、燃料集合体のプロセス量ａであるセル情報や予想ボイド率等を初期値としてサブクール沸騰特性量予測装置１１に出力する。具体的には、燃料集合体の内径、マクロセル境界位置座標、流体速度、気体速度、ボイド率、流体温度、気体温度、圧力、壁面温度等のプロセス量ａを初期値としてサブクール沸騰特性量予測装置１１の解析対象情報入力手段１３に出力する。
【００１１】
サブクール沸騰特性量予測装置１１の初期設定手段１４は、解析対象情報入力手段１３から入力された解析対象のプロセス量に基づいて解析領域を設定すると共に、その解析領域に二相流解析のための粒子を配置して粒子配置の初期設定を行う。解析領域出力手段１５は初期設定手段１４で設定された解析領域ｂをマクロコード処理装置１２に出力する。具体的には、二相流解析手段１６での解析領域の開始・終止マクロ番号を出力する。それと共に、解析領域出力手段１５は、初期設定手段１４で配置された粒子配置を二相流解析手段１６に出力する。
【００１２】
二相流解析手段１６は、解析領域の粒子配置を入力すると、マクロコード処理装置１２から解析処理情報ｃ（マクロコード処理装置１２による計算の各タイムステップの時間刻みと壁面温度）を入力し、各粒子の相変化に伴う圧力補正を加味して二相流解析を行う。二相流解析手段１６の二相流解析結果ｄは、変換手段１７によりマクロコード処理装置１２で処理可能なサブクール沸騰特性量ｅに変換され、マクロコード処理装置１２に出力される。
【００１３】
二相流解析手段１６の二相流解析結果ｄは、具体的には、熱流速、液体粒子の位置座標、速度、液体粒子の温度、圧力、気体粒子の位置座標、速度、気体粒子の温度、圧力等であり、サブクール沸騰特性量ｅは、具体的には、ボイド率、サブクール度、ペクレ数、スタントン数である。
【００１４】
このように、まず、マクロコード処理装置１２からサブクール沸騰特性量予測装置１１へセル情報とマクロコード処理装置１２によって予想されたボイド率の分布を送る。サブクール沸騰特性量予測装置１１では二相流解析手段１６で解析する領域をマクロコード処理装置１２によって予想されたボイド率分布を用いて選択する。選択されたセルのセル番号が解析領域ｂとしてサブクール沸騰特性量予測装置１１からマクロコード処理装置１２に送られる。
【００１５】
マクロコード処理装置１２による計算の各タイムステップの最初で時間刻みと壁面温度とが二相流解析手段１６に送られる。二相流解析手段１６はそのタイムステップの計算を行い、解析結果ｄを変換手段１７に出力する。変換手段１７は、二相流解析手段１６の解析結果ｄをマクロコード処理装置１２が用いることができるサブクール沸騰特性量ｅに変換し、マクロコード処理装置１２に送る。マクロコード処理装置１２は受け取ったサブクール沸騰特性量ｅを用いて全体系の解析を行い次のタイムステップに進む。
【００１６】
図２は解析対象である燃料集合体の説明図であり、図２（ａ）は燃料集合体の斜視図、図２（ｂ）は燃料集合体の一部切り欠き断面図である。沸騰水型原子炉においては燃料集合体は、図２（ａ）に示すように、正方形断面を有するチャンネルボックス（金属ダクト）１８内にスペーサ１９を介在して燃料棒２０が矩形格子で納められている。燃料集合体内の熱流動を解析する際には、この断面を図２（ｂ）のようにメッシュ（セル）に分割し、各セルにおいて流体がどのような流動様式をしているかをボイド率等に基づき決定する。これにより、流体と燃料棒との接触面積が評価でき、伝熱と摩擦とが計算できることになる。また、伝熱による相変化に伴う質量の移行と摩擦が計算できると流体の運動方程式を解くことが可能となる。
【００１７】
従来においては、ボイド率は実験相関式により求め、また、伝熱と摩擦は各々実験相関式に基づき決定している。また、熱伝達係数や摩擦係数は実験相関式に基づき決定し、セル間の流体の移動量は、圧力差や乱流効果等を考慮した実験相関式で決定している。
【００１８】
そして、燃料棒間に形成される冷却材の流路をサブチャンネルと呼んでおり、このサブチャンネルについて上述の計算を行う計算プログラムを通常サブチャンネル解析コードと呼んでいる。また、解析領域をメッシュに分割して流動様式等を実験相関式から持ち込んで気液二相流を解析する計算コードのことを本発明ではマクロコードと呼んでおり、サブチャンネル解析コードはマクロコードの一つである。
【００１９】
図３は解析対象であるサブチャンネルにおける解析領域の説明図である。ここでは、流体は垂直に上昇しており、液体（水）から気体（蒸気）に相変化をする場合を示している。すなわち、水単相領域、サブクール沸騰領域、核沸騰領域、強制対流蒸発領域、ポストドライアウト領域、過熱蒸気単相領域と相変化する場合を示している。マクロコードでは燃料集合体に形成される全てのサブチャンネルを対象としており、垂直方向においては図３の全体がマクロコードの解析対象である。これに対して、本発明のサブクール沸騰特性量予測装置１１が用いる二相流解析手段が解析対象とするのは、１つのサブチャンネルの解析領域ｂである。図３における横線はマクロコードにより垂直方向に分割されたマクロセル（マクロ計算セル）を示している。サブクール沸騰特性量予測装置１１は、複数のサブチャンネルの解析領域ｂをそれぞれ二相流解析手段を用いて解析し、それぞれのサブチャンネルのマクロセルにサブクール沸騰特性量を送っている。このように二相流解析手段により解析する解析領域を１つのサブチャンネルの解析領域ｂとしたのは、この解析領域ｂにおいてマクロコードが用いている構成式により求めたボイド率と実際のサブチャンネルのボイド率の違いが大きく、また、二相流解析手段の計算時間を短くするために解析領域を必要な範囲に小さくすることが望ましいからである。なお、複数のサブチャンネルを二相流解析手段の解析対象とすることも可能である。しかし、垂直方向の解析範囲を解析領域ｂとした場合には、この垂直方向の範囲においてサブチャンネル間を行き来する流体により垂直方向の流量の変化はわずかであり、その流量変化による二相流解析手段の結果から得られるボイド率などの変化も少ない。それに対して、二相流解析手段で１つのサブチャンネルの解析を複数行う場合に比較して、複数のサブチャンネルを解析する場合の方が著しく解析時間を要する。そのため、サブクール沸騰特性量予測装置１１は、複数のサブチャンネルについて、二相流解析手段を用いて１つのサブチャンネルの解析領域ｂの解析を複数行うことが望ましい。
【００２０】
ここで、水がサブクール状態にあるというのは、飽和状態の温度よりも低い状態にあることを意味し、サブクール沸騰とは、流路断面内に蒸気と過冷液が共存する状態における沸騰のことをいう。
【００２１】
図４は、マクロコード処理装置１２の典型的な計算フローを示すフローチャートである。マクロコード処理装置１２による従来の手法では、図４の計算フローにおける「流動様式の設定」において構成式として実験相関式を用いて正味の蒸気発生点を含むマクロセルの同定が行われる。
【００２２】
一方、本発明による手法では、この「流動様式の設定」部分において、サブクール沸騰特性量予測装置１１により二相流解析によりサブクール沸騰特性量を求めるマクロセルについては、実験相関式を適用するのではなく、サブクール沸騰特性量予測装置１１からのサブクール沸騰特性量を用いる。
【００２３】
すなわち、マクロコード処理装置１２は解析対象からプロセス量を読み込み（Ｓ１）解析を行うが、まず各時間ステップの計算を行う前にマクロコード処理装置１２がこれから行う時刻の解析処理情報ｃをサブクール沸騰特性量予測装置１１に出力し、サブクール沸騰特性量予測装置１１からこれからマクロコード処理装置１２が解析を行う時刻のサブクール沸騰特性量ｅを入力する。次に、最も上流から軸方向の１層分のマクロセルの計算を行う。まず、サブクール沸騰特性量予測装置１１からサブクール沸騰特性量ｅが入力されたマクロセルについてはサブクール沸騰特性量ｅを用いて、その他のマクロセルについては実験相関式による構成式を用いて流動様式の設定（Ｓ２）、境界面積の計算（Ｓ３）、相変化量の計算および成分間抵抗の計算を行う（Ｓ４）。そして、質量保存式、運動量保存式、エネルギー保存式の計算を行い（Ｓ５）、計算結果の印刷・保存を行う（Ｓ６）。１層分の計算が終了すると、その層のサブチャンネル間のクロスフローが計算され、その下流（上方）の入口流量が計算できる。そこで、下流に向かって、軸方向の計算ループを行い（Ｓ７）、軸方向の計算ループが完了すると、次のタイムステップに進み（Ｓ８）、時間計算ループを行う（Ｓ９）。
【００２４】
次に、サブクール沸騰特性量予測装置１１の初期設定手段１４での解析領域の設定について説明する。ここで、マクロコード処理装置１２からサブクール沸騰特性量予測装置１１へ、セル情報とマクロコード処理装置１２によって予想されたボイド率の分布とが送られるが、ここでは円管の解析を考えているためセル情報とはマクロセルの境界の垂直方向の位置のみである。
【００２５】
図５は、サブクール沸騰特性量予測装置１１の初期設定手段１４での解析領域の選択の仕方の説明図であり、図５（ａ）はマクロコード処理装置１２での解析対象の状態の説明図、図５（ｂ）はマクロコード処理装置１２による解析を行う際に設定される解析セルの説明図、図５（ｃ）は初期設定手段１４で選択された二相流解析手段１６での解析領域の説明図である。
【００２６】
図５（ｃ）に示すように、解析対象セルの選択は、予想ボイド率がゼロとなる位置を含むマクロセルより下部の２セルからボイド率が０．２となる位置を含むセルまでとする。選択されたセルのセル番号を解析領域ｂとしてサブクール沸騰特性量予測装置１１からマクロコード処理装置１２に送られる。
【００２７】
サブクール沸騰特性量予測装置１１の初期設定手段１４での初期設定においては、その解析領域に二相流解析のための粒子を初期配置する。図６は二相流解析手段１６が解析対象とする解析セルと処理粒子配置の説明図であり、図６（ａ）は二相流解析手段１６で解析対象とするために選択された領域の説明図、図６（ｂ）は粒子の初期配置と入口境界条件の設定の説明図である。図６（ａ）において、ｕｌ（Ｉ）は液体速度、ｕｇ（Ｉ）は気体速度、αｇ（Ｉ）はボイド率、Ｔｌ（Ｉ）は液体温度、Ｔｇ（Ｉ）は気体温度、Ｐ（Ｉ）は圧力、Ｉはセル番号である。図６（ｂ）において、ｕｌ，ｉは液体粒子ｉの速度ベクトル、Ｔｌ，ｉは液体粒子ｉの温度、Ｐｌ，ｉは液体粒子ｉの圧力、Ｔｓ，ｊは壁粒子ｊの温度、ｉ，ｊは粒子番号である。このように初期設定では、均一の液体粒子を壁粒子で包囲して粒子の初期配置を行う。
【００２８】
図７はサブクール沸騰特性量予測装置１１で二相流解析結果を得るまでの処理内容を示すフローチャートである。
【００２９】
図８は、二相流解析の処理内容を示すフローチャートである。このような二相流解析の処理を行う二相流解析方法は、通常、計算機上に二相流解析プログラムをインストールして実行される。
【００３０】
二相流解析の処理では、まず、解析の初期状態を設定する（Ｓｂ１）。次に表面張力や運動量の拡散による液体粒子等の位置等の計算や、液体粒子等への熱の移動等を計算する陽的計算を行う（Ｓｂ２）。次に圧力補正により液体粒子等の位置等の計算する陰的計算を行い（Ｓｂ３）、この陰的計算（Ｓｂ３）は、液体粒子等の位置が誤差範囲内に収束するまで繰り返し計算する（Ｓｂ４）。計算が収束した場合は、二相流解析結果の出力を行い（Ｓｂ５）、次の二相流解析を行う時刻までの時間刻みを計算する（Ｓｂ６）。
【００３１】
そして、液体粒子等の位置および速度を更新する（Ｓｂ７）。この更新により、二相流解析計算後の液体粒子等の位置・速度に対して、Ｓｂ６で計算された時間刻み後である次の時刻の二相流解析での陽的計算ステップ（Ｓｂ２）の計算が行われることになる。二相流解析を行う必要がある時刻の解析が行われたかの計算終了判定を行い（Ｓｂ８）、必要な時刻までの計算が終了している場合には、二相流解析の処理を終了し、必要な時刻までの計算が終了していない場合には、二相流解析の処理を継続し、陽的計算（Ｓｂ２）、陰的計算（Ｓｂ３）等の計算を繰り返し行う。二相流解析の開始時刻から終了時刻までの途中の時刻の解析結果が不要の場合には、二相流解析結果の出力（Ｓｂ５）のステップを計算終了判定（Ｓｂ７）ステップの後に行い、終了時刻での解析結果のみを出力するようにしても良い。
【００３２】
図９は、解析の初期状態を設定ステップ（Ｓｂ１）の詳細ステップである。
【００３３】
まず、パラメータの設定ステップ（Ｓｃ１）では、以下の計算を行う。
【００３４】
カーネル径係数ｒｅ０および初期粒子間隔Δを用いて、カーネル径ｒｅを次式より計算する。
【００３５】
2005050042.tif 000003
次に、初期状態設定ステップ（Ｓｃ２）においては、以下の処理を行う。
【００３６】
図６に示すように、二相流解析の処理を行う領域、解析条件および初期粒子間隔Δから液体粒子の粒子数Ｎｌ、気体粒子の粒子数Ｎｇ、構造粒子の粒子数Ｎｓを求める。この時、構造粒子は、壁などの構造体や境界条件を二相流解析で取扱うために用いられる粒子である。
【００３７】
また、図６に示すように二相流解析の処理を行う領域、解析条件および初期粒子間隔Δから、液体粒子および気体粒子の座標ベクトルｒξ，ｉ、速度ベクトルｕξ，ｉ、温度Ｔξ，ｉ、密度ρξ，ｉ、圧力Ｐξ，ｉ、比内部エネルギーの値ｅξ，ｉ、カーネル径ｒｅξ，ｉの初期値を設定する。カーネル径ｒｅ，ξ，ｉには（９−１）で計算したｒｅを入力する。添え字ξは液体粒子と気体粒子の識別子である。したがって、添え字ξの部分は液体粒子の場合にはｌ、気体粒子の場合にはｇであり、液体粒子の座標ベクトルであればｒｌ，ｉとなる。このような表記は以下の説明でも同様である。また、添え字ｊは個々の粒子の引数である。また計算機プログラムでの変数として取扱う場合には、液体粒子と気体粒子を引数で識別することができ、例えば液体粒子の引数を１、気体粒子の引数を２とすることができる。具体的には、液体粒子ｉのｘ座標の変数は、ｒｘ（１，ｉ）とすることができる。このような処理は上記した全ての変数および後に説明される変数の全てで、計算機プログラムの変数の記述として用いることができる。
【００３８】
また、図６に示すように二相流解析の処理を行う領域、解析条件および初期粒子間隔Δから、構造粒子の座標ベクトルｒｓ，ｉ、速度ベクトルｕｓ，ｉ、温度Ｔ  ｓ，ｉ、圧力Ｐ  ｓ，ｉを設定する。
【００３９】
次に初期数密度の計算ステップ（Ｓｃ３）では次の処理を行う。
【００４０】
まず、初期粒子数密度を次式により計算する。
【００４１】
2005050042.tif 000004
具体的には、次のようにして求めることができる。
【００４２】
二相流解析の液体粒子および気体粒子とは別に、粒子を例えば座標系の原点に想定し、その粒子の周りに半径がカーネル径ｒｅ長さの空間を想定し、初期粒子間隔Δで粒子を配置し、最初に想定した粒子との関係でその各粒子の座標位置を計算する。そして、最初に想定した粒子をｉとして半径ｒｅの空間に想定された全ての粒子ｊについてｗ（｜ｒｊ−ｒｉ｜）を足し合わせることにより求める。
【００４３】
粒子数密度に対する自分からの寄与ｎｓｅｌｆを設定する。本実施の形態においては、次式により計算する。
【００４４】
2005050042.tif 000005
ｎｓｅｌｆは、粒子数密度に対する自分からの寄与を考慮するものであり、（９−７）によらず、粒子数密度に対する自分からの寄与が適切に考慮されるように設定すれば良い。
【００４５】
2005050042.tif 000006
次に近傍粒子リストの作成ステップ（Ｓｃ４）において、液体粒子および気体粒子の各粒子について、全ての粒子との間の距離を計算し、カーネル径ｒｅξ，ｉ以内にある液体粒子の数およびその液体粒子の粒子番号のリスト、カーネル径ｒｅξ，ｉ以内にある気体粒子の数およびその気体粒子の粒子番号のリスト、カーネル径ｒｅξ，ｉ以内にある構造粒子の数およびその構造粒子の粒子番号のリストを作成し、記憶する。作成された液体粒子の粒子番号リスト等を近傍粒子リストＴｂｌと呼ぶことにする。
【００４６】
次に粒子数密度の計算ステップ（Ｓｃ５）では、液体粒子および気体粒子の座標ベクトル初期値ｒξ，ｉ、近傍粒子リストの作成ステップ（Ｓｃ４）で作成した近傍粒子リストＴｂｌを用いて、液体粒子および気体粒子の各粒子について液体粒子および気体粒子の粒子密度ｎξ，ｉと構造体の粒子密度ｎｓ，ｉを求める。
【００４７】
具体的には、例えば液体粒子ｉの気体粒子密度ｎｇ，ｉでは、カーネル径ｒｅξ，ｉ以内にある全ての気体粒子ｊに対して次式のように加算して、
2005050042.tif 000007
と求められる。この時（９−３）式のカーネル径ｒｅにはカーネル径ｒｅξ，ｉの値を入力して計算する。以下の計算で同様である。
【００４８】
また、例えば液体粒子ｉの液体粒子密度ｎｌ，ｉでは、カーネル径ｒｅｌ，ｉ以内にあるｉを除く全ての液体粒子ｊに対して、
2005050042.tif 000008
と求められる。
【００４９】
次に、粒子数密度の計算ステップ（Ｓｃ５）の計算が全て終了すると、初期設定ステップ（Ｓｂ１）は終了する。
【００５０】
図１０は、解析の陽的計算ステップ（Ｓｂ２）の詳細ステップである。
【００５１】
まず、仮想次元方向の位置計算ステップ（Ｓｂ１）を行う。この仮想次元方向の位置計算ステップ（Ｓｂ１）は、液体が蒸発した場合を解析する際に発生する計算不安定性を解決するためのステップである。二相流解析では図１１（ａ）に示すように液体の蒸発は、液体粒子から気体粒子（蒸気粒子）１５を新に発生させることにより模擬している。この時、新に発生させた気体粒子により。局所的に粒子数密度が非常に大きな値となる。この非常に大きな粒子数密度により、収束計算の初期のステップで気体粒子１５は非常に大きな圧力を持つ。このように収束計算の初期のステップで非常に大きな圧力を持つ粒子が発生することにより計算不安定性を生じている。
【００５２】
そこで、本実施の形態では、図１１（ｂ）に示すように、物理空間次元に新たな次元ηを付加し、その次元η方向から気体粒子（蒸発粒子）１６を挿入させる。これにより、新に発生した気体粒子（蒸発粒子）ｉと周りの粒子ｊとの距離は次式のように求めることができる。
【００５３】
2005050042.tif 000009
そして、ηｇ，ｌには予め値を与えておき、次式によりその値を速度Ｃにより減少させれば良い。
【００５４】
2005050042.tif 000010
本実施の形態では、蒸発により新に発生した気体粒子以外の全ての粒子において次元ηの位置データを持たせている。このため、新に発生した気体粒子（蒸発粒子）ｉと周りの粒子ｊとの距離は次式で求めることになる。
【００５５】
2005050042.tif 000011
新に発生した気体粒子以外の粒子は、実空間の距離と一致するようにηξ，ｉは０である。したがって、本実施の形態の仮想次元方向の位置計算（Ｓｄ１）では、ηｇ，ｉを含むηξ，ｉについて、ηξ，ｉが０でない場合に次式の計算をし、結果が負の場合には、ηｎ＋１ξ，ｉに０を入力するようにしている。
【００５６】
2005050042.tif 000012
なお、ηξ，ｉについて各時刻ごとの値を記憶しておく必要はないので、その都度、（１１−４）式によりηξ，ｉの値を更新している。
【００５７】
これにより、新たに生じた気体粒子か否かの区別をすることなく、仮想次元方向の位置計算ステップ（Ｓｄ１）の処理を行うことができ、処理が簡便化される。
【００５８】
また、計算開始時のηξ，ｉの初期値は、全ての粒子について実空間と一致するように０を入力している。
【００５９】
また、速度ｃの値は、計算不安定性が生じない最大の値以下の値を用いることができる。そして、新に発生した気体粒子１６のηｇ，ｉ（ηξ，ｉ）に入力する値としては、例えばカーネル径ｒｅを用いれば良い。
【００６０】
次に境界条件の更新計算ステップ（Ｓｄ２）は次の処理を行う。
【００６１】
まず、構造粒子の座標ベクトルｒｓ，ｉの更新計算を行う。壁などの流路の移動を構造粒子の座標ベクトルｒｓ，ｉで求めたり、液体の流入する境界条件や流出する境界条件を計算する。定常的に移動する場合には、次式により求めることができる。
【００６２】
2005050042.tif 000013
この時、液体の流入は、図１２に示すように液体の流入する境界に配置された構造粒子を移動して行う。（図１２（ａ）から図１２（ｂ）へ移動させる。）そして、この入口境界に配置された構造粒子が、液体粒子が挿入できるだけ内側（上方）へ移動した場合には、構造粒子があった場所に新に液体粒子を追加する。液体粒子の追加は、液体粒子の最大の粒子番号の次の番号を追加し、液体粒子の粒子数Ｎｌを増加させる。そして、追加された粒子番号ｊの座標ベクトルｒｌ，ｊに構造粒子の座標ベクトルｒｓ，ｉの値を入力すれば良い。その後、構造粒子の座標ベクトルｒｌ，ｊを追加した液体粒子の下の位置へ移動させる（通常は、構造粒子の座標ベクトルｒｌ，ｊを初期値に戻せば良い）。
【００６３】
液体が流出す境界においては、上記と逆の操作を行えばよい。すなわち、出口境界に配置された構造粒子が、出口側境界の端まで移動した場合、この構造粒子に最も近い液体粒子ｊを削除し、構造粒子の座標ベクトルｒｓ，ｉの値を削除した液体粒子ｊの座標ベクトルｒｌ，ｊの値とする。液体粒子の削除は、削除する液体粒子の粒子番号を記憶しておく。この時液体粒子の粒子数Ｎｌは減少させない。そして、全ての構造粒子の座標ベクトルの更新およびそれに伴う粒子の追加および削除すべき粒子番号の記憶が完了した後に、図１３に示すフローにより行う。最後に削除した粒子番号の数だけ液体粒子の粒子数Ｎｌを減少させる。
【００６４】
次に近傍粒子リストの作成ステップ（Ｓｄ３）を行いその後に粒子数密度の計算ステップ（Ｓｄ４）を行う。近傍粒子リストの作成ステップ（Ｓｄ３）および粒子数密度の計算ステップ（Ｓｄ４）は、それぞれＳｃ４ステップおよびＳｃ５ステップと同じである。
【００６５】
次に表面張力項の計算ステップ（Ｓｄ５）で液体粒子に働く表面張力の計算を行う。
【００６６】
2005050042.tif 000014
この時、初期粒子間隔Δ以内を近距離、初期粒子間隔Δより遠い場合を遠距離としている。
【００６７】
ポテンシャルパラメータを決定するための実験解析を行った結果、式（１３−１）では大きな曲率をもつ配位において曲率を過小評価することが分かった。そこで、次式のように補正係数により補正する。
【００６８】
2005050042.tif 000015
ここで、Ｎｉｌは液体粒子ｉからカーネル径ｒｅｌ，ｉ以内にある液体粒子の数、Ｎｉ０は液体粒子ｉからカーネル径ｒｅｌ，ｉ以内にある液体粒子および気体粒子を含めた粒子の数である。なお、液体粒子ｉからカーネル径ｒｅｌ，ｉ以内にある粒子の粒子番号は、近傍粒子リストＴｂｌから得る。本実施の形態では、Ｎｉ０は全ての液体粒子ｉについて、半径がカーネル径ｒｅｌ，ｉ以内に初期粒子間隔Δで粒子を配置した時の全粒子数を用いている。
【００６９】
次に、運動量拡散項の計算ステップ（Ｓｄ６）では、気体粒子および液体粒子に働く運動量拡散項を次の処理により計算している。
【００７０】
液体粒子および気体粒子の運動量拡散項Ｄξ，ｉは、液体および気体それぞれの拡散係数νξ等を用いて、近傍粒子リストＴｂｌにより得られるカーネル径ｒｅｌ，ｉ以内にある粒子番号の液体粒子および気体粒子について次式のように加算して求める。
【００７１】
2005050042.tif 000016
次に、位置速度の陽的仮更新ステップ（Ｓｄ７）では、液体粒子および気体粒子の座標ベクトルおよび速度ベクトルについての陽的な更新値を次のように計算する。
【００７２】
2005050042.tif 000017
この時、ｕξ，ｉ、ｒξ，ｉを書き換えることなく、更新値ｕ＊ξ，ｉ、ｒ＊ξ，ｉをそれぞれの粒子について別に記憶する。
【００７３】
次に、熱および相変化計算ステップ（Ｓｄ８）では、各粒子への熱の移動や蒸発や凝縮などの相変化、および各粒子の比内部エネルギーと温度、および温度により変化する拡散係数νξ、熱伝導率κξを計算している。
【００７４】
図１４は、相変化計算ステップ（Ｓｄ８）の詳細ステップである。
【００７５】
まず、液体粒子間の熱伝導計算ステップ（Ｓｈ１）では、液体粒子間の熱伝導計算について次のような処理を行う。
【００７６】
液体粒子の比内部エネルギーｅｌ，ｉを、近傍粒子リストＴｂｌにより得られるカーネル径ｒｅｌ，ｉ以内にある粒子番号の液体粒子について次式のように加算して求める。
【００７７】
2005050042.tif 000018
次に、気体粒子間の熱伝導計算ステップ（Ｓｈ２）では、気体粒子間の熱伝導計算について次のような処理を行う。
【００７８】
気体粒子の比内部エネルギーｅｇ，ｉを、近傍粒子リストＴｂｌにより得られるカーネル径ｒｅｇ，ｉ以内にある粒子番号の気体粒子について次式のように加算して求める。
【００７９】
2005050042.tif 000019
次に液体粒子ｉの界面から構造粒子ｊに移行する熱移行率ｑｌｆ→ｓと液体粒子ｉから液体粒子ｉの界面に移行する熱移行率ｑｌｆ→ｌを次の式より計算する。（Ｓｉ４）
2005050042.tif 000020
この時、ａ０は単位体積当りの伝熱面積係数、ｈｓは構造体の熱伝達率、ｈｌは液体の熱伝達率である。これら熱伝達率は熱伝導率からｈｓ  ＝κｓ／Δ、ｈｌ  ＝κｌ／Δで評価している。
【００８０】
次に、これらを足し合わせて構造粒子ｊから液体粒子ｉへの熱伝達にともない液体粒子ｉの界面から移行する熱移行率ｑｌｆ→を計算する。（Ｓｉ５）
2005050042.tif 000021
ｑｌｆ→が０以上の場合には、構造粒子ｊから液体粒子ｉへの熱伝達においては液体粒子ｉの界面に熱が蓄積されないので、構造粒子ｊから液体粒子ｊへの熱伝達量を計算する。まず、界面温度を次式で求める（Ｓｉ８）。
【００８１】
2005050042.tif 000022
次に、（１５−２）式または（１５−６）式により求められた液体粒子ｉの界面から液体粒子ｉへの熱移行率が構造粒子ｊから液体粒子ｉへ熱伝達される量であるので、次式により液体粒子ｉの比内部エネルギーを計算する（Ｓｉ１０）。
【００８２】
2005050042.tif 000023
全ての構造粒子について計算が終了するまで、Ｓｉ２からＳｉ１０までのステップを繰り返す。全ての構造粒子について計算が終了した場合には、次の粒子番号の液体粒子について計算を行うが、同じ粒子番号の液体粒子ｉについて、気体粒子から液体粒子への熱伝達量を続けて計算した後に、次の粒子番号について計算する様にしても良い。
【００８３】
2005050042.tif 000024
この時、ｈｇは気体の熱伝達率、ｈｌは液体の熱伝達率である。これら気液間の熱伝達率は熱伝導率からｈｇ＝κｇ／Δ１ｇ、ｈｌ＝κｌ／Δ１ｇで評価している。このときΔ１ｇ  ＝１．５×１０−６  （ｍ）を用いている、これは、２  ｈｌ  ｈｇ／（  ｈｌ  ＋  ｈｇ）≒４．５×（Ｗ／ｍ２／Ｋ）となるように決めた。
【００８４】
次に、これらを足し合わせて気体粒子ｊから液体粒子ｉへの熱伝達にともない液体粒子ｉの界面から移行する熱移行率ｑｌｆ→を計算する。（Ｓｊ５）
2005050042.tif 000025
全ての気体粒子について計算が終了するまで、Ｓｊ２からＳｊ８までのステップを繰り返す（Ｓｊ９）。
【００８５】
その後、全ての液体粒子について計算を繰り返して、全ての液体粒子について構造粒子から液体粒子への熱伝達量および気体粒子から液体粒子への熱伝達量が計算終了したら、液体粒子への熱伝達計算ステップ（Ｓｈ３）が終了する。
【００８６】
2005050042.tif 000026
気体粒子への熱伝達計算ステップ（Ｓｈ４）では、次のような処理を行っている。
【００８７】
液体粒子から気体粒子ｉへの熱伝達の詳細ステップを図１７に示す。
【００８８】
まず、気体粒子の粒子番号を設定する（Ｓｋ１）。次に液体粒子の粒子番号を設定する（Ｓｋ２）。
【００８９】
2005050042.tif 000027
次に気体粒子ｉから気体粒子ｉの界面に移行する熱移行率ｑｇｆ→ｇと気体粒子ｉの界面から液体粒子ｊに移行する熱移行率ｑｇｆ→ｌを次の式より計算する（Ｓｋ４）。
【００９０】
2005050042.tif 000028
次に、これらを足し合わせて液体粒子ｊから気体粒子ｉへの熱伝達にともない気体粒子ｉの界面から移行する熱移行率ｑｇｆ→を計算する。（Ｓｋ５）
2005050042.tif 000029
そして、（１７−２）式で求められる気体粒子ｉの界面から気体粒子ｉへ移行する熱移行率ｑｇｆ→ｌが気体粒子ｊから気体粒子ｉへ熱伝達される量であるので、次式により気体粒子ｉの比内部エネルギーを計算する（Ｓｋ８）。
【００９１】
2005050042.tif 000030
全ての液体粒子について計算が終了するまで、Ｓｋ２からＳｋ８までのステップを繰り返す（Ｓｋ９）。
【００９２】
その後、全ての気体粒子について計算を繰り返して、全ての気体粒子について液体粒子から気体粒子への熱伝達量が計算終了したら、気体粒子への熱伝達計算ステップ（Ｓｈ４）が終了する。
【００９３】
次に粒子発生の計算ステップ（Ｓｈ５）では、液体粒子からの新たな粒子の発生について計算する。
【００９４】
粒子発生の計算ステップ（Ｓｈ５）での液体粒子から蒸発により気体粒子の発生の詳細ステップについて図１８に示す。
【００９５】
まず、液体粒子の粒子番号を設定する（Ｓｌ１）。次に水の比エンタルピｉｌ，ｉを次式で計算する（Ｓ１２）。
【００９６】
2005050042.tif 000031
この時、ρｇは気体の密度である。Ｎｌ→ｇが１以上か否かを判断し（Ｓｌ５）、１より大きい場合には新に気体粒子を発生させる（Ｓｌ６）。この時、発生させる気体粒子の粒子番号は最も大きな気体粒子の次の番号とし、温度Ｔｇ，ｊは液体粒子ｉの圧力Ｐｌ，ｉの飽和温度とし、座標ベクトルは仮想のη次元をカーネル径ｒｅとし残りを液体粒子ｉと同じとする。
【００９７】
2005050042.tif 000032
全ての液体粒子についてＳｌ１からＳｌ８の計算が終了したら、粒子発生の計算ステップ（Ｓｈ５）は終了する。
【００９８】
次に粒子消滅の計算ステップ（Ｓｈ６）では、気体粒子の凝縮による消滅について計算する。
【００９９】
粒子消滅の計算ステップ（Ｓｈ６）の詳細ステップについて図１９に示す。
【０１００】
まず、気体粒子の粒子番号を設定する（Ｓｍ１）。
【０１０１】
2005050042.tif 000033
Ｎｇ→ｌが１より大きい場合には、気体粒子ｉを消滅させる（Ｓｍ５）。この時、各気体粒子ｉの計算時には、消滅させる気体粒子ｉの粒子番号を記憶しておき、全ての気体粒子についてＳｍ１からＳｍ６の計算が終了した後に、図１３と同様にステップ気体粒子の座標ベクトル等の更新と削除した気体粒子の数だけの気体粒子の粒子数Ｎｇを減少させる。この処理が終了することによって粒子消滅の計算ステップ（Ｓｈ６）が終了する。
【０１０２】
次に、熱的状態の計算ステップ（Ｓｈ７）では、気体粒子および液体粒子の比内部エネルギーｅξ，ｉや圧力Ｐξ，ｉから流体の状態表や、その近似式を用いて、気体粒子および液体粒子の温度Ｔξ，ｉ、気体および液体の拡散係数νξおよび熱伝導率ｈξを計算する。
【０１０３】
陰的計算ステップ（Ｓｂ３）の詳細ステップについて図２０に示す。
【０１０４】
このステップは次式の圧力方程式を解いて得られた圧力から、粒子の座標ベクトル、速度ベクトルの補正を行うものである。
【０１０５】
2005050042.tif 000034
この時、図２１（ａ）のように気体粒子および液体粒子を用いて（２０−０）をといた場合には、気体粒子と液体粒子の大きな密度の差より、計算不安定性を生じる。そこで、気体粒子の計算（図２１（ｂ）参照）と液体粒子の計算（図２１（ｃ）参照）を別々に行う。
【０１０６】
まず、全ての液体粒子を気体粒子と仮定して圧力方程式を解くために、気体圧力方程式の行列作成を作成する（Ｓｎ１）。
【０１０７】
気体圧力方程式は次式のようになる。
【０１０８】
2005050042.tif 000035
この時ｉおよびｊは全ての気体粒子、全ての液体粒子、および全ての構造粒子である。
【０１０９】
そして係数行列Ａｉｊの行列要素は次式で計算される。
【０１１０】
2005050042.tif 000036
ここで、Ｎはすべての気体、液体、構造粒子の粒子数の和である。
【０１１１】
この時、ｒ＊ｉ、ｒ＊ｊは気体粒子、液体粒子の座標ベクトルｒ＊ξ，ｉ、および構造粒子の座標ベクトルｒｓ，ｉが用いられる。ただし、圧力境界条件値を設定する粒子の粒子番号をｋとしたとき、関係する行列要素のみ（２０−２）から以下のように変更する。
【０１１２】
2005050042.tif 000037
ＩＣＣＧ法により（２０−１）式を解き、各気体粒子および気体と仮定した液体粒子の圧力を解く（Ｓｎ２）。
【０１１３】
次に、気体粒子に最も近い液体粒子は境界粒子として圧力を（２０−１）を解いて得られた値に固定し、液体圧力方程式の行列式を作成する（Ｓｎ３）
この時、（２０−２）式および（２０−３）式のｉおよびｊは全ての液体粒子と構造粒子であり、ｋは新に気体粒子の圧力を境界条件して設定された液体粒子を含む圧力が境界値として設定されている粒子である。
【０１１４】
次にＩＣＣＧ法により（２０−１）式を解き、液体粒子の圧力を解く（Ｓｎ４）。
ここで、Ｓｎ１およびＳｎ２ステップにより気体粒子の圧力方程式を先に解くのは、蒸発により新に発生した気体粒子については、圧力が求められていないからである。
【０１１５】
次に液体の圧力勾配項の計算ステップ（Ｓｎ５）では、圧力勾配項による液体粒子の速度ベクトルの変化量ｕ’ｌ，ｉを以下のような処理により計算する。
【０１１６】
液体粒子ｉの粒子番号を設定し、まず、次式により気液間圧力補正の過大評価を補正するためのパラメータｍｉを求める。
【０１１７】
2005050042.tif 000038
全ての液体粒子について計算が終了したら、液体の圧力勾配項の計算ステップ（Ｓｎ５）は終了する。
【０１１８】
次に気体の圧力勾配項の計算ステップ（Ｓｎ６）では、圧力勾配項による気体粒子の速度ベクトルの変化量ｕ’ｇ，ｉを以下のような処理により計算する。
【０１１９】
気体粒子ｉの粒子番号を設定し、まず、次式により気液間圧力補正の過大評価を補正するためのパラメータｍｉを（２０−４）式により求める。そして、次式によりｕ’ｇ，ｉを計算する。
【０１２０】
2005050042.tif 000039
座標ベクトル変化量は、液体粒子および気体粒子の速度ベクトル変化量ｕ’ｇ，ｉの計算後すぐに計算するようにして良い。
【０１２１】
次に位置速度の陰的仮更新ステップ（Ｓｎ７）では、次式により座標ベクトルｒ＊ξ，ｉ、速度ベクトルの値ｕ＊ξ，ｉを更新する。
【０１２２】
2005050042.tif 000040
次に粒子数密度の計算ステップ（Ｓｎ８）では、Ｓｄ４ステップと同様に座標ベクトルｒ＊ξ，ｉを用いて気体粒子および液体粒子の粒子密度を計算する。
【０１２３】
次に、陰的計算の収束判定ステップ（Ｓｎ９）では、液体および気体粒子の粒子数密度の初期粒子数密度ｎ０からの誤差を計算する。この陰的計算の収束判定ステップ（Ｓｎ９）が終了すると陰的計算ステップ（Ｓｂ３）は終了する。
【０１２４】
陰的計算ステップ（Ｓｂ３）が終了した後に行われる時間刻みの計算ステップ（Ｓｂ６）では、時間刻みを次式で計算する。
【０１２５】
2005050042.tif 000041
ｃｌはクーラン定数で、ｃｌ＝０．２とする。
【０１２６】
また、位置・速度の更新ステップ（Ｓｂ７）では、次式で粒子の座標ベクトルおよび速度ベクトルを更新する。
【０１２７】
2005050042.tif 000042
次に、サブクール沸騰特性量予測装置１１の変換手段１７について説明する。
図２２は変換手段１７の処理内容を示すフローチャートである。まず、二相流解析手段１６の二相流解析結果ｄを入力し（Ｓ１）、物理量の粗視化を行い（Ｓ２）、サブクール沸騰特性量の計算を行う（Ｓ３）。そして、計算したサブクール沸騰特性量、すなわち、ボイド率、サブクール度、ペクレ数、スタントン数をマクロコード処理装置１２に出力する（Ｓ４）。
【０１２８】
ここで、変換手段１７における物理量の粗視化処理（ステップＳ３）は、二相流解析手段１６による解析結果をマクロセル内部で平均化（粒子平均）する処理である。図２３は粗視化処理の具体的な処理内容を示すフローチャートである。二相流解析手段１６から二相流解析結果を入力し（Ｓ１）、各マクロセルに属する粒子の判定を行い（Ｓ２）、各マクロセル内での粒子平均操作を行い（Ｓ３）、空間平均量を蓄積する（Ｓ４）。この粗視化処理による処理結果を図２４に示す。図２４に示すように、粗視化処理はマクロセル内部で平均化（粒子平均）する処理である。
【０１２９】
また、変換手段１７におけるサブクール沸騰特性量処理（ステップＳ４）は、粗視化された物理量を基にサブクール沸騰特性量を計算する処理である。図２５は、変換手段１７におけるサブクール沸騰特性量処理の具体的な処理内容を示すフローチャートである。粗視化処理で得られた空間平均量を蓄積し（Ｓ１）、まず、各マクロセルでのボイド率を計算する（Ｓ２）。次に、各マクロセルでのサブクールどの計算を行い（Ｓ３）、各マクロセルでのペクレ数およびスタントン数を計算し（Ｓ４）、サブクール沸騰特性量を算出する（Ｓ５）。
【０１３０】
次に、サブクール沸騰特性量予測装置１２の変換手段１７からマクロコードに渡されるサブクール沸騰特性量の一例について以下説明する。まず、この一例でのサブクール沸騰特性量予測装置１２への主な入力を以下に示す。
【０１３１】
解析体系
・垂直上昇流
・内径  ０．００３  ｍ
・流路長さ  ０．２  ｍ
解析条件
・初期粒子間隔  ０．０００５  ｍ
・初期粒子数
流体（水）  １２，７０４
構造体      １１，３５６
・入口流体速度  ０．５，  １．０，  ２．０，  ４．０，  ６．０，  ８．０，  ２０．０，  ４０．０  ｍ／ｓ
・入口流体温度  ４２０  Ｋ
・流体圧力  ５．０Ｅ＋０５  Ｐａ（飽和温度は４２５  Ｋ）
・壁内面温度  ４５０  Ｋ
この一例では、流入する水の流速を変えた解析を実施しているが、これはペクレ数を変化させることに対応する。スタントン数とペクレ数の定義は次式で示される。
【０１３２】
2005050042.tif 000043
ここで、ｑｗ：熱流束、Ｇ：質量流速、Ｔｓ：飽和温度、ＴｌＢ：液体平均温度、ｃｐｌ：液体低圧比熱、ｋｌ：液体熱伝導率、Ｄ：管内径
図２６はペクレ数をパラメータとした場合の空間位置（軸方向位置）とボイド率との関係を表すマップ、図２８はペクレ数をパラメータとした場合の空間位置（軸方向位置）とスタントン数との関係を表すマップ、図２７はペクレ数をパラメータとした場合の空間位置（軸方向位置）とサブクール度との関係を表すマップである。
【０１３３】
また、変換手段１７で変換されたこれらのサブクール沸騰特性量に基づいてサブクール沸騰特性量に関する相関式を演算する相関式作成手段を設け、相関式を得るようにすることも可能である。例えば、図２６に示した各ペクレ数に対するボイド率の軸方向分布において、ボイド率が急激に上昇する軸方向位置を正味の蒸気発生点と定義し、次に、図２８から各ペクレ数に対数正味の蒸気発生点でのスタントン数の値を求め、得られるペクレ数とスタントン数との関係をプロットしてペクレ数とスタントン数との相関式を得る。
【０１３４】
図２９は、相関式作成手段で得たペクレ数とスタントン数との相関式を示している。図２９では、得られた相関式と実験相関式と比較した場合を示している。相関式作成手段は、変換手段１７から出力される結果から、上述したように、図２６のボイド率が急激に上昇する軸方向位置を正味の蒸気発生点と定義し、次に、図２８から各ペクレ数に対数正味の蒸気発生点でのスタントン数の値を求め、得られるペクレ数とスタントン数との関係をプロットしてペクレ数とスタントン数との相関式を得る。
【０１３５】
このように、本発明の実施の形態では、実験相関式を用いない二相流解析手段１６を有したサブクール沸騰特性量予測装置を作成し、このサブクール沸騰特性量予測装置により、マクロ解析コードで必要とされる領域に対して系統的にサブクール沸騰特性量間の相関式、およびサブクール沸騰領域におけるペクレ数、スタントン数と、ボイド率、サブクール度が関連付けられた特性量マップを得る。従って、計算でサブクール沸騰特性量を作成することができるため実験コストを低減できる。
【０１３６】
すなわち、サブクール沸騰特性量予測装置は、通常は計算格子を用いて離散化して扱われる流体を（カーネル関数とも呼ばれる）重み関数を用いて互いの間に相互作用を与えられた粒子を用いて離散化するＭＰＳ（Ｍｏｖｉｎｇ  Ｐａｒｔｉｃｌｅ  Ｓｅｍｉ−Ｉｍｐｌｉｃｉｔ）法に、気液間圧力補正量の高精度化のための計算アルゴリズム、相変化により発生する粒子による計算不安定性を抑制するための計算アルゴリズム、圧力補正項において浮力を導入するための計算アルゴリズム、および液体粒子の表面張力を液体粒子間のポテンシャルによる相互作用を用い、陽的な速度の更新部分に導入するための計算アルゴリズムを講じた二相流解析方法を用いる。そして、サブチャンネル解析コード等のいわゆるマクロコードで必要となるボイド率分布等を、各粒子が有する位置座標における物理量からマクロ解析コードのメッシュ構造を考慮して空間的に平均化した情報を出力する。
【０１３７】
さらに、サブクール沸騰状態を含む体系を解析しようとするマクロコード処理装置１２に対して、流体の空間位置、ペクレ数、ボイド率、スタントン数、サブクール度を出力し、サブクール沸騰に関する実験相関式を用いることなくマクロコードによる解析を可能にする。
【０１３８】
【発明の効果】
以上述べたように、本発明によれば、系統的にサブクール沸騰特性量間の相関式、およびサブクール沸騰領域におけるペクレ数、スタントン数と、ボイド率、サブクール度を得ることができ、マクロコード処理装置による解析の精度向上が可能となる。
【図面の簡単な説明】
【図１】本発明の実施の形態に係わるサブクール沸騰特性量予測装置のブロック構成図。
【図２】本発明の実施の形態における解析対象である燃料集合体の説明図。
【図３】本発明の実施の形態における解析対象である燃料集合体における解析領域の説明図。
【図４】本発明の実施の形態でのマクロコード処理装置の典型的な計算フローを示すフローチャート。
【図５】本発明の実施の形態に係わるサブクール沸騰特性量予測装置の初期設定手段での解析領域の選択の仕方の説明図。
【図６】本発明の実施の形態における二相流解析手段が解析対象とする解析セルと処理粒子配置の説明図。
【図７】本発明の実施の形態に係わるサブクール沸騰特性量予測装置で二相流解析結果を得るまでの処理内容を示すフローチャート。
【図８】本発明の実施の形態における二相流解析の処理内容を示すフローチャート。
【図９】図８における解析の初期状態を設定ステップ（Ｓｂ１）の詳細ステップを示すフローチャート。
【図１０】図８における解析の陽的計算ステップ（Ｓｂ２）の詳細ステップを示すフローチャート。
【図１１】本発明の実施の形態における二相流解析での液体の蒸発についての液体粒子から気体粒子（蒸気粒子）を新に発生させる模擬の説明図。
【図１２】本発明の実施の形態における二相流解析での液体の流入する境界に配置された構造粒子の説明図。
【図１３】本発明の実施の形態における二相流解析での液体の流出する境界での粒子数計算のフローチャート。
【図１４】図８における相変化計算ステップ（Ｓｄ８）の詳細ステップを示すフローチャート。
【図１５】構造粒子から液体粒子ｉへの熱伝達の詳細ステップを示すフローチャート。
【図１６】気体粒子から液体粒子ｉへの熱伝達の詳細ステップを示すフローチャート。
【図１７】液体粒子から気体粒子ｉへの熱伝達の詳細ステップを示すフローチャート。
【図１８】液体粒子から蒸発により気体粒子の発生の詳細ステップを示すフローチャート。
【図１９】粒子消滅の計算ステップ（Ｓｈ６）の詳細ステップを示すフローチャート。
【図２０】陰的計算ステップ（Ｓｂ３）の詳細ステップを示すフローチャート。
【図２１】本発明の実施の形態における二相流解析手段で用いるＭＰＳ法での気液の圧力方程式の取り扱いの説明図。
【図２２】本発明の実施の形態に係わるサブクール沸騰特性量予測装置の変換手段の処理内容を示すフローチャート。
【図２３】本発明の実施の形態における変換手段での粗視化処理の具体的な処理内容を示すフローチャート。
【図２４】本発明の実施の形態における変換手段での粗視化処理による処理結果の説明図。
【図２５】本発明の実施の形態における変換手段におけるサブクール沸騰特性量処理の具体的な処理内容を示すフローチャート。
【図２６】本発明の実施の形態における変換手段の出力により得られたペクレ数をパラメータとした場合の空間位置（軸方向位置）とボイド率との関係を表すマップ。
【図２７】本発明の実施の形態における変換手段の出力により得られたペクレ数をパラメータとした場合の空間位置（軸方向位置）とサブクール度との関係を表すマップ。
【図２８】本発明の実施の形態における変換手段の出力により得られたペクレ数をパラメータとした場合の空間位置（軸方向位置）とスタントン数との関係を表すマップ。
【図２９】本発明の実施の形態における相関式作成手段で得たペクレ数とスタントン数との相関式の説明図。
【符号の説明】
１１…サブクール沸騰特性量予測装置、１２…マクロコード処理装置、１３…解析対象情報入力手段、１４…初期設定手段、１５…解析領域出力手段、１６…二相流解析手段、１７…変換手段、１８…チャンネルボックス、１９…スペーサ、２０…燃料棒、２１…構造体粒子、２２…構造体粒子、２３…液体粒子、２４…気体粒子
【図１】
2005050042.tif 000044
【図２】
2005050042.tif 000045
【図３】
2005050042.tif 000046
【図４】
2005050042.tif 000047
【図５】
2005050042.tif 000048
【図６】
2005050042.tif 000049
【図７】
2005050042.tif 000050
【図８】
2005050042.tif 000051
【図９】
2005050042.tif 000052
【図１０】
2005050042.tif 000053
【図１１】
2005050042.tif 000054
【図１２】
2005050042.tif 000055
【図１３】
2005050042.tif 000056
【図１４】
2005050042.tif 000057
【図１５】
2005050042.tif 000058
【図１６】
2005050042.tif 000059
【図１７】
2005050042.tif 000060
【図１８】
2005050042.tif 000061
【図１９】
2005050042.tif 000062
【図２０】
2005050042.tif 000063
【図２１】
2005050042.tif 000064
【図２２】
2005050042.tif 000065
【図２３】
2005050042.tif 000066
【図２４】
2005050042.tif 000067
【図２５】
2005050042.tif 000068
【図２６】
2005050042.tif 000069
【図２７】
2005050042.tif 000070
【図２８】
2005050042.tif 000071
【図２９】
2005050042.tif 000072
【手続補正書】
【提出日】平成15年12月15日(2003.12.15)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】０００３
【補正方法】変更
【補正の内容】
【０００３】
サブチャンネル解析コード等のマクロコードを用いてサブクール沸騰領域に関する解析を実施する場合、ペクレ数、スタントン数と、ボイド率、サブクール度等のサブクール沸騰特性量に関する相関式が必要となる。そこで、必要となる相関式は実験結果から作成するようにしている。また、解析によりサブクール沸騰特性量に関する相関式を得るには、流体解析対象の空間に計算格子を設定し流体解析を行うことになる。
【非特許文献１】
Koshizuka, S., Tamako H. and Oka, Y.: A Particle Method for Incompressible Viscous Flow with Fluid Fragmentation, Comput. Fluid Dynamics J., Vol.4, pp.29-46, (1995). 発行国：日本,発行所：International Society of CFD
【非特許文献２】
Koshizuka, S. and Oka, Y.: Moving-Particle Semi-Implicit Method for Fragmentation of Incompressible Fluid; Nucl. Sci. and Eng., Vol.123, July,pp.421-434, (1996).発行国：米国,発行所：ANS (American Nuclear Society)