(57)【要約】 【課題】  気体と液体とが混在する気泡流の挙動を精度よく解析することが可能な気泡流の解析方法、気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置を提供する。【解決手段】  本発明は、気泡Ｇと液体Ｌとが混在した気泡流の渦度場に着目して複数の渦要素１６，１７で離散化し、気泡流内の気泡Ｇと液体Ｌとの気液体積割合（気泡体積率αｇ、液体体積率αｌ）を考慮した運動量保存関係及び質量保存関係から渦度輸送方程式を導出し、この渦度輸送方程式のラグランジェ計算により離散渦要素の挙動を求めるものである。【選択図】  図５
2006085603.tif 000002
【特許請求の範囲】
【請求項１】
気泡と液体とが混在した気泡流の解析方法であって、
解析領域内に存在する気泡及び離散渦要素の位置を時間的に変化させつつ、次の処理（ａ）〜（ｅ）を行って前記解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する気泡流の解析方法。
（ａ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記解析領域内における気泡の分布に基づき求める処理。
（ｂ）前記各解析要素内における循環強度変化を、前記気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求める処理。
（ｃ）前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、その解析要素内における前記循環強度変化と、前記各離散渦要素の位置分布とに基づきラグランジェ的に求める処理。
（ｄ）前記解析領域内の前記評価点における渦度を、前記各離散渦要素の循環強度及び位置に基づき求める処理。
（ｅ）前記評価点における渦度から前記液体のベクトルポテンシャルを求め、前記解析領域内の前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める処理。
【請求項２】
前記気泡流の過度輸送方程式は、次の数式［数１］であることを特徴とする請求項１に記載の気泡流の解析方法。
【数１】
2006085603.tif 000024
Γ  ：解析要素内における循環強度変化
αｌ：気泡と液体との気液体積割合
ｇ  ：重力加速度
ｕｌ：液体の速度
Ａ  ：解析要素の面積
【請求項３】
前記液体のスカラーポテンシャル式は、次の数式［数２］であることを特徴とする請求項１又は請求項２に記載の気泡流の解析方法。
【数２】
2006085603.tif 000025
αｌ：気泡と液体との気液体積割合
φ  ：液体のスカラーポテンシャル
ψ  ：液体のベクトルポテンシャル
【請求項４】
前記処理（ｄ）において、前記各離散渦要素のコア半径を、前記体積割合の気泡流についてコアスプレッディング法を適用して算出し、前記渦度を求めることを特徴とする請求項１〜請求項３のいずれかに記載の気泡流の解析方法。
【請求項５】
前記体積割合の気泡流についてコアスプレッディング法を適用した算出式は、次の数式［数３］であることを特徴とする請求項４に記載の気泡流の解析方法。
【数３】
2006085603.tif 000026
σα：離散渦要素のコア半径
αｌ：気泡と液体との気液体積割合
νｌ：液体の動粘度
【請求項６】
前記解析領域内に存在する気泡の形状、大きさ及び個数のうち少なくともいずれか１つを変更可能とし、前記処理（ａ）では、この変更された気泡の形状、大きさ及び個数に応じた気泡の分布に基づき前記気液体積割合を求めることを特徴とする請求項１〜請求項５のいずれかに記載の気泡流の解析方法。
【請求項７】
気泡と液体とが混在した気泡流の解析するための気泡流シミュレーションプログラムであって、コンピュータに、次のステップ（ａ）〜（ｇ）を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する気泡流シミュレーションプログラム。
（ａ）前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出するステップ。
（ｂ）前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出するステップ。
（ｃ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記ステップ（ｂ）で算出された各気泡の位置に基づき求めて記憶手段に保存するステップ。
（ｄ）前記ステップ（ｃ）によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存するステップ。
（ｅ）前記ステップ（ｄ）によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記ステップ（ａ）で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存するステップ。
（ｆ）前記ステップ（ｅ）によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求めるステップ。
（ｇ）前記ステップ（ｆ）によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求めるステップ。
【請求項８】
前記請求項７に記載の気泡流シミュレーションプログラムを記憶した記憶媒体。
【請求項９】
気泡と液体とが混在した気泡流の解析するための気泡流シミュレーション装置あって、次の手段（ａ）〜（ｇ）を備え、これらの手段を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する気泡流シミュレーション装置。
（ａ）前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出する手段。
（ｂ）前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出する手段。
（ｃ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記手段（ｂ）で算出された各気泡の位置に基づき求めて記憶手段に保存する手段。
（ｄ）前記手段（ｃ）によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存する手段。
（ｅ）前記手段（ｄ）によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記手段（ａ）で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存する手段。
（ｆ）前記手段（ｅ）によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求める手段。
（ｇ）前記手段（ｆ）によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める手段。
【請求項１０】
前記手段（ａ）で算出される前記離散渦要素の位置データが保存される渦要素情報記憶手段と、
前記手段（ｂ）で算出される前記気泡の位置データ及び前記手段（ｃ）で算出される気液体積割合データが保存される気泡情報記憶手段とを備え、
前記手段（ａ）による処理と、前記手段（ｂ）及び前記手段（ｃ）による処理とを並列的に実行することを特徴とする請求項９に記載の気泡流シミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、液中に分散した気泡が液体と相互作用を及ぼし合いながら流れる気泡流（二相流）の挙動の解析方法に関する。
【背景技術】
【０００２】
微小な気泡と液体とが混在して相互作用を及ぼし合いながら流れる気泡流は、例えば反応装置、熱交換器及び廃水処理装置など、様々な工業装置において観察される。このため、この気泡流の挙動について、合理的な解析方法の確立が望まれている。
【０００３】
従来の解析方法では、流れの変数として速度、圧力及び気泡体積率が用いられ、それらを支配する複数の偏微分方程式を様々な数値解法（差分法、有限体積法、有限要素法）によって計算するようにしていた。これらの偏微分方程式は、微分を差分で近似表現することにより、代数方程式に置き換えられるが、非線形項が数値不安定性を有するため、気泡流のレイノルズ数が高い場合には計算過程で発散が生じるおそれがある。特に、工業装置で見受けられる気泡流は、レイノルズ数が高いものが多く、精度の高い気泡流解析ができないという問題があった。
【０００４】
ところで、従来から、液体又は気体のみの単相流の解析については渦法が有効に利用されている。ここで、渦法とは、渦度をもつ渦要素を追跡して渦度場の時間変化を求めるLagrange 型解法であり，大規模渦の形成や変形など渦構造の発展過程を良好に計算できるものである。そして、下記非特許文献１には、この渦法を利用して、気泡流の挙動を解析する技術が開示されている。
【非特許文献１】Sene, K. J., et al., The role of coherent structures in bubble transport by turbulent shear flows,J. Fluid Mech., 259(1994), 219-240.
【特許文献１】特開２００４−１２６９２５公報
【特許文献２】特開２００４−１１８３９４公報
【特許文献３】特開２００１−９９７４８公報
【特許文献４】特開平９−３２９１０９号公報
【特許文献５】特開２００３−３３１２０８公報
【特許文献６】特開２００４−５０３００５公報
【特許文献７】特開２００１−１３２７００公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記非特許文献１に記載の技術は、気泡が液相に影響を及ぼさないことを前提として気泡流の挙動を解析するものであり、液中に分散した気泡が液体と相互作用を及ぼし合いながら流れる気泡流の挙動を精度良く解析できるものではなかった。
【０００６】
なお、上記特許文献１〜７には、流体の挙動をその支配方程式から計算によって予測しようとする技術が開示されているが、いずれも、本願が解析対象とする気泡流を精度よく解析できるものではない。
【０００７】
即ち、特許文献１の技術は、渦法を用いて流体の流れを解析するものであるが、水または空気などの単相流を解析対象とするものである。
特許文献２の技術は、気体と液体とが接する面、即ち気液自由界面の変形に関するものであり、気体と液体とが分離した流れを解析対象としている。
特許文献３の技術は、やはり単相流の挙動を有限差分法によって解析するものである。
【０００８】
特許文献４の技術は、流水中の固体壁面から噴き出された気泡について、その運動と分布状態とを計算し、気泡による水流の乱れ変化を求め、壁面における流体摩擦力の低減効果を解析するものである。壁面近傍での流れ現象のみを解析対象としており、水流については乱れ特性など限定された現象しか解析できず、広範囲における圧力や速度などの分布が考慮されていない。
特許文献５の技術は、飛沫や気泡の挙動を解析するものであり、その計算に有限差分法に準じた解法を用い、飛沫や気泡の存在の判定には、VOF法を利用している。VOF法は気体と飛沫、気体と気泡の界面の分解能が低い欠点を持つことが知られている。
【０００９】
特許文献６の技術は、気体と液体とが分離した状態で存在し、両者の自由界面がある流体の挙動を解析するものであり、流速が低く層流状態にある流れが解析対象となっている。
特許文献７の技術は、タンク内の水流を計算し、飽和水蒸気圧と計算で得られた圧力とを比較することにより、気泡の発生位置を予測するものであるが、気泡と液体とが相互作用を及ぼしあうことについて考慮されていない。
【００１０】
本発明は上記のような事情に基づいて完成されたものであって、気体と液体とが混在する気泡流の挙動を精度よく解析することが可能な気泡流の解析方法、気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記の目的を達成するための手段として、請求項１の発明に係る気泡流解析方法は、気泡と液体とが混在した気泡流の解析方法であって、解析領域内に存在する気泡及び離散渦要素の位置を時間的に変化させつつ、次の処理（ａ）〜（ｅ）を行って前記解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
（ａ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記解析領域内における気泡の分布に基づき求める処理。
（ｂ）前記各解析要素内における循環強度変化を、前記気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求める処理。
（ｃ）前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、その解析要素内における前記循環強度変化と、前記各離散渦要素の位置分布とに基づきラグランジェ的に求める処理。
（ｄ）前記解析領域内の前記評価点における渦度を、前記各離散渦要素の循環強度及び位置に基づき求める処理。
（ｅ）前記評価点における渦度から前記液体のベクトルポテンシャルを求め、前記解析領域内の前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める処理。
【００１２】
請求項２の発明は、請求項１に記載の気泡流の解析方法において、前記気泡流の過度輸送方程式は、次の数式［数１］であることを特徴とする。
【数１】
2006085603.tif 000003
Γ  ：解析要素内における循環強度変化
αｌ：気泡と液体との気液体積割合
ｇ  ：重力加速度
ｕｌ：液体の速度
Ａ  ：解析要素の面積
【００１３】
請求項３の発明は、請求項１又は請求項２に記載の気泡流の解析方法において、前記液体のスカラーポテンシャル式は、次の数式［数２］であることを特徴とする。
【数２】
2006085603.tif 000004
αｌ：気泡と液体との気液体積割合
φ  ：液体のスカラーポテンシャル
ψ  ：液体のベクトルポテンシャル
【００１４】
請求項４の発明は、請求項１〜請求項３のいずれかに記載の気泡流の解析方法において、前記処理（ｄ）において、前記各離散渦要素のコア半径を、前記体積割合の気泡流についてコアスプレッディング法を適用して算出し、前記渦度を求めることを特徴とする。
請求項５の発明は、請求項４に記載の気泡流の解析方法において、前記体積割合の気泡流についてコアスプレッディング法を適用した算出式は、次の数式［数３］であることを特徴とする。
【数３】
2006085603.tif 000005
σα：離散渦要素のコア半径
αｌ：気泡と液体との気液体積割合
νｌ：液体の動粘度
【００１５】
請求項６の発明は、請求項１〜請求項５のいずれかに記載の気泡流の解析方法において、前記解析領域内に存在する気泡の形状、大きさ及び個数のうち少なくともいずれか１つを変更可能とし、前記処理（ａ）では、この変更された気泡の形状、大きさ及び個数に応じた気泡の分布に基づき前記気液体積割合を求めることを特徴とする。
なお、上記請求項２，３，６の構成はいずれも、次述する請求項７〜１０の発明に対して適用することができる。また、請求項４，５の構成はいずれも、次述する請求項７，８のステップ（ｆ）、請求項９，１０の手段（ｆ）にそれぞれ適用することができる。
【００１６】
請求項７の発明に係る気泡流シミュレーションプログラムは、気泡と液体とが混在した気泡流の解析するための気泡流シミュレーションプログラムであって、コンピュータに、次のステップ（ａ）〜（ｇ）を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
（ａ）前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出するステップ。
（ｂ）前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出するステップ。
（ｃ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記ステップ（ｂ）で算出された各気泡の位置に基づき求めて記憶手段に保存するステップ。
（ｄ）前記ステップ（ｃ）によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存するステップ。
（ｅ）前記ステップ（ｄ）によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記ステップ（ａ）で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存するステップ。
（ｆ）前記ステップ（ｅ）によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求めるステップ。
（ｇ）前記ステップ（ｆ）によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求めるステップ。
【００１７】
請求項８の発明は、上記請求項７に記載の気泡流シミュレーションプログラムを記憶した記憶媒体である。なお、記憶媒体としては、フロッピーディスク（登録商標）、ＣＤ−ＲＯＭ等が好適である。
【００１８】
請求項９の発明に係る気泡流シミュレーション装置は、気泡と液体とが混在した気泡流の解析するための気泡流シミュレーション装置あって、次の手段（ａ）〜（ｇ）を備え、これらの手段を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
（ａ）前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出する手段。
（ｂ）前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出する手段。
（ｃ）前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記手段（ｂ）で算出された各気泡の位置に基づき求めて記憶手段に保存する手段。
（ｄ）前記手段（ｃ）によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存する手段。
（ｅ）前記手段（ｄ）によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記手段（ａ）で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存する手段。
（ｆ）前記手段（ｅ）によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求める手段。
（ｇ）前記手段（ｆ）によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める手段。
【００１９】
請求項１０の発明は、請求項９に記載の気泡流シミュレーション装置において、前記手段（ａ）で算出される前記離散渦要素の位置データが保存される渦要素情報記憶手段と、前記手段（ｂ）で算出される前記気泡の位置データ及び前記手段（ｃ）で算出される気液体積割合データが保存される気泡情報記憶手段とを備え、前記手段（ａ）による処理と、前記手段（ｂ）及び前記手段（ｃ）による処理とを並列的に実行することを特徴とする。
なお、より具体的な構成として、次のようなものであってもよい。
「渦要素情報保存手段と、気泡情報保存手段と、ポテンシャルデータ保存手段とを備え、
前記渦要素情報保存手段には、前記手段（ａ）で算出される前記各離散要素の位置データと、前記手段（ｄ）で算出される前記解析要素内における循環強度変化データと、前記手段（ｅ）で算出される前記各離散渦要素の循環強度及び位置データと、前記手段（ｆ）で算出される前記評価点における渦度データと、が保存され、
前記気泡情報保存手段には、前記手段（ｂ）で算出される前記各気泡の位置データと、前記手段（ｃ）で算出される前記気液体積割合データと、が保存され、
前記ポテンシャルデータ保存手段には、前記手段（ｇ）で算出される前記液体のベクトルポテンシャルデータと、前記液体のスカラーポテンシャルデータと、が保存されるよう構成される構成。」
【発明の効果】
【００２０】
本発明は、気泡と液体とが混在した気泡流の渦度場に着目して複数の渦要素で離散化し、気泡流内の気泡と液体との気液体積割合を考慮した運動量保存関係及び質量保存関係から渦度輸送方程式を導出し、この渦度輸送方程式のラグランジェ計算により離散渦要素の挙動を求めるものである。
渦度輸送方程式のラグランジェ計算には非線形項が現れないため、レイノルズ数に依存することなく、高レイノルズ数の気泡流について高精度解析が可能になる。また、渦度の高い領域には渦要素が密集する傾向があるが、このような領域に対しても気泡流について高い解析精度を確保できる。
【００２１】
また、渦度を求めるにあたって、各離散渦要素のコア半径を、体積割合の気泡流についてコアスプレッディング法を適用して算出して適用することが望ましい。このようにすれば、粘性拡散による減衰をも考慮した正確な渦度を求めることができる（請求項４，５）。
【００２２】
また、請求項６の発明によれば、解析領域内に存在する気泡の形状、大きさ及び個数のうち少なくともいずれか１つを変更可能とし、処理（ａ）では、この変更された気泡の形状、大きさ及び個数に応じた気泡の分布に基づき気液体積割合を求めるものである。従って、気泡の形状、大きさや個数を変更することにより、気泡流内に存在する気泡の変形（径の変動を含む）、分裂、合体などの各現象に対応して気泡流の解析を行うことができる。
【００２３】
更に、請求項１０の発明によれば、手段（ａ）で算出される離散渦要素の位置データが保存される渦要素情報記憶手段と、手段（ｂ）で算出される気泡の位置データ及び手段（ｃ）で算出される気液体積割合データが保存される気泡情報記憶手段とが備えられ、手段（ａ）による処理と、手段（ｂ）及び手段（ｃ）による処理とを並列的に実行する構成とした。これにより、気泡流の解析処理の高速化を図ることができる。
【００２４】
＜理論＞
次の本発明の理論について説明する。
１．記号
本発明において下記に登場する各符号は次のように定義されている。
2006085603.tif 000006
【００２５】
２．仮定
本発明は以下の仮定を前提としている。
（ａ）流れは、流れ場の代表長さに比べて微細な気泡が分散した気泡流を対象とする。
（ｂ）液体及び気泡は、非圧縮性であり、また、熱的非平衡はない。
（ｃ）気泡の質量及び運動量は、気体に比べて十分小さく無視できる。
（ｄ）気泡は球形を保って運動し、合体及び分裂はしない。
【００２６】
３．基礎式
本発明では、所定の設定データ（後述する入力データ）に基づいて，計算に必要な境界条件を設定し、流れの基礎方程式に基づいて計算を実行する。以下、境界条件の設定と計算に用いる流れの基礎方程式について説明する。
【００２７】
（１）気泡流の保存方程式
液体中に微細な気泡が分散して流れる気泡流について、上記仮定（ａ）〜（ｃ）を前提とした場合、この気泡流の質量および運動量の保存方程式は、それぞれ次式で表される。
【数１０】
2006085603.tif 000007
【数１１】
2006085603.tif 000008
但し、液体と気泡の体積率は次の関係を満たす。
【数１２】
2006085603.tif 000009
【００２８】
（２）気泡の運動方程式
気泡には仮想質量力、流動抗力、重力および揚力が作用するものとする。気泡の運動方程式は、上記仮定（ｄ）を前提とした場合、次式で与えられる。
【数１３】
2006085603.tif 000010
【００２９】
（３）液体速度場の直交分解
ヘルムホルム（Helmholts）
の定理によれば、任意のベクトル場はスカラポテンシャルφ の勾配とベクトルポテンシャルψ の回転との和として表されることから、液相速度ｕｌ
は次式のように記述できる。 
【数１４】
2006085603.tif 000011
ここで、ψはスカラ関数の勾配を加えても結果が不変である。この任意性を取り除き、ψを一意的に定めるため、ψはソレノイダルであるものとする。すなわち、
【数１５】
2006085603.tif 000012
【００３０】
数式１４ を数式１０に代入し、恒等式∇ ・ (∇×ψ)
＝０を利用して変形すれば次式が得られる。
【数１６】
2006085603.tif 000013
そして、数式１４の回転をとった後に、数式１５を代入すれば、次のベクトルPoisson
方程式が得られる。 
【数１７】
2006085603.tif 000014
ここで、ωは液体の渦度である。
【００３１】
４．数値解法
（１）渦要素と気泡とのラグランジェ（Lagrange）解析
上記数式１１の回転をとれば気泡流の渦度方程式が得られ、二次元解析の場合には次式で表される。
【数１８】
2006085603.tif 000015
【００３２】
次に、渦度場を多数の渦要素により離散化する。単相流解析に対する渦要素モデルを適用し、コア構造をもつ離散渦要素を導入する。離散渦要素α の循環強度をΓα、コア半径をσα、位置ベクトルをｘα
とすれば、離散渦要素α による位置ｘ（評価点の位置）における渦度は次式で与えられる。 
【数１９】
2006085603.tif 000016
また、流れ場に存在する離散渦要素αの移流は次のラグランジェ（Lagrange）解析から計算できる。
【数２０】
2006085603.tif 000017
【００３３】
ｘ−ｙ平面内の解析領域Ｈを、図１に示すような、寸法がΔｘとΔｙの四角形の計算格子１０（解析要素）に分割し、上記数式１６と数式１７を差分法で解析し、得られたスカラーポテンシャルφとベクトルポテンシャルψを数式１４に代入して格子点１１（評価点）上の液体の速度ｕｌを求める。この際、格子点１１で与える渦度ω（同図ω１〜ω４）の値は、各離散渦要素による渦度を数式１９から求め、全離散渦要素について加え合わせることにより計算する。各計算格子１０における液体体積率αｌは、数式１３のLagrange 解析から気泡運動を解析して気泡体積率αｇを求めたのち、数式１２を利用して算出する。
【００３４】
また、例えば、解析領域Ｈの入口境界において、ｘ−ｙ平面の垂直方向（図１で紙面奥行き方向）に気泡が間隔Δｚで一様に分布するものと仮定し、面積Ａ（＝ΔｘΔｙ）の計算格子１０におけるαｇを次式で計算する。
【数２１】
2006085603.tif 000018
ここで、ｎｇ は任意の時間における計算格子１０内に存在する直径ｄの気泡Ｇの個数（図１では４つを図示）である。
渦度は、数式１８から分かるように、粘性拡散（右辺第１項）と相分布の勾配（右辺第２項）により時間変化する。これらの変化は、数式１９における離散渦要素のコア半径σαと循環強度Γαの変化として模擬し、以下に示すように数式２０のLagrange 計算と連立させて解析する。
【００３５】
（２）粘性拡散によるコア半径の変化
渦度は粘性拡散により減衰する。ここで、液体のみ又は気体のみの単相流解析では、コア半径の大きさを時間的に変化させるコアスプレッディング（Core spreading）法により考慮される。そこで、本発明においても、これを準用し、次式をLagrange 解析してコア半径の時間変化を求める。
【数２２】
2006085603.tif 000019
なお、Core spreading 法は、Leonard, A., Vortex methods for flow simulation, J. Comput. Phys. 37, (1980), 289-335 に記載されている。
【００３６】
（３）気泡運動による循環強度の変化
任意の閉曲線まわりの循環強度Γ の時間変化率は、レイノルズ（Reynolds） の輸送定理を用いれば、次式で表される．
【数２３】
2006085603.tif 000020
そして、この数式２３の右辺に上記数式１８を代入すれば次式を得ることができる。
【数２４】
2006085603.tif 000021
ここで、粘性拡散項は、数式２２で考慮されているので、代入に際し無視している。
【００３７】
上記計算格子１０に数式２４を適用し、各計算格子１０内における循環強度Γ の時間変化率ΔΓ／Δｔを計算する。ただし、数式２４の被積分関数は格子点１１間で線形変化するものと仮定して数値積分する。各計算格子１０内にｎｖ個の離散渦要素が存在する場合、Δｔにおける離散渦要素１つ当たりの循環強度の変化量をΔΓ／ｎｖとする。離散渦要素が存在しない場合には、循環強度ΔΓの離散渦要素１つをその計算格子１０の中央から新しく発生させる。
【００３８】
（４）解析手順
時刻ｔ＝ｔでの流れ場における気泡流の挙動（各格子点１１における液体の速度に関するベクトルポテンシャル及びスカラーポテンシャル等）が既知であるならば、時刻ｔ＝ｔ＋Δｔの挙動を、例えば以下の手順で計算できる。
（ａ）気泡の運動を数式１３により計算し、気液体積割合（αｇ（＝１−αｌ））を数式２１から求める。
（ｂ）離散渦要素のコア半径σを数式２２から求める。
（ｃ）各計算格子１０内における循環強度の変化量ΔΓを数式２４から求める。
（ｄ）各離散渦要素の移流を数式２０により計算し、各格子点１１における渦度ωを数式１９から求める。
（ｅ）各格子点１１における液体の速度のベクトルポテンシャルψを数式１７から求める。
（ｆ）各格子点１１における液体の速度のスカラーポテンシャルφを数式１６から求める。
（ｇ）各格子点１１における液体の速度ｕｌを、数式１４から求める。
【００３９】
このような手順を予め定めた時間間隔毎に繰り返し実行することで、気泡流について高精度解析を行うことができる。
【発明を実施するための最良の形態】
【００４０】
以下、本発明の一実施形態を前述した理論及び図１〜図６によって説明する。
【００４１】
本実施形態に係る気泡流シミュレーションプログラムは、後述する情報記憶装置６に記録され、パーソナルコンピュータ（以下「コンピュータ１」という）に導入することにより、そのコンピュータ１に後述する気泡流の挙動の解析ステップを実行させるためのものである。
【００４２】
１．コンピュータとその周辺機器の構成
図２には、コンピュータ１とその周辺機器の構成とを示した。コンピュータ１には、演算及び周辺機器の制御等を行うための中央処理部２、この中央処理部２と周辺機器との接続を行うためのＩ／Ｏインターフェース３、及び周辺機器である表示装置４、印刷装置５、情報記憶装置６及び入力装置７が接続されている。
【００４３】
なお、表示装置４には、カソードレイチューブ（ＣＲＴ）の他に液晶装置、液晶プロジェクタ等が含まれる。また、印刷装置５は、コンピュータ１に直接に連結されている必要はなく、ＬＡＮ等のネットワークを介して連結されていてもよい。また、入力装置７は、例えばキーボードやマウスなどのコンソールであって、ここでの操作によって解析領域Ｈの境界を設定する座標位置データや流体特性を設定する入力データなどを入力することができる。
【００４４】
中央処理部２の詳細については図示しないが、中央演算処理装置（ＣＰＵ）の他に、ＲＯＭ、各種ＲＡＭ、チップセット等から構成されている。情報記憶装置６は、ハードディスクやフロッピーディスク、及びＭＯやＣＤ−ＲＯＭ等であり、コンピュータ１のプログラムや各種データを読み込みまたは記録することができる。
【００４５】
２．気泡流の解析
さて、本実施形態の作用を上記中央処理部２にて実行される気泡流シミュレーションプログラムのフローチャートを参照しつつ説明する。
【００４６】
（１）計算に必要な入力データの入力及び境界条件の設定
まず、上記入力装置７での操作によって、計算に必要な入力データを入力する。ここで、入力データとしては、液体Ｌの密度ρｌ及び動粘度νｌ、気泡Ｇの密度ρｇ及び直径ｄ、解析領域Ｈの入口における液体Ｌ及び気泡Ｇの初期速度ｕｌ(0)，ｕｇ(0)、同じく解析領域Ｈの入口における単位時間当たりの気泡Ｇの放出個数、解析領域Ｈの幾何形状及び寸法、解析時間、解析時間タイミングΔｔｖ、計算格子１０の個数及び寸法である。
【００４７】
（１）条件設定
本実施形態では、壁面で囲まれていない自由空間における、気泡Ｇを含む乱流（気泡流）に対する数値解析である。より具体的には、図３に示すような、異なる速度をもつ液体Ｌが混合して速度せん断層１５を形成し、そこに気泡Ｇが混在して流れる、気泡流を解析対象としているのである。
【００４８】
ここで、上記速度せん断層１５を形成する平板Ｗ１，箱体Ｗ２などの物体は、図４に示すように、これらの表面に循環強度Γをもつ固定渦要素１６（離散渦要素）を配置して表す。この固定渦要素１６の個数や位置は、上記入力装置７での操作によって行う。物体表面の任意の位置における液体の速度ｕｌは、固定渦要素１６と放出される放出渦要素１７（離散渦要素）の循環強度によって求められる。各固定渦要素１６の循環強度Γは、平板Ｗ１、箱体Ｗ２の表面上におけるｎ箇所で法線方向の速度成分が零となるように定められ、この計算は、コンピュータ１によって次の数式２５に示す連立方程式を解くことによって実行される。なお、この際、ケルビンの循環定理に従い、固定渦要素１６の循環強度の総和と、放出渦要素１７の循環強度の総和とが等しくなる条件が付与されている。
【数２５】
2006085603.tif 000022
ここで、固定渦要素の個数をｎとすれば、Ａはｎ×ｎの行列、Ｂはｎ列のベクトルである。
【００４９】
次いで、速度せん断層１５は、流れ場のはく離位置から放出渦要素１７を放出することにより近似的に表現される。放出渦要素１７の初期循環強度Γ０は、速度せん断層１５の両側の液体の速度Ｕ１，Ｕ２を用いた次の式で与えられる。
【数２６】
2006085603.tif 000023
ここで、Δｔｖは、放出渦要素の放出時間間隔であり、本実施形態では上記解析時間タイミングに一致させている。解析開始の操作がされると、コンピュータ１は、上記数式２６で算出される初期循環強度Γ０の放出渦要素１７を、流れ場のはく離位置即ち速度せん断層１５の始点から一定の時間間隔Δｔｖで順次放出するシミュレートを実行する。始点は、異なる速度をもつ流体が合流する位置或いは物体角部が相当し、本実施形態では、図３中、平板Ｗ１の下流側先端Ｐ１、箱体Ｗ２の両側面の下流端Ｐ２，Ｐ２である。
【００５０】
（２）解析処理
コンピュータ１は、上記条件設定後、図５のフローチャートに示す処理を上記解析時間タイミングΔｔｖで繰り返し実行し、順次算出される各格子点１１における液体Ｌの速度のベクトルポテンシャル及びスカラーポテンシャルの時間的変化に基づき気泡流の挙動を数値解析する。
【００５１】
ステップＳ１で解析時間タイミングになったかどうかを判断し、解析時間タイミングになったとき（ステップＳ１で「Ｙ」）、ステップＳ２で、上記入力装置７にて入力された気泡Ｇの放出個数分の気泡を解析領域Ｈの入口から放出する。即ち、本発明の上記「ステップ（ｂ）」及び「手段（ｂ）」でいう「解析領域内に複数の気泡を導入」に相当する処理を実行する。それとともに、速度せん断層から放出渦要素１７を放出するシミュレートを実行する。即ち、本発明の上記「ステップ（ａ）」及び「手段（ａ）」でいう「解析領域内に複数の離散渦要素を順次導入」に相当する処理を実行する。
【００５２】
次に、ステップＳ３で、気泡Ｇの運動を、前回の解析時間タイミングで算出された液体Ｌの速度ｕｌに基づき数式１３から算出する。即ち、本発明の上記「ステップ（ｂ）」及び「手段（ｂ）」でいう「解析領域内に存在する各離散渦要素の各位置を算出」に相当する処理を実行する。より具体的には、数式１３にオイラーの前進差分法を適用し、その右辺に現れる液体の速度を格子点１１における速度の補完から定め、ラグランジェ（Lagrange）計算により解析領域Ｈ内における各気泡Ｇの位置を追跡するのである。なお、補完には評価したい位置の近傍格子点１１を選択し、それらに関して内挿法が利用できる。
【００５３】
ステップＳ４では、気泡体積率αｇを算出する。上記ステップＳ３での計算結果から、各解析時間タイミング時において、流れ場を分割する各計算格子１０内に存在する気泡Ｇの個数ｎｇを知ることができる。そして、この気泡Ｇの個数ｎｇと、各格子１０の面積面積Ａ（＝ΔｘΔｙ）とを数式２１に代入する計算処理により気泡体積率αｇを算出し、この気泡体積率αｇまたは液体体積率αｌ（＝１−αｇ）を、例えば上記情報記録装置６の所定の保存領域内またはＲＡＭに保存する。即ち、本発明の「処理（ａ）」でいう「気液体積割合を、解析領域内における気泡の分布に基づき求める処理」、本発明の「ステップ（ｃ）」及び「手段（ｃ）」でいう「気液体積割合を、ステップ（ｂ）或いは手段（ｂ）で算出された各気泡の位置に基づき求めて記憶手段に保存する」に相当する処理を実行する。
【００５４】
続いて、ステップＳ５では、随時変化する放出渦要素１７のコア半径σαを数式２２から計算する。即ち、上記ステップＳ４で保存された気泡体積率αｇ等を読み出して、オイラーの前進差分法を数式２２に適用し、コア半径σαを計算するのである。即ち、請求項４でいう「各離散渦要素のコア半径を、体積割合の気泡流についてコアスプレッディング法を適用して算出」に相当する処理を実行する。
【００５５】
そして、ステップＳ６では、各計算格子１０内における総循環強度Γの時間変化率ΔΓ／Δｔを、保存された上記気泡体積率αｇ及び前回の解析時間タイミングで算出された液体Ｌの速度ｕｌを読み出して数式２４から算出し、例えば上記情報記録装置６の所定の保存領域内またはＲＡＭに保存する。即ち、本発明の処理（ｂ）でいう「各解析要素内における循環強度変化を、気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求める処理」、本発明の「ステップ（ｄ）」及び「手段（ｄ）」でいう「循環強度変化を、気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存」に相当する処理を実行する。
【００５６】
ステップＳ７では、放出渦要素１７の移流を数式２０から計算する。即ち、本発明の「ステップ（ａ）」及び「手段（ａ）」でいう「解析領域内に存在する各離散渦要素の各位置を算出」に相当する処理を実行する。より具体的には、上記気泡Ｇの運動計算時と同様、オイラーの前進差分法により数式２０を計算し、解析領域Ｈ内における各放出渦要素１７の位置を追跡するのである。
【００５７】
ステップＳ８では、上記ステップＳ６で保存された各計算格子１０内における総循環強度の時間変化率ΔΓ／Δｔデータと、上記ステップＳ７にて計算される各放出渦要素１７の位置とから、各計算格子１０内に存在する放出渦要素１７の循環強度Γα及び位置ｘαを求めて例えば上記情報記録装置６の所定の保存領域内またはＲＡＭに保存する。ここで、計算格子１０内に放出渦要素１７がｎｖ個存在する場合、Δｔｖ間における放出渦要素１７の１つ当たりの循環強度の変化量は、ΔΓ／ｎｖとする。また、放出渦要素１７が存在しない場合には、その計算格子１０の中央位置に１つについては循環強度の変化量ΔΓの放出渦要素１７をい１つ新しく発生させる。即ち、本発明の「処理（ｃ）」でいう「各離散渦要素の循環強度及び位置を、その解析要素内における循環強度変化と、各離散渦要素の位置分布とに基づきラグランジェ的に求める処理」、本発明の「ステップ（ｅ）」及び「手段（ｅ）」でいう「各解析要素内に存在する各離散渦要素の循環強度及び位置を、循環強度変化と、ステップ（ａ）で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存」に相当する処理を実行する。
【００５８】
ステップＳ９では、上記ステップＳ５で保存された放出渦要素１７のコア半径σαデータ、及び、ステップＳ８で保存された各放出渦要素１７の循環強度Γα及び位置ｘαデータを読み出して、各放出渦要素１７が与える渦度ωαの分布を、数式１９から算出する。また、上記固定渦要素１６の循環強度及び位置データから各固定渦要素１６が与える渦度ωの分布をも読み出す。そして、各格子点１１について、解析領域Ｈ内に存在する全ての離散渦要素（固定渦要素１６及び放出渦要素１７）によって与えられる渦度分布を重ね合わせることにより、各格子点１１における渦度ωを計算し、例えば上記情報記録装置６の所定の保存領域内またはＲＡＭに保存する。即ち、本発明の「処理（ｄ）」でいう「評価点における渦度を、各離散渦要素の循環強度及び位置に基づき求める処理」、本発明の「ステップ（ｆ）」及び「手段（ｆ）」でいう「評価点における渦度を、各離散渦要素の循環強度及び位置に基づき求める」に相当する処理を実行する。
【００５９】
そして、ステップ１０では、ステップＳ９で保存された各格子点１１における渦度ωを読み出して、数式１７の右辺に代入し、この数式１７を差分法により解析して液体Ｇの速度のベクトルポテンシャルψを計算する。なお、解析にはＳＯＲ法が利用できる。即ち、本発明の「処理（ｅ）」「ステップ（ｇ）」及び「手段（ｇ）」でいう「渦度から液体のベクトルポテンシャルを求め」に相当する処理を実行する。
【００６０】
ステップＳ１１では、ステップＳ４で求められた気泡体積率αｇ及び上記ステップＳ１０で求められたベクトルポテンシャルψとを数式１６に代入し、この数式１６を差分法により解析してスカラーポテンシャルφを計算する。なお、解析もＳＯＲ法が利用できる。即ち、本発明の「処理（ｅ）」「ステップ（ｇ）」及び「手段（ｇ）」でいう「液体のスカラーポテンシャルを、液体のベクトルポテンシャルと、体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める」に相当する処理を実行する。
【００６１】
ステップＳ１２では、ステップ１０でのベクトルポテンシャルψと、ステップＳ１１でのスカラーポテンシャルφとを用いて、全ての格子点１１における液体の速度を数式１４から計算し、例えば上記情報記録装置６の所定の保存領域内またはＲＡＭに保存する。ステップＳ１２で、解析継続の場合（ステップＳ１２で「Ｙ」）、ステップＳ１に戻る。
【００６２】
以上のステップを繰り返し実行することにより、気泡Ｇと液体Ｌの流れ、すなわち気泡流の時間変化をシミュレートすることができる。そして、以上の演算処理が終了すると、コンピュータ１は、その解析データを例えば表示装置４や印刷装置５に引き渡す。解析データは、気泡Ｇと液体Ｌの速度ｕｇ，ｕｌ、気泡体積率αｇ、気泡Ｇおよび離散渦要素１６，１７の位置座標であり、一般的に利用できる可視化ソフトを用いることにより、ベクトル表示や等高線表示できる。
【００６３】
なお、図６は計算結果の一例であり、平板Ｗ１の両側を異なる速度で鉛直上向きに流れる水が平板Ｗ１下流で混合し、その先端から気泡Ｇが放出された場合の結果であり、上記解析データを基にして、気泡Ｇ、離散渦要素１６，１７、水の速度ベクトルの分布を表示してある。
【００６４】
３．本実施形態の効果
本実施形態による気泡流の解析方法では、微分項を差分近似する必要がない。このため、従来の解析方法では数値不安定性に起因して計算が破綻するような、高レイノルズ数の流れの解析が可能であり、渦度の高い領域に渦要素が能動的に移流するため、高い解析精度が確保できる。
また、流体の流速や圧力を未知変数とする従来の解析方法では乱流モデルの利用が不可欠であったが、本実施形態の解析方法では乱流モデルを利用することなく行うことができ、気泡流の精密な解析が可能である。
【００６５】
更に、気泡流の解析精度は、単位時間当たりに放出される離散渦要素の個数に依存する。即ち、離散渦要素の個数が増すほど計算負荷の増加をもたらすが、高精度の解析結果を得ることができる。そこで、本実施形態では、上記離散渦要素１１６，１７の導入時に入力装置７（離散渦要素の導入個数を調整する調整手段として機能）での操作によって離散渦要素１６，１７の導入個数を変更可能に設定できる構成とすることができる。このような構成であれば、使用する計算機の能力と要求する解析精度に応じて、離散渦要素１６，１７の導入個数を調整することができ、使用者の要望に応じた自由度の高い解析を行うことができる。
＜他の実施形態＞
本発明は上記記述及び図面によって説明した実施形態に限定されるものではなく、例えば次のような実施形態も本発明の技術的範囲に含まれ、さらに、下記以外にも要旨を逸脱しない範囲内で種々変更して実施することができる。
（１）上記実施形態では、解析要素として矩形状の計算格子１０としたがこれに限らず、矩形以外の形状であってもよい。
【００６６】
（２）また、評価点としては計算格子の四隅の格子点１１としたが、これに限らず、計算格子とは別に設けた格子の格子点であっても勿論よい。また、評価点は１つであってもよい。
【００６７】
（３）上記実施形態では、放出渦要素１７の放出時間間隔と、解析時間タイミングとを一致させていたが、これに限らず、一致させない構成であってもよい。また、解析時間タイミングは、一定時間間隔に限らず、例えば入力装置７から指令信号が入力される毎のタイミングであってもよい。
【００６８】
（４）上記実施形態では、気泡流シミュレーションプログラムを例に挙げて説明したが、この気泡流シミュレーションプログラムと同様の動作をする各手段を備えた気泡流シミュレーション装置であってもよい。
【００６９】
（５）上記実施形態では、解析領域Ｈについて二次元解析を行った例を説明したが、本発明は三次元解析についても勿論適用することができる。
【００７０】
（６）上記実施形態では、気泡Ｇは円形形状としたが、これに限らず、高機能化を図る観点で、例えば楕円形状等であってもよい。この場合、例えば入力装置７にて気泡Ｇについての複数種の形状を選択可能とし、それら複数種の形状のそれぞれに対応した気泡体積率αｇの算出式（上記数式２１に準じた式）情報、或いは、複数種の形状のそれぞれに対応した気泡体積情報を記憶した記憶手段から、選択された種類の形状に対応した上記算出式または気泡体積情報を読み出してそれに基づき気泡体積率αｇを算出する構成としてもよい。
【００７１】
（７）上記実施形態では、平板Ｗ１や箱体Ｗ２の周りを流れる気泡流について解析した例を説明したが、これに限らず、曲面形状など任意の形状をもつ物体の周りを流れる気泡流についても解析することができるので、高機能のシミュレーションプログラムとなる。この場合、上記固定渦要素１６の配置作業において、解析領域Ｈ内に配置される物体の形状に沿って固定渦要素１６を配置するように構成すれば簡単に対応できる。
【００７２】
（８）上記実施形態では、気泡の運動及び気泡体積率の計算を行った後に、各放出渦要素の循環強度変化及び位置の算出を行う構成としたが、これに限られない。即ち、例えば、図１中の情報記録装置６に、渦要素情報保存領域と、気泡情報保存領域と、ポテンシャルデータ保存領域とを確保する。そして、渦要素情報保存領域に、上記ステップＳ６で算出される計算格子１０内における循環強度変化データと、上記ステップＳ７で算出される前記各離散要素１６，１７の位置及び循環強度データと、ステップＳ８で算出される格子点１１における渦度データとを保存する。気泡情報保存領域には、ステップＳ３で算出される各気泡Ｇの位置データと、ステップＳ４で算出される気泡体積率データを保存する。ポテンシャルデータ保存領域に、ステップＳ９，１０で算出される液体のベクトルポテンシャルデータと、液体のスカラーポテンシャルデータを保存する。そして、気泡Ｇに関する計算処理（ステップＳ３，４）と、離散渦要素１６，１７に関する計算処理（ステップＳ７）とを並列的に処理する構成であってもよい。
【図面の簡単な説明】
【００７３】
【図１】本発明の解析要素（計算格子）を説明するための図
【図２】本発明の一実施形態に係る気泡流シミュレーションプログラムが導入されるコンピュータ及び周辺機器を示す図
【図３】本発明の計算の対象とする気泡流を説明するための図
【図４】速度せん断層を誘起する平板および箱体の固定渦要素による表現を説明するための図
【図５】コンピュータによる処理内容を示すフローチャート
【図６】本実施形態による解析結果を示す図
【符号の説明】
【００７４】
１…コンピュータ
２…中央処理部
６…情報記録装置（記憶手段）
１０…計算格子（解析要素）
１１…格子点（評価点）
１５…速度せん断層
１６…固定渦要素（離散渦要素）
１７…放出渦要素（離散渦要素）
Ｇ…気泡
Ｈ…解析領域
Ｌ…液体
【図１】
2006085603.tif 000027
【図２】
2006085603.tif 000028
【図３】
2006085603.tif 000029
【図４】
2006085603.tif 000030
【図５】
2006085603.tif 000031
【図６】
2006085603.tif 000032