JP4370372B2 - キャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置 - Google Patents

Description

本発明は、管内、流体機械、物体まわりの流体に発生するキャビテーションを気泡径分布モデルを用いて、精度良く安定的に解析するのに適したキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置に関する。

ロケットエンジン等に用いられる液体酸素、液体水素等のポンプにおいて、キャビテーション、即ち、飽和蒸気圧よりも局所的に圧力が下がってしまい部分的に気化する現象が発生している。

液体中に気体が混合することにより、音速が著しく低下してしまい、圧力情報が伝わりにくくなる。このことにより、気体が混合していない部分はポンプにより吸い込みが可能であり、気体が混合している部分は吸い込めなくなる、いわゆる脈動現象が発生してしまう。この脈動現象により、圧力振動が羽根を振動させ、破壊させるおそれがある。脈動現象による破壊を防止するためには、キャビテーションの挙動を正確に捉える必要がある。

従来、管内、流体機械及び物体回りに発生するキャビテーションにおいて、どの程度圧力が下がり、それにより、どの程度気化し、どの程度音速が低下し、どの程度圧力が伝わりにくいのかを解析することは非常に困難とされていた。

そこで、液体中に発生するキャビテーションを球形気泡から構成される気泡群として取り扱い、その気泡群の挙動を捉えることでキャビテーションの特徴を把握しようとする技術が提案されている（非特許文献１）。

しかし、一般的な流れ場では、キャビーテーションは、一点から気泡が発生するものではなく、様々な地点から発生した気泡が混在する複雑な現象を伴うものであり、キャビテーションの解析は、精度、安定性に乏しくなってしまったり、計算機の負荷が大きくなる等の問題を抱えている。

Ito, Y., Wakamatsu, H. and Nagasaki, T., Numerical Simulation of Sub-cooled Cavitating Flow by using Bubble Size Distribution, 2003, Proceedings of the Sixth International Symposium on Experimental and Computational Aerothermodynamics of Internal Flows, pp. 280-286 Robert J. Simoneau and Robert C. Hendricks, Two-Phase Chocked Flow of Cryogenic Fluids in Converging-Diverging Nozzles, NASA TP-1484, 1979

本発明の目的は、流体中に発生するキャビテーションを、球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を取り扱う気泡径分布モデルを採用したキャビテーション解析において、その精度、安定性を高めることができ、計算負荷を低減することができるキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置を提供することにある。

この目的を達成するため、本発明に係るキャビテーション解析方法は、少なくとも、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析方法において、ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出する。

また、この目的を達成するため、本発明に係るキャビテーション解析プログラムは、少なくとも、所定のステップをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析プログラムにおいて、上記所定のステップは、ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップである。

また、この目的を達成するため、本発明に係るプログラム記録媒体は、少なくとも所定のステップを有し、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するコンピュータによって読み取り可能にプログラムを記憶したプログラム記録媒体において、上記所定のステップは、ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップである。

また、この目的を達成するため、本発明に係る解析装置は、少なくとも、所定のステップをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析プログラムを記憶した記録部と、上記記録部から上記解析プログラムを読み出して解析処理を実行する演算部とを備える解析装置において、上記所定のステップは、ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップである。

以上のように、本発明は、液体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するものであり、各階級内の気泡数密度分布を算出し、この気泡数密度分布から各階級の境界での気泡数密度を算出し、この気泡数密度から各階級を通過する蒸発量・凝縮量を算出するので、数値拡散を防ぐことができ、解析の精度、安定性を高めることができると共に、高い汎用性を有するものである。また、本発明は、気泡の成長・収縮運動と液体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、解析における時間刻みが小さくなることを防ぐとともに、計算安定性を増すことができる。

本発明によれば、液体圧力、液体温度、各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を考慮することにより、数値拡散を防ぐことができ、気泡の蒸発・凝縮計算の精度を高めることができるとともに、階級間隔を不等間隔にした場合の計算安定性を高めることができる。よって、本発明は、精度及び安定性の高いキャビテーション解析を可能とするものであり、汎用性を高めることを可能とするものである。

また、本発明によれば、気泡の成長・収縮計算と流体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、時間刻みを大きく取ることができ、キャビテーション解析の計算安定性を増すことができ、よって、計算負荷を低減することを可能とするものである。

以下、本発明が適用された解析装置の実施の形態を図面を参照して説明する。

本発明を適用した解析装置１は、図１に示すように、入力部２と、演算部３と、メモリ４と、表示部５と、出力部６と、記憶部７とを備える。この解析装置１では、入力部２を介して入力が行われ、記憶部７に記憶されたキャビテーション解析プログラムからメモリ４に必要な係数及び式を読み出し、演算部３で演算することによりキャビテーション解析を行う。解析結果は、表示部５に表示したり、出力部６から出力することができる。

解析装置１の記憶部７に記憶されたキャビテーション解析プログラムは、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムであって、ステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有するものである。

このキャビテーション解析プログラムを記憶部７に有する解析装置１は、液体中に発生するキャビテーションを球形気泡から構成される気泡群として取り扱い、その気泡群の挙動を捉えることでキャビテーションの特徴と掴もうとするものであり、キャビテーション気泡をその質量に応じて階級化し、各階級ごとに気泡半径、圧力、温度、蒸発・凝縮量、液体とのスリップ速度などを計算することで、キャビテーション気泡を厳密に捉え、解析しようとするものである。

以下に、本発明を適用した解析装置１を用いて、図１３に示すような先細末広ノズルにおける２次元計算の例により、キャビテーション解析を行うキャビテーション解析方法の例を示す。

まず、全体の計算手順に先立ち、「１．基礎式」、「２．物性値」、「３．気泡を含む液体の音速」、「４．気泡新規発生」、「５．気泡質量に基づいた階級の離散化」、「６．気泡の並進運動」、「７．液体の圧力計算と気泡の半径計算のカップリング」及び「８．時間刻み」について説明する。尚、条件としては、Ｓｉｍｏｎｅａｕ ａｎｄ Ｈｅｎｄｒｉｃｋｓの実験条件と同一とする。以下の説明に用いる記号の対応表を表１、表２に示す。

