TOP > 国内特許検索 > 気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置 > 明細書

明細書 :気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4032123号 (P4032123)
公開番号 特開2006-085603 (P2006-085603A)
登録日 平成19年11月2日(2007.11.2)
発行日 平成20年1月16日(2008.1.16)
公開日 平成18年3月30日(2006.3.30)
発明の名称または考案の名称 気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置
国際特許分類 G06F  19/00        (2006.01)
FI G06F 19/00 110
請求項の数または発明の数 4
全頁数 20
出願番号 特願2004-272011 (P2004-272011)
出願日 平成16年9月17日(2004.9.17)
権利譲渡・実施許諾 特許権者において、実施許諾の用意がある。
審査請求日 平成18年8月9日(2006.8.9)
特許権者または実用新案権者 【識別番号】504139662
【氏名又は名称】国立大学法人名古屋大学
発明者または考案者 【氏名】内山 知実
【氏名】出川 智啓
早期審査対象出願または早期審理対象出願 早期審査対象出願
審査官 【審査官】松田 直也
参考文献・文献 特開2004-126925(JP,A)
特開2004-118394(JP,A)
特開2001-99748(JP,A)
特開平9-329109(JP,A)
特開2003-331208(JP,A)
特開2004-503005(JP,A)
特開2001-132700(JP,A)
特開平10-11420(JP,A)
調査した分野 G06F 19/00
JSTPlus(JDream2)
特許請求の範囲 【請求項1】
気泡と液体とが混在した気泡流の解析するための気泡流シミュレーションプログラムであって、コンピュータに、次のステップ(a)~(h)を含むステップを繰り返し実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて、時間的に変化する前記気泡流の挙動を解析する気泡流シミュレーションプログラム。
前記コンピュータは、前記ステップ(c)、(d)、(e)、(f)用の記憶手段を備え、前記気泡流シミュレーションプログラムは、
(a)前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置データを算出するステップ。
(b)前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出するステップ。
(c)前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記ステップ(b)で算出された各気泡の位置に基づき求めて記憶手段に保存するステップ。
(d)前記ステップ(c)によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存するステップ。
前記気泡流の渦度輸送方程式は、次の数式[数1]である。
【数1】
JP0004032123B2_000023t.gif

Γ :解析要素内における循環強度変化
α:気泡と液体との気液体積割合
g :重力加速度
:液体の速度
A :解析要素の面積
(e)前記ステップ(d)によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記ステップ(a)で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存するステップ。
(f)前記ステップ(e)によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求めて記憶手段に保存するステップ。
(g)前記ステップ(f)によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求めるステップ。
前記液体のスカラーポテンシャル式は、次の数式[数2]である。
【数2】
JP0004032123B2_000024t.gif

α:気泡と液体との気液体積割合
φ :液体のスカラーポテンシャル
ψ :液体のベクトルポテンシャル
(h) 前記ステップ(g)によって求められた前記評価点における前記液体の速度ポテンシャルおよびスカラーポテンシャルを用いて、前記評価点における前記液体の速度を計算するステップ。
【請求項2】
前記請求項1に記載の気泡流シミュレーションプログラムを記憶した記憶媒体。
【請求項3】
気泡と液体とが混在した気泡流の解析するための気泡流シミュレーション装置であって、次の手段(a)~(h)を備えるとともに、
前記手段(a)で算出される前記離散渦要素の位置データが保存される渦要素情報記憶手段と、 前記手段(b)で算出される前記気泡の位置データが保存される気泡情報記憶手段と、前記手段(c)、(d)、(e)、(f)用の記憶手段とを備えており、
前記手段(a)~(h)を繰り返し実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて、時間的に変化する前記気泡流の挙動を解析する気泡流シミュレーション装置。
(a)前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置データを算出する手段。
(b)前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出する手段。
(c)前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記手段(b)で算出された各気泡の位置に基づき求めて記憶手段に保存する手段。
(d)前記手段(c)によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存する手段。
前記気泡流の渦度輸送方程式は、次の数式[数1]である。
【数1】
JP0004032123B2_000025t.gif

Γ :解析要素内における循環強度変化
α:気泡と液体との気液体積割合
g :重力加速度
:液体の速度
A :解析要素の面積

(e)前記手段(d)によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記手段(a)で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存する手段。
(f)前記手段(e)によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求めて記憶手段に保存する手段。
(g)前記手段(f)によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める手段。
液体のスカラーポテンシャル式は、次の数式[数2]である。
【数2】
JP0004032123B2_000026t.gif

