(57)【要約】【課題】  流体中に発生するキャビテーションを、気泡径分布モデルを用いて気泡の挙動を取り扱うキャビテーション解析において、その精度、安定性及び汎用性を高める。【解決手段】    流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析する。液体圧力、液体温度、各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を考慮し、この気泡数密度分布から各階級の境界での気泡数密度を算出し、気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出すると共に、気泡の成長・収縮運動の計算を繰り返し計算により算出する。【選択図】  図１１
2005222339.tif 000002
【特許請求の範囲】
【請求項１】
初期条件及び境界条件を入力するステップ１と、
上記入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、
繰り返し計算の初期値として液体圧力を与えるステップ３と、
上記液体保存量並びに上記各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、
上記液体圧力及び上記液体温度から各階級のボイド率を算出すると共に、上記各階級の気泡数密度及び気泡質量を補正するステップ５と、
上記補正された各階級の気泡数密度及び気泡質量から上記液体保存量を補正するステップ６と、
上記液体圧力及び上記液体温度から算出される仮の液体密度と、上記各階級のボイド率から算出される液体の体積率及び上記液体保存量から算出される液体密度とを比較することにより、上記液体圧力を補正してステップ４〜６の処理を行うステップ７と、
上記ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、
上記仮の液体密度と、上記液体の体積率及び液体保存量から算出される液体密度とが等しくなったとき上記ステップ２〜８において算出された各計算値を出力するステップ９とを備え、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析方法において、
ステップ５は、上記液体圧力と上記液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、上記新規気泡数密度及び上記新規気泡質量から最小階級の上記気泡数密度及び上記気泡質量を補正するステップ５−１と、
繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、
上記各階級の気泡半径、上記液体圧力及び上記液体保存量から各階級の気泡圧力を算出するステップ５−３と、
上記各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、
上記液体圧力、上記液体温度、上記補正された各階級の気泡数密度及び気泡質量、上記各階級の気泡圧力、上記各階級の気泡温度、並びに、上記各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、
上記各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、上記各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、
上記蒸発量及び上記凝縮量から上記各階級の気泡数密度及び気泡質量を補正するステップ５−７と、
上記補正された各階級の気泡数密度及び気泡質量並びに上記各階級の気泡密度から算出される各階級の気泡半径を上記仮定した各階級の気泡半径とを比較し、判定結果に基づいて上記仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、
上記仮定した各階級の気泡半径と上記算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とするキャビテーション解析方法。
【請求項２】
上記ステップ５−５では、上記各階級の気泡数密度及び気泡質量に基づいて各階級の平均気泡質量を算出し、上記各階級の平均気泡質量の値と、各階級の境界における質量の値とに基づいて上記各階級内の気泡数密度分布を算出することを特徴とする請求項１記載のキャビテーション解析方法。
【請求項３】
上記ステップ５−５では、上記各階級の平均気泡質量の値ｍ１，ｋと、この階級の隣接する質量の軽い側の境界における気泡質量の値ｗｌｉｇｈｔ，ｋと、この階級の隣接する質量の重い側の境界における気泡質量の値ｗｈｅａｖｙ，ｋとから、下記の式（１）〜式（３）により各階級を分布１型、分布２型又は分布３型に分類し、
上記分布１型、分布２型、及び分布３型に分類された各階級の境界値ＮＧＷ及び端点ｗｆ，ｗｂの値を下記の式（４）〜（６）により算出することで各階級内の気泡数密度分布を算出することを特徴とする請求項１記載のキャビテーション解析方法。
【数１】
2005222339.tif 000060
【請求項４】
初期条件及び境界条件を入力するステップ１と、
上記入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、
繰り返し計算の初期値として液体圧力を与えるステップ３と、
上記液体保存量並びに上記各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、
上記液体圧力及び上記液体温度から各階級のボイド率を算出すると共に、上記各階級の気泡数密度及び気泡質量を補正するステップ５と、
上記補正された各階級の気泡数密度及び気泡質量から上記液体保存量を補正するステップ６と、
上記液体圧力及び上記液体温度から算出される仮の液体密度と、上記各階級のボイド率から算出される液体の体積率及び上記液体保存量から算出される液体密度とを比較することにより、上記液体圧力を補正してステップ４〜６の処理を行うステップ７と、
上記ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、
上記仮の液体密度と、上記液体の体積率及び液体保存量から算出される液体密度とが等しくなったとき上記ステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムにおいて、
ステップ５は、上記液体圧力と上記液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、上記新規気泡数密度及び上記新規気泡質量から最小階級の上記気泡数密度及び上記気泡質量を補正するステップ５−１と、
繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、
上記各階級の気泡半径、上記液体圧力及び上記液体保存量から各階級の気泡圧力を算出するステップ５−３と、
上記各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、
上記液体圧力、上記液体温度、上記補正された各階級の気泡数密度及び気泡質量、上記各階級の気泡圧力、上記各階級の気泡温度、並びに、上記各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、
上記各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、上記各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、
上記蒸発量及び上記凝縮量から上記各階級の気泡数密度及び気泡質量を補正するステップ５−７と、
上記補正された各階級の気泡数密度及び気泡質量並びに上記各階級の気泡密度から算出される各階級の気泡半径を上記仮定した各階級の気泡半径とを比較し、判定結果に基づいて上記仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、
上記仮定した各階級の気泡半径と上記算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とするキャビテーション解析プログラム。
【請求項５】
上記ステップ５−５では、上記各階級の気泡数密度及び気泡質量に基づいて各階級の平均気泡質量を算出し、上記各階級の平均気泡質量の値と、各階級の境界における質量の値とに基づいて上記各階級内の気泡数密度分布を算出することを特徴とする請求項４記載のキャビテーション解析プログラム。
【請求項６】
上記ステップ５−５では、上記各階級の平均気泡質量の値ｍ１，ｋと、この階級の隣接する質量の軽い側の境界における気泡質量の値ｗｌｉｇｈｔ，ｋと、この階級の隣接する質量の重い側の境界における気泡質量の値ｗｈｅａｖｙ，ｋとから、下記の式（７）〜式（９）により各階級を分布１型、分布２型又は分布３型に分類し、
上記分布１型、分布２型、及び分布３型に分類された各階級の境界値ＮＧＷ及び端点ｗｆ，ｗｂの値を下記の式（１０）〜（１２）により算出することで各階級内の気泡数密度分布を算出することを特徴とする請求項４記載のキャビテーション解析プログラム。
【数２】
2005222339.tif 000061
【請求項７】
初期条件及び境界条件を入力するステップ１と、
上記入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、
繰り返し計算の初期値として液体圧力を与えるステップ３と、
上記液体保存量並びに上記各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、
上記液体圧力及び上記液体温度から各階級のボイド率を算出すると共に、上記各階級の気泡数密度及び気泡質量を補正するステップ５と、
上記補正された各階級の気泡数密度及び気泡質量から上記液体保存量を補正するステップ６と、
上記液体圧力及び上記液体温度から算出される仮の液体密度と、上記各階級のボイド率から算出される液体の体積率及び上記液体保存量から算出される液体密度とを比較することにより、上記液体圧力を補正してステップ４〜６の処理を行うステップ７と、
上記ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、
上記仮の液体密度と、上記液体の体積率及び液体保存量から算出される液体密度とが等しくなったとき上記ステップ２〜８において算出された各計算値を出力するステップ９とを有し、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するコンピュータによって読み取り可能にプログラムを記憶したプログラム記録媒体において、
ステップ５は、上記液体圧力と上記液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、上記新規気泡数密度及び上記新規気泡質量から最小階級の上記気泡数密度及び上記気泡質量を補正するステップ５−１と、
繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、
上記各階級の気泡半径、上記液体圧力及び上記液体保存量から各階級の気泡圧力を算出するステップ５−３と、
上記各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、
上記液体圧力、上記液体温度、上記補正された各階級の気泡数密度及び気泡質量、上記各階級の気泡圧力、上記各階級の気泡温度、並びに、上記各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、
上記各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、上記各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、
上記蒸発量及び上記凝縮量から上記各階級の気泡数密度及び気泡質量を補正するステップ５−７と、
上記補正された各階級の気泡数密度及び気泡質量並びに上記各階級の気泡密度から算出される各階級の気泡半径を上記仮定した各階級の気泡半径とを比較し、判定結果に基づいて上記仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、
上記仮定した各階級の気泡半径と上記算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とするプログラム記録媒体。
【請求項８】
上記ステップ５−５では、上記各階級の気泡数密度及び気泡質量に基づいて各階級の平均気泡質量を算出し、上記各階級の平均気泡質量の値と、各階級の境界における質量の値とに基づいて上記各階級内の気泡数密度分布を算出することを特徴とする請求項７記載のプログラム記録媒体。
【請求項９】
上記ステップ５−５では、上記各階級の平均気泡質量の値ｍ１，ｋと、この階級の隣接する質量の軽い側の境界における気泡質量の値ｗｌｉｇｈｔ，ｋと、この階級の隣接する質量の重い側の境界における気泡質量の値ｗｈｅａｖｙ，ｋとから、下記の式（１３）〜式（１５）により各階級を分布１型、分布２型又は分布３型に分類し、
上記分布１型、分布２型、及び分布３型に分類された各階級の境界値ＮＧＷ及び端点ｗｆ，ｗｂの値を下記の式（１６）〜（１８）により算出することで各階級内の気泡数密度分布を算出することを特徴とする請求項７記載のプログラム記録媒体。
【数３】
2005222339.tif 000062
【請求項１０】
初期条件及び境界条件を入力するステップ１と、上記入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、上記液体保存量並びに上記各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、上記液体圧力及び上記液体温度から各階級のボイド率を算出すると共に、上記各階級の気泡数密度及び気泡質量を補正するステップ５と、上記補正された各階級の気泡数密度及び気泡質量から上記液体保存量を補正するステップ６と、上記液体圧力及び上記液体温度から算出される仮の液体密度と、上記各階級のボイド率から算出される液体の体積率及び上記液体保存量から算出される液体密度とを比較することにより、上記液体圧力を補正してステップ４〜６の処理を行うステップ７と、上記ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、上記仮の液体密度と、上記液体の体積率及び液体保存量から算出される液体密度とが等しくなったとき上記ステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムを記憶した記録部と、
上記記録部から上記解析プログラムを読み出して解析処理を実行する演算部とを備える解析装置において、
上記解析プログラムのステップ５は、上記液体圧力と上記液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、上記新規気泡数密度及び上記新規気泡質量から最小階級の上記気泡数密度及び上記気泡質量を補正するステップ５−１と、
繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、
上記各階級の気泡半径、上記液体圧力及び上記液体保存量から各階級の気泡圧力を算出するステップ５−３と、
上記各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、
上記液体圧力、上記液体温度、上記補正された各階級の気泡数密度及び気泡質量、上記各階級の気泡圧力、上記各階級の気泡温度、並びに、上記各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、
上記各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、上記各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、
上記蒸発量及び上記凝縮量から上記各階級の気泡数密度及び気泡質量を補正するステップ５−７と、
上記補正された各階級の気泡数密度及び気泡質量並びに上記各階級の気泡密度から算出される各階級の気泡半径を上記仮定した各階級の気泡半径とを比較し、判定結果に基づいて上記仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、
上記仮定した各階級の気泡半径と上記算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とする解析装置。
【請求項１１】
上記ステップ５−５では、上記各階級の気泡数密度及び気泡質量に基づいて各階級の平均気泡質量を算出し、上記各階級の平均気泡質量の値と、各階級の境界における質量の値とに基づいて上記各階級内の気泡数密度分布を算出することを特徴とする請求項１０記載の解析装置。
【請求項１２】
上記ステップ５−５では、上記各階級の平均気泡質量の値ｍ１，ｋと、この階級の隣接する質量の軽い側の境界における気泡質量の値ｗｌｉｇｈｔ，ｋと、この階級の隣接する質量の重い側の境界における気泡質量の値ｗｈｅａｖｙ，ｋとから、下記の式（１９）〜式（２１）により各階級を分布１型、分布２型又は分布３型に分類し、
上記分布１型、分布２型、及び分布３型に分類された各階級の境界値ＮＧＷ及び端点ｗｆ，ｗｂの値を下記の式（２２）〜（２４）により算出することで各階級内の気泡数密度分布を算出することを特徴とする請求項１０記載の解析装置。
【数４】
2005222339.tif 000063
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、管内、流体機械、物体まわりの流体に発生するキャビテーションを気泡径分布モデルを用いて、精度良く安定的に解析するのに適したキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置に関する。
【背景技術】
【０００２】
ロケットエンジン等に用いられる液体酸素、液体水素等のポンプにおいて、キャビテーション、即ち、飽和蒸気圧よりも局所的に圧力が下がってしまい部分的に気化する現象が発生している。
【０００３】
液体中に気体が混合することにより、音速が著しく低下してしまい、圧力情報が伝わりにくくなる。このことにより、気体が混合していない部分はポンプにより吸い込みが可能であり、気体が混合している部分は吸い込めなくなる、いわゆる脈動現象が発生してしまう。この脈動現象により、圧力振動が羽根を振動させ、破壊させるおそれがある。脈動現象による破壊を防止するためには、キャビテーションの挙動を正確に捉える必要がある。
【０００４】
従来、管内、流体機械及び物体回りに発生するキャビテーションにおいて、どの程度圧力が下がり、それにより、どの程度気化し、どの程度音速が低下し、どの程度圧力が伝わりにくいのかを解析することは非常に困難とされていた。
【０００５】
そこで、液体中に発生するキャビテーションを球形気泡から構成される気泡群として取り扱い、その気泡群の挙動を捉えることでキャビテーションの特徴を把握しようとする技術が提案されている（非特許文献１）。
【０００６】
しかし、一般的な流れ場では、キャビーテーションは、一点から気泡が発生するものではなく、様々な地点から発生した気泡が混在する複雑な現象を伴うものであり、キャビテーションの解析は、精度、安定性に乏しくなってしまったり、計算機の負荷が大きくなる等の問題を抱えている。
【０００７】
【非特許文献１】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
【発明の開示】
【発明が解決しようとする課題】
【０００８】
本発明の目的は、流体中に発生するキャビテーションを、球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、階級ごとに気泡の挙動を取り扱う気泡径分布モデルを採用したキャビテーション解析において、その精度、安定性を高めることができ、計算負荷を低減することができるキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置を提供することにある。
【課題を解決するための手段】
【０００９】
この目的を達成するため、本発明に係るキャビテーション解析方法は、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とを備え、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析方法において、ステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とする。
【００１０】
また、この目的を達成するため、本発明に係るキャビテーション解析プログラムは、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムにおいて、ステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とする。
【００１１】
また、この目的を達成するため、本発明に係るプログラム記録媒体は、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とを有し、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するコンピュータによって読み取り可能にプログラムを記憶したプログラム記録媒体において、ステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とする。
【００１２】
また、この目的を達成するため、本発明に係る解析装置は、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムを記憶した記録部と、記録部から解析プログラムを読み出して解析処理を実行する演算部とを備える解析装置において、解析プログラムのステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有することを特徴とする。
【００１３】
以上のように、本発明は、液体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するものであり、各階級内の気泡数密度分布を算出し、この気泡数密度分布から各階級の境界での気泡数密度を算出し、この気泡数密度から各階級を通過する蒸発量・凝縮量を算出するので、数値拡散を防ぐことができ、解析の精度、安定性を高めることができると共に、高い汎用性を有するものである。また、本発明は、気泡の成長・収縮運動と液体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、解析における時間刻みが小さくなることを防ぐとともに、計算安定性を増すことができる。
【発明の効果】
【００１４】
本発明によれば、液体圧力、液体温度、各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を考慮することにより、数値拡散を防ぐことができ、気泡の蒸発・凝縮計算の精度を高めることができるとともに、階級間隔を不等間隔にした場合の計算安定性を高めることができる。よって、本発明は、精度及び安定性の高いキャビテーション解析を可能とするものであり、汎用性を高めることを可能とするものである。
【００１５】
また、本発明によれば、気泡の成長・収縮計算と流体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、時間刻みを大きく取ることができ、キャビテーション解析の計算安定性を増すことができ、よって、計算負荷を低減することを可能とするものである。
【発明を実施するための最良の形態】
【００１６】
以下、本発明が適用された解析装置の実施の形態を図面を参照して説明する。
【００１７】
本発明を適用した解析装置１は、図１に示すように、入力部２と、演算部３と、メモリ４と、表示部５と、出力部６と、記憶部７とを備える。この解析装置１では、入力部２を介して入力が行われ、記憶部７に記憶されたキャビテーション解析プログラムからメモリ４に必要な係数及び式を読み出し、演算部３で演算することによりキャビテーション解析を行う。解析結果は、表示部５に表示したり、出力部６から出力することができる。
【００１８】
解析装置１の記憶部７に記憶されたキャビテーション解析プログラムは、初期条件及び境界条件を入力するステップ１と、入力された初期条件及び境界条件から、時間刻み、仮の液体保存量、各階級の気泡数密度及び気泡質量、並びに、気泡の発生からの経過時間を算出するステップ２と、繰り返し計算の初期値として液体圧力を与えるステップ３と、液体保存量並びに各階級の気泡数密度及び気泡質量から液体温度を算出するステップ４と、液体圧力及び液体温度から各階級のボイド率を算出すると共に、各階級の気泡数密度及び気泡質量を補正するステップ５と、補正された各階級の気泡数密度及び気泡質量から液体保存量を補正するステップ６と、液体圧力及び液体温度から算出される仮の液体密度と、各階級のボイド率から算出される液体の体積率及び液体保存量から算出される液体密度とを比較することにより、液体圧力を補正してステップ４〜６の処理を行うステップ７と、ステップ２〜７において算出された各計算値の境界値を上記境界条件により補正するステップ８と、仮の液体密度と、液体の体積率及び液体保存量から算出される液体密度とが等しくなったときステップ２〜８において算出された各計算値を出力するステップ９とをコンピュータによって実行可能とし、流体中に発生するキャビテーションを球形気泡からなる気泡群として取り扱い、気泡をその質量により階級分けし、各階級ごとに気泡の挙動を解析するキャビテーション解析プログラムであって、ステップ５は、液体圧力と液体温度から新規発生する新規気泡半径、新規気泡数密度及び新規気泡質量を算出するとともに、新規気泡数密度及び新規気泡質量から最小階級の気泡数密度及び気泡質量を補正するステップ５−１と、繰り返し計算の初期値として各階級の気泡半径を仮定するステップ５−２と、各階級の気泡半径、液体圧力及び液体保存量から各階級の気泡圧力を算出するステップ５−３と、各階級の気泡圧力から各階級の気泡温度及び気泡密度を算出するステップ５−４と、液体圧力、液体温度、補正された各階級の気泡数密度及び気泡質量、各階級の気泡圧力、各階級の気泡温度、並びに、各階級の気泡密度に基づいて各階級内の気泡数密度分布を算出するステップ５−５と、各階級内の気泡数密度分布から各階級の境界での気泡数密度を算出し、各階級の境界での気泡数密度から各階級境界を通過する蒸発量及び凝縮量を算出するステップ５−６と、蒸発量及び凝縮量から各階級の気泡数密度及び気泡質量を補正するステップ５−７と、補正された各階級の気泡数密度及び気泡質量並びに各階級の気泡密度から算出される各階級の気泡半径を仮定した各階級の気泡半径とを比較し、判定結果に基づいて仮定した各階級の気泡半径を補正してステップ５−３〜５−７の処理を行うステップ５−８と、仮定した各階級の気泡半径と算出された各階級の気泡半径とが等しくなったとき、この各階級の気泡半径から各階級のボイド率を算出するステップ５−９とを有するものである。
【００１９】
このキャビテーション解析プログラムを記憶部７に有する解析装置１は、液体中に発生するキャビテーションを球形気泡から構成される気泡群として取り扱い、その気泡群の挙動を捉えることでキャビテーションの特徴と掴もうとするものであり、キャビテーション気泡をその質量に応じて階級化し、各階級ごとに気泡半径、圧力、温度、蒸発・凝縮量、液体とのスリップ速度などを計算することで、キャビテーション気泡を厳密に捉え、解析しようとするものである。
【００２０】
以下に、本発明を適用した解析装置１を用いて、図１３に示すような先細末広ノズルにおける２次元計算の例により、キャビテーション解析を行うキャビテーション解析方法の例を示す。
【００２１】
まず、全体の計算手順に先立ち、「１．基礎式」、「２．物性値」、「３．気泡を含む液体の音速」、「４．気泡新規発生」、「５．気泡質量に基づいた階級の離散化」、「６．気泡の並進運動」、「７．液体の圧力計算と気泡の半径計算のカップリング」及び「８．時間刻み」について説明する。尚、条件としては、Ｓｉｍｏｎｅａｕ  ａｎｄ  Ｈｅｎｄｒｉｃｋｓの実験条件と同一とする。以下の説明に用いる記号の対応表を表１、表２に示す。
【００２２】
【表１】
2005222339.tif 000003
【００２３】
【表２】
2005222339.tif 000004
【００２４】
１．基礎式
このキャビテーションモデルを採用した計算コードを開発する。基礎式は２次元の例を示す。一般座標系（ξ，η）を用いて以下のように表せる。単位体積あたりの液相部質量αＬρＬ、運動量ｘ方向成分αＬρＬｕＬ、運動量ｙ方向成分αＬρＬνＬ、全エネルギーαＬρＬｅＬの４式に加え、単位体積あたりの気泡数密度の確率密度関数ＮＧｗ、気泡発生からの平均経過時間ｔｅｌａｐｓｅも気泡移流速度で移流するとして計算した。また、気泡関係諸量は気泡の成長・収縮に伴い気泡質量軸方向（ｗ）にも移流するとして計算した。
【００２５】
【数１】
2005222339.tif 000005
【００２６】
２．物性値
キャビテーションは液体温度に敏感なため、正確な計算には温度に対する影響を厳密に扱うことが欠かせない。そこで、温度による物性値の変化を考慮した。以下に示す物性値を温度の関数として取り扱った。これらは全て物性値表のデータを連続な関数として近似したものである。これらのデータは基本的には想定温度範囲内（三重点温度〜臨界点温度）のデータを基にしているが、計算中にプログラムにより想定温度範囲外の温度を一時的に参照してしまう場合がある。このときにエラーが発生すると計算の安定性が著しく損なわれるため、想定温度範囲内の物性値と連続な適当な物性値を温度範囲外に対しても与えるようにした。ここで与えられる適当な物性値は計算過程で必要な値であって最終的な収束値には成り得ないため計算の信頼性への影響は少ないと思われる。
【００２７】
気液界面の物性値を表３に、液体の物性値を表４に、蒸気の物性値を表５に示す。また、計算上必要な逆関数を表６に示す。逆関数は安定かつ正確な計算のため元の値に正確に戻す必要がある。そのためこれらの逆関数を必要とする関数に関しては誤差の管理を厳密に行った。
【００２８】
【表３】
2005222339.tif 000006
【００２９】
【表４】
2005222339.tif 000007
【００３０】
【表５】
2005222339.tif 000008
【００３１】
【表６】
2005222339.tif 000009
【００３２】
３．気泡を含む液体の音速
液体は気相を内部に含むと圧縮性が単相時に比べて著しく大きくなり、図２に示すように音速が大きく低下する。ここで、図２は、２相流体の平均音速とボイド率の関係を示す図であり、「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」に記載されたものである。高速キャビテーションではこの効果が無視できないため、液体が気泡を内部に含んで音速が低下する挙動をモデルを導入して評価した。そのモデルを以下に示す。
【００３３】
気相は、液相中に均質に分散した気泡であると考える。このとき２相流体の平均音速は次式となる。
【００３４】
【数２】
2005222339.tif 000010
【００３５】
ここで、Ｐは２相流体の圧力（ここでは気液間の圧力差は考慮しない）で、ρＶは２相流の平均密度であり、式（３０）となる。
【００３６】
【数３】
2005222339.tif 000011
【００３７】
αはボイド率（気相部分の体積分率）である。一方、クウォリティ（気相部分の質量分率）は次式となる。
【００３８】
【数４】
2005222339.tif 000012
【００３９】
式（３０）、式（３１）よりクウォリティχは、式（３２）となる。
【００４０】
【数５】
2005222339.tif 000013
【００４１】
ここで、式（３１）をＰについて微分すると式（３３）となる。
【００４２】
【数６】
2005222339.tif 000014
【００４３】
よって、式（２９）、式（３３）より式（３４）が得られる。
【００４４】
【数７】
2005222339.tif 000015
【００４５】
この式（３４）を再度式（２９）を用いて整理すると次式となる。
【００４６】
【数８】
2005222339.tif 000016
【００４７】
本モデルの採用により図２に示すボイド率と音速の関係を得ることができる。
【００４８】
４．気泡新規発生
新規発泡は主として不純物上に発生するものと壁面上に発生するものがある。本モデルでは両初期気泡径とも次式（３６）で示される過飽和度に反比例するものとして与えた。
【００４９】
【数９】
2005222339.tif 000017
【００５０】
発泡頻度は不純物上が次式（３７）で示されるｎ０であり、壁面上においては、次式（３８）で示されるｎ０ｗａｌｌである。
【００５１】
【数１０】
2005222339.tif 000018
【００５２】
【数１１】
2005222339.tif 000019
【００５３】
ここで、式中で用いられる壁面からの単位体積あたりの発泡数ｎＳｗａｌｌを求める次式（３９）の面積Ａｗａｌｌと壁面近傍の空間Ｖｏｌｗａｌｌの与え方である。
【００５４】
【数１２】
2005222339.tif 000020
【００５５】
図３に示すような領域で気泡が発生する場合を考える。面積Ａｗａｌｌは壁面隣接セル上の壁面積（２Ｄでは壁長さ）を与える。壁面近傍の空間Ｖｏｌｗａｌｌは式（３６）で求められる初期半径Ｒｏが発生し十分離脱できる領域である必要がある。壁面隣接セルの体積Ｖｏｌｃｅｌｌ，１、その１つ壁面から遠いセルの体積Ｖｏｌｃｅｌｌ，２、壁面からｊ番目のセル体積Ｖｏｌｃｅｌｌ，ｊとする。また、壁面隣接セルの壁面から一番遠い点までの距離ｄｊ＝１、壁面からｊ番目のセルの壁面から一番遠い点までの距離ｄｊとする。このときにｄｊ−１≦２Ｒｏ≦ｄｊとなる壁面から１番目からｊ番目のセルの中で発泡するとする。すなわち、式（４０）とし、壁面から１番目からｊ番目の全てのセルの中でｎＳｗａｌｌの発泡数で発生するとした。
【００５６】
【数１３】
2005222339.tif 000021
【００５７】
尚、図３に示すような領域で気泡が発生する場合のボイド率の取り方について説明する。ｄｊ≧２Ｒ０、かつ、ｄｊ−１＜２Ｒ０であるｊ以下のセルでは、ｊ以下のセルの平均的ボイド率を使用する。すなわち、図４に示すように、ｊ＝２では、実際のボイド率はαｊ＝２であるが、本発明の計算では、ボイド率α＝αｊ＝１：２が用いられる。また、図５に示すように、ｊ＝１では、実際のボイド率はαｊ＝１であるが、本発明の計算では、ボイド率α＝αｊ＝１：２が用いられる。ここで、ｊ＝１：２の平均ボイド率：αｊ＝１：２は、図６に示すように、ｊ＝１及びｊ＝２の両方の領域を１つの領域とした場合のボイド率を算出することにより得られる。
【００５８】
この取り扱いにより壁面近傍のグリッドを適度に小さくしておけば、壁面近傍の格子点の細かさに対する気泡発泡のグリッド依存性はほとんどなくなる。
【００５９】
気泡発生数は次式（４１）で示されるΠＧ，気泡発生質量は次式（４２）で示されるΛＧである。
【００６０】
【数１４】
2005222339.tif 000022
【００６１】
【数１５】
2005222339.tif 000023
【００６２】
新規発泡気泡は一番小さな気泡質量階級に属するため既存の気泡の保存量に付加される。発泡前の領域に存在する単位体積あたりの気泡数をｎＧ＿ｏｌｄ，質量をｍＧ＿ｏｌｄとすると、発泡後の単位体積あたりの気泡数ｎＧ＿ｎｅｗ，質量ｍＧ＿ｎｅｗは次式となる。
【００６３】
【数１６】
2005222339.tif 000024
【００６４】
５．気泡質量に基づいた階級の離散化
気泡径分布モデルでは、空間の次元に加えて１個あたりの気泡質量の次元を追加する。以下その次元の向きをｗ方向と呼ぶことにする。ｗ軸の正方向が質量の増加を意味するようにその向きを定義する。ある空間におけるｗ軸方向の気泡数密度の確率密度関数ＮＧｗの分布を図７に示す。気泡群が蒸発に伴い成長すると気泡質量が増えるためw軸の正の方向に移動し、気泡が凝縮すると気泡質量が減るためw軸の負の方向に移動する。すなわち、空間的移動とは独立に気泡の蒸発凝縮に伴い気泡質量軸方向に気泡が運動するという新しい概念を導入する。ここで空間的に静止した気泡1個あたりの蒸発量γ１とするときの気泡数密度関数ＮＧｗは、次式となる。
【００６５】
【数１７】
2005222339.tif 000025
【００６６】
この式（４５）は、気泡の成長がｗ軸方向に移流速度γ１で運動していることを表している。過熱液中（気泡成長の場合）または過冷液中（気泡凝縮の場合）に存在する気泡周りに生成される温度境界層から供給される熱量によって駆動される蒸発・凝縮量がγ１である。過熱・過冷液中に投入されてからｔｅｌａｐｓｅ後の球形気泡周りに形成される温度境界層の気泡界面上での温度勾配は次式となる。
【００６７】
【数１８】
2005222339.tif 000026
【００６８】
よって、気泡界面を通過する熱量はフーリエの法則より次式となる。すなわち、熱伝導により伝えられる熱流束は，温度勾配ｄＴ／ｄｒに比例するという法則による。尚、比例定数は、熱伝導率ｋＬとなる。
【００６９】
【数１９】
2005222339.tif 000027
【００７０】
この熱量が蒸発・凝縮に用いられるので蒸発・凝縮量は、式（４８）となる。
【００７１】
【数２０】
2005222339.tif 000028
【００７２】
以上は静止気泡について述べたが空間移動する気泡においても一般性は失われない。すなわち、蒸発を伴いながら空間移動する気泡数密度関数の時間変化は次式となる。
【００７３】
【数２１】
2005222339.tif 000029
【００７４】
ここで、δ（ｗ−ｗ０）は、式（５０）で表されるデルタ関数を示す。
【００７５】
【数２２】
2005222339.tif 000030
【００７６】
デルタ関数とは、カッコ内の数が０のときを含んで積分すると１となり、カッコ内の数が０でないとき０となるパルス関数である。ここで、ｗ０は、新規発生気泡の質量を表す。
【００７７】
【数２３】
2005222339.tif 000031
【００７８】
２次元の場合∇，ＶＧは、次式となる。
【００７９】
【数２４】
2005222339.tif 000032
【００８０】
【数２５】
2005222339.tif 000033
【００８１】
同様に、気泡発生からの経過時間であるｔｅｌａｐｓｅも気泡の移動に伴って移流すると考えられるため、式（５４）のように表せる。
【００８２】
【数２６】
2005222339.tif 000034
【００８３】
ここで右辺の１は、１タイムステップごとに時間刻み分だけ経過時間が増えることを表している。
【００８４】
実際の計算では、気泡数密度関数ＮＧｗをｗで積分して求めた階級内の単位体積あたりの気泡数密度ｎＧ，ｋと、気泡数密度関数ＮＧｗにｗをかけてからｗで積分して求めた階級内の単位体積あたりの気泡質量ｍＧ，ｋとを用いて行う。ここで、ｗは１個あたりの質量を表すものである。気泡数密度ｎＧ，ｋ及び気泡質量ｍＧ，ｋは、式（５５）、式（５６）のように定義される。
【００８５】
【数２７】
2005222339.tif 000035
【００８６】
この気泡数密度、気泡質量を気泡の質量軸をある質量範囲で定義された階級に分割する。ある階級ｋについて考える。階級の軽い方の境界をｌｉｇｈｔ境界，重い方をｈｅａｖｙ境界と呼び、各々の気泡質量をｗｌｉｇｈｔ，ｋ，ｗｈｅａｖｙ，ｋとすると、階級内に存在する単位質量あたりの気泡数密度ｎＧ，ｋ、気泡質量ｍＧ，ｋは次式で定義される。
【００８７】
【数２８】
2005222339.tif 000036
【００８８】
ここで、ｍ１，ｋは、階級内の気泡１つあたりの平均質量を表し、次式（５９）となる。
【００８９】
【数２９】
2005222339.tif 000037
【００９０】
この平均質量ｍ１，ｋは、階級の境界質量ｗｌｉｇｈｔ，ｋ，ｗｈｅａｖｙ，ｋと以下の関係を満たさなければならない。
【００９１】
【数３０】
2005222339.tif 000038
【００９２】
本式（６０）は、実際に計算する際に、気泡の階級分けが適切に行われているのかをチェックするのに有効な条件式となる。
【００９３】
基礎式（４９）を単位体積あたりの気泡数密度ｎＧ，ｋと単位体積あたりの気泡質量ｍＧ，ｋの定義式（５７）、（５８）を用いて変形すると次式のように表せる。
【００９４】
【数３１】
2005222339.tif 000039
【００９５】
さらに、階級ｋ内の気泡数密度関数ＮＧｗの分布について考える。ＮＧｗの分布で一番単純なものは図８に示すように、階級内で全て同一の値を採る均一な分布であると仮定したものである。このとき、ＮＧｗは階級内でどこでも同じ値で、かつ、平均値を取るため、次式（６３）となる。
【００９６】
【数３２】
2005222339.tif 000040
【００９７】
本式（６３）は、式（５７）を満たす。しかし、式（５８）を満たすのは、次式（６４）のときのみである。
【００９８】
【数３３】
2005222339.tif 000041
【００９９】
階級ｋ内では常に式（５７）、式（５８）を満たさなければならないという条件を満たさない。これは気泡数密度関数ＮＧｗの分布が図８に示すように階級内で全て同一の値を採るという仮定が誤っていることを示す。そこで、図９に示すように重心をｍ１，ｋとする３つの１次関数分布を持つと仮定し、分類した。最小階級用の分布１型１１、中間階級用の分布２型１２、最大階級用の分布３型１３である。どの条件でどの分布を持つかを式（６５）〜（６７）に示す。
【０１００】
【数３４】
2005222339.tif 000042
【０１０１】
また、それぞれの境界値ＮＧｗ（ｗｌｉｇｈｔ，ｋ），ＮＧｗ（ｗｈｅａｖｙ，ｋ）及び端点ｗｆ，ｗｂの値は、以下の式（６８）〜（７０）の通りとなる。
【０１０２】
【数３５】
2005222339.tif 000043
【０１０３】
この分布は、式（５７）だけでなく式（５８）も満たす。この分布は、さらに数値計算上の大きな効果を持っている。それは、階級内の気泡がある程度成長しないとそれより大きな階級に気泡が全く移動しないという点である。最大階級の分布である分布３型を見ると重い側の境界において気泡数密度関数ＮＧｗ（ｗｈｅａｖｙ，ｋ）＝０であることから分かる。不要な数値拡散を抑制する効果をもち、精度の向上をもたらす。
【０１０４】
６．気泡の並進運動
次式を用いて時間発展的にスリップ速度ＶＳ＝ＶＧ−ＶＬを解く。
【０１０５】
【数３６】
2005222339.tif 000044
【０１０６】
流れ場の計算から求められたＶＬを基準に気相速度ＶＧを求める。尚、式（７１）において、右辺第１項は、スリップ速度ＶＳの気泡と液体の間に働く力、すなわち、気液間抗力を表すものである。気液間抗力は、抗力係数ＣＤ、気泡前面投影面積πＲ２、及び、スリップ速度の動圧ρＬ（｜ＶＳ｜ＶＳ）／２の積である。また、式（７１）の右辺第２項は、圧力勾配中に存在する気泡が受ける力で、いわゆる浮力にあたるものである。
【０１０７】
７．液体の圧力計算と気泡の半径計算のカップリング
その計算点においてあらゆる階級に気泡の存在しない純液体領域（液体体積率αＬ＝１）の温度ＴＬ、圧力ＰＬは流れ場の計算から求められた密度ρＬ、内部エネルギーｅＳＬを用いて、式（７２）、式（７３）から一義に定まる。
【０１０８】
【数３７】
2005222339.tif 000045
【０１０９】
一方、気相の存在する２相領域の温度ＴＬは液相部分の各保存量αＬρＬ，αＬρＬｕＬ，αＬρＬνＬ，αＬρＬｅＬより内部エネルギーｅＳＬ＝ｅＬ−０．５（ｕ２Ｌ＋ν２Ｌ）が求まるので式（７２）を用いて定まる。しかし、圧力ＰＬは液体体積率αＬが定まらないと密度ρＬが定まらず、加えて液体体積率αＬ自身が圧力ＰＬの関数となっており収束計算を必要とする。以下に圧力ＰＬと液体体積率αＬの収束計算ループの過程について記す。この液体圧力計算のフローチャートを図１１に示す。はじめに仮の液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙを設定する。本計算ではＰＬ＿ｔｅｍｐｏｒａｒｙ＝０を初期値とし、徐々にＰＬ＿ｔｅｍｐｏｒａｒｙを大きくして行き、収束したＰＬ＿ｔｅｍｐｏｒａｒｙを圧力ＰＬとする。
【０１１０】
ここで、仮の液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙに対する液体体積率αＬを求める手順を示す。液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙに平衡する気泡半径Ｒｋを求める。この気泡半径Ｒｋも新規気泡発生、蒸発・凝縮、消滅などの現象が複雑に絡み合っているため収束計算を必要とする。この気泡半径計算のフローチャートを図１２に示す。気泡半径Ｒｋを求める小ループでは着目しているコントロールボリューム全てが階級ｋの気泡で占められた場合の気泡半径Ｒｋの最大値Ｒｍａｘ，ｋを初期値として徐々にＲｋを減少させて収束させる。各階級ごとの気泡数密度ｎＧ，ｋ，気泡質量ｍＧ，ｋと求める。この計算には液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと液体飽和圧力Ｐｓａｔ（ＴＬ）とから求められる新規気泡発生の効果、液体温度ＴＬと気泡温度ＴＧ，ｋから求められる蒸発・凝縮量の効果、気泡温度ＴＧ，ｋが臨界温度ＴＣ以上、つまり、超臨界状態になったとき気泡が液体に溶解する効果を陰的に含んでいる。次に既存の気泡では、一般に次式（７４）で表されるレイリープレセット式（Ｒａｙｌｅｉｇｈ  Ｐｌｅｓｓｅｔ  Ｅｑｕａｔｉｏｎ）が適用されることが知られている。
【０１１１】
【数３８】
2005222339.tif 000046
【０１１２】
発生直後や崩壊直前の小さな気泡には、界面半径方向の慣性力の効果が大きいため右辺第４〜６項の時間微分項の効果が大きいが、大きな気泡では界面の半径方向の慣性力の効果が相対的に小さくなるため右辺第４〜６項の時間微分項の効果が小さくなる。そこで、１タイムステップ前の気泡半径Ｒｏｌｄ，ｋ、２タイムステップ前の気泡半径Ｒｏｌｄｏｌｄ，ｋとするとき、Ｒｏｌｄ，ｋが一定の気泡径Ｒ１以下（実際の計算ではＲ１=１μｍ）では、時間について離散化した次式（７５）に示すように時間微分項を１−Ｒｏｌｄ，ｋ／Ｒ１倍して気泡圧力ＰＧ，ｋを計算し、一定の気泡径Ｒ１以上では次式（７５）に示すように時間微分項を無視して気泡圧力ＰＧ，ｋを計算する。
【０１１３】
【数３９】
2005222339.tif 000047
【０１１４】
この取り扱いにより計算負荷を低減しつつ、気泡径Ｒ１前後で連続的に時間微分項の効果を導入できるようになった。
【０１１５】
さらに既存気泡のない計算セルに新規に気泡が発生した場合にはスリップが無視できるので、次式（７６）を用いて気泡圧力ＰＧ，ｋを計算する。
【０１１６】
【数４０】
2005222339.tif 000048
【０１１７】
これらの気泡圧力ＰＧ，ｋに対応する気泡温度ＴＧ，ｋ、気泡密度ρＧ，ｋが飽和条件より求められる。
【０１１８】
【数４１】
2005222339.tif 000049
【０１１９】
ここで求められたρＧ，ｋを用いて、先ほど求めた気泡数密度ｎＧ，ｋ、気泡質量ｍＧ，ｋが次式（７９）を満たすまで、気泡半径Ｒｋを求める小ループを繰り返す。
【０１２０】
【数４２】
2005222339.tif 000050
【０１２１】
よって、階級ｋのボイド率αＧ，ｋは次式（８０）と表せる。
【０１２２】
【数４３】
2005222339.tif 000051
【０１２３】
次に液体温度ＴＬ，仮の液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙに対して求めた階級ｋのボイド率αＧ，ｋを用いて真の液体圧力ＰＬを求める過程を示す。仮の液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙに対する液体体積率αＬ＿ｔｅｍｐｏｒａｒｙと液体密度ρＬ＿ｔｅｍｐｏｒａｒｙとが、既知の保存量αＬρＬとの関係が式（８３）となるまで、液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙを求めるループを繰り返す。
【０１２４】
【数４４】
2005222339.tif 000052
【０１２５】
【数４５】
2005222339.tif 000053
【０１２６】
【数４６】
2005222339.tif 000054
【０１２７】
以上のような繰り返し計算をしているとある階級ｋの気泡が消滅することがある。そのときの体積あたりの気泡質量はｍＧ，ｋなので各保存式より、次式（８４）〜（８７）となる。
【０１２８】
【数４７】
2005222339.tif 000055
【０１２９】
この過程を液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙを求めるループ中の最後に行う。
【０１３０】
８．時間刻み
時間刻みは移流項に対してＣＦＬ条件を満たす必要がある。また、圧力勾配による加速項、ｓｔｉｆｆな生成項と粘性項に対する時間制限も満足しなくてはならない。ここで、ｓｔｉｆｆな生成項とは、計算を硬直化させる大きな生成項である。すなわち、生成項、即ち式（１）のＳが大きく、かつ、時間刻みが大きい場合、計算が不安定となる。このことを防ぐために、時間刻みを小さくする必要があり、計算を硬直化（ｓｔｉｆｆｎｅｓｓ）させる。このような生成項をｓｔｉｆｆな生成項である。
【０１３１】
【数４８】
2005222339.tif 000056
【０１３２】
ここで，ＡＬｘ，ＡＧｘ，ＡＬｙ，ＡＧｙは、各相各方向の圧力勾配項、粘性項および生成項などの時間微分値である。
【０１３３】
【数４９】
2005222339.tif 000057
【０１３４】
ここでＤｉｍは次元数。１次元ならＤｉｍ＝１，２次元ならＤｉｍ＝２、３次元ならＤｉｍ＝３となる。
【０１３５】
次に、本発明を適用した解析装置１を用いたキャビテーション解析の全体の計算手順について説明する。表７に、以下で用いられる計算条件を示す。尚、本条件は、Ｓｉｍｏｎｅａｕ  ａｎｄ  Ｈｅｎｄｒｉｃｋｓの実験条件と同一とした。
【０１３６】
【表７】
2005222339.tif 000058
【０１３７】
まず、図１０に示すように、ステップＳ１〜ステップＳ２の流れ場の計算について説明する。ステップＳ１−１では、計算対象の格子を図１３に示すように読み込む。ここでは、図１３に示すよう２次元で、流れ方向１４３点、高さ方向１２点とした。
【０１３８】
ステップＳ１−２では、初期条件を表７のように、また、境界条件を表８のように設定する。本解析においては、表７に示すように、壁面を断熱壁として設定される。尚、本条件は、Ｓｉｍｏｎｅａｕ  ａｎｄ  Ｈｅｎｄｒｉｃｋｓの実験条件と同一とした。
【０１３９】
【表８】
2005222339.tif 000059
【０１４０】
ステップＳ２−１では、液体温度、気泡平均温度、液体体積率、ボイド率から上述の式（３５）を用いて、２相平均音速ｃＶを計算した。
【０１４１】
ステップＳ２−２では、式（１）を用いて、液体の各方向速度より粘性応力項ＦＶを計算する。
【０１４２】
ステップＳ２−３では、式（１）の各基礎式の空間移流成分の各計算点での時間微分値を計算し、空間移流項Ｅ，Ｆを算出する。尚、空間移流成分はＭＵＳＣＬ法（Ｍｏｎｏｔｏｎｉｃ  Ｕｐｓｔｒｅａｍ−ｃｅｎｔｅｒｅｄ  Ｓｃｈｅｍｅ  ｆｏｒ  Ｃｏｎｓｅｒｖａｔｉｏｎ  Ｌａｗｓ）により精度を高めたＳＨＵＳ（Ｓｉｍｐｌｅ  Ｈｉｇｈ−ｒｅｓｏｌｕｔｉｏｎ  Ｕｐｗｉｎｄ  Ｓｃｈｅｍｅ）を利用し最大で空間３次精度を達成しつつＴＶＤ（Ｔｏｔａｌ  Ｖａｒｉａｔｉｏｎ  Ｄｉｍｉｎｉｓｈｉｎｇ）性を保った風上差分を採用することで、衝撃波等の不連続面も安定的に扱うことが可能となった。
【０１４３】
ステップＳ２−４では、式（７１）に示すように、液体の各方向成分、気泡半径Ｒから、気液間抗力及び各階級の気泡と液相のスリップ速度ＶＳを計算する。ここで、液体の各方向成分とは、スリップ速度ＶＳ＝（ｕＳ，ｖＳ）、気相速度ＶＧ＝（ｕＧ，ｖＧ）、液層速度ＶＬ＝（ｕＬ，ｖＬ）である。また、ｘ方向成分は、ｕＳ，ｕＧ，ｕＬであり、ｙ方向成分は、ｖＳ，ｖＧ，ｖＬである。尚、スリップ速度、すなわち、空間微分∇ＶＧ，∇ＶＬの計算には、風上１次精度の差分法を用いた。
【０１４４】
ステップＳ２−５では、式（１）に示すように、上述のステップＳ２−２〜Ｓ２−４で算出した、粘性応力項ＦＶ、空間移流項Ｅ，Ｆ、気液間抗力を用いて液体保存量（液相部質量、Ｘ方向成分運動量、Ｙ方向成分運動量及び全エネルギー量）及び各階級の気泡の時間微分値を計算する。
【０１４５】
ステップＳ２−６では、式（８８）〜（９１）に示すように、上述のステップ７の結果からＣＦＬ（Ｃｏｕｒａｎｔ  Ｆｒｉｅｄｒｉｃｈ  Ｌｅｗｙ）条件を満たす時間刻みを決定する。
【０１４６】
ステップＳ２−７では、ステップＳ２−５で算出した液体保存量及び各階級の気泡の時間微分値とステップＳ２−６で算出した時間刻みを用いて、時間進行させ、７つの保存量を計算する。すなわち、液体の保存量である単位体積当たりの液相部質量αＬρＬ、Ｘ方向成分運動量αＬρＬｕＬ、Ｙ方向成分運動量αＬρＬνＬ及び全エネルギーαＬρＬｅＬ、並びに、気泡数密度関数ＮＧｗから式（５７）、式（５８）で求められる各階級の仮の気泡保存量である各階級の仮の気泡数密度ｎＧ，ｋ及び気泡質量ｍＧ，ｋが算出されると共に気泡の発生からの経過時間ｔｅｌａｐｓｅ，ｋが決定する。
【０１４７】
次に、図１１に示すように、ステップＳ３〜ステップＳ７の上述のステップＳ２で算出した７つの保存量を用いて液体圧力および体積率を繰り返し計算する方法について説明する。ここでは、気泡の成長・収縮計算と液体の圧力・温度計算の計算安定性を高めるため両者をカップリングして繰り返し計算、いわゆる陰的計算を行う。
【０１４８】
ステップＳ３では、繰り返し計算の初期値として、液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙを適当に仮定する。
【０１４９】
ステップＳ４では、式（７２）を用いて、単位体積当たりの液相部質量αＬρＬ、Ｘ方向成分運動量αＬρＬｕＬ、Ｙ方向成分運動量αＬρＬνＬ及び全エネルギーαＬρＬｅＬから液体の内部エネルギーｅＳＬを求め、この値と物性値から液体温度ＴＬを算出する。
【０１５０】
ステップＳ５では、ステップＳ２−７で算出した７つの保存量と、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと、ステップＳ４で算出した液体温度ＴＬとを用いて気泡半径の繰り返し計算を行い、各階級のボイド率αＧ，ｋを算出すると共に、各階級の気泡数密度ｎＧ，ｋ及び気泡質量ｍＧ，ｋを補正する。
【０１５１】
ステップＳ６では、式（８４）〜（８７）に示すように、後述するステップ５−１で算出した新規気泡発生量並びにステップＳ５−６で算出された蒸発量・凝縮量ΓＩ及び気泡崩壊量により液体保存量（単位体積当たりの液相部質量αＬρＬ、Ｘ方向成分運動量αＬρＬｕＬ、Ｙ方向成分運動量αＬρＬνＬ及び全エネルギーαＬρＬｅＬ）を補正する。
【０１５２】
ステップＳ７では、まず、式（８２）に示すように、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと、ステップＳ４で算出された液体温度ＴＬから仮の液体密度ρＬ＿ｔｅｍｐｏｒａｒｙを算出する。次に、ステップＳ５で算出した階級ｋのボイド率αＧ，ｋから式（８１）を用いて液体体積率αＬ＿ｔｅｍｐｏｒａｒｙを算出する。そして、この仮の液体密度ρＬ＿ｔｅｍｐｏｒａｒｙと液体体積率αＬ＿ｔｅｍｐｏｒａｒｙとから計算される仮の液体保存量αＬ＿ｔｅｍｐｏｒａｒｙρＬ＿ｔｅｍｐｏｒａｒｙと液体保存量αＬρＬとが式（８３）で示される関係を満たすかどうかを比較して確かめる。
【０１５３】
尚、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと、ステップＳ４で算出された液体温度ＴＬから仮の液体密度ρＬ＿ｔｅｍｐｏｒａｒｙを算出すると共に、ステップＳ５で算出した階級ｋのボイド率αＧ，ｋから液体体積率αＬ＿ｔｅｍｐｏｒａｒｙを算出し、この液体体積率αＬ＿ｔｅｍｐｏｒａｒｙと仮の液体保存量αＬρＬとから算出される液体密度ρＬとを比較してもよい。
【０１５４】
ステップＳ５により算出されるαＬ＿ｔｅｍｐｏｒａｒｙρＬ＿ｔｅｍｐｏｒａｒｙとステップＳ２−７で算出した仮の液体保存量αρが異なる場合は、ステップＳ３に戻り、液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙを再度仮定して、繰り返し計算を行う。１又は複数回ステップＳ３〜Ｓ７を繰り返した後に、ステップＳ７で算出されるαＬ＿ｔｅｍｐｏｒａｒｙρＬ＿ｔｅｍｐｏｒａｒｙとステップＳ２−７で算出した液体保存量αρが等しくなったときに、この繰り返し計算は終了し、このときの液体圧力ＰＬが収束値となる。
【０１５５】
ここで、図１２を用いて、ステップＳ５で、上述のステップＳ２で算出した７つの保存量と、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと、ステップＳ４で算出した液体温度ＴＬとを用いて気泡半径の繰り返し計算を行い、液体の体積率を計算する方法について更に詳しく説明する。ここでは、気泡の成長・収縮計算から液体ボイド率を求める方法を説明する。
【０１５６】
すなわち、ステップＳ５−１では、式（４１）〜（４２）に示すように、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと、ステップＳ４で算出した液体温度ＴＬとを用いて過飽和度、新規発生する新規気泡半径、新規気泡数密度ΠＧ、新規気泡質量ΛＧを計算する。ここで、算出された新規発生した気泡は一番小さな気泡質量階級に属するので、式（４３）〜（４４）に示すように最小階級の気泡数密度ｎＧ，ｋ、気泡質量ｍＧ，ｋの値に加算される。尚、新規気泡数密度は、式（３７）に示すように、過飽和度である（ＰＧ−ＰＬ）を用いて算出する。ここで、式（７７）の関係、及び、新規発生気泡においてはＴＧ＝ＴＬであることを仮定していることから、ＰＧ−ＰＬ＝ＰＧ（ＴＬ）−ＰＬとすることができ、過飽和度（液体過熱度）を圧力表記により間接的に表記したものである。また、新規気泡質量は、新規気泡数密度及び新規気泡半径を用いて算出する。
【０１５７】
ステップＳ５−２では、各階級の既存の気泡半径ＲＫを繰り返し計算の初期値として仮定する。
【０１５８】
ステップＳ５−３では、式（７５）〜（７６）に示すように、液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙと液体保存量を用いて、気泡界面の運動方程式から各階級の気泡圧力ＰＧ，ｋを計算する。気泡圧力は、式（７４）で表せるレイリープレセット式（Ｒａｙｌｅｉｇｈ  Ｐｌｅｓｓｅｔ  Ｅｑｕａｔｉｏｎ）が適用されることが知られている。一般に、発生直後及び崩壊直前の小さな気泡（１μｍ以下）には、界面半径方向の慣性力の効果が大きいため右辺第４〜６項の時間微分項の効果が大きいが、比較的大きな気泡（１μｍ以上）では界面の半径方向の慣性力の効果が相対的に小さくなるため右辺第４〜６項の時間微分項の効果が小さくなる。そこで気泡径に応じて右辺第４〜６項の時間微分項に重み付けをした式（７５）を用いる。新規発生時には気泡と液体のスリップが無視できるのでスリップ項を無視した式（７６）を用いる。
【０１５９】
ステップＳ５−４では、式（７７）、（７８）に示すように、ステップＳ５−３で算出した気泡圧力ＰＧ，ｋを用いて、飽和条件から各階級の気泡温度ＴＧ，ｋ、気泡密度ρＧ，ｋを計算する。
【０１６０】
ステップＳ５−５では、ステップＳ３で仮定した液体圧力ＰＬ＿ｔｅｍｐｏｒａｒｙ、ステップＳ４で算出した液体温度ＴＬ、ステップＳ５−１で算出した気泡数密度ｎＧ，ｋ、気泡質量ｍＧ，ｋ、ステップＳ５−３で算出した気泡圧力ＰＧ，ｋ、ステップＳ５−４で算出した気泡温度ＴＧ，ｋ、気泡密度ρＧ，ｋを用いて、式（５７）〜式（７０）に示すように、各階級内の気泡数密度分布を仮定する。すなわち、ここでは、階級の境界値と、流れ場および階級相互の計算から求められる気泡数密度ｎＧ，ｋと体積当たりの気泡質量ｍＧ，ｋを満たすように、階級内で１次関数を用いた気泡数密度の分布を算出する。
【０１６１】
気泡数密度の分布を算出するには、まず、式（５９）に示すように、階級ｋ内の気泡１つあたりの平均気泡質量ｍ１，ｋを算出する。次に、この平均気泡質量ｍ１，ｋと、階級ｋに隣接する質量の軽い側の階級との境界での質量値ｗｌｉｇｈｔ，ｋと、階級ｋに隣接する質量の重い側の階級との境界での質量値ｗｈｅａｖｙ，ｋとに基づいて、式（６５）〜（６７）に示すように、最小階級用の分布１型１１、中間階級用の分布２型１２、最大階級用の分布３型１３に分類する。分類された各階級の境界値ＮＧＷ及び端点Ｗｆ，Ｗｂは、式（６８）〜（７０）により算出される。各階級内の気泡数密度分布は、例えば図９に示すように、この境界値及び端点から１次関数分布が得られる。各々の分布の１次式より階級の境界での気泡数密度がより厳密に計算できる。
【０１６２】
階級内で１次関数を用いた気泡数密度の分布を仮定できることは、数値計算上大きな効果を有する。例えば、階級内の気泡が有る程度成長しないとそれより大きな階級に気泡が全く移動しないので、最大階級の分布３型においては、重い側の境界において数密度の質量微分ＮＧｗ（ｗｈｅａｖｙ，ｋ）＝０である。このように不要な数値拡散を抑制することができ、計算精度を向上させることができる。
【０１６３】
ステップＳ５−６では、階級境界での蒸発量・凝縮量を計算する。階級内での蒸発・凝縮量ΓＩは、式（６２）の中辺２第２項と右辺第２項の関係で示されるように気泡１つあたりの蒸発・凝縮量γＩと気泡数密度関数ＮＧｗの積のｗに関する積分により算出される。
【０１６４】
ステップＳ５−７では、この各階級の気泡１個あたりの蒸発量γＩと階級内の蒸発・凝縮量ΓＩから、式（６１）〜（６２）により、階級ごとの気泡数密度ｎＧ，ｋと気泡質量ｍＧ，ｋが変更（修正）される。気泡数密度ｎＧ，ｋと気泡質量ｍＧ，ｋから、式（７９）に示すように、気泡１個あたりの体積を算出し、気泡半径Ｒ’Ｋを求める。
【０１６５】
これらの計算により気泡質量、および、液体密度、各方向液体運動量、全液体エネルギーに変化が生じる。さらに気泡質量が０となった場合には気泡が消滅する。これにより、再度、気泡数密度、気泡質量に変化が生じる。
【０１６６】
ステップＳ５−８では、ステップＳ５−７で算出された気泡半径Ｒ’ＫとステップＳ５−２で仮定した気泡半径Ｒｋを比較する。図１２に示すように、ステップＳ５−７で算出された気泡半径Ｒ’ＫとステップＳ５−２で仮定した気泡半径Ｒｋとが等しくならない場合は、ステップＳ５−２に戻り、気泡半径を再度仮定して、繰り返し計算を行う。
【０１６７】
１又は複数回ステップＳ５−２〜ステップＳ５−７を繰り返した後に、ステップＳ５−７で算出された気泡半径Ｒ’ＫとステップＳ５−２で仮定した気泡半径Ｒｋが等しくなったときに、この繰り返し計算は終了し、この気泡半径Ｒｋが収束値となる。
【０１６８】
ステップＳ５−９では、式（８０）に示すように、この気泡半径Ｒｋから各階級のボイド率αＧ，ｋを算出する。
【０１６９】
以上のように、本発明を適用した解析装置１は、仮に決定された液相部分の保存量を満たすように、新規発生・消滅気泡数、気泡の蒸発・凝縮量等を考慮した繰り返し計算により、各計算点での圧力、温度、液体各方向運動量、各階級ごとのボイド率、気泡数密度、気泡質量、及び、発生からの経過時間を決定することができる。
【０１７０】
ステップＳ８では、図１０に示すように、Ｓ２〜Ｓ７において算出された各計算値の境界値をＳ１−２で設定した表８の境界条件により補正、すなわち、再設定する。設定される境界条件は、流入境界、流出境界及び壁面境界の条件である。
【０１７１】
流入境界では、温度は一定とし、液単相であるサブクール流入とした。サブクール流入とは、液体温度が飽和温度（沸点）より低い状態で流入することを表すことを意味する。キャビテーションは圧力が下がることにより飽和温度が下がり、飽和温度が液体温度（一定）より低くなることで沸騰し気泡が発生する現象である。そのため、サブクールの状態では、一般的に気泡が発生せず、気泡の存在しない液単相の状態になる。流入境界の速度は、勾配０で外挿し、１点内側の計算格子点の速度と等しくする。ここで外挿とは、内挿の反対のことを意味する。内挿とはある点を挟む２点の値から、その内側にある任意の点を類推することであり、すなわち、外挿とはある点に隣接する２点から、その外側にある任意の点を類推することである。ここで算出する境界点は、計算点の外側に存在するために外挿によりこの値を算出する。全圧Ｐｔが一定となるように、流入境界の圧力Ｐは、Ｐ＝Ｐｔ−０．５Ｖ２と設定される。
【０１７２】
流出境界では、温度、速度、気泡諸量ともに勾配０で外挿し、１点内側の計算格子点の温度、速度、気泡諸量と等しく設定される。流出境界上の速度がその点での音速よりも小さい亜音速流出の場合は、流出境界上の出口圧力を一定にした。一方、流出境界上の速度がその点での音速よりも大きい超音速流出の場合は、流出境界上の出口圧力は、勾配０で外挿し１点内側の計算格子点の圧力と等しく設定される。
【０１７３】
壁面境界では、壁面に流体及び気泡が流入しないことから、圧力、気泡諸量は勾配０で外挿し１点内側の計算格子点の圧力、気泡諸量と等しくし、速度は壁面上で静止されるようにした。壁面から熱が流入しない断熱壁の場合には、壁面上の温度勾配が０とされ外挿し、１点内側の計算格子点の温度と等しくされ、壁面温度が一定の等温壁の場合は、壁面上の温度を壁温度で一定とされる。
【０１７４】
このステップＳ８では、壁面近傍における数値等が現実的な値からかけ離れた値となることを防止し、この境界条件により調整することで、よりキャビテーション現象を正確に再現することができる。ここで一定値を取る境界条件は「ディリクレ条件」を意味するものである。また、勾配一定値を取る条件は「ノイマン条件」を意味するものである。
【０１７５】
ステップＳ９では、図１０に示すように、あらかじめ設定された出力条件を満たさなかった場合、若しくは、ユーザ等が手動により出力しないことを決定した場合に、ステップ３に戻り再計算を行うことができる。ここで、あらかじめ設定される出力条件には、例えば、一定時間おきに出力をすることがある。この場合、指定された一定タイムステップに達していないときは、ステップ３に戻り再計算を行う。このとき、計算結果を旧時間のデータとしてメモリに保存し、次の計算を行うことができる。また、あらかじめ設定された出力条件を満たした場合、若しくは、ユーザー等が手動により出力を決定した場合に、計算結果を出力することができる。このとき、指定された割合で画像や保存データを出力することができる。計算結果は、例えば、図１４に示すように、圧力、マッハ数、温度、ボイド率、気泡半径（階級１）及び気泡半径（階級２）を、数値の大きさを色の分布で示すことができるとともに、図１５に示すように、壁面上圧力分布を線図上に得ることができる。図１４において、面塗りは液体圧力，マッハ数，液体温度，ボイド率，階級１と階級２の気泡半径を示す．図１５において、線図は壁面上圧力分布の計算結果を示し、白丸はＳｉｍｏｎｅａｕ  ａｎｄ  Ｈｅｎｄｒｉｃｋｓが行った壁面上圧力分布の実験結果である。図１４、図１５の結果により、本解析装置１において、キャビテーション現象が再現されていることが分かる。
【０１７６】
ステップＳ１０では、あらかじめ設定された終了条件を満たさなかった場合、若しくは、ユーザー等が手動により終了しないことを決定した場合に、ステップ３に戻り再計算を行うことができる。また、あらかじめ設定された終了条件を満たした場合、若しくは、ユーザー等が手動により終了する決定した場合に、解析装置１による解析を終了する。尚、ここで終了条件としては、非定常計算における、指定した時間または、指定したタイムステップまで計算することを設定することができ、さらに、定常計算における、定常状態になるまで計算することを設定することができる。この終了条件を設定自在としたことにより、様々な用途の解析に用いることができる。
【０１７７】
本発明を適用した解析装置１は、気泡密度分布を考慮することにより、数値拡散を防ぐことができたので、従来の解析装置の計算プログラムにおいて気泡の隣接階級の間隔比を１のみでしか行われていなかった計算を、気泡の隣接階級の間隔比を２．１５（１０１／３）程度まで高めても計算を行えるようになった。このことにより、本発明は、例えば、１ｎｇ〜１００ｎｇの気泡を同時に扱う場合、従来では、１００階級以上の階級を必要であったのに対し、本解析装置では、６階級で計算することができるようになった。
【０１７８】
また、本発明を適用した解析装置１は、気泡の成長・収縮計算と液体の圧力・温度計算をカップリングして繰り返し計算を行うことにより、時間刻みが極端に小さくなる現象を回避できたので、従来の解析装置における計算プログラムにおいてキャビテーション気泡無し流れで１０−８（ｓ）程度の時間刻みで計算可能であった流れ場において、気泡が成長してくると１０−１２（ｓ）以下の時間刻みが必要であるのに対し、本解析装置では、キャビテーション気泡の有無に拘わらず１０−８（ｓ）程度の時間刻みで安定的に計算できるようになった。
【０１７９】
上述のように、本発明は、各階級内の気泡数密度分布を考慮することにより、数値拡散を防ぐことができ、気泡の蒸発・凝縮計算の精度を高めることができるとともに、階級間隔を不等間隔にした場合の計算安定性を高めることができる。よって、本発明は、精度及び安定性の高いキャビテーション解析を可能とするものであり、汎用性を高めることを可能とする。
【０１８０】
また、本発明によれば、気泡の成長・収縮計算と流体の圧力・温度計算とを組合せ、繰り返し計算を行うことにより、時間刻みを大きく取ることができ、キャビテーション解析の計算安定性を増すことができ、よって、計算負荷を低減することが可能となる。
【０１８１】
尚、上述した解析装置は、２次元計算を行うものであるが、本発明は、３次元計算を行うこともできる。
【産業上の利用可能性】
【０１８２】
本発明を適用したキャビテーション解析方法、キャビテーション解析プログラム、プログラム記録媒体及び解析装置は、ロケットエンジン等に用いられる液体酸素、液体水素用のポンプ及びＬＰＧ等の液体メタンに用いられるポンプのキャビテーション解析に用いることができるのみならず、液体窒素、液体酸素、液体水素及び液体メタン等の各種低温流体中に発生するキャビテーション解析に用いることができる。
【図面の簡単な説明】
【０１８３】
【図１】本発明を適用した解析装置の構成を示す図である。
【図２】本発明を適用した解析装置により得られるボイド率と音速の関係を示す図である。
【図３】壁面近傍に発生する新規気泡の壁面の面積と計算空間との関係を示す図である。
【図４】図３に示すｊ＝２のセルにおけるボイド率を説明する図である。
【図５】図３に示すｊ＝１のセルにおけるボイド率を説明する図である。
【図６】壁面近傍に発生する新規気泡のボイド率の与え方を説明するための図である。
【図７】気泡数密度の確率密度関数の気泡１つあたりの質量に対する分布を示す図である。
【図８】従来技術で用いられる各階級の気泡の気泡数密度分布を示す図である。
【図９】本発明を適用した解析装置により算出される各階級の気泡の数密度分布を示す図である。
【図１０】本発明を適用した解析装置のフローチャートである。
【図１１】本発明を適用した解析装置の液体圧力計算のフローチャートである。
【図１２】本発明を適用した解析装置の気泡半径計算のフローチャートである。
【図１３】本発明を適用した解析装置により解析されるノズル形状を示す図である。
【図１４】本発明を適用した解析装置により得られる計算結果の一例を示す図である。
【図１５】本発明を適用した解析装置により得られる壁面上圧力分布の計算結果の一例を示す図である。
【符号の説明】
【０１８４】
１  解析装置、２  入力部、３  演算部、４  メモリ、５  表示部、６  出力部、７  記憶部
【図１】
2005222339.tif 000064
【図２】
2005222339.tif 000065
【図３】
2005222339.tif 000066
【図４】
2005222339.tif 000067
【図５】
2005222339.tif 000068
【図６】
2005222339.tif 000069
【図７】
2005222339.tif 000070
【図８】
2005222339.tif 000071
【図９】
2005222339.tif 000072
【図１０】
2005222339.tif 000073
【図１１】
2005222339.tif 000074
【図１２】
2005222339.tif 000075
【図１３】
2005222339.tif 000076
【図１４】
2005222339.tif 000077
【図１５】
2005222339.tif 000078