１．基礎式
このキャビテーションモデルを採用した計算コードを開発する。基礎式は２次元の例を示す。一般座標系（ξ，η）を用いて以下のように表せる。単位体積あたりの液相部質量α_Ｌρ_Ｌ、運動量ｘ方向成分α_Ｌρ_Ｌｕ_Ｌ、運動量ｙ方向成分α_Ｌρ_Ｌν_Ｌ、全エネルギーα_Ｌρ_Ｌｅ_Ｌの４式に加え、単位体積あたりの気泡数密度の確率密度関数Ｎ_Ｇｗ、気泡発生からの平均経過時間ｔ_{ｅｌａｐｓｅ}も気泡移流速度で移流するとして計算した。また、気泡関係諸量は気泡の成長・収縮に伴い気泡質量軸方向（ｗ）にも移流するとして計算した。

２．物性値
キャビテーションは液体温度に敏感なため、正確な計算には温度に対する影響を厳密に扱うことが欠かせない。そこで、温度による物性値の変化を考慮した。以下に示す物性値を温度の関数として取り扱った。これらは全て物性値表のデータを連続な関数として近似したものである。これらのデータは基本的には想定温度範囲内（三重点温度〜臨界点温度）のデータを基にしているが、計算中にプログラムにより想定温度範囲外の温度を一時的に参照してしまう場合がある。このときにエラーが発生すると計算の安定性が著しく損なわれるため、想定温度範囲内の物性値と連続な適当な物性値を温度範囲外に対しても与えるようにした。ここで与えられる適当な物性値は計算過程で必要な値であって最終的な収束値には成り得ないため計算の信頼性への影響は少ないと思われる。

気液界面の物性値を表３に、液体の物性値を表４に、蒸気の物性値を表５に示す。また、計算上必要な逆関数を表６に示す。逆関数は安定かつ正確な計算のため元の値に正確に戻す必要がある。そのためこれらの逆関数を必要とする関数に関しては誤差の管理を厳密に行った。

３．気泡を含む液体の音速
液体は気相を内部に含むと圧縮性が単相時に比べて著しく大きくなり、図２に示すように音速が大きく低下する。ここで、図２は、２相流体の平均音速とボイド率の関係を示す図であり、「Akmandor, I. S. and Nagashima, T., Predictions for Cryogenic Homogeneous Two-Phase Flows in a Choked Laval Nozzle, Journal of Thermophysics and Heat Transfer, Vol. 13, No. 3, 355-363」に記載されたものである。高速キャビテーションではこの効果が無視できないため、液体が気泡を内部に含んで音速が低下する挙動をモデルを導入して評価した。そのモデルを以下に示す。

気相は、液相中に均質に分散した気泡であると考える。このとき２相流体の平均音速は次式となる。

ここで、Ｐは２相流体の圧力（ここでは気液間の圧力差は考慮しない）で、ρ_Ｖは２相流の平均密度であり、式（３０）となる。

αはボイド率（気相部分の体積分率）である。一方、クウォリティ（気相部分の質量分率）は次式となる。

式（３０）、式（３１）よりクウォリティχは、式（３２）となる。

ここで、式（３１）をＰについて微分すると式（３３）となる。

よって、式（２９）、式（３３）より式（３４）が得られる。

この式（３４）を再度式（２９）を用いて整理すると次式となる。

本モデルの採用により図２に示すボイド率と音速の関係を得ることができる。

４．気泡新規発生
新規発泡は主として不純物上に発生するものと壁面上に発生するものがある。本モデルでは両初期気泡径とも次式（３６）で示される過飽和度に反比例するものとして与えた。

発泡頻度は不純物上が次式（３７）で示されるｎ_０であり、壁面上においては、次式（３８）で示されるｎ_{０ｗａｌｌ}である。

ここで、式中で用いられる壁面からの単位体積あたりの発泡数ｎ_{Ｓｗａｌｌ}を求める次式（３９）の面積Ａ_ｗａｌｌと壁面近傍の空間Ｖｏｌ_ｗａｌｌの与え方である。

図３に示すような領域で気泡が発生する場合を考える。面積Ａ_ｗａｌｌは壁面隣接セル上の壁面積（２Ｄでは壁長さ）を与える。壁面近傍の空間Ｖｏｌ_ｗａｌｌは式（３６）で求められる初期半径Ｒ_ｏが発生し十分離脱できる領域である必要がある。壁面隣接セルの体積Ｖｏｌ_{ｃｅｌｌ，１}、その１つ壁面から遠いセルの体積Ｖｏｌ_{ｃｅｌｌ，２}、壁面からｊ番目のセル体積Ｖｏｌ_{ｃｅｌｌ，ｊ}とする。また、壁面隣接セルの壁面から一番遠い点までの距離ｄ_ｊ＝１、壁面からｊ番目のセルの壁面から一番遠い点までの距離ｄ_ｊとする。このときにｄ_ｊ−１≦２Ｒ_ｏ≦ｄ_ｊとなる壁面から１番目からｊ番目のセルの中で発泡するとする。すなわち、式（４０）とし、壁面から１番目からｊ番目の全てのセルの中でｎ_{Ｓｗａｌｌ}の発泡数で発生するとした。

尚、図３に示すような領域で気泡が発生する場合のボイド率の取り方について説明する。ｄ_ｊ≧２Ｒ_０、かつ、ｄ_ｊ−１＜２Ｒ_０であるｊ以下のセルでは、ｊ以下のセルの平均的ボイド率を使用する。すなわち、図４に示すように、ｊ＝２では、実際のボイド率はα_ｊ＝２であるが、本発明の計算では、ボイド率α＝α_{ｊ＝１：２}が用いられる。また、図５に示すように、ｊ＝１では、実際のボイド率はα_ｊ＝１であるが、本発明の計算では、ボイド率α＝α_{ｊ＝１：２}が用いられる。ここで、ｊ＝１：２の平均ボイド率：α_{ｊ＝１：２}は、図６に示すように、ｊ＝１及びｊ＝２の両方の領域を１つの領域とした場合のボイド率を算出することにより得られる。

この取り扱いにより壁面近傍のグリッドを適度に小さくしておけば、壁面近傍の格子点の細かさに対する気泡発泡のグリッド依存性はほとんどなくなる。

気泡発生数は次式（４１）で示されるΠ_Ｇ，気泡発生質量は次式（４２）で示されるΛ_Ｇである。

新規発泡気泡は一番小さな気泡質量階級に属するため既存の気泡の保存量に付加される。発泡前の領域に存在する単位体積あたりの気泡数をｎ_{Ｇ＿ｏｌｄ}，質量をｍ_{Ｇ＿ｏｌｄ}とすると、発泡後の単位体積あたりの気泡数ｎ_{Ｇ＿ｎｅｗ}，質量ｍ_{Ｇ＿ｎｅｗ}は次式となる。

