(57)【要約】移流方程式を用いて、界面を形成する各流体の既知物性値に基づいて、流体間の界面近傍における第１の未知物性値を算出し、流体の支配方程式を用いて、算出された第１の未知物性値に基づいて、流体間の界面近傍における挙動をあらわす第１の未知流体挙動値を算出し、デジタイザ法を用いて、算出された第１の未知流体挙動値に基づいて、流体間の界面近傍における界面関数を算出し、平滑化手法を用いて、算出された界面関数を変換する。また、移流方程式を用いて、界面を形成する各流体の既知物性値と、変換された界面関数と、に基づいて、流体間の界面近傍における第２の未知物性値を算出し、流体の支配方程式を用いて、第２の未知物性値と、第１の未知流体挙動値と、に基づいて、流体間の界面近傍における挙動をあらわす第２の未知流体挙動値を算出する。
2004088563.tif 000002
【特許請求の範囲】
【請求項１】
移流方程式を用いて、界面を形成する各流体の既知物性値に基づいて、前記流体間の界面近傍における第１の未知物性値を算出する第１の未知物性値算出工程と、
流体の支配方程式を用いて、前記第１の未知物性値算出工程によって算出された第１の未知物性値に基づいて、前記流体間の界面近傍における挙動をあらわす第１の未知流体挙動値を算出する第１の未知流体挙動値算出工程と、
デジタイザ法を用いて、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出工程と、
平滑化手法を用いて、前記界面関数算出工程によって算出された界面関数を変換する界面関数変換工程と、
を含んだことを特徴とする流体解析方法。
【請求項２】
前記移流方程式を用いて、前記界面を形成する前記各流体の既知物性値と、前記界面関数変換工程によって変換された界面関数と、に基づいて、前記流体間の界面近傍における第２の未知物性値を算出する第２の未知物性値算出工程と、
前記流体の支配方程式を用いて、前記第２の未知物性値算出工程によって算出された第２の未知物性値と、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値と、に基づいて、前記流体間の界面近傍における挙動をあらわす第２の未知流体挙動値を算出する第２の未知流体挙動値算出工程と、
を含んだことを特徴とする請求の範囲第１項に記載の流体解析方法。
【請求項３】
移流方程式を用いて、界面を形成する各流体の既知物性値に基づいて、前記流体間の界面近傍における第１の未知物性値を算出する第１の未知物性値算出工程と、
流体の支配方程式を用いて、前記第１の未知物性値算出工程によって算出された第１の未知物性値と、前記いずれか一方の流体の表面張力係数と、前記界面の曲率をあらわす曲率関数と、に基づいて、前記流体間の界面近傍における挙動をあらわす第１の未知流体挙動値を算出する第１の未知流体挙動値算出工程と、
デジタイザ法を用いて、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出工程と、
平滑化手法を用いて、前記界面関数算出工程によって算出された界面関数を変換する界面関数変換工程と、
を含んだことを特徴とする流体解析方法。
【請求項４】
前記移流方程式を用いて、前記界面を形成する前記各流体の既知物性値と、前記界面関数変換工程によって変換された界面関数と、に基づいて、前記流体間の界面近傍における第２の未知物性値を算出する第２の未知物性値算出工程と、
前記流体の支配方程式を用いて、前記第２の未知物性値算出工程によって算出された第２の未知物性値と、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値と、前記いずれか一方の流体の表面張力係数と、前記変換された界面関数から得られる前記流体間の界面における曲率と、前記界面の曲率をあらわす曲率関数と、に基づいて、前記流体間の界面近傍における挙動をあらわす第２の未知流体挙動値を算出する第２の未知流体挙動値算出工程と、
を含んだことを特徴とする請求の範囲第３項に記載の流体解析方法。
【請求項５】
請求の範囲第１項〜第４項のいずれか一つに記載の流体解析方法を、コンピュータに実行させることを特徴とする流体解析プログラム。
【請求項６】
流体ごとに前記流体の性質をあらわす既知物性値を格納するデータ格納手段と、
解析対象となる流体を２つ選択し、この選択された流体の既知物性値を前記データ格納手段から抽出するデータ抽出手段と、
流体の支配方程式を用いて、前記流体の界面近傍における物性値に基づいて、前記流体間の界面近傍における挙動をあらわす流体挙動値を算出する流体挙動値算出手段と、
デジタイザ法を用いて、前記流体挙動値算出手段によって算出された流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出手段と、
平滑化手法を用いて、前記界面関数算出手段によって算出された界面関数を変換する界面関数変換手段と、
移流方程式を用いて、前記データ抽出手段によって抽出された流体の物性値と、前記界面関数変換手段によって変換された界面関数と、に基づいて、前記流体間の界面近傍における物性値を算出する物性値算出手段と、
を備えることを特徴とする流体解析装置。
【請求項７】
前記データ格納手段は、さらに、前記流体ごとに、表面張力係数を格納しており、
前記データ抽出手段は、さらに、前記流体選択手段によって選択された流体のうちいずれか一方の表面張力係数を、前記データ格納手段から抽出し、
前記流体挙動値算出手段は、前記流体の支配方程式を用いて、さらに、前記データ抽出手段から抽出されたいずれか一方の流体の表面張力係数と、前記界面関数変換手段によって変換された界面関数と、前記変換された界面関数から得られる前記流体間の界面における曲率と、に基づいて、前記流体間の界面近傍における挙動をあらわす流体挙動値を算出することを特徴とする請求の範囲第６項に記載の流体解析装置。
【発明の詳細な説明】
【技術分野】  この発明は、自由表面問題における数値解法において、界面捕捉法に基づいた自由表面流れの数値シミュレーションをおこなうための流体解析方法、流体解析プログラムおよび流体解析装置に関する。
【背景技術】  自由表面を持つ流れの問題は、より一般的な分類として移動境界問題と称されており、その移動境界面を扱う方法として、界面捕捉法と界面追跡法が開発されている。界面追跡法には、たとえば、１９７２年に、Ｃ．Ｗ．Ｈｉｒｔらによって開発されたＡＬＥ法や、１９７７年にＲ．Ｂｏｎｎｅｒｏｔらによって開発されたＳｐａｃｅ−Ｔｉｍｅ法が存在する。
この界面追跡法は、界面をより正確に表現できるため有効であるが、自由表面（界面）上の計算メッシュを流体粒子とともに移動させるため、移動境界を表現する場合には、自由表面（界面）が大きく移動すると、計算メッシュのゆがみが大きくなり、計算が困難になる。さらに、計算メッシュは毎時間ステップごとに変更しなければならないため、そのアルゴリズムが複雑となる（たとえば、数値流体力学編集委員会編、４移動境界流れ解析、財団法人東京大学出版会、１９９５年３月１０日参照。）。
一方、界面捕捉法は、計算メッシュを固定しており、移動境界を表現する場合には特別な処理を施す必要があるが、計算メッシュは毎時間ステップごとに変化しないため、アルゴリズムを単純化することができ、界面が大きく変動する現象に有利である。この界面捕捉法の代表的な解法として、たとえば、１９８１年にＣ．Ｗ．ＨｉｒｔらによってＶＯＦ（Ｖｏｌｕｍｅ  ｏｆ  ｆｌｕｉｄ）法が開発されている（たとえば、Ｃ．Ｗ．Ｈｉｒｔ  ａｎｄ  Ｂ．Ｄ．Ｎｉｃｈｏｌｓ、Ｖｏｌｕｍｅ  ｏｆ  Ｆｌｕｉｄ（ＶＯＦ）  Ｍｅｔｈｏｄ  ｆｏｒ  ｔｈｅ  Ｄｙｎａｍｉｃｓ  ｏｆ  Ｆｒｅｅ  Ｂｏｕｎｄａｒｉｅｓ、Ｊｏｕｒｎａｌ  ｏｆ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ、１９８１、３９、ｐｐ．２０１−ｐｐ．２０５参照。）。
このＶＯＦ法は、流体の支配方程式であるＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式（下記式（１）および下記式（２））において、界面関数φによって物性値（密度ρ、粘性係数μ）の異なる流体（流体Ａ、流体Ｂ）を定義し、界面関数φを未知量にもつ移流方程式（下記式（３））において界面関数φを定義して、定義された界面関数φを用いて上記物性値を、下記式（４）および下記式（５）によって算出し、その界面を確定するという界面捕捉法である。
2004088563.tif 000003
上記式（１）〜式（５）において、ｔは時間をあらわし、ｕｉはｘｉ空間方向（ｉ＝１，２のときは２次元、ｉ＝１，２，３のときは３次元）の速度成分、ｐは圧力、ｆｉは外力である。また、式（４）および式（５）のρＡは流体Ａの密度、ρＢは流体Ｂの密度、μＡは流体Ａの粘性係数、μＢは流体Ｂの粘性係数である。
なお、界面関数φは０≦φ≦１の値をとり、流体Ａおよび流体Ｂは界面関数φを用いて、第４０図に示すように、流体Ａの界面関数φをφ＝１、流体Ｂの界面関数φをφ＝０として定義し、界面関数φ＝０．５となる部分を、異なる流体Ａ、Ｂ間の界面として評価する。
ここで、流体Ａをφ＝１、流体Ｂをφ＝０とした移流方程式は、上記式（３）により、第４１図に示すような界面が不連続な関数を解くこととなる。これらＶＯＦ法をベースにした解法においては、界面関数をいかに精度よく求めるかが問題となる。
なお、第４１図および第４２図に示すような不連続な分布を持つ関数を用いて上記式（３）の移流方程式を数値計算によって正確に解析できる解法は存在せず、現時点で精度のよい安定化有限要素法（たとえば安定化気泡関数法）などの解法を用いた場合、第４３図に示す計算結果となる（たとえば、松本純一，梅津剛，川原睦人、浅水長波方程式に対する安定化有限要素法と安定化気泡関数法、応用力学論文集、土木学会、２００２年、５巻、２３５頁〜２４２頁参照。）。
また、ＶＯＦ法に基づく界面捕捉法の一つとして、２０００年に、Ｓ．Ａｌｉａｂａｄｉらによって、ＳＦＥ／ＩＣ法が開発されている。このＳＦＥ／ＩＣ法は、数値計算によって得たれた崩れた不連続な分布を持つ関数を、シャープな分布にリモデリングする方法である（たとえば、Ｓ．Ａｌｉａｂａｄｉ  ａｎｄ  Ｔ．Ｅ．Ｔｅｚｄｕｙａｒ    、Ｓｔａｂｉｌｉｚｅｄ−Ｆｉｎｉｔｅ−Ｅｌｅｍｅｎｔ／Ｉｎｔｅｒｆａｃｅ−Ｃａｐｔｕｒｉｎｇ  Ｔｅｃｈｎｉｑｕｅ  ｆｏｒ  Ｐａｒａｌｌｅｌ  Ｃｏｍｐｕｔａｔｉｏｎ  ｏｆ  Ｕｎｓｔｅａｄｙ  Ｆｌｏｗｓ  ｗｉｔｈ  Ｉｎｔｅｒｆａｃｅ、Ｃｏｍｐｕｔｅｒ  Ｍｅｔｈｏｄｓ  ｉｎ  Ａｐｐｌｉｅｄ  Ｍｅｃｈａｎｉｃｓ  ａｎｄ  Ｅｎｇｉｎｅｅｒｉｎｇ、２０００、１９０、ｐｐ．２４３−ｐｐ．２６１参照。）。
さらに、ＶＯＦ法に基づく界面捕捉法の一つとして、１９９４年に、Ｓｕｓｓｍａｎらによって、Ｌｅｖｅｌ  Ｓｅｔ法が開発されている。Ｌｅｖｅｌ  Ｓｅｔ法は、不連続な分布を持つ界面関数を解く代わりに、滑らかな分布をもつ距離関数を導入して、界面を評価する方法である。この解法は、滑らかな関数を計算するため、界面の評価を精度よくおこなうことが可能である（たとえば、Ｍ．Ｓｕｓｓｍａｎ，  Ｐ．Ｓｍｅｒｅｋａ，  ａｎｄ  Ｓ．Ｏｓｈｅｒ、Ａ  Ｌｅｖｅｌ  Ｓｅｔ  Ａｐｐｒｏａｃｈ  ｆｏｒ  Ｃｏｍｐｕｔｉｎｇ  Ｓｏｌｕｔｉｏｎｓ  ｔｏ  Ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ  Ｔｗｏ−Ｐｈａｓｅ  Ｆｌｏｗ、Ｊｏｕｒｎａｌ  ｏｆ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ、１９９４、１１４、ｐｐ．１４６−ｐｐ．１５９参照。）。
また、流体解析法の一つとして、ＣＩＰ（Ｃｕｂｉｃ−Ｉｎｔｅｒｐｏｌａｔｅｄ  Ｐｒｏｐａｇａｔｉｏｎ）法が開発されている。１９９３年に矢部らによって、このＣＩＰ法を用いた解法の中で提案されたデジタイザ法は、実際の不連続な界面関数を解く変わりにタンジェント変換を用いた関数を経由して不連続な分布を持つ関数を求めることができる（たとえば、Ｔ．Ｙａｂｅ  ａｎｄ  Ｆ．Ｘｉａｏ、Ｄｅｓｃｒｉｐｔｉｏｎ  ｏｆ  Ｃｏｍｐｌｅｘ  ａｎｄ  Ｓｈａｒｐ  Ｉｎｔｅｒｆａｃｅ  ｄｕｒｉｎｇ  Ｓｈｏｃｋ  Ｗａｖｅ  Ｉｎｔｅｒａｃｔｉｏｎ  ｗｉｔｈ  Ｌｉｑｕｉｄ  Ｄｒｏｐ、Ｊｏｕｒｎａｌ  ｏｆ  Ｐｈｙｓｉｃａｌ  Ｓｏｃｉｅｔｙ  ｏｆ  Ｊａｐａｎ、１９９３、６２、ｐｐ．２５３７−ｐｐ．２５４０参照。）。
しかしながら、第４３図（ａ）に示すような初期条件から、界面関数φを原点まわりに１周期計算させると、第４３図（ｂ）に示すように、界面関数φの分布が徐々にくずれてしまい、５周期後には、第４３図（ｃ）に示すように、不連続な界面関数φの分布となっている部分が細くなり、略円錐状となる。このように、現時点で有効であると考えられていた安定化有限要素法を用いても、時間変化があると、界面関数φの界面状態を維持することができず、界面関数φの不連続面の評価が不十分であり、界面関数φを精度よく算出することができなかった。
なお、通常の有限要素法を用いた場合には、第４４図（ａ）に示すように、１／４周期後には、界面関数φの分布が、第４４図（ｂ）に示す安定化有限要素法を用いた場合よりも、顕著に不安定となり、界面関数φの不連続面の評価がより一層不十分となり、界面関数φを精度よく算出することができなかった。
ここで、ＶＯＦ法を適用したときの界面関数φの不連続面での評価が不十分な場合の計算結果の影響について、第４５図に示すモデルを用いて説明する。このモデルは水中の空気が自由表面（水面）に向かって上昇するモデルである。この第４５図の解析モデルの計算結果は第４６図に示すような結果となる。すなわち、第４６図では、不連続な分布をもつ界面関数の精度が悪いために、計算が進むにしたがって、２つの流体、すなわち空気と水の界面がぼやけてしまい、正確に２つの流体の界面が計算できていないこととなる。
また、ＳＦＥ／ＩＣ法は、第４７図に示すように、崩れた界面関数φの分布からリモデリングを行うために、崩れた分布の影響を少しずつ受けてしまうため、計算を進めていくにつれて、第４７図（ａ）に示す初期条件から、第４７図（ｂ）および第４７図（ｃ）に示すように、分布形状が変形してしまうこととなり、界面関数φの形状を維持することができず、界面関数φを精度よく算出することができなかった。
また、Ｌｅｖｅｌ  Ｓｅｔ法では、距離関数の定義（｜▽Ｆ｜＝１）を満たす関数を反復計算によって求めるため、流体の挙動が複雑な場合、反復を繰り返しても計算結果が収束するとは限らない。また、解を算出することが可能な場合であっても、反復回数が非常に多くなった場合は大きな労力（計算コスト）を要する。
一方、上述したＣＩＰ法で使用されていたデジタイザ法を有限要素法の解法に取り入れると、第４８図（ａ）に示す初期条件の界面関数φは、第４８図（ｂ）に示すように１周期後であっても、第４８図（ｃ）に示すように５周期後であっても、界面関数φの形状は保たれており、不連続関数について正確な計算をおこなうことができる。
しかしながら、第４５図に示す解析モデルについて、上述したデジタイザ法を経由して得られた界面関数は、流体Ａと流体Ｂに対応している界面関数φの値が、１と０であるため、極端に不連続な関数となっている。有限要素法では三角形、四面体などの有限な形状で任意領域をメッシュにして分割しているために、この関数をそのまま界面として評価してしまうと、空気と水の２つの流体の界面上で、第４９図のようにガタガタな界面となってしまうこととなる。
第４９図に示すように、実際の２つの流体における界面を考えてみると、表面張力の影響によって、界面はある曲率を持ったなめらかな分布をしている。この表面張力により、極端に不連続な界面関数を使用した場合には、界面がガタガタなために、界面における表面張力の効果を考慮した評価がなされず、その結果、計算がかえって不安定になり、得られる結果も実際の現象を表現していない結果となる。また、この手法により反復計算を実行すると、計算不能に陥ることとなり、流体解析をおこなうことができなくなってしまうという問題があった。
この発明は、上述した従来技術による問題点を解消するため、簡便で汎用性のある手法により、不連続な分布となる界面関数を正確に算出し、さらに、物性値の異なる流体間における表面張力がもたらす界面上の滑らかさを界面関数に取り入れることができる流体解析方法、流体解析プログラムおよび流体解析装置を提供することを目的とする。
【発明の開示】
上述した課題を解決し、目的を達成するため、この発明にかかる流体解析方法は、移流方程式を用いて、界面を形成する各流体の既知物性値に基づいて、前記流体間の界面近傍における第１の未知物性値を算出する第１の未知物性値算出工程と、流体の支配方程式を用いて、前記第１の未知物性値算出工程によって算出された第１の未知物性値に基づいて、前記流体間の界面近傍における挙動をあらわす第１の未知流体挙動値を算出する第１の未知流体挙動値算出工程と、デジタイザ法を用いて、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出工程と、平滑化手法を用いて、前記界面関数算出工程によって算出された界面関数を変換する界面関数変換工程と、を含んだことを特徴とする。
この発明によれば、不連続な分布となる界面関数を、流体間の実際の界面状態に対応する界面関数に変換することができる。
また、上述した流体解析方法は、前記移流方程式を用いて、前記界面を形成する前記各流体の既知物性値と、前記界面関数変換工程によって変換された界面関数と、に基づいて、前記流体間の界面近傍における第２の未知物性値を算出する第２の未知物性値算出工程と、前記流体の支配方程式を用いて、前記第２の未知物性値算出工程によって算出された第２の未知物性値と、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値と、に基づいて、前記流体間の界面近傍における挙動をあらわす第２の未知流体挙動値を算出する第２の未知流体挙動値算出工程と、を含んだこととしてもよい。
この発明によれば、経時的に変化する流体間の界面近傍の挙動を再現することができる。
また、上述した流体解析方法は、移流方程式を用いて、界面を形成する各流体の既知物性値に基づいて、前記流体間の界面近傍における第１の未知物性値を算出する第１の未知物性値算出工程と、流体の支配方程式を用いて、前記第１の未知物性値算出工程によって算出された第１の未知物性値と、前記いずれか一方の流体の表面張力係数と、前記界面の曲率をあらわす曲率関数と、に基づいて、前記流体間の界面近傍における挙動をあらわす第１の未知流体挙動値を算出する第１の未知流体挙動値算出工程と、デジタイザ法を用いて、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出工程と、平滑化手法を用いて、前記界面関数算出工程によって算出された界面関数を変換する界面関数変換工程と、を含んだこととしてもよい。
この発明によれば、不連続な分布となる界面関数を、表面張力が影響する流体間の実際の界面状態に対応する界面関数に変換することができる。
また、上述した流体解析方法は、請求の範囲第３項に記載の発明において、前記移流方程式を用いて、前記界面を形成する前記各流体の既知物性値と、前記界面関数変換工程によって変換された界面関数と、に基づいて、前記流体間の界面近傍における第２の未知物性値を算出する第２の未知物性値算出工程と、前記流体の支配方程式を用いて、前記第２の未知物性値算出工程によって算出された第２の未知物性値と、前記第１の未知流体挙動値算出工程によって算出された第１の未知流体挙動値と、前記いずれか一方の流体の表面張力係数と、前記変換された界面関数から得られる前記流体間の界面における曲率と、前記界面の曲率をあらわす曲率関数と、に基づいて、前記流体間の界面近傍における挙動をあらわす第２の未知流体挙動値を算出する第２の未知流体挙動値算出工程と、を含んだこととしてもよい。
この発明によれば、経時的に変化する流体間の界面近傍における表面張力の影響を取り入れた挙動を再現することができる。
また、流体解析プログラムは、上記いずれか一つに記載の流体解析方法を、コンピュータに実行させることを特徴とする。
この発明によれば、上記いずれか一つに記載の流体解析方法を、コンピュータ処理によって実現することができる。
また、流体解析装置は、流体ごとに前記流体の性質をあらわす既知物性値を格納するデータ格納手段と、解析対象となる流体を２つ選択し、この選択された流体の既知物性値を前記データ格納手段から抽出するデータ抽出手段と、流体の支配方程式を用いて、前記流体の界面近傍における物性値に基づいて、前記流体間の界面近傍における挙動をあらわす流体挙動値を算出する流体挙動値算出手段と、デジタイザ法を用いて、前記流体挙動値算出手段によって算出された流体挙動値に基づいて、前記流体間の界面近傍における界面関数を算出する界面関数算出手段と、平滑化手法を用いて、前記界面関数算出手段によって算出された界面関数を変換する界面関数変換手段と、移流方程式を用いて、前記データ抽出手段によって抽出された流体の物性値と、前記界面関数変換手段によって変換された界面関数と、に基づいて、前記流体間の界面近傍における物性値を算出する物性値算出手段と、を備えることを特徴とする。
この発明によれば、不連続な分布となる界面関数を、流体間の実際の界面状態に対応する界面関数に変換することができる。
また、上述した流体解析装置は、前記データ格納手段は、さらに、前記流体ごとに、表面張力係数を格納しており、前記データ抽出手段は、さらに、前記流体選択手段によって選択された流体のうちいずれか一方の表面張力係数を、前記データ格納手段から抽出し、前記流体挙動値算出手段は、前記流体の支配方程式を用いて、さらに、前記データ抽出手段から抽出されたいずれか一方の流体の表面張力係数と、前記界面関数変換手段によって変換された界面関数と、前記変換された界面関数から得られる前記流体間の界面における曲率と、に基づいて、前記流体間の界面近傍における挙動をあらわす流体挙動値を算出することを特徴とする。
この発明によれば、不連続な分布となる界面関数を、表面張力が影響する流体間の実際の界面状態に対応する界面関数に変換することができる。
【図面の簡単な説明】  第１図は、この発明の本実施の形態にかかる流体解析装置のハードウェア構成を示すブロック図であり、第２図は、この発明の本実施の形態にかかる流体解析装置の機能構成を示すブロック図であり、第３図は、この発明の本実施の形態にかかる流体解析装置におけるデータ格納部に格納されている流体情報テーブルを示す説明図であり、第４図は、この発明の本実施の形態にかかる流体解析装置における前処理手順を示すフローチャートであり、第５図は、この発明の本実施の形態にかかる流体解析装置における流体解析処理手順を示すフローチャートであり、第６図は、この発明の本実施の形態にかかる流体解析装置における流体解析処理手順を示すフローチャートであり、第７図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第８図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第９図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１０図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１１図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１２図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１３図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１４図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１５図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１６図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１７図は、第４５図の解析モデルを例にした、ＶＯＦ法とＶＯＦ法＋デジタイザ法と本実施の形態にかかる流体解析との比較図であり、第１８図は、水中の気泡流れの解析モデルを示す説明図であり、第１９図は、第１８図の解析モデルにおいて、気泡形状（Ｓｐｈｅｒｉｃａｌ）を例にした場合の流体の挙動を示す説明図であり、第２０図は、第１８図の解析モデルにおいて、気泡形状（Ｓｋｉｒｔｅｄ）を例にした場合の流体の挙動を示す説明図であり、第２１図は、ランスから純酸素ガスを高速流入して鉄や銅などの金属を精錬する転炉内部における流体の挙動を解析するモデルを示す説明図であり、第２２図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２３図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２４図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２５図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトル���示す説明図であり、第２６図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２７図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２８図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第２９図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第３０図は、第２１図の解析モデルにおいて、転炉内部の流速ベクトルを示す説明図であり、第３１図は、ノズルからインクを噴射する流体の挙動を解析するモデルを示す説明図であり、第３２図は、第３１図の解析モデルにおいて、固定メッシュに分割したメッシュ図である。第３３図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３４図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３５図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３６図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３７図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３８図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第３９図は、第３１図の解析モデルにおいて、気液界面を示す説明図であり、第４０図は、流体および流体間における界面の解析モデルを示す説明図であり、第４１図は、第４０図の解析例となる関数を示す説明図であり、第４２図は、第４１図の関数をｘｙ座標に展開したグラフであり、第４３図は、第４０図の解析モデルにおいて、安定化気泡関数法を適用した場合の解析例となる関数を示す説明図であり、第４４図は、第４０図の解析モデルにおいて、有限要素法を適用した場合の解析例となる関数を示す説明図であり、第４５図は、水中の気泡の流況を示す解析モデルを示す説明図であり、第４６図は、第４５図の解析モデルにおいて、ＶＯＦ法＋安定化気泡関数法を適用した場合の気液界面と流速ベクトルを示す説明図であり、第４７図は、第４０図の解析モデルにおいて、ＳＦＥ／ＩＣ法を適用した場合の解析例となる関数を示す説明図であり、第４８図は、第４０図の解析モデルにおいて、ＶＯＦ法＋デジタイザ法を適用した場合の解析例となる関数を示す説明図であり、第４９図は、第４５図の解析モデルにおいて、ＶＯＦ法＋デジタイザ法を適用した場合の流速ベクトルを示す説明図である。
【発明を実施するための最良の形態】  以下に添付図面を参照して、この発明にかかる流体解析方法、流体解析プログラムおよび流体解析装置の好適な実施の形態を詳細に説明する。本実施の形態では、複雑な界面を扱うためにＶＯＦ法を採用し、精度の良い自由表面を計算するためにＣＩＰ法で提案・採用されているデジタイザ法を有限要素法（安定化気泡関数法）の解法に取り込んでいる。また、有限要素法でデジタイザ法を取り込むにあたって、界面の評価のしかたを有限要素法で計算をするうえで正確に表現できるような平滑化手法を採用している。
（流体解析装置のハードウェア構成）
まず、第１図を用いて、本実施の形態にかかる流体解析装置１のハードウェア構成について説明する。流体解析装置１は、ＣＰＵ１０１と、ＲＯＭ１０２とＲＡＭ１０３と、ＨＤＤ（ハードディスクドライブ）１０４と、ＨＤ（ハードディスク）１０５と、ＦＤＤ（フレキシブルディスクドライブ）１０６と、着脱可能な記録媒体の一例としてＦＤ（フレキシブルディスク）１０７と、ディスプレイ１０８と、Ｉ／Ｆ（インターフェース）１０９と、キーボード１１０と、マウス１１１と、スキャナ１１２と、プリンタ１１３と、を備えている。また、各構成部は、バス１００によってそれぞれ接続されている。
ここで、ＣＰＵ１０１は、流体解析装置１の全体の制御を司る。ＲＯＭ１０２は、ブートプログラムなどのプログラムを記憶している。ＲＡＭ１０３は、ＣＰＵ１０１のワークエリアとして使用される。ＨＤＤ１０４は、ＣＰＵ１０１の制御にしたがってＨＤ１０５に対するデータのリード／ライトを制御する。ＨＤ１０５は、ＨＤＤ１０４の制御で書き込まれたデータを記憶する。
ＦＤＤ１０６は、ＣＰＵ１０１の制御にしたがってＦＤ１０７に対するデータのリード／ライトを制御する。ＦＤ１０７は、ＦＤＤ１０６の制御で書き込まれたデータを記憶したり、ＦＤ１０７に記憶されたデータを流体解析装置１に読み取らせたりする。着脱可能な記録媒体として、ＦＤ１０７のほか、ＣＤ−ＲＯＭ（ＣＤ−Ｒ、ＣＤ−ＲＷ）、ＭＯ、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）、メモリーカードなどであってもよい。ディスプレイ１０８は、カーソル、アイコンあるいはツールボックスをはじめ、文書、画像、機能情報などのデータを表示する。このディスプレイ１０８は、たとえば、ＣＲＴ、ＴＦＴ液晶ディスプレイ、プラズマディスプレイ等を採用することができる。
Ｉ／Ｆ１０９は、通信回線を通じてインターネットなどのネットワークに接続され、このネットワークを介して他の装置に接続される。そして、Ｉ／Ｆ１０９は、ネットワークと内部のインターフェースを司り、外部装置からのデータの入出力を制御する。Ｉ／Ｆ１０９には、たとえばモデムやＬＡＮアダプタなどを採用することができる。
キーボード１１０は、文字、数字、各種指示などの入力のためのキーを備え、データの入力をおこなう。また、タッチパネル式の入力パッドやテンキーなどであってもよい。マウス１１１は、カーソルの移動や範囲選択、あるいはウィンドウの移動やサイズの変更などをおこなう。ポインティングデバイスとして同様に機能を備えるものであれば、トラックボールやジョイスティックなどであってもよい。
スキャナ１１２は、画像を光学的に読み取り、流体解析装置１内に画像データを取り込む。また、プリンタ１１３は、画像データや文書データを印刷する。プリンタ１１３には、たとえば、レーザプリンタやインクジェットプリンタを採用することができる。
（流体解析装置１の機能的構成）
つぎに、本実施の形態にかかる流体解析装置１の機能的構成について説明する。第２図に示すように、流体解析装置１は、データ格納部２と、データ抽出部３と、初期設定部４と、流体解析部５と、表示制御部６と、から構成されている。
データ格納部２には、第３図に示すように、流体ごとに、流体ＩＤ、流体名、流体物性値（密度、粘性係数、表面張力係数）、初期流体挙動値（解析開始流速および解析開始圧力）からなる流体情報テーブル１２が格納されている。また、データ格納部２には、図示はしないが、初期設定情報として、時間方向に対して、逐次解析計算をおこなうために時間増分量、計算ステップ数などの情報も格納されている。
さらに、解析対象範囲に分割形成される固定メッシュの分割幅、節点数およびメッシュ数などのメッシュ情報も格納されている。この固定メッシュは、解析対象範囲が２次元空間である場合は、三角形または格子状、解析対象範囲が３次元空間である場合は、四面体または立方体などの有限な形状によって構成される。メッシュ情報は、スキャナ１１２により読み取ることによっても取り込むことができる。このデータ格納部２は、具体的には、例えば第１図に示したＲＡＭ１０３、ＨＤ１０５またはＦＤ１０７によりその機能を実現する。
データ抽出部３は、解析対象となる流体の流体ＩＤを選択し、データ格納部２から、その選択された流体ＩＤに関連付けられている流体名、既知物性値（密度、粘性係数、表面張力係数）、初期流体挙動値（解析開始流速および解析開始圧力）を抽出する。このデータ抽出部３は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
初期設定部４は、データ抽出部３による各データの抽出後、選択された一方の流体に解析開始界面関数φ０（時刻ｔ＝０における界面関数φ）としてφ＝１、選択された他方の流体に解析開始界面関数φ０としてφ＝０を設定する。また、初期設定部４は、選択された２つの流体が気体と液体である場合、液体の表面張力係数を流体挙動値算出部８に出力する。この初期設定部４は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
流体解析部５は、物性値算出部７と、流体挙動値算出部８と、界面関数算出部９と、界面関数変換部１０と、から構成されている。流体解析部５は、データ格納部２に格納された固定メッシュごとに流体の挙動を解析する。この流体解析部５は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
物性値算出部７は、下記式（６）および下記式（７）の移流方程式を用いて、界面関数変換部１０によって変換された界面関数φと、データ抽出部３から抽出された流体の密度および粘性係数と、に基づいて、流体間の界面における未知密度ρおよび未知粘性係数μを算出する（以下、算出された未知密度ρを「算出密度ρ」と称し、算出された未知粘性係数μを「算出粘性係数μ」と称す。）。
2004088563.tif 000004
ここで、ρＡは流体Ａの密度、ρＢは流体Ｂの密度である。また、μＡは流体Ａの粘性係数、μＢは流体Ｂの粘性係数である。なお、この物性値算出部７は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
流体挙動値算出部８は、下記式（８）および下記式（９）に示すＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式を用いて、流体間の界面における算出密度および算出粘性係数に基づいて、流体間の界面における未知流速ｕ，ｖ，ｗを算出する。また、同時に未知圧力ｐも算出する（以下、算出された未知流速ｕ，ｖ，ｗを「算出流速ｕ，ｖ，ｗ」と称し、算出された未知圧力ｐを「算出圧力ｐ」と称す。）。
2004088563.tif 000005
ここで、ｉおよびｊは、１〜３から選ばれる自然数である。このｉおよびｊは、解析対象範囲の空間方向における次元をあらわす数字であり、ｕｉ＝ｕ１のときは、Ｘ軸方向における流体の算出流速ｕをあらわし、ｕｉ＝ｕ２のときは、Ｙ軸方向における流体の算出流速ｖをあらわし、ｕｉ＝ｕ３のときは、Ｚ軸方向における流体の算出流速ｗをあらわす。
また、ｆｉは外力をあらわし、流体解析をおこなう解析対象範囲の環境によって設定される。すなわち、外力ｆｉは、下記式（１０）によってあらわすことができる。
2004088563.tif 000006
ここで、外力ｆｇｉは、重力加速度ｇによる外力であり、２次元の場合は、下記式（１１）によりあらわされる。なお、３次元の場合は、下記式（１２）によってあらわされる。
2004088563.tif 000007
また、外力ｆｓｖｉは、１９９２年にＪ．Ｕ．Ｂｒａｃｋｂｉｌｌらによってモデル化（ＣＳＦモデル）された表面張力による外力であり、下記式（１３）によってあらわすことができる。
2004088563.tif 000008
ここで、κ（φ）は界面の曲率をあらわす曲率関数であり、下記式（１４）によってあらわされる。
2004088563.tif 000009
この流体挙動値算出部８は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
界面関数算出部９は、デジタイザ法を用いて、流体挙動値算出部８によって算出された算出流速に基づいて、解析対象範囲における界面関数φを算出する。具体的には、移流方程式にデジタイザ法を適用した処理であり、流体挙動値算出部８から入力されてくる算出流速ｕ，ｖ（２次元の場合。３次元の場合は、算出流速ｕ，ｖ，ｗ）を用い、下記式（１５）によって、デジタイザ処理をおこなう。このデジタイザ法は、実際の不連続な界面関数φを解く代わりにタンジェント変換を用いた関数Ｆ（φ）を経由して不連続な分布を持つ関数を求める方法である。
2004088563.tif 000010
すなわち、式（１５）を解くと、下記式（１６）のタンジェント変換式が導き出される。そして、この式（１６）を下記式（１７）に変換して、界面関数φを算出する。
2004088563.tif 000011
2004088563.tif 000012
なお、上記式（１７）のεは、たとえばε＝１０−５０程度の微小な値であり、式（１７）の無限大への発散を抑制する調整値である。また、この界面関数算出部９は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
界面関数変換部１０は、下記式（１８）を用いて、界面関数算出部９によって算出された界面関数φを、平滑化された界面関数φに変換して、界面の近傍における平滑化処理を実行する。具体的には、界面関数算出部９から算出された界面関数φを、下記式（１８）に示す拡散方程式を用いて変換する。
2004088563.tif 000013
ここで、ｓは仮想的な時間、Ｋは仮想的な拡散係数である。ｓとＫの選択の指標は以下の通りである。
（ｉ）拡散方程式はデジタイザ法によって算出された界面関数φによりギザギザとなった形状をスムーズにするために用いられる。
（ｉｉ）仮想的な時間ｓと仮想的な拡散係数Ｋは、界面関数φによりギザギザとなった形状が初期形状と同等の形状にするように選択し、この仮想的な時間ｓと仮想的な拡散係数Ｋの選択は、初期解析の１回のみ行われる。
（ｉｉｉ）ギザギザの度合いは分割幅に依存するが、仮想的な拡散係数Ｋは要素ごとに異なった値を指定でき、要素分割が非均質な場合も適用可能である。界面捕捉法は、一般的に時間進展に対して固定メッシュを使用するので、初期に選択された仮想的な時間ｓと仮想的な拡散係数Ｋは、以後の時間ステップに有効である。
また、平滑化した界面の具体的な表現方法としては、平滑化処理のときに（拡散方程式を解くときに）用いた補間関数を用いている。たとえば、拡散方程式を解くときに１次要素を用いた場合には、１次の補間によって界面の近似を行っている。なお、この界面関数変換部１０は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
また、表示制御部６は、流体解析部５における解析結果を第１図に示すディスプレイ１０８に出力する。この表示制御は、ＶＯＦ法による流体解析に一般的に用いられるためその説明を省略する。この表示制御部６は、具体的には、たとえば第１図に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また第１図に示したＩ／Ｆ１０９によってその機能を実現する。
つぎに、本実施の形態にかかる処理手順について、第４図〜第６図のフローチャートを用いて説明する。なお、ここでは流体Ａは気体、流体Ｂは液体とする。まず、第４図に示すように、入力装置（たとえば第１図に示すキーボード１１０またはマウス１１１）を介して、解析対象となる流体ＩＤが選択された場合（ステップＳ４０１：Ｙｅｓ）、データ格納部２から、選択された流体ＩＤを参照して、選択された流体Ａ，Ｂの流体物性値である密度ρＡ，ρＢ、粘性係数μＡ，μＢおよび表面張力係数σＡ，σＢと、選択された流体Ａ，Ｂの初期流体挙動値である解析開始流速ｕＡ，ｖＡ，ｗＡ，ｕＢ，ｖＢ，ｗＢおよび解析開始圧力ｐＡ，ｐＢを抽出する（ステップＳ４０２）。
このデータ抽出のあと、選択された一方の流体に解析開始界面関数φ０（界面関数φ＝１）、他方の流体に解析開始界面関数φ０（界面関数φ＝０）を設定する（ステップＳ４０３）。つぎに、液体である流体Ｂの表面張力を選択する（Ｓ４０４）。そして、入力装置から解析開始入力を受け付けた場合（ステップＳ４０５：Ｙｅｓ）、データ格納部２から、あらかじめ格納された解析対象範囲の固定メッシュを抽出する（ステップＳ４０６）。そして、式（１８）の拡散方程式の仮想的な時間ｓと仮想的な拡散係数Ｋの値が入力された場合（ステップＳ４０７：Ｙｅｓ）、抽出された固定メッシュごとに、流体解析処理を実行する（ステップＳ４０８）。
第５図および第６図は、流体解析処理の処理手順を示すフローチャートである。まず、第５図に示すように、データ抽出した各流体Ａ，Ｂの密度ρＡ，ρＢおよび粘性係数μＡ，μＢと、初期設定された解析開始界面関数φ０を、式（６）および式（７）の移流方程式に代入する（ステップＳ５０１）。そして、この移流方程式を解いて、流体Ａ，Ｂ間の界面における未知密度ρおよび未知粘性係数μを算出する（ステップＳ５０２）。
つぎに、算出密度ρおよび算出粘性係数μと、データ抽出された流体の解析開始流速ｕＡ，ｖＡ，ｗＡ，ｕＢ，ｖＢ，ｗＢ、解析開始圧力ｐＡ，ｐＢおよびステップＳ４０４において選択された表面張力係数σＢと、解析開始界面関数φ０とを、式（８）および式（９）のＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式（式（１０）〜式（１４）も含む）に代入する（ステップＳ５０３）。そして、このＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式を解いて、所定時間経過後の流体Ａ，Ｂ間の界面における未知流速ｕ，ｖ，ｗおよび未知圧力ｐを算出する（ステップＳ５０４）。
この算出された算出流速ｕ，ｖ，ｗを、上記式（１５）に代入し（ステップＳ５０５）、上記式（１５）を上記式（１６）に展開する（ステップＳ５０６）。この展開された上記式（１６）に、解析開始界面関数φ０（界面関数φ＝１または界面関数φ＝０）を代入して、タンジェント変換値Ｆ（φ）を算出する（ステップＳ５０７）。そして、上記式（１６）を式（１７）に変形し（ステップＳ５０８）、この変形した式（１７）に、算出されたＦ（φ）の値を代入して、界面関数φを算出する（ステップＳ５０９）。
つぎに、算出された界面関数φを、上記式（１８）の拡散方程式に代入して、この拡散方程式を解き、ステップＳ５０９で算出された界面関数φを、界面を平滑にする界面関数φに変換する（ステップＳ５１０）。
つぎに、変換された界面関数φと、各流体Ａ，Ｂの密度ρＡ，ρＢおよび粘性係数μＡ，μＢを、上記式（６）、式（７）の移流方程式に代入する（ステップＳ６０１）。そして、この移流方程式を解き、未知密度ρおよび未知粘性係数μを算出する（ステップＳ６０２）。
つぎに、ステップＳ６０２で算出された算出密度ρおよび算出粘性係数μと、ステップＳ５０４で算出された算出流速ｕ，ｖ，ｗおよび算出圧力ｐと、ステップＳ５０９で変換された界面関数φと、ステップＳ４０４で選択された表面張力係数σＢと、を上記式（８）および式（９）のＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式（式（１０）〜式（１４）も含む）に代入する（ステップＳ６０３）。そして、このＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式を解いて、さらに所定時間経過後の流体Ａ，Ｂ間の界面における算出流速ｕ，ｖ，ｗおよび算出圧力ｐを算出する（ステップＳ６０４）。
つぎに、ステップＳ６０４で算出された算出流速ｕ，ｖ，ｗを上記式（１５）に代入し（ステップＳ６０５）、上記式（１６）に展開する（ステップＳ６０６）。この展開された式（１６）に、ステップＳ５１０で算出された、変換された界面関数φを代入して、タンジェント変換値Ｆ（φ）を算出する（ステップＳ６０７）。
そして、上記式（１６）を式（１７）に変形し（ステップＳ６０８）、この変形した式（１７）に、ステップＳ６０７で算出されたＦ（φ）の値を代入して、あらたな界面関数φを算出する（ステップＳ６０９）。そして、このステップＳ６０９で算出した界面関数φを、上記式（１８）の拡散方程式に代入して、この拡散方程式を解き、ステップＳ６０９で算出された界面関数φを、界面近傍を平滑にする界面関数φに変換する（ステップＳ６１０）。以後、ステップＳ６０１〜ステップＳ６１０を繰り返し実行する。
このように、ＶＯＦ法のみによる流体解析では、非現実的な流体解析となる一方、本実施の形態にかかる流体解析装置１は、ＶＯＦ法にデジタイザ法を組み合わせて不連続な関数を正確に計算し、さらに拡散方程式により、平滑化処理を施すことにより、不連続すぎる関数によって再現されている界面を、実際の現象で生じている滑らかな（もっともらしい）界面に、簡単な解析方法により評価することができる。
なお、上述した実施の形態にかかる流体解析装置１においては、流体の支配方程式として、式（４）および式（５）に示すＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式を適用することとしたが、流体の支配方程式であれば、この式に限定されることはなく、他の流体の支配方程式も用いることができる。この場合、データ格納部２に格納されるデータは、適用された支配方程式に適用可能（代入可能）なデータを格納することとなる。
また、上述した実施の形態にかかる流体解析装置１においては、平滑化手法として、式（１８）に示す拡散方程式によって界面関数を変換して、界面近傍を平滑化することとしておるが、この拡散方程式に限定されることはなく、他に、たとえば、川原らの提案したセレクティブランピングパラメータ法などの界面関数に数値拡散（粘性）の作用を与える手法も適用することができる。
このセレクティブランピングパラメータ法は、下記式（１９）を用いることによって、対象とする界面関数φｎを、数値拡散の作用を与えた界面関数φｎ＋１に変換する手法である。
2004088563.tif 000014
この式（１９）中のＭＬｕｍｐｅｄは、時間微分項を有限要素法によって定式化したときに導出される質量行列Ｍを集中化した行列であり集中質量行列と呼ばれている。ここで、ｎ（ｎ＝１，２，・・・・・，ｍ）は反復回数であり、ｎ＝１のときの界面関数φ１が、上記式（１７）によって算出される界面関数である。
この式（１９）によって算出された界面関数φｎ＋１は、つぎの反復計算のときには、界面関数φｎとして用いられ、あらたな界面関数φｎ＋１を算出する。そして、ｍ回目に反復された界面関数φｍ＋１が、界面近傍を平滑にするために界面関数φ１から変換された界面関数φとなり、上記式（６）および式（７）の移流方程式に用いることとなる。
また、ＭＭｉｘｅｄは、下記式（２０）に示すように、集中質量行列ＭＬｕｍｐｅｄと質量行列Ｍを混合した行列であり混合質量行列と呼ばれている。
2004088563.tif 000015
ここで、ｅ  はランピングパラメータと呼ばれている。ｅは、式（１９）を使用するときの人工粘性の度合いを調節する値で、０．０≦ｅ＜１．０であり、その性質は  ｅの値が１．０に近づくほど人工的な粘性が少なくなるというものである。
拡散方程式に代えて、このセレクティブランピングパラメータ法を用いる場合、第４図に示すステップＳ４０６において、反復回数ｍと、ランピングパラメータｅを入力する。これにより、ステップＳ４０７において、対象とする界面関数φｎに対し、式（１９）の計算をｍ回反復させることができ、界面関数φの平滑化を実行することができる。
セレクティブランピングパラメータ法は、上記式（１９）を求めるのみでよいので、数値拡散の作用を与えるための計算が簡便である。この手法は有限要素法の定式化において質量行列の集中化が可能な場合に適用できる。また、入力する数値は、拡散方程式にあっては、仮想的な時間ｓおよび拡散係数Ｋをメッシュ幅に基づいて正確に規定している。
セレクティブランピングパラメータ法では、反復回数（１以上の整数）とランピングパラメータｅ（０．０≦ｅ＜１．０）であるため、入力する数値は単純な値でよく、拡散方程式における仮想的な時間ｓおよび拡散係数Ｋの設定よりも簡単に数値設定をおこなうことができる。
なお、上述した実施の形態にかかる流体解析装置１においては、選択した流体が気体と液体の場合について説明したが、液体同士、または気体同士の場合には、これらの界面における表面張力係数の実験値をあらかじめデータ格納部２に格納しておき、選択された流体の組み合わせによって、式（１３）に用いる表面張力係数を選択することとしてもよい。
つぎに、第４５図の解析モデルを用いて、ＶＯＦ法のみによる流体解析、ＶＯＦ法とデジタイザ法とを組み合わせた流体解析および本実施の形態にかかる流体解析装置１による流体解析と、を比較した例を第７図〜第１７図に示す。この第７図〜第１７図は、０．５秒ごとに流体解析した例である。
第７図に示すように、ＶＯＦ法とデジタイザ法との組み合わせでは、界面がガタガタに形成されているが、本実施の形態にかかる流体解析では、デジタイザ法によって得られた不連続すぎる界面関数φを拡散方程式によって平滑化処理を施している。このため、第７図〜第１７図に示すように、本実施の形態にかかる流体解析では、実際の現象における（表面張力の影響によって生じる）界面の滑らかさをあらわすことができる。
なお、ＶＯＦ法では、第７図〜第１７図に示すように、不連続な分布をもつ界面関数φの精度が悪いために、計算が進むにしたがって、２つの流体、すなわち空気と水の界面がぼやけてしまい、正確に２つの流体の界面が計算できていないこととなる。また、ＶＯＦ法とデジタイザ法との組み合わせでは、界面関数φの値が０と１のみであるため、極端に不連続となるため、第９図〜第１７図に示すように、０．７秒後には、計算不能に陥っていることがわかる。
つぎに、第１８図の解析モデルを用いた流体解析例を、第１９図および第２０図に示す。このシミュレーションは、液体中の気体（気泡）が水面に向かって上昇していく現象を計算したものである。液体中を上昇する気泡の形状は、液体と気体の物性値および気泡の大きさのバランスによって様々に変化する。
気泡流に関する実験は１９７３年にＪ．Ｒ．Ｇｒａｃｅらによって、液体と気体の物性値からなる３つの無次元された物性指標として、数種類の代表的な気泡の形状に分類できることが示されている。第１９図は、計算条件として気泡がＳｐｈｅｒｉｃａｌ（球状）と呼ばれる形状になる物性指標を採用した計算結果である。一方、第２０図は気泡がＳｋｉｒｔｅｄ（スカート形状）と呼ばれる形状になる物性指標を採用した計算結果である。第１９図および第２０図のいずれも無次元された物性指標に対応した特徴的な形状が再現されており、本実施の形態にかかる流体解析の信頼性を評価することができる。
つぎに、第２１図の解析モデルを用いた流体解析例を、第２２図〜第３０図に示す。この解析モデルは、ランスから高温高圧な純酸素ガスを流入して鉄や銅などの金属を精錬する転炉内部における流れ現象を解析するモデルである。第２２図〜第３０図に示すように、高速な酸素流入および両流体の密度が大きく異なる（純酸素ガスの密度３．０ｋｇ／ｍ３に対し、金属の密度は７．０×１０３ｋｇ／ｍ３）という過酷な条件下であっても、安定した解析結果を得ることができ、正確な界面評価をおこなうことができる。
このように、転炉内では、高温のため内部の状態を監視することが不可能であるが、本実施の形態の流体解析によれば、転炉内部の流体の挙動を捕捉することができる。これにより、転炉の大きさなど転炉の設計変更にも活用することができる。
つぎに、第３１図の解析モデルを用いた流体解析例を、第３２図〜第３９図に示す。第３２図は、固定メッシュに分割したメッシュ図である。この解析モデルは、ノズルからインクを噴射する流体の挙動を解析するモデルである。第３３図〜第３９図に示すように、インクがノズルから噴射されて底面に落ちる状態がわかる。この例で示すように、ミクロンオーダーの流体解析も忠実に再現することができ、肉眼では捉えることができない流体の挙動を捕捉することができる。
また、図示はしないが、本実施の形態にかかる流体解析装置１では、このほか、河川の氾濫、海上の波の流体挙動、容器内の液体の挙動、運転中におけるタンクローリー車が搭載するタンク内部の液体の挙動またはこれらをコンピュータグラフィックに応用した場合など、各種分野において適用することができる。
以上説明したように、本実施の形態にかかる流体解析装置１では、複雑な界面を扱うためにＶＯＦ法を採用し、精度の良い自由表面を計算するためにＣＩＰ法で提案・採用されているデジタイザ法を有限要素法（安定化気泡関数法）の解法に取り込んだ。
ＣＩＰ法は時間方向の定式化において陽解法を使用する必要があり、非常に過酷な計算条件（特に高速な流入速度の条件）のもとでは、数値的な安定性に限界があるが（ＣＦＬ条件）、安定化気泡関数法は陰解法を採用できるので時間方向に関する数値的な安定性に対して非常に有利である。さらに、有限要素法でデジタイザ法を取り込むにあたって、界面の評価のしかたを有限要素法で計算をするうえで正確に表現できるような平滑化手法を採用したことにより、不連続すぎる界面変数の不連続面を平滑化する処理をおこなうことができる。
なお、本実施の形態で説明した流体解析方法は、予め用意されたプログラムをパーソナル・コンピュータやワークステーション等のコンピュータで実行することにより実現することができる。このプログラムは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行される。またこのプログラムは、インターネット等のネットワークを介して配布することが可能な伝送媒体であってもよい。
【産業上の利用可能性】  以上説明したように、この発明の流体解析方法、流体解析プログラムおよび流体解析装置によれば、簡便で汎用的な解析手法により、実際の自由表面問題を持つ流れの挙動を正確に解析することができるという効果を奏する。さらに、物性値の異なる流体間における表面張力がもたらす界面上の滑らかさを界面関数に取り入れることができるという効果を奏する。
【図１】
2004088563.tif 000016
【図２】
2004088563.tif 000017
【図３】
2004088563.tif 000018
【図４】
2004088563.tif 000019
【図５】
2004088563.tif 000020
【図６】
2004088563.tif 000021
【図７】
2004088563.tif 000022
【図８】
2004088563.tif 000023
【図９】
2004088563.tif 000024
【図１０】
2004088563.tif 000025
【図１１】
2004088563.tif 000026
【図１２】
2004088563.tif 000027
【図１３】
2004088563.tif 000028
【図１４】
2004088563.tif 000029
【図１５】
2004088563.tif 000030
【図１６】
2004088563.tif 000031
【図１７】
2004088563.tif 000032
【図１８】
2004088563.tif 000033
【図１９】
2004088563.tif 000034
【図２０】
2004088563.tif 000035
【図２１】
2004088563.tif 000036
【図２２】
2004088563.tif 000037
【図２３】
2004088563.tif 000038
【図２４】
2004088563.tif 000039
【図２５】
2004088563.tif 000040
【図２６】
2004088563.tif 000041
【図２７】
2004088563.tif 000042
【図２８】
2004088563.tif 000043
【図２９】
2004088563.tif 000044
【図３０】
2004088563.tif 000045
【図３１】
2004088563.tif 000046
【図３２】
2004088563.tif 000047
【図３３】
2004088563.tif 000048
【図３４】
2004088563.tif 000049
【図３５】
2004088563.tif 000050
【図３６】
2004088563.tif 000051
【図３７】
2004088563.tif 000052
【図３８】
2004088563.tif 000053
【図３９】
2004088563.tif 000054
【図４０】
2004088563.tif 000055
【図４１】
2004088563.tif 000056
【図４２】
2004088563.tif 000057
【図４３】
2004088563.tif 000058
【図４４】
2004088563.tif 000059
【図４５】
2004088563.tif 000060
【図４６】
2004088563.tif 000061
【図４７】
2004088563.tif 000062
【図４８】
2004088563.tif 000063
【図４９】
2004088563.tif 000064
【国際調査報告】
2004088563.tif 000065
2004088563.tif 000066
2004088563.tif 000067