α:気泡と液体との気液体積割合
φ :液体のスカラーポテンシャル
ψ :液体のベクトルポテンシャル
(h) 前記手段(g)によって求められた前記評価点における前記液体の速度ポテンシャルおよびスカラーポテンシャルを用いて、前記評価点における前記液体の速度を計算する手段。
【請求項4】
前記手段(a)で算出される前記離散渦要素の位置データが保存される渦要素情報記憶手段と、
前記手段(b)で算出される前記気泡の位置データ及び前記手段(c)で算出される気液体積割合データが保存される気泡情報記憶手段とを備え、
前記手段(a)による処理と、前記手段(b)及び前記手段(c)による処理とを並列的に実行することを特徴とする請求項3に記載の気泡流シミュレーション装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、液中に分散した気泡が液体と相互作用を及ぼし合いながら流れる気泡流(二相流)の挙動の解析方法に関する。
【背景技術】
【0002】
微小な気泡と液体とが混在して相互作用を及ぼし合いながら流れる気泡流は、例えば反応装置、熱交換器及び廃水処理装置など、様々な工業装置において観察される。このため、この気泡流の挙動について、合理的な解析方法の確立が望まれている。
【0003】
従来の解析方法では、流れの変数として速度、圧力及び気泡体積率が用いられ、それらを支配する複数の偏微分方程式を様々な数値解法(差分法、有限体積法、有限要素法)によって計算するようにしていた。これらの偏微分方程式は、微分を差分で近似表現することにより、代数方程式に置き換えられるが、非線形項が数値不安定性を有するため、気泡流のレイノルズ数が高い場合には計算過程で発散が生じるおそれがある。特に、工業装置で見受けられる気泡流は、レイノルズ数が高いものが多く、精度の高い気泡流解析ができないという問題があった。
【0004】
ところで、従来から、液体又は気体のみの単相流の解析については渦法が有効に利用されている。ここで、渦法とは、渦度をもつ渦要素を追跡して渦度場の時間変化を求めるLagrange 型解法であり,大規模渦の形成や変形など渦構造の発展過程を良好に計算できるものである。そして、下記非特許文献1には、この渦法を利用して、気泡流の挙動を解析する技術が開示されている。