５．気泡質量に基づいた階級の離散化
気泡径分布モデルでは、空間の次元に加えて１個あたりの気泡質量の次元を追加する。以下その次元の向きをｗ方向と呼ぶことにする。ｗ軸の正方向が質量の増加を意味するようにその向きを定義する。ある空間におけるｗ軸方向の気泡数密度の確率密度関数Ｎ_Ｇｗの分布を図７に示す。気泡群が蒸発に伴い成長すると気泡質量が増えるためw軸の正の方向に移動し、気泡が凝縮すると気泡質量が減るためw軸の負の方向に移動する。すなわち、空間的移動とは独立に気泡の蒸発凝縮に伴い気泡質量軸方向に気泡が運動するという新しい概念を導入する。ここで空間的に静止した気泡1個あたりの蒸発量γ_１とするときの気泡数密度関数Ｎ_Ｇｗは、次式となる。

この式（４５）は、気泡の成長がｗ軸方向に移流速度γ_１で運動していることを表している。過熱液中（気泡成長の場合）または過冷液中（気泡凝縮の場合）に存在する気泡周りに生成される温度境界層から供給される熱量によって駆動される蒸発・凝縮量がγ_１である。過熱・過冷液中に投入されてからｔ_{ｅｌａｐｓｅ}後の球形気泡周りに形成される温度境界層の気泡界面上での温度勾配は次式となる。

よって、気泡界面を通過する熱量はフーリエの法則より次式となる。すなわち、熱伝導により伝えられる熱流束は，温度勾配ｄＴ／ｄｒに比例するという法則による。尚、比例定数は、熱伝導率ｋ_Ｌとなる。

この熱量が蒸発・凝縮に用いられるので蒸発・凝縮量は、式（４８）となる。

以上は静止気泡について述べたが空間移動する気泡においても一般性は失われない。すなわち、蒸発を伴いながら空間移動する気泡数密度関数の時間変化は次式となる。

ここで、δ（ｗ−ｗ_０）は、式（５０）で表されるデルタ関数を示す。

デルタ関数とは、カッコ内の数が０のときを含んで積分すると１となり、カッコ内の数が０でないとき０となるパルス関数である。ここで、ｗ_０は、新規発生気泡の質量を表す。

２次元の場合∇，Ｖ_Ｇは、次式となる。

同様に、気泡発生からの経過時間であるｔ_{ｅｌａｐｓｅ}も気泡の移動に伴って移流すると考えられるため、式（５４）のように表せる。

ここで右辺の１は、１タイムステップごとに時間刻み分だけ経過時間が増えることを表している。

実際の計算では、気泡数密度関数Ｎ_Ｇｗをｗで積分して求めた階級内の単位体積あたりの気泡数密度ｎ_Ｇ，ｋと、気泡数密度関数Ｎ_Ｇｗにｗをかけてからｗで積分して求めた階級内の単位体積あたりの気泡質量ｍ_Ｇ，ｋとを用いて行う。ここで、ｗは１個あたりの質量を表すものである。気泡数密度ｎ_Ｇ，ｋ及び気泡質量ｍ_Ｇ，ｋは、式（５５）、式（５６）のように定義される。

この気泡数密度、気泡質量を気泡の質量軸をある質量範囲で定義された階級に分割する。ある階級ｋについて考える。階級の軽い方の境界をｌｉｇｈｔ境界，重い方をｈｅａｖｙ境界と呼び、各々の気泡質量をｗ_{ｌｉｇｈｔ，ｋ}，ｗ_{ｈｅａｖｙ，ｋ}とすると、階級内に存在する単位質量あたりの気泡数密度ｎ_Ｇ，ｋ、気泡質量ｍ_Ｇ，ｋは次式で定義される。

ここで、ｍ_１，ｋは、階級内の気泡１つあたりの平均質量を表し、次式（５９）となる。

この平均質量ｍ_１，ｋは、階級の境界質量ｗ_{ｌｉｇｈｔ，ｋ}，ｗ_{ｈｅａｖｙ，ｋ}と以下の関係を満たさなければならない。

本式（６０）は、実際に計算する際に、気泡の階級分けが適切に行われているのかをチェックするのに有効な条件式となる。

基礎式（４９）を単位体積あたりの気泡数密度ｎ_Ｇ，ｋと単位体積あたりの気泡質量ｍ_Ｇ，ｋの定義式（５７）、（５８）を用いて変形すると次式のように表せる。

さらに、階級ｋ内の気泡数密度関数Ｎ_Ｇｗの分布について考える。Ｎ_Ｇｗの分布で一番単純なものは図８に示すように、階級内で全て同一の値を採る均一な分布であると仮定したものである。このとき、Ｎ_Ｇｗは階級内でどこでも同じ値で、かつ、平均値を取るため、次式（６３）となる。

本式（６３）は、式（５７）を満たす。しかし、式（５８）を満たすのは、次式（６４）のときのみである。

階級ｋ内では常に式（５７）、式（５８）を満たさなければならないという条件を満たさない。これは気泡数密度関数Ｎ_Ｇｗの分布が図８に示すように階級内で全て同一の値を採るという仮定が誤っていることを示す。そこで、図９に示すように重心をｍ_１，ｋとする３つの１次関数分布を持つと仮定し、分類した。最小階級用の分布１型１１、中間階級用の分布２型１２、最大階級用の分布３型１３である。どの条件でどの分布を持つかを式（６５）〜（６７）に示す。

また、それぞれの境界値Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}），Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）及び端点ｗ_ｆ，ｗ_ｂの値は、以下の式（６８）〜（７０）の通りとなる。

この分布は、式（５７）だけでなく式（５８）も満たす。この分布は、さらに数値計算上の大きな効果を持っている。それは、階級内の気泡がある程度成長しないとそれより大きな階級に気泡が全く移動しないという点である。最大階級の分布である分布３型を見ると重い側の境界において気泡数密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）＝０であることから分かる。不要な数値拡散を抑制する効果をもち、精度の向上をもたらす。

６．気泡の並進運動
次式を用いて時間発展的にスリップ速度Ｖ_Ｓ＝Ｖ_Ｇ−Ｖ_Ｌを解く。

流れ場の計算から求められたＶ_Ｌを基準に気相速度Ｖ_Ｇを求める。尚、式（７１）において、右辺第１項は、スリップ速度Ｖ_Ｓの気泡と液体の間に働く力、すなわち、気液間抗力を表すものである。気液間抗力は、抗力係数Ｃ_Ｄ、気泡前面投影面積πＲ^２、及び、スリップ速度の動圧ρ_Ｌ（｜Ｖ_Ｓ｜Ｖ_Ｓ）／２の積である。また、式（７１）の右辺第２項は、圧力勾配中に存在する気泡が受ける力で、いわゆる浮力にあたるものである。

