(57)【要約】【課題】  本発明に関する課題としては、数値不安定性を回避し多成分非圧縮性流体と移動壁とを取り扱う数値流体計算において、安定した時間発展を行うことにある。【解決手段】  本発明によると、多成分非圧縮性流体と移動壁とを取り扱う数値流体計算において、ナビエ・ストークス方程式と体積占有率関数の移流方程式とに同一の移流速度を用いて時間発展を行い、かつ質量保存則、非圧縮拘束条件を同時に満足させる数値流体計算方法が提供される。【選択図】  図３
2010225088.tif 000002
【特許請求の範囲】
【請求項１】
少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算において、ナビエ・ストークス方程式と体積占有率関数の移流方程式との時間発展を同一の移流速度を用いて行い、かつ質量保存則と非圧縮拘束条件とを同時に満足させる多成分非圧縮性流体計算方法。
【請求項２】
請求項１の非圧縮性流体計算方法において最初に圧力計算を行い、該圧力を用いた運動量の修正と流速の決定を行った後にナビエ・ストークス方程式と体積占有率関数の移流方程式とに含まれる移流項の計算を行うことにより両式に含まれる移流項の計算が同一の速度で行われる請求項１に記載の多成分非圧縮性流体計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算において、安定した時間発展を行うための多成分非圧縮性流体計算方法に関する。
【背景技術】
【０００２】
数値流体計算における支配方程式はナビエ・ストークス方程式である。
【０００３】
この微分方程式を初期値、境界条件を設定して解析的、数値的に解くことにより、流体の運動を予測することが可能となる。
【０００４】
古くは様々な仮定、近似を導入して解析的な解を求めることが主流であった。
【０００５】
しかし近年は、コンピュータの発展により流体の支配方程式であるナビエ・ストークス方程式を精度良く数値的に解くことが可能となり、産業界では様々な流体現象の予測技術として数値流体計算が用いられている。
【０００６】
例えば航空機翼の最大揚力を数値流体計算により推定し、風洞実験の回数を削減することにより航空機の開発期間、開発コストの削減に成功している(特許文献1参照)。
【０００７】
また、自動車ボディーの流動抵抗を数値流体計算を用いて予測することにより、実際の実験を行うことなく車体の最適形状の決定を行っている(特許文献2参照)。
【０００８】
他にも地球シミュレータのような大規模計算機を用いた気象予測が行われている。
【０００９】
多成分流体を対象とした数値流体計算は製造プロセスにおける基本原理解明のツールとしても期待されており、特殊ガラス、特殊合金、多孔質プラスチック等の製造プロセスの解析に用いられている(特許文献3参照)。
【００１０】
半導体開発においては、固体表面に均質な液体膜を形成する技術開発の一助として移動境界(自由表面)を含む数値流体計算を用いた現象解析が行われている(特許文献4、特許文献5参照)。
【００１１】
このように数値流体計算技術はその進歩により様々な現象を取り扱えるよう発展してきているが、未だ数値流体計算における問題が残されている。
【００１２】
例えば、気体と液体のように密度比の大きな複数成分の流体を対象とした数値流体計算には数値的な不安定が伴い、計算が困難であることが知られている。
【００１３】
不安定性が誘起される原因の一つは、ナビエ・ストークス方程式と体積占有率関数の移流方程式との時間発展において、両式に含まれる移流項の計算に異なる移流速度が用いられることである。
【００１４】
他にも、このような系においては質量保存則と非圧縮拘束条件とを同時にかつ厳密に満たすことが数値計算上の誤差により一般に困難であり、これが数値不安定性の原因となり得ることが知られている。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】特開2004-12248号公報
【特許文献２】特開2003-178099号公報
【特許文献３】特開2001-9974号公報
【特許文献４】特開平10-015479号公報
【特許文献５】特開平10-185755号公報
【発明の概要】
【発明が解決しようとする課題】
【００１６】
本発明に関する課題としては、前述のような数値不安定性を回避し少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算において、安定した時間発展を行うことにある。
【００１７】
少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算において、ナビエ・ストークス方程式の移流項と体積占有率関数の移流方程式の移流項とが異なる移流速度を用いて計算されることに起因する数値不安定性が存在する。
【００１８】
また、付帯される非圧縮拘束条件が数値誤差により厳密に満たされないことに起因する数値不安定性も存在する。
【課題を解決するための手段】
【００１９】
上記課題を解決するために本発明は、少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算において、ナビエ・ストークス方程式と体積占有率関数の移流方程式との時間発展を同一の移流速度を用いて行い、かつ質量保存則と非圧縮拘束条件とを同時に満足させる多成分非圧縮性流体計算方法である。
【発明の効果】
【００２０】
本発明によると、少なくとも2つ以上の非圧縮性流体と移動壁とを取り扱う数値流体計算におけるナビエ・ストークス方程式と体積占有率関数の移流方程式との時間発展において両式の移流項の計算に同一の移流速度を用いることにより、質量保存則と非圧縮拘束条件とを同時に満足させる数値的に安定な数値流体計算方法が提供される。
【図面の簡単な説明】
【００２１】
【図１】本発明の数値流体計算方法を搭載した装置の説明図１。
【図２】本発明の数値流体計算方法を搭載した装置の説明図２。
【図３】本発明による多成分数値流体計算方法の流れ図 (第一の実施例)。
【図４】1次元定常流れ場(時刻ｔ)。
【図５】1次元定常流れ場(時刻ｔ＋Δｔ)。
【図６】移動壁と空隙率(時刻ｔ)。
【図７】移動壁と空隙率(時刻ｔ＋Δｔ)。
【図８】本発明による多成分数値流体計算方法の流れ図(第二の実施例)。
【発明を実施するための形態】
【００２２】
本発明に示すアルゴリズムは図１及び図２に示す計算機内のハードディスクに収められている。
【００２３】
ここでいう計算機とは演算処理部、記憶装置部、入力部、出力表示部からなるシステムである。
【００２４】
演算処理部は図１及び図２に示すようにCPU1やメモリ2等からなる。
【００２５】
記憶装置部はハードディスク3やフレキシブルディスク4からなる。
【００２６】
入力部はキーボード5やマウス6等からなる。
【００２７】
そして出力表示部はディスプレイ7等からなる。
【００２８】
各演算部等は、本体8に収まっている。
【００２９】
計算機はメモリ上に図3に示したフローチャートに沿って設計、コード化されたプログラムをロードし、必要な計算領域を確保する。
【００３０】
その後、適当な方法で入力された画像やパラメータに基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク等に書き込んで記憶させ、ディスプレイ上に表示させる。
【実施例１】
【００３１】
先述のとおり、ナビエ・ストークス方程式の移流項と体積占有率関数の移流方程式の移流項とが同一の移流速度を用いて計算されることは本実施例の特徴の1つである。
【００３２】
この2つの移流速度が同一でない場合には、運動量と流速との関係に不合理が生じ、流速に物理的に不合理な加減速が生じる。
【００３３】
以下、この不合理について示す。
【００３４】
例として、図４、５に示すように胞に流体1が流入し、胞から流体2が流出する1次元問題を考える。
【００３５】
本例において流体は非粘性であり外力は働かないものと仮定している。
【００３６】
また流路面積が一定であることと非圧縮拘束条件とを仮定すると、時刻ｔにおける流体1と流体2の流速ｕ1（ｔ）、ｕ2（ｔ）は同一となる。
【００３７】
この速度をｃ（ｔ）とする。
【００３８】
理論上は流体1、流体2の体積占有率関数の移流速度ｕF1（ｔ）、ｕF2（ｔ）運動量の移流速度ｕJ1（ｔ）、ｕJ2（ｔ）はすべて同一の値となり、［数１］が成立する。
【００３９】
【数１】
2010225088.tif 000003
【００４０】
さらに、運動量と速度との関係から［数２］、［数３］が成立することは明らかである。
【００４１】
ここでＪ1（ｔ）、Ｊ2（ｔ）は該胞内の流体1、流体2の運動量ρ1、ρ2は流体1、流体2の密度とする。
【００４２】
Ｆ1（ｔ）、Ｆ2（ｔ）は流体1、流体2の体積占有率関数、Ｖは胞の体積である。
【００４３】
【数２】
2010225088.tif 000004
【００４４】
【数３】
2010225088.tif 000005
【００４５】
Δｔ秒後の該胞における体積占有率関数Ｆ1（ｔ＋Δｔ）、Ｆ2（ｔ＋Δｔ）および胞内の流体1および流体2の運動量Ｊ1（ｔ＋Δｔ）、Ｊ2（ｔ＋Δｔ）、は［数４］〜［数７］で与えられる。
【００４６】
Ｓは流入出境界の面積とする。
【００４７】
【数４】
2010225088.tif 000006
【００４８】
【数５】
2010225088.tif 000007
【００４９】
【数６】
2010225088.tif 000008
【００５０】
【数７】
2010225088.tif 000009
【００５１】
この時［数１］の関係が保持されているならば、Δｔ秒後の状態量に対して、［数８］、［数９］が成立することは明らかであり、流速ｕ1（ｔ＋Δｔ）とｕ2（ｔ＋Δｔ）は一致する。
【００５２】
これより、仮定した非圧縮拘束条件が保持されていることが分かる。
【００５３】
【数８】
2010225088.tif 000010
【００５４】
【数９】
2010225088.tif 000011
【００５５】
しかしながら、数値計算上の理由等によりｕF1（ｔ）とｕJ1（ｔ）あるいはｕF2（ｔ）とｕJ2（ｔ）とが一致しない場合には、流速ｕ1（ｔ＋Δｔ）とｕ2（ｔ＋Δｔ）とが一致することは保証されない。
【００５６】
また、それらが一致しない場合には仮定された非圧縮拘束条件に反する結果となる。
【００５７】
このように2つの移流方程式に用いる移流速度が同一で無いことは計算結果に不合理を生じさせる。
【００５８】
そのため、2つの移流方程式に同一の移流速度を用いることは重要である。
【００５９】
上述のとおり体積占有率関数の移流方程式と運動量の移流方程式とに同一の移流速度を用いることが重要である。
【００６０】
それに加えて、この移流速度を用いた移流方程式によって時間発展された体積占有率関数は移動壁を考慮した非圧縮拘束条件を満足する必要がある。
【００６１】
例えば図６、７に示すように、ある時刻ｔにおける空隙率がＨ（ｔ）である胞の空隙率が、移動壁の影響により時刻ｔ＋ΔｔにおいてＨ（ｔ＋Δｔ）に変化した場合を考える。
【００６２】
このような場合には、時刻ｔにおいて非圧縮拘束条件に関する方程式［数１０］を満たしていた体積占有率関数を時刻ｔ＋Δｔにおける非圧縮拘束条件に関する方程式［数１１］を満たすように体積占有率関数の時間発展を行う必要がある。
【００６３】
【数１０】
2010225088.tif 000012
【００６４】
【数１１】
2010225088.tif 000013
【００６５】
以下、ナビエ・ストークス方程式に含まれる移流項の計算と体積占有率関数の移流方程式の移流項とに同一の移流速度を用いるための理論を示す。
【００６６】
多成分非圧縮性流体計算における支配方程式は［数１２］の体積占有率関数の移流方程式、［数１３］のナビエ・ストークス方程式である。
【００６７】
［数１４］は非圧縮拘束条件に関する方程式である。
【００６８】
［数１３］において、対象とする系によっては表面張力項を省略することも可能である。
【００６９】
【数１２】
2010225088.tif 000014
【００７０】
【数１３】
2010225088.tif 000015
【００７１】
【数１４】
2010225088.tif 000016
【００７２】
【数１５】
2010225088.tif 000017
【００７３】
【数１６】
2010225088.tif 000018
【００７４】
【数１７】
2010225088.tif 000019
【００７５】
これらの支配方程式を時間方向に離散化することにより、［数２１］〜［数２４］が導かれる。
【００７６】
［数２１］〜［数２４］においてhatが上に付された変数は次ステップの量であることを表し、付されないものは現ステップの量であることを表している。
【００７７】
【数１８】
2010225088.tif 000020
【００７８】
は［数２２］に従ってナビエ・ストークス方程式の圧力項部分のみを用いて運動量を更新することにより得た運動量の予測子である。
【００７９】
［数２１］で用いられる
【００８０】
【数１９】
2010225088.tif 000021
【００８１】
は［数２２］で求められた
【００８２】
【数２０】
2010225088.tif 000022
【００８３】
を用いて［数２５］で定義される。
【００８４】
［数２２］と［数２３］との両式を実行することによりナビエ・ストークス方程式を用いた時間更新が行われる。
【００８５】
【数２１】
2010225088.tif 000023
【００８６】
【数２２】
2010225088.tif 000024
【００８７】
【数２３】
2010225088.tif 000025
【００８８】
【数２４】
2010225088.tif 000026
【００８９】
【数２５】
2010225088.tif 000027
【００９０】
圧力に対する方程式は［数２１］〜［数２５］を用いることにより導かれる。
【００９１】
以下、その手順を示す。
【００９２】
［数２１］をすべての成分に対して足し合わせることにより［数２６］が導かれる。
【００９３】
【数２６】
2010225088.tif 000028
【００９４】
ここで、［数２４］の関係を代入することにより［数２７］が導かれる。
【００９５】
また、移項操作により［数２８］を得る
【００９６】
【数２７】
2010225088.tif 000029
【００９７】
【数２８】
2010225088.tif 000030
【００９８】
［数２２］［数２５］［数２８］より、［数２９］を得る。
【００９９】
また、展開して［数３０］が得られる。
【０１００】
【数２９】
2010225088.tif 000031
【０１０１】
【数３０】
2010225088.tif 000032
【０１０２】
以上に示した［数１２］〜［数１７］、［数２１］〜［数３０］を用いて数値流体計算を行うための具体的な手順を以下に示す。
【０１０３】
以下の手順においては空間的な離散化が必要となる。
【０１０４】
その方法は有限差分法、有限要素法、メッシュフリー法等、手法は問わない。
【０１０５】
また、以上に示した［数１２］〜［数１７］、［数２１］〜［数３０］においては時間的な離散化にオイラー陽解法を用いることを前提としているが、他の時間的な離散化手法を適用することも可能である。
【０１０６】
第1のステップでは、初期値の設定を行う。
【０１０７】
ここでいう初期値の設定とは格子の初期形状と該格子における各胞の空隙率、流体に関連した物理定数、流体に関連する状態変数(流速、圧力)、体積占有率関数等を含む数値計算に必要なすべての変数、定数を設定することを指すものとする。
【０１０８】
任意の初期値を設定することが可能であるが、付帯される拘束条件式を満足する初期値を設定する必要がある。
【０１０９】
また、初期ステップにおいて急激な流体の運動を誘発するような初期値の設定は避けることが望ましい。
【０１１０】
そのような初期値を設定することも可能であるが、時間刻みを小さく設定する等により計算の安定性を確保することが必要となる。
【０１１１】
第2のステップでは、圧力の算出を行う。
【０１１２】
圧力の算出は［数３０］を用いて行う。
【０１１３】
［数３０］右辺第2項は理論上は零となる。
【０１１４】
そのため、［数３０］右辺第2項の計算を省略することも可能である。
【０１１５】
しかしながら実際には数値計算上の誤差によって一般に非零の値となるため、［数３０］右辺第2項の計算を行うことが望ましい。
【０１１６】
この計算により、数値誤差を補正する効果が得られ計算をより安定に行うことが可能となる。
【０１１７】
第3のステップでは、第2のステップで求まった圧力を用いて運動量の修正を行う。
【０１１８】
運動量の修正は［数２２］を用いて行う。
【０１１９】
運動量と速度との間には［数１５］の関係があり、運動量の修正を行うことは速度の修正を行うことに等価である。
【０１２０】
第4のステップでは、各成分の体積占有率関数の更新を行う。
【０１２１】
体積占有率関数の更新は［数２１］を用いて行われる。
【０１２２】
［数２１］は典型的な移流方程式であり、移流速度には第3のステップで修正された運動量を用いて［数２５］の関係から求められた値を用いる。
【０１２３】
これにより更新された体積占有率関数が非圧縮拘束条件に関する方程式［数２４］を満たすことが理論上は保証される。
【０１２４】
ただし、実際には数値計算上の誤差により厳密に満たされることは期待できない。
【０１２５】
第5のステップでは、各成分の運動量の更新を行う。
【０１２６】
運動量の更新は［数２３］を用いて行われる。
【０１２７】
［数２３］は移流拡散方程式であり、移流速度には第3のステップで修正された運動量を用いて［数２５］の関係から求められた値を用いる。
【０１２８】
これは第4のステップで体積占有率関数の更新を行う時に用いられた移流速度と同一であり、この点は本発明の特徴の1つである。
【０１２９】
第6のステップでは、数値計算の終了条件を満足しているか否かの判定を行う。
【０１３０】
満足していない場合には第2のステップに戻る。
【０１３１】
満足している場合には計算を終了する。
【０１３２】
上記に示した通り、第4のステップの体積占有率関数の更新と第5のステップの運動量の更新において両者の更新に用いる移流速度は同一のものであり、かつ質量保存則、非圧縮拘束条件を理論的に満足させるものである。
【０１３３】
また、この速度を決定するために必要な［数３０］により算出される圧力は［数３０］右辺第2項により、数値計算上の誤差により発生する質量保存則の不満足を修正する効果を含んでいる。
【０１３４】
そのため、拘束条件に対する補正効果が期待できる。
【０１３５】
このように付帯される拘束条件を満足させかつ、すべての状態変数の移流に同一の移流速度を用いることにより安定した計算を行うことを可能にすることが本発明の特徴である。
【実施例２】
【０１３６】
以下、本発明による多成分数値流体計算方法のその他の実施例を図8に示すフローチャートに沿って具体的に述べる。
【０１３７】
以後は、実施例1と同一の点については説明を省略する。
【０１３８】
第7のステップから第9のステップは実施例1における第1のステップから第3のステップに同じである。
【０１３９】
第10のステップにおいては、第9のステップで計算された修正された運動量
【０１４０】
【数３１】
2010225088.tif 000033
【０１４１】
を別のメモリ上に退避する。
【０１４２】
この変数を
【０１４３】
【数３２】
2010225088.tif 000034
【０１４４】
とする。
【０１４５】
第11のステップにおいては、運動量の更新を行う。
【０１４６】
本ステップは実施例1における第5のステップに相当し、［数２３］を用いて更新が行われる。
【０１４７】
第12のステップにおいては、体積占有率関数の更新を行う。
【０１４８】
本ステップは実施例1における第4のステップに相当し、［数２１］を用いて更新が行われる。
【０１４９】
この時、第9のステップで求めた運動量
【０１５０】
【数３３】
2010225088.tif 000035
【０１５１】
は第11のステップにおいて
【０１５２】
【数３４】
2010225088.tif 000036
【０１５３】
へと更新されてしまっているため、これを用いて計算された速度を移流速度として用いると第11のステップにおける運動量の移流速度と異なる移流速度を用いる結果となってしまう。
【０１５４】
そのため、第10のステップで別のメモリ上に保存した運動量
【０１５５】
【数３５】
2010225088.tif 000037
【０１５６】
を用いて［数３７］により求めた移流速度
【０１５７】
【数３６】
2010225088.tif 000038
【０１５８】
を用いることにより体積占有率関数と運動量との移流速度を同一とすることが可能となる。
【０１５９】
【数３７】
2010225088.tif 000039
【０１６０】
第13のステップは第一の実施例の第6のステップと同一であり計算の終了判定を行う。
【０１６１】
本実施例においては、第10のステップにおいては流速
【０１６２】
【数３８】
2010225088.tif 000040
【０１６３】
を
【０１６４】
【数３９】
2010225088.tif 000041
【０１６５】
として保存し第12のステップにおける体積占有率関数の移流速度として
【０１６６】
【数４０】
2010225088.tif 000042
【０１６７】
を用いることも可能である。
【産業上の利用可能性】
【０１６８】
本発明により、多成分非圧縮性流体と移動壁を取り扱う数値流体計算において、数値不安定性を回避して安定した時間発展が可能となる。
【符号の説明】
【０１６９】
１  CPU
２  メモリ
３  ハードディスク
４  フレキシブルディスク
５  キーボード
６  マウス
７  ディスプレイ
８  本体
９  フレキシブルディスクドライブ
１０  流体1の領域
１１  流体2の領域
１２  流体1の領域
１３  流体2の領域
１４  移動壁の領域
１５  胞の空隙
１６  移動壁の領域
１７  胞の空隙
Step １  初期値の設定を行うステップ
Step ２  圧力を計算するステップ
Step ３  運動量を修正するステップ
Step ４  体積占有率関数を更新するステップ
Step ５  運動量を更新するステップ
Step ６  終了判定を行うステップ
Step ７  初期値の設定を行うステップ
Step ８  圧力を計算するステップ
Step ９  運動量を修正するステップ
Step １０  運動量をメモリ上に退避するステップ
Step １１  運動量を更新するステップ
Step １２  体積占有率関数を更新するステップ
Step １３  終了判定を行うステップ
【図１】
2010225088.tif 000043
【図２】
2010225088.tif 000044
【図３】
2010225088.tif 000045
【図４】
2010225088.tif 000046
【図５】
2010225088.tif 000047
【図６】
2010225088.tif 000048
【図７】
2010225088.tif 000049
【図８】
2010225088.tif 000050