【非特許文献1】Sene, K. J., et al., The role of coherent structures in bubble transport by turbulent shear flows,J. Fluid Mech., 259(1994), 219-240.
【特許文献1】特開2004-126925公報
【特許文献2】特開2004-118394公報
【特許文献3】特開2001-99748公報
【特許文献4】特開平9-329109号公報
【特許文献5】特開2003-331208公報
【特許文献6】特開2004-503005公報
【特許文献7】特開2001-132700公報
【発明の開示】
【発明が解決しようとする課題】
【0005】
しかしながら、上記非特許文献1に記載の技術は、気泡が液相に影響を及ぼさないことを前提として気泡流の挙動を解析するものであり、液中に分散した気泡が液体と相互作用を及ぼし合いながら流れる気泡流の挙動を精度良く解析できるものではなかった。
【0006】
なお、上記特許文献1~7には、流体の挙動をその支配方程式から計算によって予測しようとする技術が開示されているが、いずれも、本願が解析対象とする気泡流を精度よく解析できるものではない。
【0007】
即ち、特許文献1の技術は、渦法を用いて流体の流れを解析するものであるが、水または空気などの単相流を解析対象とするものである。
特許文献2の技術は、気体と液体とが接する面、即ち気液自由界面の変形に関するものであり、気体と液体とが分離した流れを解析対象としている。
特許文献3の技術は、やはり単相流の挙動を有限差分法によって解析するものである。
【0008】
特許文献4の技術は、流水中の固体壁面から噴き出された気泡について、その運動と分布状態とを計算し、気泡による水流の乱れ変化を求め、壁面における流体摩擦力の低減効果を解析するものである。壁面近傍での流れ現象のみを解析対象としており、水流については乱れ特性など限定された現象しか解析できず、広範囲における圧力や速度などの分布が考慮されていない。
特許文献5の技術は、飛沫や気泡の挙動を解析するものであり、その計算に有限差分法に準じた解法を用い、飛沫や気泡の存在の判定には、VOF法を利用している。VOF法は気体と飛沫、気体と気泡の界面の分解能が低い欠点を持つことが知られている。
【0009】
特許文献6の技術は、気体と液体とが分離した状態で存在し、両者の自由界面がある流体の挙動を解析するものであり、流速が低く層流状態にある流れが解析対象となっている。
特許文献7の技術は、タンク内の水流を計算し、飽和水蒸気圧と計算で得られた圧力とを比較することにより、気泡の発生位置を予測するものであるが、気泡と液体とが相互作用を及ぼしあうことについて考慮されていない。
【0010】
本発明は上記のような事情に基づいて完成されたものであって、気体と液体とが混在する気泡流の挙動を精度よく解析することが可能な気泡流の解析方法、気泡流シミュレーションプログラム及びそれを記憶した記憶媒体並びに気泡流シミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【0011】
上記の目的を達成するための手段として、請求項1の発明に係る気泡流解析方法は、気泡と液体とが混在した気泡流の解析方法であって、解析領域内に存在する気泡及び離散渦要素の位置を時間的に変化させつつ、次の処理(a)~(e)を行って前記解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
(a)前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記解析領域内における気泡の分布に基づき求める処理。
(b)前記各解析要素内における循環強度変化を、前記気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求める処理。
(c)前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、その解析要素内における前記循環強度変化と、前記各離散渦要素の位置分布とに基づきラグランジェ的に求める処理。
(d)前記解析領域内の前記評価点における渦度を、前記各離散渦要素の循環強度及び位置に基づき求める処理。
(e)前記評価点における渦度から前記液体のベクトルポテンシャルを求め、前記解析領域内の前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める処理。
【0012】
請求項2の発明は、請求項1に記載の気泡流の解析方法において、前記気泡流の過度輸送方程式は、次の数式[数1]であることを特徴とする。
【数1】
JP0004032123B2_000002t.gif
Γ :解析要素内における循環強度変化
α:気泡と液体との気液体積割合
g :重力加速度
:液体の速度
A :解析要素の面積
【0013】
請求項3の発明は、請求項1又は請求項2に記載の気泡流の解析方法において、前記液体のスカラーポテンシャル式は、次の数式[数2]であることを特徴とする。
【数2】
JP0004032123B2_000003t.gif
α:気泡と液体との気液体積割合
φ :液体のスカラーポテンシャル
ψ :液体のベクトルポテンシャル
【0014】
請求項4の発明は、請求項1~請求項3のいずれかに記載の気泡流の解析方法において、前記処理(d)において、前記各離散渦要素のコア半径を、前記体積割合の気泡流についてコアスプレッディング法を適用して算出し、前記渦度を求めることを特徴とする。
請求項5の発明は、請求項4に記載の気泡流の解析方法において、前記体積割合の気泡流についてコアスプレッディング法を適用した算出式は、次の数式[数3]であることを特徴とする。
【数3】
JP0004032123B2_000004t.gif
σα:離散渦要素のコア半径
α:気泡と液体との気液体積割合
ν:液体の動粘度
【0015】
請求項6の発明は、請求項1~請求項5のいずれかに記載の気泡流の解析方法において、前記解析領域内に存在する気泡の形状、大きさ及び個数のうち少なくともいずれか1つを変更可能とし、前記処理(a)では、この変更された気泡の形状、大きさ及び個数に応じた気泡の分布に基づき前記気液体積割合を求めることを特徴とする。
なお、上記請求項2,3,6の構成はいずれも、次述する請求項7~10の発明に対して適用することができる。また、請求項4,5の構成はいずれも、次述する請求項7,8のステップ(f)、請求項9,10の手段(f)にそれぞれ適用することができる。
【0016】
請求項7の発明に係る気泡流シミュレーションプログラムは、気泡と液体とが混在した気泡流の解析するための気泡流シミュレーションプログラムであって、コンピュータに、次のステップ(a)~(g)を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
(a)前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出するステップ。
(b)前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出するステップ。
(c)前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記ステップ(b)で算出された各気泡の位置に基づき求めて記憶手段に保存するステップ。
(d)前記ステップ(c)によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存するステップ。
(e)前記ステップ(d)によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記ステップ(a)で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存するステップ。
(f)前記ステップ(e)によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求めるステップ。
(g)前記ステップ(f)によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求めるステップ。
【0017】
請求項8の発明は、上記請求項7に記載の気泡流シミュレーションプログラムを記憶した記憶媒体である。なお、記憶媒体としては、フロッピーディスク(登録商標)、CD-ROM等が好適である。
【0018】
請求項9の発明に係る気泡流シミュレーション装置は、気泡と液体とが混在した気泡流の解析するための気泡流シミュレーション装置あって、次の手段(a)~(g)を備え、これらの手段を実行させて、解析領域内に想定した評価点について求まる前記液体のスカラーポテンシャル及びベクトルポテンシャルの時間的変化に基づいて前記気泡流の挙動を解析する。
(a)前記解析領域内に複数の離散渦要素を順次導入し、そのときに前記解析領域内に存在する各離散渦要素の各位置を算出する手段。
(b)前記解析領域内に複数の気泡を導入し、そのときに前記解析領域内に存在する各気泡の位置を算出する手段。
(c)前記解析領域を複数に区画した各解析要素内における気泡と液体との気液体積割合を、前記手段(b)で算出された各気泡の位置に基づき求めて記憶手段に保存する手段。
(d)前記手段(c)によって保存された気液体積割合データを読み出して、前記各解析要素内における循環強度変化を、当該気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存する手段。
(e)前記手段(d)によって保存された前記循環強度変化データを読み出して、前記各解析要素内に存在する各離散渦要素の循環強度及び位置を、当該循環強度変化と、前記手段(a)で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存する手段。
(f)前記手段(e)によって保存された各離散渦要素の循環強度及び位置データを読み出して、前記解析領域内の前記評価点における渦度を、当該各離散渦要素の循環強度及び位置に基づき求める手段。
(g)前記手段(f)によって保存された前記評価点における渦度データを読み出して当該渦度から前記液体のベクトルポテンシャルを求め、前記評価点における前記液体のスカラーポテンシャルを、当該液体のベクトルポテンシャルと、前記体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める手段。
【0019】
請求項10の発明は、請求項9に記載の気泡流シミュレーション装置において、前記手段(a)で算出される前記離散渦要素の位置データが保存される渦要素情報記憶手段と、前記手段(b)で算出される前記気泡の位置データ及び前記手段(c)で算出される気液体積割合データが保存される気泡情報記憶手段とを備え、前記手段(a)による処理と、前記手段(b)及び前記手段(c)による処理とを並列的に実行することを特徴とする。
なお、より具体的な構成として、次のようなものであってもよい。
「渦要素情報保存手段と、気泡情報保存手段と、ポテンシャルデータ保存手段とを備え、
前記渦要素情報保存手段には、前記手段(a)で算出される前記各離散要素の位置データと、前記手段(d)で算出される前記解析要素内における循環強度変化データと、前記手段(e)で算出される前記各離散渦要素の循環強度及び位置データと、前記手段(f)で算出される前記評価点における渦度データと、が保存され、
前記気泡情報保存手段には、前記手段(b)で算出される前記各気泡の位置データと、前記手段(c)で算出される前記気液体積割合データと、が保存され、
前記ポテンシャルデータ保存手段には、前記手段(g)で算出される前記液体のベクトルポテンシャルデータと、前記液体のスカラーポテンシャルデータと、が保存されるよう構成される構成。」
【発明の効果】
【0020】
本発明は、気泡と液体とが混在した気泡流の渦度場に着目して複数の渦要素で離散化し、気泡流内の気泡と液体との気液体積割合を考慮した運動量保存関係及び質量保存関係から渦度輸送方程式を導出し、この渦度輸送方程式のラグランジェ計算により離散渦要素の挙動を求めるものである。
渦度輸送方程式のラグランジェ計算には非線形項が現れないため、レイノルズ数に依存することなく、高レイノルズ数の気泡流について高精度解析が可能になる。また、渦度の高い領域には渦要素が密集する傾向があるが、このような領域に対しても気泡流について高い解析精度を確保できる。
【0021】
また、渦度を求めるにあたって、各離散渦要素のコア半径を、体積割合の気泡流についてコアスプレッディング法を適用して算出して適用することが望ましい。このようにすれば、粘性拡散による減衰をも考慮した正確な渦度を求めることができる(請求項4,5)。
【0022】
また、請求項6の発明によれば、解析領域内に存在する気泡の形状、大きさ及び個数のうち少なくともいずれか1つを変更可能とし、処理(a)では、この変更された気泡の形状、大きさ及び個数に応じた気泡の分布に基づき気液体積割合を求めるものである。従って、気泡の形状、大きさや個数を変更することにより、気泡流内に存在する気泡の変形(径の変動を含む)、分裂、合体などの各現象に対応して気泡流の解析を行うことができる。
【0023】
更に、請求項10の発明によれば、手段(a)で算出される離散渦要素の位置データが保存される渦要素情報記憶手段と、手段(b)で算出される気泡の位置データ及び手段(c)で算出される気液体積割合データが保存される気泡情報記憶手段とが備えられ、手段(a)による処理と、手段(b)及び手段(c)による処理とを並列的に実行する構成とした。これにより、気泡流の解析処理の高速化を図ることができる。
【0024】
<理論>
次の本発明の理論について説明する。
1.記号
本発明において下記に登場する各符号は次のように定義されている。
JP0004032123B2_000005t.gif
【0025】
2.仮定
本発明は以下の仮定を前提としている。
(a)流れは、流れ場の代表長さに比べて微細な気泡が分散した気泡流を対象とする。
(b)液体及び気泡は、非圧縮性であり、また、熱的非平衡はない。
(c)気泡の質量及び運動量は、気体に比べて十分小さく無視できる。
(d)気泡は球形を保って運動し、合体及び分裂はしない。
【0026】
3.基礎式
本発明では、所定の設定データ(後述する入力データ)に基づいて,計算に必要な境界条件を設定し、流れの基礎方程式に基づいて計算を実行する。以下、境界条件の設定と計算に用いる流れの基礎方程式について説明する。
【0027】
(1)気泡流の保存方程式
液体中に微細な気泡が分散して流れる気泡流について、上記仮定(a)~(c)を前提とした場合、この気泡流の質量および運動量の保存方程式は、それぞれ次式で表される。
【数10】
JP0004032123B2_000006t.gif
【数11】
JP0004032123B2_000007t.gif
但し、液体と気泡の体積率は次の関係を満たす。
【数12】
JP0004032123B2_000008t.gif