７．液体の圧力計算と気泡の半径計算のカップリング
その計算点においてあらゆる階級に気泡の存在しない純液体領域（液体体積率α_Ｌ＝１）の温度Ｔ_Ｌ、圧力Ｐ_Ｌは流れ場の計算から求められた密度ρ_Ｌ、内部エネルギーｅ_ＳＬを用いて、式（７２）、式（７３）から一義に定まる。

一方、気相の存在する２相領域の温度Ｔ_Ｌは液相部分の各保存量α_Ｌρ_Ｌ，α_Ｌρ_Ｌｕ_Ｌ，α_Ｌρ_Ｌν_Ｌ，α_Ｌρ_Ｌｅ_Ｌより内部エネルギーｅ_ＳＬ＝ｅ_Ｌ−０．５（ｕ^２ _Ｌ＋ν^２ _Ｌ）が求まるので式（７２）を用いて定まる。しかし、圧力Ｐ_Ｌは液体体積率α_Ｌが定まらないと密度ρ_Ｌが定まらず、加えて液体体積率α_Ｌ自身が圧力Ｐ_Ｌの関数となっており収束計算を必要とする。以下に圧力Ｐ_Ｌと液体体積率α_Ｌの収束計算ループの過程について記す。この液体圧力計算のフローチャートを図１１に示す。はじめに仮の液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を設定する。本計算ではＰ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}＝０を初期値とし、徐々にＰ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を大きくして行き、収束したＰ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を圧力Ｐ_Ｌとする。

ここで、仮の液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}に対する液体体積率α_Ｌを求める手順を示す。液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}に平衡する気泡半径Ｒ_ｋを求める。この気泡半径Ｒ_ｋも新規気泡発生、蒸発・凝縮、消滅などの現象が複雑に絡み合っているため収束計算を必要とする。この気泡半径計算のフローチャートを図１２に示す。気泡半径Ｒ_ｋを求める小ループでは着目しているコントロールボリューム全てが階級ｋの気泡で占められた場合の気泡半径Ｒ_ｋの最大値Ｒ_{ｍａｘ，ｋ}を初期値として徐々にＲ_ｋを減少させて収束させる。各階級ごとの気泡数密度ｎ_Ｇ，ｋ，気泡質量ｍ_Ｇ，ｋと求める。この計算には液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と液体飽和圧力Ｐ_ｓａｔ（Ｔ_Ｌ）とから求められる新規気泡発生の効果、液体温度Ｔ_Ｌと気泡温度Ｔ_Ｇ，ｋから求められる蒸発・凝縮量の効果、気泡温度Ｔ_Ｇ，ｋが臨界温度Ｔ_Ｃ以上、つまり、超臨界状態になったとき気泡が液体に溶解する効果を陰的に含んでいる。次に既存の気泡では、一般に次式（７４）で表されるレイリープレセット式（Ｒａｙｌｅｉｇｈ Ｐｌｅｓｓｅｔ Ｅｑｕａｔｉｏｎ）が適用されることが知られている。

発生直後や崩壊直前の小さな気泡には、界面半径方向の慣性力の効果が大きいため右辺第４〜６項の時間微分項の効果が大きいが、大きな気泡では界面の半径方向の慣性力の効果が相対的に小さくなるため右辺第４〜６項の時間微分項の効果が小さくなる。そこで、１タイムステップ前の気泡半径Ｒ_{ｏｌｄ，ｋ}、２タイムステップ前の気泡半径Ｒ_{ｏｌｄｏｌｄ，ｋ}とするとき、Ｒ_{ｏｌｄ，ｋ}が一定の気泡径Ｒ_１以下（実際の計算ではＲ_１=１μｍ）では、時間について離散化した次式（７５）に示すように時間微分項を１−Ｒ_{ｏｌｄ，ｋ}／Ｒ_１倍して気泡圧力Ｐ_Ｇ，ｋを計算し、一定の気泡径Ｒ_１以上では次式（７５）に示すように時間微分項を無視して気泡圧力Ｐ_Ｇ，ｋを計算する。

この取り扱いにより計算負荷を低減しつつ、気泡径Ｒ_１前後で連続的に時間微分項の効果を導入できるようになった。

さらに既存気泡のない計算セルに新規に気泡が発生した場合にはスリップが無視できるので、次式（７６）を用いて気泡圧力Ｐ_Ｇ，ｋを計算する。

これらの気泡圧力Ｐ_Ｇ，ｋに対応する気泡温度Ｔ_Ｇ，ｋ、気泡密度ρ_Ｇ，ｋが飽和条件より求められる。

ここで求められたρ_Ｇ，ｋを用いて、先ほど求めた気泡数密度ｎ_Ｇ，ｋ、気泡質量ｍ_Ｇ，ｋが次式（７９）を満たすまで、気泡半径Ｒ_ｋを求める小ループを繰り返す。

よって、階級ｋのボイド率α_Ｇ，ｋは次式（８０）と表せる。

次に液体温度Ｔ_Ｌ，仮の液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}に対して求めた階級ｋのボイド率α_Ｇ，ｋを用いて真の液体圧力Ｐ_Ｌを求める過程を示す。仮の液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}に対する液体体積率α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と液体密度ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}とが、既知の保存量α_Ｌρ_Ｌとの関係が式（８３）となるまで、液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を求めるループを繰り返す。

以上のような繰り返し計算をしているとある階級ｋの気泡が消滅することがある。そのときの体積あたりの気泡質量はｍ_Ｇ，ｋなので各保存式より、次式（８４）〜（８７）となる。

この過程を液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を求めるループ中の最後に行う。

８．時間刻み
時間刻みは移流項に対してＣＦＬ条件を満たす必要がある。また、圧力勾配による加速項、ｓｔｉｆｆな生成項と粘性項に対する時間制限も満足しなくてはならない。ここで、ｓｔｉｆｆな生成項とは、計算を硬直化させる大きな生成項である。すなわち、生成項、即ち式（１）のＳが大きく、かつ、時間刻みが大きい場合、計算が不安定となる。このことを防ぐために、時間刻みを小さくする必要があり、計算を硬直化（ｓｔｉｆｆｎｅｓｓ）させる。このような生成項をｓｔｉｆｆな生成項である。

ここで，Ａ_Ｌｘ，Ａ_Ｇｘ，Ａ_Ｌｙ，Ａ_Ｇｙは、各相各方向の圧力勾配項、粘性項および生成項などの時間微分値である。

