(57)【要約】【課題】圧縮性流体の２相分離流の体積率非導入計算方法において高精度な数値計算を行う。【解決手段】圧縮性流体の２相分離流の数値計算方法であって、コンピュータが、体積率を用いない１流体モデルにおいて、各相の混合の程度を表わす変数のラグランジュ微分を計算する第１の工程と、各相の混合の程度について仮定されたある混合則から仮想的な各相の体積率および密度を含む物理量を計算する第２の工程と、各相の圧力が平衡であること及び等エントロピー流れであることを仮定して前記第１の工程における計算結果を補正する第３の工程と、を有することを特徴とする圧縮性流体の２相分離流の数値計算方法。【選択図】図２
2014021904.tif 000002
【特許請求の範囲】
【請求項１】
圧縮性流体の２相分離流の数値計算方法であって、コンピュータが、
体積率を用いない１流体モデルにおいて、各相の混合の程度を表わす変数のラグランジュ微分を計算する第１の工程と、
各相の混合の程度についてある混合則を仮定し、それに基づき仮想的な各相の体積率および密度を含む物理量を計算する第２の工程と、
各相の圧力が平衡であること及びエントロピー変化に基づき前記第１の工程における計算結果を補正する第３の工程と、
を有することを特徴とする圧縮性流体の２相分離流の数値計算方法。
【請求項２】
コンピュータに、請求項１に記載の数値計算方法を実行させるプログラム。
【請求項３】
請求項２に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項４】
圧縮性流体の２相分離流を数値計算する計算装置であって、
体積率を用いない１流体モデルにおいて、各相の混合の程度を表わす変数のラグランジュ微分を計算する第１の工程と、
各相の混合の程度についてある混合則を仮定し、それに基づき仮想的な各相の体積率および密度を含む物理量を計算する第２の工程と、
各相の圧力が平衡であること及びエントロピー変化に基づき前記第１の工程における計算結果を補正する第３の工程と、
を実行するコンピュータを備える計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は数値計算方法、プログラム、記録媒体、及び計算装置に関する。
【背景技術】
【０００２】
圧縮性流体の２相流は工学の広範な分野において課題となる重要な流動現象である。２相流の形態としては、気液、固気および固液混相流がその代表であるが、中でも気液混相流は産業上広範に見られる現象であり特に重要である。例えば、翼近傍に発生するキャビテーション現象、気泡崩壊現象を用いた医療技術、インクジェットの吐出過程における発泡現象などが挙げられ、製品性能に影響を与える場合も多い。このような状況に対し、近年性能向上が著しい計算機を用いた数値流体計算の活用が盛んに行われている（非特許文献１参照）。
【０００３】
気液２相流の数値計算では、空間サイズが細かい気泡もしくは液滴の全てを解像しながら計算することは膨大な計算資源を要するために難しい。その回避策として、或る領域内に存在する液体もしくは気体を空間平均して体積率として表現するという考え方がある。現在行われている数値計算の計算方法は、この体積率を支配方程式に導入する方法、しない方法の２種類に大別できる。
【０００４】
前者には例えばSaurelらによる方法（非特許文献２参照）が挙げられ、圧縮性流体計算の代表的なベンチマーク問題である衝撃波管問題や気泡への衝撃波干渉問題において計算精度の高さが示されている。
【０００５】
後者は所謂、体積率を導入しない１流体モデルと呼ばれる分類に属し、例えばγ-based
model（非特許文献３）が挙げられる。この方法は体積率を導入しない代わりに混合状態の物理量を計算対象とし、相の識別は状態方程式のパラメータ方法により行うため、前者よりも簡便なモデルとなる。そのため工学分野での応用が期待されるが、我々の検討では少なくとも低マッハ数において前者よりも計算精度が劣る場合があることがわかっている。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】L.Liu, J.Li and Z.Feng, A numerical method for simulation of attached cavitation flows, Int. J. Numer. Meth. Fluid 2006; 52; 639-658
【非特許文献２】R.Saurel, F.Petitpas and R.Abgrall, Modelling phase transition in metastable liquid: application to cavitating and flashing flows, J.Fluid. Mech. 607, 2008
【非特許文献３】K.M.Shyue, An Efficisent Shock-Capturing Algorithm for Compressible Multicomponent Problems, Journal of Computationla Physics, 142, 208-242, 1998
【発明の概要】
【発明が解決しようとする課題】
【０００７】
従来の体積率を用いない１流体モデルは混合状態の物理量を対象としていたため、各相の圧力平衡条件とエントロピー変化を適切に計算に取り込むことができなかった。
【０００８】
本発明は、圧縮性流体の２相分離流の体積率非導入計算方法において高精度な数値計算を行う方法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
本発明は、圧縮性流体の２相分離流の数値計算方法であって、コンピュータが、
体積率を用いない１流体モデルにおいて、各相の混合の程度を表わす変数のラグランジュ微分を計算する第１の工程と、
各相の混合の程度についてある混合則を仮定し、それに基づき仮想的な各相の体積率および密度を含む物理量を計算する第２の工程と、
各相の圧力が平衡であること及びエントロピー変化に基づき前記第１の工程における計算結果を補正する第３の工程と、
を有することを特徴とする圧縮性流体の２相分離流の数値計算方法である。
【００１０】
本発明は、圧縮性流体の２相分離流を数値計算する計算装置であって、
体積率を用いない１流体モデルにおいて、各相の混合の程度を表わす変数のラグランジュ微分を計算する第１の工程と、
各相の混合の程度についてある混合則を仮定し、それに基づき仮想的な各相の体積率および密度を含む物理量を計算する第２の工程と、
各相の圧力が平衡であること及びエントロピー変化に基づき前記第１の工程における計算結果を補正する第３の工程と、
を実行するコンピュータを備える計算装置である。
【発明の効果】
【００１１】
本発明によれば、圧縮性流体の２相分離流の体積率非導入計算方法において高精度な数値計算を行うことができる。
【図面の簡単な説明】
【００１２】
【図１】ハードウェア構成例を示すブロック図
【図２】実施例に係る２相分離流の計算方法における処理を示すフローチャート
【図３】１次元テスト問題の条件を説明するための図
【図４】１次元テスト問題の結果を示す図
【発明を実施するための形態】
【００１３】
図１〜３を参照して本発明の実施形態について詳細に説明する。
図１は、本発明の一実施形態に係る流体解析計算装置のハードウェア構成例を示すブロック図である。図１において、１０００は２相分離流の計算を実現するコンピュータである。ＲＯＭ１０３やハードディスク（ＨＤ）や着脱可能な可搬記憶媒体（例えば、フレキシブルディスク（ＦＤ））等の外部メモリ１１１には、各種プログラム等が記憶されている。各種プログラムとは、例えば、ＣＰＵ１０１の制御プログラムであるＢＩＯＳ（Basic Input/Output System）やオペレーティングシステム（ＯＳ）や、ＰＣの実行する各機
能を実現するためのプログラムである。
【００１４】
また、図１に示すコンピュータ１０００において、ＲＡＭ１０２は、ＣＰＵ１０１の主メモリやワークエリア等として機能する。尚、ＣＰＵ１０１は、処理の実行に際して必要なプログラム等をＲＡＭ１０２にロードして、プログラムを実行することで各種動作を実現する。
【００１５】
また、入力コントローラ（入力Ｃ）１０５は、キーボード（ＫＢ）１０９や不図示のマウス等のポインティングデバイスからの入力を制御する。さらに、ビデオコントローラ（ＶＣ）１０６は、ＣＲＴディスプレイ（ＣＲＴ）１１０等の表示器への表示を制御する。尚、表示器はＣＲＴに限らず、液晶ディスプレイであっても構わない。これらは必要に応じて管理者が使用するものであって、本実施形態には直接関係があるものではない。
【００１６】
メモリコントローラ（ＭＣ）１０７は、ハードディスク（ＨＤ）やフレキシブルディスク（ＦＤ）、或いはアダプタやスロットを介して接続されるコンパクトフラッシュ（登録商標）メモリ等の外部メモリ１１１へのアクセスを制御する。これらの外部メモリには、ブートプログラム、ブラウザソフトウエア、各種のアプリケーション、フォントデータ、ユーザファイル、編集ファイル、各種データ等が記憶される。
【００１７】
通信Ｉ／Ｆコントローラ（通信Ｉ／ＦＣ）１０８は、ネットワークを介して外部機器と接続・通信するものであり、ネットワークでの通信制御処理を実行する。例えば、ＴＣＰ／ＩＰを用いたインターネット通信や他装置との間でデータ送受信等が可能である。尚、ＣＰＵ１０１は、例えばＲＡＭ１０２内の表示情報用領域へアウトラインフォントの展開（ラスタライズ）処理を実行することにより、ＣＲＴ１１０上での表示を可能としている。また、ＣＰＵ１０１は、ＣＲＴ１１０上の不図示のマウスカーソル等のＧＵＩを表示させることで、ユーザがＣＰＵ１０１に指示を入力することを可能とする。本発明は、以下で説明する数値計算方法をコンピュータに実行させるプログラムや、そのプログラムを記録したコンピュータ読み取り可能な記録媒体や、数値計算方法を実行するコンピュータを備える計算装置を含む。
【００１８】
以下では体積率を用いない１流体モデルの具体的な計算式をγ-based modelを基にして説明する。
γ-based modelでは混合状態の物理量を計算の対象とし、質量保存、運動量保存、エネルギー保存式について下式を考える。
【数１】
2014021904.tif 000003
【００１９】
ここでρは密度、ｕは速度、ｐは圧力、ｅｔは全エネルギー、Ｉは単位行列を意味する。この基本式は不可逆過程を含まず、等エントロピー流れを対象としている。
【００２０】
なお、上式だけでは問題が閉じないため更に状態方程式が必要となる。
圧縮性流体の２相分離流の計算において一般的であるstiffened gas EOS(equation of state)は、
【数２】
2014021904.tif 000004
ここでγは比熱比、ｐ∞は物質の特性を表わすためのパラメータである。
【００２１】
γ-based modelでは相の識別のために状態方程式のパラメータの時間発展式、
【数３】
2014021904.tif 000005
を解く。式（５）および式（６）の変数の取り方は圧力場および速度場が空間一様な初期条件下で界面が移動したときに場に非物理な振動が生じないように決める。
【００２２】
この１／（γ―１）は２相の混合の程度を表す変数であり、２相の混合状態における値であるが、本提案では次式で表わされる規則によって相が混合されたとする混合則の仮定を行う。
【数４】
2014021904.tif 000006
ここでαは体積率、添字ｃは或る相ｃを意味する。
【００２３】
体積率の定義、
【数５】
2014021904.tif 000007
を用い、式（７）を整理すると例えば２相の場合には、添字１を相１、添字２を相２として、
【数６】
2014021904.tif 000008
と書くことができる。上式から明らかなように、
【数７】
2014021904.tif 000009
である。
【００２４】
一方、体積率を導入する方法では各相の速度に差がないとすると、各相の質量保存、運動量保存、エネルギー保存に下式を考えることができる（非特許文献２参照）。
【数８】
2014021904.tif 000010
【００２５】
ここでμは各相の圧力の緩和係数、添字Ｉは界面での値を意味する。なお、上式は各相の圧力に差がないとき、式（１）〜式（３）と同様に等エントロピー流れを対象としたものになる。
【００２６】
各相の圧力に差がない場合は、その値をｐとして次式、
【数９】
2014021904.tif 000011
が成り立つ。ここでｓｃは相ｃのエントロピーである。
【００２７】
式（１５）から明らかなように、例えば相１と相２に関しては、
【数１０】
2014021904.tif 000012
が成り立つが、等エントロピー流れを対象とするならば、
【数１１】
2014021904.tif 000013
と変形できる。上式に式（１２）を代入し、ｃｃを相ｃの音速と表わすと、
【数１２】
2014021904.tif 000014
という圧力平衡と等エントロピーを考慮した体積率の時間発展式が得られる。
【００２８】
体積率の時間発展は式（１１）の関係式があるから、
【数１３】
2014021904.tif 000015
として１／（γ―１）の時間微分に繋げることができる。
【００２９】
式（１９）の右辺は式（５）と比較すれば明らかなように本来γ-based modelには存在しない項である。この項は導出過程から判るように、混合則の仮定から求めた仮想的体積率を用いて、体積率を導入する方法の基本式で成り立つ圧力平衡条件とエントロピー変化の関係式を導入した結果として得られている。このことから、体積率を導入しない方法へ、より詳細な物理現象が反映されたモデルであると言える。
【００３０】
図２に本発明の一実施形態に係る２相分離流の計算方法における処理のフローチャートを示す。
【００３１】
まず、ステップＳ１において、ＣＰＵ１０１は、速度、圧力など場の変数とする物理量をメモリ上に確保する。
【００３２】
続いてステップＳ２において、ＣＰＵ１０１は、流路形状や初期条件等の計算条件の入力を受け付ける。
【００３３】
続いてステップＳ３において、ＣＰＵ１０１は、時間ループに入り、計算中の時間を更新する。
【００３４】
続いてステップＳ４において、ＣＰＵ１０１は、一般的なγ-based modelに該当する計算を行う。これにより、各相の混合程度を表わす変数の全微分（ラグランジュ微分）、つまり式（１９）の左辺に含まれる移流項が計算される（第１の工程）。
【００３５】
続いてステップＳ５において、ＣＰＵ１０１は、仮定された混合則である式（７）を用いて仮想的な各相の体積率および密度などの物理量を計算する（第２の工程）。
【００３６】
続いてステップＳ６において、ＣＰＵ１０１は、ステップＳ５で求めた体積率と物理量を用いて、各相の圧力が平衡であること及びエントロピー変化を考慮してステップＳ４の計算結果である各相の混合程度を表わす変数を補正する（第３の工程）。より具体的には式（１９）の右辺を計算し、それによってステップＳ４で求めていた式（１９）左辺を補正する。
【００３７】
最後にステップＳ７において、ステップＳ１で入力した所定の時間を計算したかを判断し、達していない場合は、ステップＳ３に戻って更に計算を行う。一方、達した場合には計算終了となる。
【００３８】
本モデルの精度を確認するためにテスト計算を実施した。
テスト計算は我々の検討によって通常のγ-based modelでは理論解から大きく差が生じてしまう１次元の気泡の膨張、収縮問題により行った。問題の設定を図３に示す。
【００３９】
計算空間の中心に気体１を配置し、その周囲に液体２、さらにその外側に気体３を設定する。初期圧力を中心気体は周囲よりも高く設定し、計算の進行に応じて気体が膨張する。或る大きさまで膨らむと収縮に移行し、もとの大きさに戻るというテスト問題である。この系は液体が非圧縮性流体と見なせる条件では、例えば中野らによる文献（第２３回数値流体力学シンポジウム、界面運動を含む圧縮性気液二相流の半陰解法）によれば、
【数１４】
2014021904.tif 000016
により表わされる。ここでＶは中心気体の体積、ｕＡは中心気体界面の速度、ρ（０）およびＶ（０）は初期密度および初期体積、ＸＡおよびＸＢは液体の中心側座標および外側座標、ｐ∞は外側気体の圧力で大気圧、ｐｂは中心気体圧力である。上式をニュートン法によって解き、理論解を得ることができる。
【００４０】
中心気体の圧力の時間変化について、理論解とγ-based modelおよび本提案を比較した結果を図４に示す。図４から明らかなように、本提案は中心気体の圧力についてγ-based
modelよりも精度が大きく向上し、ほぼ理論解を再現している。この結果からも本提案が高精度な計算方法であると言える。
【００４１】
なお、本明細書では簡単のために散逸のない基本式を用いて説明を行ったが、熱伝導などの散逸を導入する場合も上述の考え方に基づいて計算式を構成すればよい。同様に状態方程式についてもstiffened gas EOSを用いて説明したが、別の状態方程式を用いても上
述の考え方に基づいて計算式を構成すればよい。
【符号の説明】
【００４２】
１０１：ＣＰＵ
１０００：コンピュータ
【図１】
2014021904.tif 000017
【図２】
2014021904.tif 000018
【図３】
2014021904.tif 000019
【図４】
2014021904.tif 000020