【0028】
(2)気泡の運動方程式
気泡には仮想質量力、流動抗力、重力および揚力が作用するものとする。気泡の運動方程式は、上記仮定(d)を前提とした場合、次式で与えられる。
【数13】
JP0004032123B2_000009t.gif

【0029】
(3)液体速度場の直交分解
ヘルムホルム(Helmholts)
の定理によれば、任意のベクトル場はスカラポテンシャルφ の勾配とベクトルポテンシャルψ の回転との和として表されることから、液相速度u
は次式のように記述できる。
【数14】
JP0004032123B2_000010t.gif
ここで、ψはスカラ関数の勾配を加えても結果が不変である。この任意性を取り除き、ψを一意的に定めるため、ψはソレノイダルであるものとする。すなわち、
【数15】
JP0004032123B2_000011t.gif

【0030】
数式14 を数式10に代入し、恒等式∇ ・ (∇×ψ)
=0を利用して変形すれば次式が得られる。
【数16】
JP0004032123B2_000012t.gif
そして、数式14の回転をとった後に、数式15を代入すれば、次のベクトルPoisson
方程式が得られる。
【数17】
JP0004032123B2_000013t.gif
ここで、ωは液体の渦度である。
【0031】
4.数値解法
(1)渦要素と気泡とのラグランジェ(Lagrange)解析
上記数式11の回転をとれば気泡流の渦度方程式が得られ、二次元解析の場合には次式で表される。
【数18】
JP0004032123B2_000014t.gif

【0032】
次に、渦度場を多数の渦要素により離散化する。単相流解析に対する渦要素モデルを適用し、コア構造をもつ離散渦要素を導入する。離散渦要素α の循環強度をΓα、コア半径をσα、位置ベクトルをxα
とすれば、離散渦要素α による位置x(評価点の位置)における渦度は次式で与えられる。
【数19】
JP0004032123B2_000015t.gif
また、流れ場に存在する離散渦要素αの移流は次のラグランジェ(Lagrange)解析から計算できる。
【数20】
JP0004032123B2_000016t.gif