ここでＤｉｍは次元数。１次元ならＤｉｍ＝１，２次元ならＤｉｍ＝２、３次元ならＤｉｍ＝３となる。

次に、本発明を適用した解析装置１を用いたキャビテーション解析の全体の計算手順について説明する。表７に、以下で用いられる計算条件を示す。尚、本条件は、Ｓｉｍｏｎｅａｕ ａｎｄ Ｈｅｎｄｒｉｃｋｓの実験条件と同一とした。

まず、図１０に示すように、ステップＳ１〜ステップＳ２の流れ場の計算について説明する。ステップＳ１−１では、計算対象の格子を図１３に示すように読み込む。ここでは、図１３に示すよう２次元で、流れ方向１４３点、高さ方向１２点とした。

ステップＳ１−２では、初期条件を表７のように、また、境界条件を表８のように設定する。本解析においては、表７に示すように、壁面を断熱壁として設定される。尚、本条件は、Ｓｉｍｏｎｅａｕ ａｎｄ Ｈｅｎｄｒｉｃｋｓの実験条件と同一とした。

ステップＳ２−１では、液体温度、気泡平均温度、液体体積率、ボイド率から上述の式（３５）を用いて、２相平均音速ｃ_Ｖを計算した。

ステップＳ２−２では、式（１）を用いて、液体の各方向速度より粘性応力項Ｆ_Ｖを計算する。

ステップＳ２−３では、式（１）の各基礎式の空間移流成分の各計算点での時間微分値を計算し、空間移流項Ｅ，Ｆを算出する。尚、空間移流成分はＭＵＳＣＬ法（Ｍｏｎｏｔｏｎｉｃ Ｕｐｓｔｒｅａｍ−ｃｅｎｔｅｒｅｄ Ｓｃｈｅｍｅ ｆｏｒ Ｃｏｎｓｅｒｖａｔｉｏｎ Ｌａｗｓ）により精度を高めたＳＨＵＳ（Ｓｉｍｐｌｅ Ｈｉｇｈ−ｒｅｓｏｌｕｔｉｏｎ Ｕｐｗｉｎｄ Ｓｃｈｅｍｅ）を利用し最大で空間３次精度を達成しつつＴＶＤ（Ｔｏｔａｌ Ｖａｒｉａｔｉｏｎ Ｄｉｍｉｎｉｓｈｉｎｇ）性を保った風上差分を採用することで、衝撃波等の不連続面も安定的に扱うことが可能となった。

ステップＳ２−４では、式（７１）に示すように、液体の各方向成分、気泡半径Ｒから、気液間抗力及び各階級の気泡と液相のスリップ速度Ｖ_Ｓを計算する。ここで、液体の各方向成分とは、スリップ速度Ｖ_Ｓ＝（ｕ_Ｓ，ｖ_Ｓ）、気相速度Ｖ_Ｇ＝（ｕ_Ｇ，ｖ_Ｇ）、液層速度Ｖ_Ｌ＝（ｕ_Ｌ，ｖ_Ｌ）である。また、ｘ方向成分は、ｕ_Ｓ，ｕ_Ｇ，ｕ_Ｌであり、ｙ方向成分は、ｖ_Ｓ，ｖ_Ｇ，ｖ_Ｌである。尚、スリップ速度、すなわち、空間微分∇Ｖ_Ｇ，∇Ｖ_Ｌの計算には、風上１次精度の差分法を用いた。

ステップＳ２−５では、式（１）に示すように、上述のステップＳ２−２〜Ｓ２−４で算出した、粘性応力項Ｆ_Ｖ、空間移流項Ｅ，Ｆ、気液間抗力を用いて液体保存量（液相部質量、Ｘ方向成分運動量、Ｙ方向成分運動量及び全エネルギー量）及び各階級の気泡の時間微分値を計算する。

ステップＳ２−６では、式（８８）〜（９１）に示すように、上述のステップ７の結果からＣＦＬ（Ｃｏｕｒａｎｔ Ｆｒｉｅｄｒｉｃｈ Ｌｅｗｙ）条件を満たす時間刻みを決定する。

ステップＳ２−７では、ステップＳ２−５で算出した液体保存量及び各階級の気泡の時間微分値とステップＳ２−６で算出した時間刻みを用いて、時間進行させ、７つの保存量を計算する。すなわち、液体の保存量である単位体積当たりの液相部質量α_Ｌρ_Ｌ、Ｘ方向成分運動量α_Ｌρ_Ｌｕ_Ｌ、Ｙ方向成分運動量α_Ｌρ_Ｌν_Ｌ及び全エネルギーα_Ｌρ_Ｌｅ_Ｌ、並びに、気泡数密度関数Ｎ_Ｇｗから式（５７）、式（５８）で求められる各階級の仮の気泡保存量である各階級の仮の気泡数密度ｎ_Ｇ，ｋ及び気泡質量ｍ_Ｇ，ｋが算出されると共に気泡の発生からの経過時間ｔ_{ｅｌａｐｓｅ，ｋ}が決定する。

次に、図１１に示すように、ステップＳ３〜ステップＳ７の上述のステップＳ２で算出した７つの保存量を用いて液体圧力および体積率を繰り返し計算する方法について説明する。ここでは、気泡の成長・収縮計算と液体の圧力・温度計算の計算安定性を高めるため両者をカップリングして繰り返し計算、いわゆる陰的計算を行う。

ステップＳ３では、繰り返し計算の初期値として、液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を適当に仮定する。

ステップＳ４では、式（７２）を用いて、単位体積当たりの液相部質量α_Ｌρ_Ｌ、Ｘ方向成分運動量α_Ｌρ_Ｌｕ_Ｌ、Ｙ方向成分運動量α_Ｌρ_Ｌν_Ｌ及び全エネルギーα_Ｌρ_Ｌｅ_Ｌから液体の内部エネルギーｅ_ＳＬを求め、この値と物性値から液体温度Ｔ_Ｌを算出する。

ステップＳ５では、ステップＳ２−７で算出した７つの保存量と、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と、ステップＳ４で算出した液体温度Ｔ_Ｌとを用いて気泡半径の繰り返し計算を行い、各階級のボイド率α_Ｇ，ｋを算出すると共に、各階級の気泡数密度ｎ_Ｇ，ｋ及び気泡質量ｍ_Ｇ，ｋを補正する。

ステップＳ６では、式（８４）〜（８７）に示すように、後述するステップ５−１で算出した新規気泡発生量並びにステップＳ５−６で算出された蒸発量・凝縮量Γ_Ｉ及び気泡崩壊量により液体保存量（単位体積当たりの液相部質量α_Ｌρ_Ｌ、Ｘ方向成分運動量α_Ｌρ_Ｌｕ_Ｌ、Ｙ方向成分運動量α_Ｌρ_Ｌν_Ｌ及び全エネルギーα_Ｌρ_Ｌｅ_Ｌ）を補正する。

ステップＳ７では、まず、式（８２）に示すように、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と、ステップＳ４で算出された液体温度Ｔ_Ｌから仮の液体密度ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を算出する。次に、ステップＳ５で算出した階級ｋのボイド率α_Ｇ，ｋから式（８１）を用いて液体体積率α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を算出する。そして、この仮の液体密度ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と液体体積率α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}とから計算される仮の液体保存量α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と液体保存量α_Ｌρ_Ｌとが式（８３）で示される関係を満たすかどうかを比較して確かめる。

尚、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と、ステップＳ４で算出された液体温度Ｔ_Ｌから仮の液体密度ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を算出すると共に、ステップＳ５で算出した階級ｋのボイド率α_Ｇ，ｋから液体体積率α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を算出し、この液体体積率α_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と仮の液体保存量α_Ｌρ_Ｌとから算出される液体密度ρ_Ｌとを比較してもよい。

ステップＳ５により算出されるα_{Ｌ＿ｔｅｍｐｏｒａｒｙ}ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}とステップＳ２−７で算出した仮の液体保存量αρが異なる場合は、ステップＳ３に戻り、液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}を再度仮定して、繰り返し計算を行う。１又は複数回ステップＳ３〜Ｓ７を繰り返した後に、ステップＳ７で算出されるα_{Ｌ＿ｔｅｍｐｏｒａｒｙ}ρ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}とステップＳ２−７で算出した液体保存量αρが等しくなったときに、この繰り返し計算は終了し、このときの液体圧力Ｐ_Ｌが収束値となる。

ここで、図１２を用いて、ステップＳ５で、上述のステップＳ２で算出した７つの保存量と、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と、ステップＳ４で算出した液体温度Ｔ_Ｌとを用いて気泡半径の繰り返し計算を行い、液体の体積率を計算する方法について更に詳しく説明する。ここでは、気泡の成長・収縮計算から液体ボイド率を求める方法を説明する。

すなわち、ステップＳ５−１では、式（４１）〜（４２）に示すように、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と、ステップＳ４で算出した液体温度Ｔ_Ｌとを用いて過飽和度、新規発生する新規気泡半径、新規気泡数密度Π_Ｇ、新規気泡質量Λ_Ｇを計算する。ここで、算出された新規発生した気泡は一番小さな気泡質量階級に属するので、式（４３）〜（４４）に示すように最小階級の気泡数密度ｎ_Ｇ，ｋ、気泡質量ｍ_Ｇ，ｋの値に加算される。尚、新規気泡数密度は、式（３７）に示すように、過飽和度である（Ｐ_Ｇ−Ｐ_Ｌ）を用いて算出する。ここで、式（７７）の関係、及び、新規発生気泡においてはＴ_Ｇ＝Ｔ_Ｌであることを仮定していることから、Ｐ_Ｇ−Ｐ_Ｌ＝Ｐ_Ｇ（Ｔ_Ｌ）−Ｐ_Ｌとすることができ、過飽和度（液体過熱度）を圧力表記により間接的に表記したものである。また、新規気泡質量は、新規気泡数密度及び新規気泡半径を用いて算出する。

ステップＳ５−２では、各階級の既存の気泡半径Ｒ_Ｋを繰り返し計算の初期値として仮定する。

ステップＳ５−３では、式（７５）〜（７６）に示すように、液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}と液体保存量を用いて、気泡界面の運動方程式から各階級の気泡圧力Ｐ_Ｇ，ｋを計算する。気泡圧力は、式（７４）で表せるレイリープレセット式（Ｒａｙｌｅｉｇｈ Ｐｌｅｓｓｅｔ Ｅｑｕａｔｉｏｎ）が適用されることが知られている。一般に、発生直後及び崩壊直前の小さな気泡（１μｍ以下）には、界面半径方向の慣性力の効果が大きいため右辺第４〜６項の時間微分項の効果が大きいが、比較的大きな気泡（１μｍ以上）では界面の半径方向の慣性力の効果が相対的に小さくなるため右辺第４〜６項の時間微分項の効果が小さくなる。そこで気泡径に応じて右辺第４〜６項の時間微分項に重み付けをした式（７５）を用いる。新規発生時には気泡と液体のスリップが無視できるのでスリップ項を無視した式（７６）を用いる。

ステップＳ５−４では、式（７７）、（７８）に示すように、ステップＳ５−３で算出した気泡圧力Ｐ_Ｇ，ｋを用いて、飽和条件から各階級の気泡温度Ｔ_Ｇ，ｋ、気泡密度ρ_Ｇ，ｋを計算する。

ステップＳ５−５では、ステップＳ３で仮定した液体圧力Ｐ_{Ｌ＿ｔｅｍｐｏｒａｒｙ}、ステップＳ４で算出した液体温度Ｔ_Ｌ、ステップＳ５−１で算出した気泡数密度ｎ_Ｇ，ｋ、気泡質量ｍ_Ｇ，ｋ、ステップＳ５−３で算出した気泡圧力Ｐ_Ｇ，ｋ、ステップＳ５−４で算出した気泡温度Ｔ_Ｇ，ｋ、気泡密度ρ_Ｇ，ｋを用いて、式（５７）〜式（７０）に示すように、各階級内の気泡数密度分布を仮定する。すなわち、ここでは、階級の境界値と、流れ場および階級相互の計算から求められる気泡数密度ｎ_Ｇ，ｋと体積当たりの気泡質量ｍ_Ｇ，ｋを満たすように、階級内で１次関数を用いた気泡数密度の分布を算出する。

気泡数密度の分布を算出するには、まず、式（５９）に示すように、階級ｋ内の気泡１つあたりの平均気泡質量ｍ_１，ｋを算出する。次に、この平均気泡質量ｍ_１，ｋと、階級ｋに隣接する質量の軽い側の階級との境界での質量値ｗ_{ｌｉｇｈｔ，ｋ}と、階級ｋに隣接する質量の重い側の階級との境界での質量値ｗ_{ｈｅａｖｙ，ｋ}とに基づいて、式（６５）〜（６７）に示すように、最小階級用の分布１型１１、中間階級用の分布２型１２、最大階級用の分布３型１３に分類する。分類された各階級の境界値Ｎ_ＧＷ及び端点Ｗ_ｆ，Ｗ_ｂは、式（６８）〜（７０）により算出される。各階級内の気泡数密度分布は、例えば図９に示すように、この境界値及び端点から１次関数分布が得られる。各々の分布の１次式より階級の境界での気泡数密度がより厳密に計算できる。