【0033】
x-y平面内の解析領域Hを、図1に示すような、寸法がΔxとΔyの四角形の計算格子10(解析要素)に分割し、上記数式16と数式17を差分法で解析し、得られたスカラーポテンシャルφとベクトルポテンシャルψを数式14に代入して格子点11(評価点)上の液体の速度uを求める。この際、格子点11で与える渦度ω(同図ω1~ω4)の値は、各離散渦要素による渦度を数式19から求め、全離散渦要素について加え合わせることにより計算する。各計算格子10における液体体積率αは、数式13のLagrange 解析から気泡運動を解析して気泡体積率αを求めたのち、数式12を利用して算出する。
【0034】
また、例えば、解析領域Hの入口境界において、x-y平面の垂直方向(図1で紙面奥行き方向)に気泡が間隔Δzで一様に分布するものと仮定し、面積A(=ΔxΔy)の計算格子10におけるαを次式で計算する。
【数21】
JP0004032123B2_000017t.gif
ここで、n は任意の時間における計算格子10内に存在する直径dの気泡Gの個数(図1では4つを図示)である。
渦度は、数式18から分かるように、粘性拡散(右辺第1項)と相分布の勾配(右辺第2項)により時間変化する。これらの変化は、数式19における離散渦要素のコア半径σαと循環強度Γαの変化として模擬し、以下に示すように数式20のLagrange 計算と連立させて解析する。
【0035】
(2)粘性拡散によるコア半径の変化
渦度は粘性拡散により減衰する。ここで、液体のみ又は気体のみの単相流解析では、コア半径の大きさを時間的に変化させるコアスプレッディング(Core spreading)法により考慮される。そこで、本発明においても、これを準用し、次式をLagrange 解析してコア半径の時間変化を求める。
【数22】
JP0004032123B2_000018t.gif
なお、Core spreading 法は、Leonard, A., Vortex methods for flow simulation, J. Comput. Phys. 37, (1980), 289-335 に記載されている。
【0036】
(3)気泡運動による循環強度の変化
任意の閉曲線まわりの循環強度Γ の時間変化率は、レイノルズ(Reynolds) の輸送定理を用いれば、次式で表される.
【数23】
JP0004032123B2_000019t.gif
そして、この数式23の右辺に上記数式18を代入すれば次式を得ることができる。
【数24】
JP0004032123B2_000020t.gif
ここで、粘性拡散項は、数式22で考慮されているので、代入に際し無視している。
【0037】
上記計算格子10に数式24を適用し、各計算格子10内における循環強度Γ の時間変化率ΔΓ/Δtを計算する。ただし、数式24の被積分関数は格子点11間で線形変化するものと仮定して数値積分する。各計算格子10内にnv個の離散渦要素が存在する場合、Δtにおける離散渦要素1つ当たりの循環強度の変化量をΔΓ/nvとする。離散渦要素が存在しない場合には、循環強度ΔΓの離散渦要素1つをその計算格子10の中央から新しく発生させる。
【0038】
(4)解析手順
時刻t=tでの流れ場における気泡流の挙動(各格子点11における液体の速度に関するベクトルポテンシャル及びスカラーポテンシャル等)が既知であるならば、時刻t=t+Δtの挙動を、例えば以下の手順で計算できる。
(a)気泡の運動を数式13により計算し、気液体積割合(α(=1-α))を数式21から求める。
(b)離散渦要素のコア半径σを数式22から求める。
(c)各計算格子10内における循環強度の変化量ΔΓを数式24から求める。
(d)各離散渦要素の移流を数式20により計算し、各格子点11における渦度ωを数式19から求める。
(e)各格子点11における液体の速度のベクトルポテンシャルψを数式17から求める。
(f)各格子点11における液体の速度のスカラーポテンシャルφを数式16から求める。
(g)各格子点11における液体の速度uを、数式14から求める。
【0039】
このような手順を予め定めた時間間隔毎に繰り返し実行することで、気泡流について高精度解析を行うことができる。
【発明を実施するための最良の形態】
【0040】
以下、本発明の一実施形態を前述した理論及び図1~図6によって説明する。
【0041】
本実施形態に係る気泡流シミュレーションプログラムは、後述する情報記憶装置6に記録され、パーソナルコンピュータ(以下「コンピュータ1」という)に導入することにより、そのコンピュータ1に後述する気泡流の挙動の解析ステップを実行させるためのものである。
【0042】
1.コンピュータとその周辺機器の構成
図2には、コンピュータ1とその周辺機器の構成とを示した。コンピュータ1には、演算及び周辺機器の制御等を行うための中央処理部2、この中央処理部2と周辺機器との接続を行うためのI/Oインターフェース3、及び周辺機器である表示装置4、印刷装置5、情報記憶装置6及び入力装置7が接続されている。
【0043】
なお、表示装置4には、カソードレイチューブ(CRT)の他に液晶装置、液晶プロジェクタ等が含まれる。また、印刷装置5は、コンピュータ1に直接に連結されている必要はなく、LAN等のネットワークを介して連結されていてもよい。また、入力装置7は、例えばキーボードやマウスなどのコンソールであって、ここでの操作によって解析領域Hの境界を設定する座標位置データや流体特性を設定する入力データなどを入力することができる。
【0044】
中央処理部2の詳細については図示しないが、中央演算処理装置(CPU)の他に、ROM、各種RAM、チップセット等から構成されている。情報記憶装置6は、ハードディスクやフロッピーディスク、及びMOやCD-ROM等であり、コンピュータ1のプログラムや各種データを読み込みまたは記録することができる。
【0045】
2.気泡流の解析
さて、本実施形態の作用を上記中央処理部2にて実行される気泡流シミュレーションプログラムのフローチャートを参照しつつ説明する。
【0046】
(1)計算に必要な入力データの入力及び境界条件の設定
まず、上記入力装置7での操作によって、計算に必要な入力データを入力する。ここで、入力データとしては、液体Lの密度ρ及び動粘度ν、気泡Gの密度ρ及び直径d、解析領域Hの入口における液体L及び気泡Gの初期速度u(0),u(0)、同じく解析領域Hの入口における単位時間当たりの気泡Gの放出個数、解析領域Hの幾何形状及び寸法、解析時間、解析時間タイミングΔtv、計算格子10の個数及び寸法である。
【0047】
(1)条件設定
本実施形態では、壁面で囲まれていない自由空間における、気泡Gを含む乱流(気泡流)に対する数値解析である。より具体的には、図3に示すような、異なる速度をもつ液体Lが混合して速度せん断層15を形成し、そこに気泡Gが混在して流れる、気泡流を解析対象としているのである。
【0048】
ここで、上記速度せん断層15を形成する平板W1,箱体W2などの物体は、図4に示すように、これらの表面に循環強度Γをもつ固定渦要素16(離散渦要素)を配置して表す。この固定渦要素16の個数や位置は、上記入力装置7での操作によって行う。物体表面の任意の位置における液体の速度uは、固定渦要素16と放出される放出渦要素17(離散渦要素)の循環強度によって求められる。各固定渦要素16の循環強度Γは、平板W1、箱体W2の表面上におけるn箇所で法線方向の速度成分が零となるように定められ、この計算は、コンピュータ1によって次の数式25に示す連立方程式を解くことによって実行される。なお、この際、ケルビンの循環定理に従い、固定渦要素16の循環強度の総和と、放出渦要素17の循環強度の総和とが等しくなる条件が付与されている。
【数25】
JP0004032123B2_000021t.gif
ここで、固定渦要素の個数をnとすれば、Aはn×nの行列、Bはn列のベクトルである。
【0049】
次いで、速度せん断層15は、流れ場のはく離位置から放出渦要素17を放出することにより近似的に表現される。放出渦要素17の初期循環強度Γ0は、速度せん断層15の両側の液体の速度U1,U2を用いた次の式で与えられる。
【数26】
JP0004032123B2_000022t.gif
ここで、Δtvは、放出渦要素の放出時間間隔であり、本実施形態では上記解析時間タイミングに一致させている。解析開始の操作がされると、コンピュータ1は、上記数式26で算出される初期循環強度Γ0の放出渦要素17を、流れ場のはく離位置即ち速度せん断層15の始点から一定の時間間隔Δtvで順次放出するシミュレートを実行する。始点は、異なる速度をもつ流体が合流する位置或いは物体角部が相当し、本実施形態では、図3中、平板W1の下流側先端P1、箱体W2の両側面の下流端P2,P2である。
【0050】
(2)解析処理
コンピュータ1は、上記条件設定後、図5のフローチャートに示す処理を上記解析時間タイミングΔtvで繰り返し実行し、順次算出される各格子点11における液体Lの速度のベクトルポテンシャル及びスカラーポテンシャルの時間的変化に基づき気泡流の挙動を数値解析する。
【0051】
ステップS1で解析時間タイミングになったかどうかを判断し、解析時間タイミングになったとき(ステップS1で「Y」)、ステップS2で、上記入力装置7にて入力された気泡Gの放出個数分の気泡を解析領域Hの入口から放出する。即ち、本発明の上記「ステップ(b)」及び「手段(b)」でいう「解析領域内に複数の気泡を導入」に相当する処理を実行する。それとともに、速度せん断層から放出渦要素17を放出するシミュレートを実行する。即ち、本発明の上記「ステップ(a)」及び「手段(a)」でいう「解析領域内に複数の離散渦要素を順次導入」に相当する処理を実行する。
【0052】
次に、ステップS3で、気泡Gの運動を、前回の解析時間タイミングで算出された液体Lの速度uに基づき数式13から算出する。即ち、本発明の上記「ステップ(b)」及び「手段(b)」でいう「解析領域内に存在する各気泡の各位置を算出」に相当する処理を実行する。より具体的には、数式13にオイラーの前進差分法を適用し、その右辺に現れる液体の速度を格子点11における速度の補完から定め、ラグランジェ(Lagrange)計算により解析領域H内における各気泡Gの位置を追跡するのである。なお、補完には評価したい位置の近傍格子点11を選択し、それらに関して内挿法が利用できる。
【0053】
ステップS4では、気泡体積率αを算出する。上記ステップS3での計算結果から、各解析時間タイミング時において、流れ場を分割する各計算格子10内に存在する気泡Gの個数nを知ることができる。そして、この気泡Gの個数nと、各格子10の面積面積A(=ΔxΔy)とを数式21に代入する計算処理により気泡体積率αを算出し、この気泡体積率αまたは液体体積率α(=1-α)を、例えば上記情報記録装置6の所定の保存領域内またはRAMに保存する。即ち、本発明の「処理(a)」でいう「気液体積割合を、解析領域内における気泡の分布に基づき求める処理」、本発明の「ステップ(c)」及び「手段(c)」でいう「気液体積割合を、ステップ(b)或いは手段(b)で算出された各気泡の位置に基づき求めて記憶手段に保存する」に相当する処理を実行する。
【0054】
続いて、ステップS5では、随時変化する放出渦要素17のコア半径σαを数式22から計算する。即ち、上記ステップS4で保存された気泡体積率α等を読み出して、オイラーの前進差分法を数式22に適用し、コア半径σαを計算するのである。即ち、請求項4でいう「各離散渦要素のコア半径を、体積割合の気泡流についてコアスプレッディング法を適用して算出」に相当する処理を実行する。
【0055】
そして、ステップS6では、各計算格子10内における総循環強度Γの時間変化率ΔΓ/Δtを、保存された上記気泡体積率α及び前回の解析時間タイミングで算出された液体Lの速度uを読み出して数式24から算出し、例えば上記情報記録装置6の所定の保存領域内またはRAMに保存する。即ち、本発明の処理(b)でいう「各解析要素内における循環強度変化を、気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求める処理」、本発明の「ステップ(d)」及び「手段(d)」でいう「循環強度変化を、気液体積割合の気泡流についての運動量保存関係から導出された気泡流の渦度輸送方程式と、レイノルズの輸送定理とに基づき求めて記憶手段に保存」に相当する処理を実行する。
【0056】
ステップS7では、放出渦要素17の移流を数式20から計算する。即ち、本発明の「ステップ(a)」及び「手段(a)」でいう「解析領域内に存在する各離散渦要素の各位置を算出」に相当する処理を実行する。より具体的には、上記気泡Gの運動計算時と同様、オイラーの前進差分法により数式20を計算し、解析領域H内における各放出渦要素17の位置を追跡するのである。
【0057】
ステップSでは、上記ステップS6で保存された各計算格子10内における総循環強度の時間変化率ΔΓ/Δtデータと、上記ステップS7にて計算される各放出渦要素17の位置とから、各計算格子10内に存在する放出渦要素17の循環強度Γα及び位置xαを求めて例えば上記情報記録装置6の所定の保存領域内またはRAMに保存する。ここで、計算格子10内に放出渦要素17がnv個存在する場合、Δtv間における放出渦要素17の1つ当たりの循環強度の変化量は、ΔΓ/nvとする。また、放出渦要素17が存在しない場合には、その計算格子10の中央位置に1つについては循環強度の変化量ΔΓの放出渦要素17をつ新しく発生させる。即ち、発明の「ステップ(e)」及び「手段(e)」でいう「各解析要素内に存在する各離散渦要素の循環強度及び位置を、循環強度変化と、ステップ(a)で算出された各離散渦要素の位置に基づきラグランジェ的に求めて記憶手段に保存」に相当する処理を実行する。
【0058】
ステップSでは、上記ステップS5で保存された放出渦要素17のコア半径σαデータ、及び、ステップSで保存された各放出渦要素17の循環強度Γα及び位置xαデータを読み出して、各放出渦要素17が与える渦度ωαの分布を、数式19から算出する。また、上記固定渦要素16の循環強度及び位置データから各固定渦要素16が与える渦度ωの分布をも読み出す。そして、各格子点11について、解析領域H内に存在する全ての離散渦要素(固定渦要素16及び放出渦要素17)によって与えられる渦度分布を重ね合わせることにより、各格子点11における渦度ωを計算し、例えば上記情報記録装置6の所定の保存領域内またはRAMに保存する。即ち本発明の「ステップ(f)」及び「手段(f)」でいう「評価点における渦度を、各離散渦要素の循環強度及び位置に基づき求める」に相当する処理を実行する。
【0059】
そして、ステップでは、ステップSで保存された各格子点11における渦度ωを読み出して、数式17の右辺に代入し、この数式17を差分法により解析して液体Gの速度のベクトルポテンシャルψを計算する。なお、解析にはSOR法が利用できる。即ち、本発明のステップ(g)」及び「手段(g)」でいう「渦度から液体のベクトルポテンシャルを求め」に相当する処理を実行する。
【0060】
ステップS10では、ステップS4で求められた気泡体積率α及び上記ステップSで求められたベクトルポテンシャルψとを数式16に代入し、この数式16を差分法により解析してスカラーポテンシャルφを計算する。なお、解析もSOR法が利用できる。即ち、本発明のステップ(g)」及び「手段(g)」でいう「液体のスカラーポテンシャルを、液体のベクトルポテンシャルと、体積割合の気泡流の質量保存関係及びヘルムホルツの定理によって導出された液体のスカラーポテンシャル式とに基づき求める」に相当する処理を実行する。
【0061】
ステップS1では、ステップ9でのベクトルポテンシャルψと、ステップS1でのスカラーポテンシャルφとを用いて、全ての格子点11における液体の速度を数式14から計算し(即ち、本発明の「ステップ(h)」及び「手段(h)」でいう「・・評価点における前記液体の速度ポテンシャルおよびスカラーポテンシャルを用いて、前記評価点における液体の速度を計算する」に相当する処理を実行し)、例えば上記情報記録装置6の所定の保存領域内またはRAMに保存する。ステップS12で、解析継続の場合(ステップS12で「Y」)、ステップS1に戻る。
【0062】
以上のステップを繰り返し実行することにより、気泡Gと液体Lの流れ、すなわち気泡流の時間変化をシミュレートすることができる。そして、以上の演算処理が終了すると、コンピュータ1は、その解析データを例えば表示装置4や印刷装置5に引き渡す。解析データは、気泡Gと液体Lの速度u,u、気泡体積率α、気泡Gおよび離散渦要素16,17の位置座標であり、一般的に利用できる可視化ソフトを用いることにより、ベクトル表示や等高線表示できる。
【0063】
なお、図6は計算結果の一例であり、平板W1の両側を異なる速度で鉛直上向きに流れる水が平板W1下流で混合し、その先端から気泡Gが放出された場合の結果であり、上記解析データを基にして、気泡G、離散渦要素16,17、水の速度ベクトルの分布を表示してある。
【0064】
3.本実施形態の効果
本実施形態による気泡流の解析方法では、微分項を差分近似する必要がない。このため、従来の解析方法では数値不安定性に起因して計算が破綻するような、高レイノルズ数の流れの解析が可能であり、渦度の高い領域に渦要素が能動的に移流するため、高い解析精度が確保できる。
また、流体の流速や圧力を未知変数とする従来の解析方法では乱流モデルの利用が不可欠であったが、本実施形態の解析方法では乱流モデルを利用することなく行うことができ、気泡流の精密な解析が可能である。
【0065】
更に、気泡流の解析精度は、単位時間当たりに放出される離散渦要素の個数に依存する。即ち、離散渦要素の個数が増すほど計算負荷の増加をもたらすが、高精度の解析結果を得ることができる。そこで、本実施形態では、上記離散渦要素116,17の導入時に入力装置7(離散渦要素の導入個数を調整する調整手段として機能)での操作によって離散渦要素16,17の導入個数を変更可能に設定できる構成とすることができる。このような構成であれば、使用する計算機の能力と要求する解析精度に応じて、離散渦要素16,17の導入個数を調整することができ、使用者の要望に応じた自由度の高い解析を行うことができる。
<他の実施形態>
本発明は上記記述及び図面によって説明した実施形態に限定されるものではなく、例えば次のような実施形態も本発明の技術的範囲に含まれ、さらに、下記以外にも要旨を逸脱しない範囲内で種々変更して実施することができる。
(1)上記実施形態では、解析要素として矩形状の計算格子10としたがこれに限らず、矩形以外の形状であってもよい。
【0066】
(2)また、評価点としては計算格子の四隅の格子点11としたが、これに限らず、計算格子とは別に設けた格子の格子点であっても勿論よい。また、評価点は1つであってもよい。
【0067】
(3)上記実施形態では、放出渦要素17の放出時間間隔と、解析時間タイミングとを一致させていたが、これに限らず、一致させない構成であってもよい。また、解析時間タイミングは、一定時間間隔に限らず、例えば入力装置7から指令信号が入力される毎のタイミングであってもよい。
【0068】
(4)上記実施形態では、気泡流シミュレーションプログラムを例に挙げて説明したが、この気泡流シミュレーションプログラムと同様の動作をする各手段を備えた気泡流シミュレーション装置であってもよい。
【0069】
(5)上記実施形態では、解析領域Hについて二次元解析を行った例を説明したが、本発明は三次元解析についても勿論適用することができる。
【0070】
(6)上記実施形態では、気泡Gは円形形状としたが、これに限らず、高機能化を図る観点で、例えば楕円形状等であってもよい。この場合、例えば入力装置7にて気泡Gについての複数種の形状を選択可能とし、それら複数種の形状のそれぞれに対応した気泡体積率αの算出式(上記数式21に準じた式)情報、或いは、複数種の形状のそれぞれに対応した気泡体積情報を記憶した記憶手段から、選択された種類の形状に対応した上記算出式または気泡体積情報を読み出してそれに基づき気泡体積率αを算出する構成としてもよい。
【0071】
(7)上記実施形態では、平板W1や箱体W2の周りを流れる気泡流について解析した例を説明したが、これに限らず、曲面形状など任意の形状をもつ物体の周りを流れる気泡流についても解析することができるので、高機能のシミュレーションプログラムとなる。この場合、上記固定渦要素16の配置作業において、解析領域H内に配置される物体の形状に沿って固定渦要素16を配置するように構成すれば簡単に対応できる。
【0072】
(8)上記実施形態では、気泡の運動及び気泡体積率の計算を行った後に、各放出渦要素の循環強度変化及び位置の算出を行う構成としたが、これに限られない。即ち、例えば、図1中の情報記録装置6に、渦要素情報保存領域と、気泡情報保存領域と、ポテンシャルデータ保存領域とを確保する。そして、渦要素情報保存領域に、上記ステップS6で算出される計算格子10内における循環強度変化データと、上記ステップS7で算出される前記各離散要素16,17の位置及び循環強度データと、ステップS8で算出される格子点11における渦度データとを保存する。気泡情報保存領域には、ステップS3で算出される各気泡Gの位置データと、ステップS4で算出される気泡体積率データを保存する。ポテンシャルデータ保存領域に、ステップS9,10で算出される液体のベクトルポテンシャルデータと、液体のスカラーポテンシャルデータを保存する。そして、気泡Gに関する計算処理(ステップS3,4)と、離散渦要素16,17に関する計算処理(ステップS7)とを並列的に処理する構成であってもよい。
【図面の簡単な説明】
【0073】
【図1】本発明の解析要素(計算格子)を説明するための図
【図2】本発明の一実施形態に係る気泡流シミュレーションプログラムが導入されるコンピュータ及び周辺機器を示す図
【図3】本発明の計算の対象とする気泡流を説明するための図
【図4】速度せん断層を誘起する平板および箱体の固定渦要素による表現を説明するための図
【図5】コンピュータによる処理内容を示すフローチャート
【図6】本実施形態による解析結果を示す図
【符号の説明】
【0074】
1…コンピュータ
2…中央処理部
6…情報記録装置(記憶手段)
10…計算格子(解析要素)
11…格子点(評価点)
15…速度せん断層
16…固定渦要素(離散渦要素)
17…放出渦要素(離散渦要素)
G…気泡
H…解析領域
L…液体
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5