階級内で１次関数を用いた気泡数密度の分布を仮定できることは、数値計算上大きな効果を有する。例えば、階級内の気泡が有る程度成長しないとそれより大きな階級に気泡が全く移動しないので、最大階級の分布３型においては、重い側の境界において数密度の質量微分Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）＝０である。このように不要な数値拡散を抑制することができ、計算精度を向上させることができる。

ステップＳ５−６では、階級境界での蒸発量・凝縮量を計算する。階級内での蒸発・凝縮量Γ_Ｉは、式（６２）の中辺２第２項と右辺第２項の関係で示されるように気泡１つあたりの蒸発・凝縮量γ_Ｉと気泡数密度関数Ｎ_Ｇｗの積のｗに関する積分により算出される。

ステップＳ５−７では、この各階級の気泡１個あたりの蒸発量γ_Ｉと階級内の蒸発・凝縮量Γ_Ｉから、式（６１）〜（６２）により、階級ごとの気泡数密度ｎ_Ｇ，ｋと気泡質量ｍ_Ｇ，ｋが変更（修正）される。気泡数密度ｎ_Ｇ，ｋと気泡質量ｍ_Ｇ，ｋから、式（７９）に示すように、気泡１個あたりの体積を算出し、気泡半径Ｒ’_Ｋを求める。

これらの計算により気泡質量、および、液体密度、各方向液体運動量、全液体エネルギーに変化が生じる。さらに気泡質量が０となった場合には気泡が消滅する。これにより、再度、気泡数密度、気泡質量に変化が生じる。

ステップＳ５−８では、ステップＳ５−７で算出された気泡半径Ｒ’_ＫとステップＳ５−２で仮定した気泡半径Ｒ_ｋを比較する。図１２に示すように、ステップＳ５−７で算出された気泡半径Ｒ’_ＫとステップＳ５−２で仮定した気泡半径Ｒ_ｋとが等しくならない場合は、ステップＳ５−２に戻り、気泡半径を再度仮定して、繰り返し計算を行う。

１又は複数回ステップＳ５−２〜ステップＳ５−７を繰り返した後に、ステップＳ５−７で算出された気泡半径Ｒ’_ＫとステップＳ５−２で仮定した気泡半径Ｒ_ｋが等しくなったときに、この繰り返し計算は終了し、この気泡半径Ｒ_ｋが収束値となる。

ステップＳ５−９では、式（８０）に示すように、この気泡半径Ｒ_ｋから各階級のボイド率α_Ｇ，ｋを算出する。

以上のように、本発明を適用した解析装置１は、仮に決定された液相部分の保存量を満たすように、新規発生・消滅気泡数、気泡の蒸発・凝縮量等を考慮した繰り返し計算により、各計算点での圧力、温度、液体各方向運動量、各階級ごとのボイド率、気泡数密度、気泡質量、及び、発生からの経過時間を決定することができる。

ステップＳ８では、図１０に示すように、Ｓ２〜Ｓ７において算出された各計算値の境界値をＳ１−２で設定した表８の境界条件により補正、すなわち、再設定する。設定される境界条件は、流入境界、流出境界及び壁面境界の条件である。

流入境界では、温度は一定とし、液単相であるサブクール流入とした。サブクール流入とは、液体温度が飽和温度（沸点）より低い状態で流入することを表すことを意味する。キャビテーションは圧力が下がることにより飽和温度が下がり、飽和温度が液体温度（一定）より低くなることで沸騰し気泡が発生する現象である。そのため、サブクールの状態では、一般的に気泡が発生せず、気泡の存在しない液単相の状態になる。流入境界の速度は、勾配０で外挿し、１点内側の計算格子点の速度と等しくする。ここで外挿とは、内挿の反対のことを意味する。内挿とはある点を挟む２点の値から、その内側にある任意の点を類推することであり、すなわち、外挿とはある点に隣接する２点から、その外側にある任意の点を類推することである。ここで算出する境界点は、計算点の外側に存在するために外挿によりこの値を算出する。全圧Ｐ_ｔが一定となるように、流入境界の圧力Ｐは、Ｐ＝Ｐ_ｔ−０．５Ｖ^２と設定される。

流出境界では、温度、速度、気泡諸量ともに勾配０で外挿し、１点内側の計算格子点の温度、速度、気泡諸量と等しく設定される。流出境界上の速度がその点での音速よりも小さい亜音速流出の場合は、流出境界上の出口圧力を一定にした。一方、流出境界上の速度がその点での音速よりも大きい超音速流出の場合は、流出境界上の出口圧力は、勾配０で外挿し１点内側の計算格子点の圧力と等しく設定される。

壁面境界では、壁面に流体及び気泡が流入しないことから、圧力、気泡諸量は勾配０で外挿し１点内側の計算格子点の圧力、気泡諸量と等しくし、速度は壁面上で静止されるようにした。壁面から熱が流入しない断熱壁の場合には、壁面上の温度勾配が０とされ外挿し、１点内側の計算格子点の温度と等しくされ、壁面温度が一定の等温壁の場合は、壁面上の温度を壁温度で一定とされる。

このステップＳ８では、壁面近傍における数値等が現実的な値からかけ離れた値となることを防止し、この境界条件により調整することで、よりキャビテーション現象を正確に再現することができる。ここで一定値を取る境界条件は「ディリクレ条件」を意味するものである。また、勾配一定値を取る条件は「ノイマン条件」を意味するものである。

ステップＳ９では、図１０に示すように、あらかじめ設定された出力条件を満たさなかった場合、若しくは、ユーザ等が手動により出力しないことを決定した場合に、ステップ３に戻り再計算を行うことができる。ここで、あらかじめ設定される出力条件には、例えば、一定時間おきに出力をすることがある。この場合、指定された一定タイムステップに達していないときは、ステップ３に戻り再計算を行う。このとき、計算結果を旧時間のデータとしてメモリに保存し、次の計算を行うことができる。また、あらかじめ設定された出力条件を満たした場合、若しくは、ユーザー等が手動により出力を決定した場合に、計算結果を出力することができる。このとき、指定された割合で画像や保存データを出力することができる。計算結果は、例えば、図１４に示すように、圧力、マッハ数、温度、ボイド率、気泡半径（階級１）及び気泡半径（階級２）を、数値の大きさを色の分布で示すことができるとともに、図１５に示すように、壁面上圧力分布を線図上に得ることができる。図１４において、面塗りは液体圧力，マッハ数，液体温度，ボイド率，階級１と階級２の気泡半径を示す．図１５において、線図は壁面上圧力分布の計算結果を示し、白丸はＳｉｍｏｎｅａｕ ａｎｄ Ｈｅｎｄｒｉｃｋｓが行った壁面上圧力分布の実験結果である。図１４、図１５の結果により、本解析装置１において、キャビテーション現象が再現されていることが分かる。

ステップＳ１０では、あらかじめ設定された終了条件を満たさなかった場合、若しくは、ユーザー等が手動により終了しないことを決定した場合に、ステップ３に戻り再計算を行うことができる。また、あらかじめ設定された終了条件を満たした場合、若しくは、ユーザー等が手動により終了する決定した場合に、解析装置１による解析を終了する。尚、ここで終了条件としては、非定常計算における、指定した時間または、指定したタイムステップまで計算することを設定することができ、さらに、定常計算における、定常状態になるまで計算することを設定することができる。この終了条件を設定自在としたことにより、様々な用途の解析に用いることができる。

本発明を適用した解析装置１は、気泡密度分布を考慮することにより、数値拡散を防ぐことができたので、従来の解析装置の計算プログラムにおいて気泡の隣接階級の間隔比を１のみでしか行われていなかった計算を、気泡の隣接階級の間隔比を２．１５（１０^１／３）程度まで高めても計算を行えるようになった。このことにより、本発明は、例えば、１ｎｇ〜１００ｎｇの気泡を同時に扱う場合、従来では、１００階級以上の階級を必要であったのに対し、本解析装置では、６階級で計算することができるようになった。

また、本発明を適用した解析装置１は、気泡の成長・収縮計算と液体の圧力・温度計算をカップリングして繰り返し計算を行うことにより、時間刻みが極端に小さくなる現象を回避できたので、従来の解析装置における計算プログラムにおいてキャビテーション気泡無し流れで１０^−８（ｓ）程度の時間刻みで計算可能であった流れ場において、気泡が成長してくると１０^−１２（ｓ）以下の時間刻みが必要であるのに対し、本解析装置では、キャビテーション気泡の有無に拘わらず１０^−８（ｓ）程度の時間刻みで安定的に計算できるようになった。

上述のように、本発明は、各階級内の気泡数密度分布を考慮することにより、数値拡散を防ぐことができ、気泡の蒸発・凝縮計算の精度を高めることができるとともに、階級間隔を不等間隔にした場合の計算安定性を高めることができる。よって、本発明は、精度及び安定性の高いキャビテーション解析を可能とするものであり、汎用性を高めることを可能とする。

また、本発明によれば、気泡の成長・収縮計算と流体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、時間刻みを大きく取ることができ、キャビテーション解析の計算安定性を増すことができ、よって、計算負荷を低減することが可能となる。

尚、上述した解析装置は、２次元計算を行うものであるが、本発明は、３次元計算を行うこともできる。

本発明を適用したキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置は、ロケットエンジン等に用いられる液体酸素、液体水素用のポンプ及びＬＰＧ等の液体メタンに用いられるポンプのキャビテーション解析に用いることができるのみならず、液体窒素、液体酸素、液体水素及び液体メタン等の各種低温流体中に発生するキャビテーション解析に用いることができる。

本発明を適用した解析装置の構成を示す図である。 本発明を適用した解析装置により得られるボイド率と音速の関係を示す図である。 壁面近傍に発生する新規気泡の壁面の面積と計算空間との関係を示す図である。 図３に示すｊ＝２のセルにおけるボイド率を説明する図である。 図３に示すｊ＝１のセルにおけるボイド率を説明する図である。 壁面近傍に発生する新規気泡のボイド率の与え方を説明するための図である。 気泡数密度の確率密度関数の気泡１つあたりの質量に対する分布を示す図である。 従来技術で用いられる各階級の気泡の気泡数密度分布を示す図である。 本発明を適用した解析装置により算出される各階級の気泡の数密度分布を示す図である。 本発明を適用した解析装置のフローチャートである。 本発明を適用した解析装置の液体圧力計算のフローチャートである。 本発明を適用した解析装置の気泡半径計算のフローチャートである。 本発明を適用した解析装置により解析されるノズル形状を示す図である。 本発明を適用した解析装置により得られる計算結果の一例を示す図である。 本発明を適用した解析装置により得られる壁面上圧力分布の計算結果の一例を示す図である。

符号の説明

１ 解析装置、２ 入力部、３ 演算部、４ メモリ、５ 表示部、６ 出力部、７ 記憶部

Claims (4)

1. 少なくとも、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析方法において、
   ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するキャビテーション解析方法。
   

   
2. 少なくとも、所定のステップをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析プログラムにおいて、
   上記所定のステップは、
   ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップであるキャビテーション解析プログラム。
   

   
3. 少なくとも、所定のステップを有し、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するコンピュータによって読み取り可能にプログラムを記憶したプログラム記録媒体において、
   上記所定のステップは、
   ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップであるプログラム記録媒体。
   

   
4. 少なくとも、所定のステップをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を解析するキャビテーション解析プログラムを記憶した記録部と、
   上記記録部から上記解析プログラムを読み出して解析処理を実行する演算部とを備える解析装置において、
   上記所定のステップは、
   ある階級ｋにおける気泡数密度ｎ_Ｇ，ｋ、単位体積当たりの気泡質量ｍ_Ｇ，ｋに基づいて式（１）により本階級の平均気泡質量ｍ_１，ｋを算出し、隣接する階級との軽い側の境界における気泡質量の値ｗ_{ｌｉｇｈｔ，ｋ}、隣接する階級との重い側の境界における気泡質量の値ｗ_{ｈｅａｖｙ，ｋ}を用いて、式（２）〜式（４）により本階級を分布１型、分布２型又は分布３型に分類し、上記分布１型、分布２型、及び分布３型に分類された本階級の軽い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｌｉｇｈｔ，ｋ}）、重い側の境界の単位体積当たりの気泡数密度の確率密度関数Ｎ_Ｇｗ（ｗ_{ｈｅａｖｙ，ｋ}）、軽い側の端点ｗ_ｂ及び重い側の端点ｗ_ｆの値を式（５）〜（７）により算出するステップである解析装置。
   

   

Images