(57)【要約】【課題】安定化有限要素法に比べて定式化が飛躍的に簡便な流体解析をおこなうこと。【解決手段】安定化気泡関数有限要素流体解析装置１は、データ格納部２と、データ抽出部３と、流体挙動値算出部５と、表示制御部６と、から構成されている。データ格納部２には、流体ごとに、流体ＩＤ、流体名、既知物性値（密度、粘性係数）、初期流体挙動値（解析開始流速および解析開始圧力）からなる流体情報テーブルが格納されている。データ抽出部３は、入力装置を操作することにより、解析対象となる流体の流体ＩＤを選択し、データ格納部２から、その選択された流体ＩＤに関連付けられている流体名、既知物性値（密度、粘性係数）、初期流体挙動値（解析開始流速および解析開始圧力）を抽出する。流体挙動値算出部５は、Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ方程式に基づいて、流体の挙動値を算出する。表示制御部６は、流体挙動値算出部５における解析結果を出力する。【選択図】        図２
2005115934.tif 000002
【特許請求の範囲】
【請求項１】
所定の流体を選択する流体選択工程と、
前記流体選択工程によって選択された流体に関する物性値を抽出する物性値抽出工程と、
前記流体選択工程によって選択された流体に関する、予め設定された初期流体挙動値を抽出する初期流体挙動値抽出工程と、
前記物性値抽出工程によって抽出された物性値と、前記初期流体挙動値抽出工程によって抽出された初期流体挙動値と、に基づいて、流体の支配方程式を用いて、前記流体選択工程によって選択された流体の未知の挙動をあらわす未知流体挙動値を算出する未知流体挙動値算出工程と、
を含んだことを特徴とする安定化気泡関数有限要素流体解析方法。
【請求項２】
さらに、前記未知流体挙動値算出工程によって算出された未知流体挙動値と、前記物性値抽出工程によって抽出された物性値と、に基づいて、前記流体の支配方程式を用いて、前記流体選択工程によって選択された流体の未知の挙動をあらわす第２の未知流体挙動値を算出する第２の未知流体挙動値算出工程を含んだことを特徴とする請求項１に記載の安定化気泡関数有限要素流体解析方法。
【請求項３】
前記流体の支配方程式の粘性項における行列に、安定化項を用いたことを特徴とする請求項１または２に記載の安定化気泡関数有限要素流体解析方法。
【請求項４】
所定の流体を選択させる流体選択工程と、
前記流体選択工程によって選択された流体に関する物性値を抽出させる物性値抽出工程と、
前記流体選択工程によって選択された流体に関する、予め設定された初期流体挙動値を抽出する初期流体挙動値抽出工程と、
前記物性値抽出工程によって抽出された物性値と、前記初期流体挙動値抽出工程によって抽出された初期流体挙動値と、に基づいて、流体の支配方程式を用いて、前記流体選択工程によって選択された流体の未知の挙動をあらわす未知流体挙動値を算出させる未知流体挙動値算出工程と、
をコンピュータに実行させることを特徴とする安定化気泡関数有限要素流体解析プログラム。
【請求項５】
さらに、前記未知流体挙動値算出工程によって算出された未知流体挙動値と、前記物性値抽出工程によって抽出された物性値と、に基づいて、前記流体の支配方程式を用いて、前記流体選択工程によって選択された流体の未知の挙動をあらわす第２の未知流体挙動値を算出させる第２の未知流体挙動値算出工程を含んだことを特徴とする請求項４に記載の安定化気泡関数有限要素流体解析プログラム。
【請求項６】
前記流体の支配方程式の粘性項における行列に、安定化項を用いたことを特徴とする請求項４または５に記載の安定化気泡関数有限要素流体解析プログラム。
【請求項７】
所定の流体を選択する流体選択手段と、
前記流体選択手段によって選択された流体に関する物性値を抽出する物性値抽出手段と、
前記物性値抽出手段によって抽出された物性値と、前記流体選択手段によって選択された流体の現在の既知挙動値と、に基づいて、流体の支配方程式を用いて、前記流体選択手段によって選択された流体の未知の挙動をあらわす未知流体挙動値を算出する未知流体挙動値算出手段と、
を備えることを特徴とする安定化気泡関数有限要素流体解析装置。
【請求項８】
前記流体の支配方程式の粘性項における行列に、安定化項を用いたことを特徴とする請求項７に記載の安定化気泡関数有限要素流体解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、数値流体力学における非圧縮粘性流れの解析について、数値流体の応用問題に有効かつ計算精度の高い数値シミュレーションを行うための安定化気泡関数有限要素流体解析方法、安定化気泡関数有限要素流体解析プログラムおよび安定化気泡関数有限要素流体解析装置に関する。
【背景技術】
【０００２】
非圧縮粘性流体の支配方程式（運動方程式、連続式）を通常の有限要素法で定式化すると下記式（１）、（２）のように表すことができる（たとえば、下記非特許文献１参照。）。
【０００３】
【数１】
2005115934.tif 000003
【０００４】
【数２】
2005115934.tif 000004
【０００５】
式（１）、（２）のｕi、ｐは流速、圧力であり求めるべき未知量である。Ｍ、Ｋ（ｕj）、Ｂ、Ｓ、Ｎ、ＭΓは質量行列、移流行列、圧力行列、粘性行列、外力行列、境界積分行列である。低次要素（三角形一次要素、四面体一次要素）を使用した場合には、流速ｕi、圧力ｐの未知量は要素レベルにて図３１、図３２のように保持される。
【０００６】
非定常問題における一般的な有限要素法は、時間方向に差分近似を用いて解析をおこなう。圧力項と連続式を完全陰的に扱い、移流速度にAdams-Bashforth近似を適用し、他の項はCrank-Nicolson近似を用いた場合の式（１）および式（２）は、下記式（３）および式（４）になる（たとえば、下記非特許文献２参照）。
【０００７】
【数３】
2005115934.tif 000005
【０００８】
【数４】
2005115934.tif 000006
【０００９】
式（３）中の上付き添え字ｎは現在時刻ｎにおける既知物理量を表しており、ｎ＋１は時刻ｎから微小時間Δｔ経過後における未知物理量をあらわしている。
【００１０】
非圧縮粘性流体の解析では、ある物体を流れによって運ぼうとする力（移流力）を表現する項を数値的に離散化した場合、計算上、非常に不安定な振る舞いをする。この数値的な不安定性のために通常の有限要素法を用いた場合には、計算が不能になることが多く、また、計算が行えたとしても得られた解は、非常に精度が悪く実際の現象とはかけ離れた結果を得る。
【００１１】
さらに、流速と圧力が未知量となるが、この未知量を解くべき行列においてその対角成分に一部ゼロの値が現れてしまうことから（式（３）参照）、流速と圧力の補間関数の組み合わせによっては、逆行列が存在しなくなるような状態になり結果として非常に激しい圧力振動が現れる。
【００１２】
例えば、流速と圧力の補間関数の組み合わせとして同じ補間関数を用いると、解くべき行列の状態が非常に悪くなり、図３３のように圧力振動が原因で実際の現象とは違った計算結果を得る。
【００１３】
この数値的な不安定性、および行列の状態の悪さを取り除くための代表的な解法として、スタンフォード大学の研究グループで開発された安定化有限要素法がある。この解法は、数値的な不安定性を回避するために数学的に意味のある安定化項を考慮して計算を行うものである（たとえば、下記非特許文献３参照。）。下記式（５）および式（６）に低次要素にて代表的な安定化有限要素法の一つであるSUPG/PSPG法を適用した有限要素方程式とCavity流れの計算結果を図３４に示す。
【００１４】
【数５】
2005115934.tif 000007
【００１５】
【数６】
2005115934.tif 000008
【００１６】
式（５）および式（６）のそれぞれ第２行目の項が安定化項である。圧力項と連続式を完全陰的に扱い、移流速度にAdams-Bashforth近似を適用し、他の項はCrank-Nicolson近似を用いた場合の式（５）および式（６）は、下記式（７）になる。
【００１７】
【数７】
2005115934.tif 000009
【００１８】
安定化有限要素法の問題点としては、定式化が最も簡単になる低次要素を使用した場合でも安定化項が非常に複雑であり、また、この項数が多い（約５０〜１００項）。これは、実問題、流体―構造連成問題、逆解析など応用問題への適用を考えると、更にアルゴリズムが複雑になり、具体的な定式化、およびそのプログラミングが困難になる。
【００１９】
定式化の簡便な有限要素解析手法として、安定化有限要素法とは異なった定式化で流体解析を行うことのできる気泡関数要素を用いた有限要素法が提案されている。気泡関数要素とは、最も簡単な補間関数である低次要素に図３５、図３６に示す重心点の自由度を加えた要素である。重心点の自由度は、安定化有限要素法で用いられている安定化項と密接な関係があり重心点の付加によって複雑な安定化項の考慮を還元している（たとえば、下記非特許文献４参照。）。
【００２０】
【非特許文献１】日本数値流体力学会有限要素法研究委員会編、有限要素法による流れのシミュレーション、シュプリンガー・フェアラーク東京、１９９８年１１月９日
【非特許文献２】ジェイ・シー・シモ、エフ・アーメロ（J.C.Simo and F.Armero）著、「アンコンディショナル  スタビリティ  アンド  ロングターム  ビヘイバー  オブ  トランジェント  アルゴリズムズ  フォー  ザ  インコンプレッシブル  ナビエ−ストークス  アンド  オイラー  イクエイジョンズ（Unconditional Stability and Long-Term Behavior of Transient Algorithms for the Incompressible Navier-Stokes and Euler Equations）」、コンピュータ  メソッズ  イン  アプライド  メカニクス  アンド  エンジニアリング（Computer Methods in Applied Mechanics and Engineering）、１９９４年、１１１号、１１１頁〜１５４頁
【非特許文献３】ティー・イー・テジュヤール、エス・ミタル、エス・イー・レイ、アール・シン（T.E.Tezduyar, S.Mittal, S.E.Ray and R.Shin）著、「インコンプレッシブル  フロー  コンピューテイションズ  ウィズ  スタビライズド  ビリニアー  アンド  リニアー  イコール−オーダー−インターポレイション  ベロシティ−プレッシャ  エレメント（Incompressible Flow Computations with Stabilized Bilinear and Linear Equal-Order-Interpolation Velocity-Pressure Element）」、コンピュータ  メソッズ  イン  アプライド  メカニクス  アンド  エンジニアリング（Computer Methods in Applied Mechanics and Engineering）、１９９２年、９５号、２２１頁〜２４２頁
【非特許文献４】ジェイ・シー・シモ、エフ・アーメロ、シー・エー・テイラー（J.C.Simo, F.Armero and C.A.Taylor）著、「ステイブル  アンド  タイム−ディシパティブ  ファイナイト  エレメント  メソッズ  フォー  ジ  インコンプレッシブル  ナビエ−ストークス  イクエイジョンズ  イン  アドベクション  ドミネイテッド  フローズ（Stable and Time-Dissipative Finite Element Methods for the Incompressible Navier-Stokes Equations in Advection Dominated Flows）  インターナショナル  ジャーナル  フォー  ニューメリカル  メソッズ  イン  エンジニアリング（International Journal for Numerical Methods in Engineering）、１９９５年、３８号、１４７５頁〜１５０６頁
【発明の開示】
【発明が解決しようとする課題】
【００２１】
しかしながら、気泡関数要素の大きな問題として、重心点の付加による計算上の安定性には限界があり、その結果、圧力振動が発生してしまうという問題があり、流体の挙動を正確に表現することができなかった（図３７参照）。
【００２２】
上述した従来技術による問題点を解消するため、定式化の簡便性には気泡関数要素を適用し、気泡関数のもつ数値計算上の不安定性には、有限要素方程式（下記式（８）および式（９））において、式（８）の左辺第５項の重心点のみの粘性項を用いる。この発明は、有限要素法の初級者レベルの知識があれば、だれもが簡単に作成することのできる低次要素（三角形・四面体一次要素）のプログラムから、簡単なプログラムの変更で高精度な安定化気泡関数有限要素流体解析方法、安定化気泡関数有限要素流体解析プログラムおよび安定化気泡関数有限要素流体解析装置を提供することを目的とする。
【００２３】
【数８】
2005115934.tif 000010
【００２４】
【数９】
2005115934.tif 000011
【課題を解決するための手段】
【００２５】
上記従来技術の問題点を解消するため、この発明の安定化気泡関数有限要素流体解析方法、安定化気泡関数有限要素流体解析プログラムおよび安定化気泡関数有限要素流体解析装置は、所定の流体を選択し、選択された流体に関する物性値を抽出し、選択された流体に関する、予め設定された初期流体挙動値を抽出し、抽出された物性値と初期流体挙動値とに基づいて、流体の支配方程式を用いて、選択された流体の未知の挙動をあらわす未知流体挙動値を算出することを特徴とする。
【発明の効果】
【００２６】
本発明にかかる流体解析方法、流体解析プログラムおよび流体解析装置によれば、気泡関数要素の安定化作用の限界について、非定常問題においても重心点のみの粘性項を、2レベル-3レベル有限要素近似（松本純一, 気泡関数を用いた非圧縮性粘性流れ解析のための2レベル-3レベル有限要素法, 応用力学論文集（土木学会）, 7, pp.339-346, 2004.）に基づいた方法（安定化気泡関数有限要素法）により導出し、その計算上の数値不安定性を回避することができるという効果を奏する。
【００２７】
また、定式化が安定化有限要素法に比べて飛躍的に簡便である（各方程式に対して重心点のみに一つの安定化項）という効果を奏する。このことから、実問題や応用問題適用への具体的な定式化、およびそのプログラミングにおいて非常に有効な方法である。さらに、安定化有限要素法と同等の計算精度を有するという効果を奏する。
【発明を実施するための最良の形態】
【００２８】
以下に添付図面を参照して、この発明にかかる安定化気泡関数有限要素流体解析方法、安定化気泡関数有限要素流体解析プログラムおよび安定化気泡関数有限要素流体解析装置の好適な実施の形態を詳細に説明する。
【００２９】
（安定化気泡関数有限要素流体解析装置のハードウェア構成）
まず、図１を用いて、この実施の形態にかかる安定化気泡関数有限要素流体解析装置のハードウェア構成について説明する。安定化気泡関数有限要素流体解析装置は、ＣＰＵ１０１と、ＲＯＭ１０２とＲＡＭ１０３と、ＨＤＤ（ハードディスクドライブ）１０４と、ＨＤ（ハードディスク）１０５と、ＦＤＤ（フレキシブルディスクドライブ）１０６と、着脱可能な記録媒体の一例としてＦＤ（フレキシブルディスク）１０７と、ディスプレイ１０８と、Ｉ／Ｆ（インターフェース）１０９と、キーボード１１０と、マウス１１１と、スキャナ１１２と、プリンタ１１３と、を備えている。また、各構成部は、バス１００によってそれぞれ接続されている。
【００３０】
ここで、ＣＰＵ１０１は、安定化気泡関数有限要素流体解析装置の全体の制御を司る。ＲＯＭ１０２は、ブートプログラムなどのプログラムを記憶している。ＲＡＭ１０３は、ＣＰＵ１０１のワークエリアとして使用される。ＨＤＤ１０４は、ＣＰＵ１０１の制御にしたがってＨＤ１０５に対するデータのリード／ライトを制御する。ＨＤ１０５は、ＨＤＤ１０４の制御で書き込まれたデータを記憶する。
【００３１】
ＦＤＤ１０６は、ＣＰＵ１０１の制御にしたがってＦＤ１０７に対するデータのリード／ライトを制御する。ＦＤ１０７は、ＦＤＤ１０６の制御で書き込まれたデータを記憶したり、ＦＤ１０７に記憶されたデータを安定化気泡関数有限要素流体解析装置に読み取らせたりする。着脱可能な記録媒体として、ＦＤ１０７のほか、ＣＤ−ＲＯＭ（ＣＤ−Ｒ、ＣＤ−ＲＷ）、ＭＯ、ＤＶＤ（Digital Versatile Disk）、メモリーカードなどであってもよい。ディスプレイ１０８は、カーソル、アイコンあるいはツールボックスをはじめ、文書、画像、機能情報などのデータを表示する。このディスプレイ１０８は、たとえば、ＣＲＴ、ＴＦＴ液晶ディスプレイ、プラズマディスプレイ等を採用することができる。
【００３２】
Ｉ／Ｆ１０９は、通信回線を通じてインターネットなどのネットワークに接続され、このネットワークを介して他の装置に接続される。そして、Ｉ／Ｆ１０９は、ネットワークと内部のインターフェースを司り、外部装置からのデータの入出力を制御する。Ｉ／Ｆ１０９には、たとえばモデムやＬＡＮアダプタなどを採用することができる。
【００３３】
キーボード１１０は、文字、数字、各種指示などの入力のためのキーを備え、データの入力をおこなう。また、タッチパネル式の入力パッドやテンキーなどであってもよい。マウス１１１は、カーソルの移動や範囲選択、あるいはウィンドウの移動やサイズの変更などをおこなう。ポインティングデバイスとして同様に機能を備えるものであれば、トラックボールやジョイスティックなどであってもよい。
【００３４】
スキャナ１１２は、画像を光学的に読み取り、安定化気泡関数有限要素流体解析装置内に画像データを取り込む。また、プリンタ１１３は、画像データや文書データを印刷する。プリンタ１１３には、たとえば、レーザプリンタやインクジェットプリンタを採用することができる。
【００３５】
（安定化気泡関数有限要素流体解析装置の機能的構成）
つぎに、本実施の形態にかかる安定化気泡関数有限要素流体解析装置の機能的構成について説明する。図２に示すように、安定化気泡関数有限要素流体解析装置１は、データ格納部２と、データ抽出部３と、流体挙動値算出部５と、表示制御部６と、から構成されている。
【００３６】
データ格納部２には、図３に示すように、流体ごとに、流体ＩＤ、流体名、流体の既知物性値（密度、粘性係数、（無次元化した場合はレイノルズ数））、初期流体挙動値（解析開始流速および解析開始圧力）からなる流体情報テーブルが格納されている。また、データ格納部２には、図示はしないが、初期設定情報として、時間方向に対して、逐次解析計算をおこなうために時間増分量、計算ステップ数などの情報も格納されている。
【００３７】
さらに、解析対象範囲に分割形成される固定メッシュの分割幅、節点数およびメッシュ数などのメッシュ情報も格納されている。この固定メッシュは、解析対象範囲が２次元空間である場合は、三角形、解析対象範囲が３次元空間である場合は、四面体の有限な形状によって構成される。メッシュ情報は、スキャナ１１２により読み取ることによっても取り込むことができる。このデータ格納部２は、具体的には、例えば図１に示したＲＡＭ１０３、ＨＤ１０５またはＦＤ１０７によりその機能を実現する。
【００３８】
データ抽出部３は、入力装置（たとえば図１に示すキーボード１１０またはマウス１１１）を操作することにより、解析対象となる流体の流体ＩＤを選択し、データ格納部２から、その選択された流体ＩＤに関連付けられている流体名、流体の既知物性値（密度、粘性係数、（無次元化した場合はレイノルズ数））、初期流体挙動値（解析開始流速および解析開始圧力）を抽出する。このデータ抽出部３は、具体的には、たとえば図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。
【００３９】
流体挙動値算出部５は、Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ方程式に基づいて、流体の挙動値を算出する。この流体挙動値算出部５は、具体的には、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。また、流体挙動値算出部５における算出処理の具体的説明については後述する。
【００４０】
また、表示制御部６は、流体挙動値算出部５における解析結果を図１に示すディスプレイ１０８に出力する。この表示制御は、流体解析に一般的に用いられるためその説明を省略する。この表示制御部６は、具体的には、たとえば図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。
【００４１】
つぎに、本実施の形態にかかる処理手順について、図４のフローチャートを用いて説明する。まず、入力装置（たとえば図１に示すキーボード１１０またはマウス１１１）を介して、解析対象となる流体ＩＤが選択された場合（ステップＳ４０１：Ｙｅｓ）、データ格納部２から、選択された流体ＩＤを参照して、選択された流体の流体物性値である密度ρ、粘性係数μ、（無次元化した場合はレイノルズ数Ｒｅ）、選択された流体の初期流体挙動値である解析開始流速ｕ，ｖ，ｗおよび解析開始圧力ｐを抽出する（ステップＳ４０２）。
【００４２】
そして、入力装置から解析開始入力を受け付けた場合（ステップＳ４０３：Ｙｅｓ）、データ格納部２から、あらかじめ格納された解析対象範囲の固定メッシュを抽出する（ステップＳ４０４）。そして、流体挙動値算出部５による算出処理を実行する（ステップＳ４０５）。この算出処理は、所定の微小時間（たとえば１秒）間隔で算出処理を実行され、終了入力があるまでおこなう（ステップＳ４０６：Ｎｏ）。そして、終了入力があった場合（ステップＳ４０６：Ｙｅｓ）は、終了する。これにより、時々刻々と変化する流体の挙動値を算出することができる。
【００４３】
つぎに、流体挙動値算出部５およびその算出処理内容を具体的に説明する。
＜１．Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ方程式＞
非圧縮粘性流れにおける基礎方程式は、以下の運動方程式（１０）と連続式（１１）によって表される。ここで、ｕi，ｐ，ｆiは、ｘi(ｉ=１，２，３)方向の流速，圧力，外力である。また、νは、動粘性係数μ／ρ、無次元化した場合は１／Ｒｅであり、ρ、μ、Ｒｅは、密度、粘性係数、レイノルズ数である。
【００４４】
【数１０】
2005115934.tif 000012
【００４５】
上付き添え字「・」、下付き添え字「,ｉ」をそれぞれ時間に関する偏微分、空間に関する偏微分とすると、上式（１０）、（１１）は、以下の式（１２）、（１３）のように書き換えられる。
【００４６】
【数１１】
2005115934.tif 000013
【００４７】
境界Γは，Γ1とΓ2に分けられ、以下の式（１４）、（１５）のように境界条件が規定される。
【００４８】
【数１２】
2005115934.tif 000014
【００４９】
ここで，δijはクロネッカーのデルタ関数、ｎjは境界Γ2の外向き法線ベクトルである。
【００５０】
＜２．有限要素方程式＞
上記に示したＮａｖｉｅｒ‐Ｓｔｏｋｅｓ方程式を有限要素法によって定式化をおこなった方程式（有限要素方程式（１６）〜（１８））は次のように表すことができる。
【００５１】
【数１３】
2005115934.tif 000015
【００５２】
有限要素法は、任意形状の解析領域を簡単な幾何学形状（要素と呼ぶ）を用いて分割し表現する。実際に、上式の有限要素方程式を作成するためには、各要素において方程式を作成し、要素レベルの有限要素方程式を重ね合せるといった手段を用いて全体形（任意形状）の有限要素方程式を作成する。上記に示した有限要素方程式（１６）〜（１８）を、要素レベルの有限要素方程式を用いて表現すると次式（１９）、式（２０）のようになる。
【００５３】
【数１４】
2005115934.tif 000016
【００５４】
Ｎkは総要素数、ｋはｋ番目の要素である。また、各要素領域における積分記号を下記式（２１）のように定義する。
【００５５】
【数１５】
2005115934.tif 000017
【００５６】
上式（１９）〜式（２１）により、以下の各要素レベルの有限要素方程式（２２）、（２３）を得る。
【数１６】
2005115934.tif 000018
【００５７】
要素の最も簡単な幾何学形状としては、２次元では三角形形状、３次元では四面体形状がある。以下に要素として三角形形状（２次元）、四面体形状（３次元）を用いて幾何学形状内の値を補間するための関数を示す。
【００５８】
＜３．三角形要素と四面体要素の補間関数＞
三角形一次要素（圧力の補間関数）
圧力の補間は一次関数を用いた補間関数を使用する（下記式（２４）〜式（２８）を参照。）。
【００５９】
【数１７】
2005115934.tif 000019
【００６０】
三角形一次要素（流速の補間関数）については、下記式（２９）〜式（３４）を用いる。
【数１８】
2005115934.tif 000020
【００６１】
有限要素法では、ｕ(k)i、ｖ(k)i、ｐ(k)、ｑ(k)を補間関数と呼び、Ψα（α＝１，２，３）を形状関数と呼んでいる。また、補間関数内に出てくる番号は、図５〜図８の番号に対応している。
【００６２】
三角形気泡関数要素（流速の補間関数、表現形式１）（下記式（３５）〜式（４１）を参照。）
【数１９】
2005115934.tif 000021
【００６３】
三角形気泡関数要素（流速の補間関数、表現形式２）
三角形気泡関数要素の補間関数は次式のような記述でも表現することができる。通常は、補間関数として表現形式１を使用するが、時間微分項に出てくる行列（質量行列）を集中化して時間微分項の近似を行う場合には、表現形式２が使用される（下記式（４２）〜式（４７）を参照。）。
【００６４】
【数２０】
2005115934.tif 000022
【００６５】
四面体一次要素（圧力の補間関数）
圧力の補間は一次関数を用いた補間関数を使用する（下記式（４８）〜式（５２）を参照。）。
【００６６】
【数２１】
2005115934.tif 000023
【００６７】
四面体一次要素（流速の補間関数）（下記式（５３）〜式（５８）を参照。）
【数２２】
2005115934.tif 000024
【００６８】
四面体気泡関数要素（流速の補間関数、表現形式１）（下記式（５９）〜式（６５）を参照。）
【数２３】
2005115934.tif 000025
【００６９】
四面体気泡関数要素（流速の補間関数、表現形式２）
四面体気泡関数要素の補間関数は次式のような記述でも表現することができる。通常は、補間関数として表現形式１を使用するが、時間微分項に出てくる行列（質量行列）を集中化して時間微分項の近似を行う場合には、表現形式２が使用される（下記式（６６）〜式（７１）を参照。）。
【００７０】
【数２４】
2005115934.tif 000026
【００７１】
＜４．要素レベルの有限要素方程式＞
以上に示した、補間関数を導入して、移流項の移動速度ｕ(k)jを要素内平均値で近似し、外力項のｆiを要素内定数値と仮定すると、要素レベルの有限要素方程式は次式のようになる（下記式（７２）、式（７３）を参照。）。
【００７２】
【数２５】
2005115934.tif 000027
【００７３】
ｖTi、ｑTiは、任意の値をとるベクトルであるので、要素レベルの有限要素方程式は以下の式（７４）および式（７５）のように表すことができる。
【００７４】
【数２６】
2005115934.tif 000028
【００７５】
＜５．気泡関数について＞
重心点での形状関数φBは気泡関数と呼ばれており、形状関数の定義（重心点が１、その他の点では０となる関数）を満たせば任意の関数を選択することができる。気泡関数の例として良く使われている気泡関数を以下に示す。図９および図１０に示したアイソパラメトリック座標系｛ｒ,ｓ｝（｛ｒ,ｓ,ｔ｝）を用いて次式（７６）〜式（７９）のように定義される。
【００７６】
【数２７】
2005115934.tif 000029
【００７７】
【数２８】
2005115934.tif 000030
【００７８】
つぎに、式（７４）に代入されるパラメータを列挙する。
Ａ．流速の補間関数に三角形気泡関数要素（表現形式１）、圧力の補間関数に三角形１次要素を用い、粘性項に重心点のみの安定化項を使用した場合、粘性項は以下の式（８０）の通りである。
【００７９】
【数２９】
2005115934.tif 000031
【００８０】
上式中のτkSは安定化パラメータと呼ばれており、計算を安定にかつ、高精度に行うためのパラメータである。τkSの形式は幾つか提案されている。例として良く使われている安定化パラメータτkSを以下の式（８１）および式（８２）に示す。
【００８１】
【数３０】
2005115934.tif 000032
【００８２】
Ｂ．流速の補間関数に三角形気泡関数要素（表現形式２）、圧力の補間関数に三角形１次要素を用い、粘性項に安定化項を使用した場合、粘性項は以下の式（８３）の通りである。
【００８３】
【数３１】
2005115934.tif 000033
【００８４】
Ｃ．流速の補間関数に四面体１次要素、圧力の補間関数に四面体１次要素を用いた場合、粘性項は以下の式（８４）の通りである。
【００８５】
【数３２】
2005115934.tif 000034
【００８６】
Ｄ．流速の補間関数に四面体気泡関数要素（表現形式２）、圧力の補間関数に四面体１次要素を用い、粘性項に安定化項を使用した場合、粘性項は以下の式（８５）の通りである。
【００８７】
【数３３】
2005115934.tif 000035
【００８８】
＜６．Ｃａｖｉｔｙ流れ＞
発明手法の数値的な安定性を評価するために、検証としてＣａｖｉｔｙ内の流れを取り上げる。解析領域と境界条件を図１１および図１２に示す。レイノルズ数をＲｅ＝１、Ｒｅ＝４００として検証の計算を実行した。図１３および図１４に示すように３２×３２、３２×３２×３２の均等分割を行ったメッシュを用いている。
【００８９】
一般的にＣａｖｉｔｙ流れは、定常解を持つ問題であるので、計算は定常状態に収束するまで行っている。時間増分はΔｔ＝０．０００５（Ｒｅ＝１）、０．２（Ｒｅ＝４００）を採用した。検証の計算は、「ＳＵＰＧ／ＰＳＰＧ法：一次要素(流速)/一次要素(圧力)」（時間方向には直接法を使用：非特許文献２参照。）、「気泡関数要素(流速)/一次要素(圧力)」（時間方向には分離型解法を使用：J.Matsumoto and M.Kawahara,  Shape Identification for Fluid-Structure Interaction Problem Using Improved Bubble Element,  International Journal of Computational Fluid Dynamics,  2001, 15, pp.33-pp.45を参照。）、「発明解法」（時間方向には分離型解法を使用：J.Matsumoto and M.Kawahara,  Shape Identification for Fluid-Structure Interaction Problem Using Improved Bubble Element,  International Journal of Computational Fluid Dynamics,  2001, 15, pp.33-pp.45を参照。）の三つの解法で比較した。
【００９０】
図１５〜図２６に圧力の計算結果を示す。図１５〜図２６から解るように、２次元（三角形）、３次元（四面体）ともに、「気泡関数要素(流速)/一次要素(圧力)」を用いた場合の計算結果は、数値的な安定性に限界があり圧力振動が著しく発生している。これに対して、「発明解法」を用いた場合の計算結果は、適切な数値粘性が与えられていることから圧力振動が発生していないことが解る。また、「発明解法」は「ＳＵＰＧ／ＰＳＰＧ法：一次要素(流速)/一次要素(圧力)」と同様の計算結果を得ている。
【００９１】
＜７．Ｓｔａｎｄｉｎｇ  Ｖｏｒｔｅｘ問題＞
この問題の目的は，計算手法がどの程度の数値的な散逸誤差をもつかを調べる問題である。非粘性流れを仮定し，初期条件として人工的な渦を発生させ計算を進める。厳密解は初期条件そのものであるが、安定な計算ならかならず解は減衰する。初期の運動エネルギーがどの程度、減衰をしたかを調べることによって数値的な散逸誤差を評価する。解析領域(図２７を参照。)は、１×１であり、２０×２０の均等な有限要素分割を用いている（図２８を参照。）。初期条件は領域の中心からの極座標（ｒ，θ）による式によって表される（図２９を参照。）。
【００９２】
検証の計算は、「ＳＵＰＧ／ＰＳＰＧ法：一次要素(流速)/一次要素(圧力)」（時間方向には直接法を使用：非特許文献２参照。）、「気泡関数要素(流速)/一次要素(圧力)」（時間方向には分離型解法を使用：J.Matsumoto and M.Kawahara,  Shape Identification for Fluid-Structure Interaction Problem Using Improved Bubble Element,  International Journal of Computational Fluid Dynamics,  2001, 15, pp.33-pp.45を参照。）、「発明解法」（時間方向には分離型解法を使用：J.Matsumoto and M.Kawahara,  Shape Identification for Fluid-Structure Interaction Problem Using Improved Bubble Element,  International Journal of Computational Fluid Dynamics,  2001, 15, pp.33-pp.45を参照。）の三つの解法で比較を行っている。無次元時間ｔ＝３．０（Ｒｅ＝０．０、時間増分量Δｔ＝０．０５を使用）のときの運動エネルギーの保存率を図３０に示す。図３０をみると、通常の気泡関数（気泡関数）に比べて、発明解法のほうが、運動エネルギーの減衰が非常に少ないことが解る。また、計算精度の高いＳＵＰＧ／ＰＳＰＧ法（ＳＵＰＧ／ＰＳＰＧ）と比べてもほぼ同等の計算精度を有しており、発明解法の精度の高さがわかる。
【００９３】
以上説明したように、本実施の形態によれば、非定常解析においても気泡関数要素の安定化作用の限界を重心点のみのたった１項の粘性項を、2レベル-3レベル有限要素近似に基づいた方法（安定化気泡関数有限要素法）により導出し、その計算上の不安定性を回避することができる。
【００９４】
また、定式化が安定化有限要素法に比べて飛躍的に簡便であることから、実問題や応用問題適用への具体的な定式化、およびそのプログラミングにおいて非常に有効な方法である。さらに、安定化有限要素法と同等の計算精度をもっている。
【００９５】
なお、本実施の形態で説明した安定化気泡関数有限要素流体解析方法は、予め用意されたプログラムをパーソナル・コンピュータやワークステーション等のコンピュータで実行することにより実現することができる。このプログラムは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行される。またこのプログラムは、インターネット等のネットワークを介して配布することが可能な伝送媒体であってもよい。
【産業上の利用可能性】
【００９６】
以上のように、本発明にかかる安定化気泡関数有限要素流体解析方法、安定化気泡関数有限要素流体解析プログラムおよび安定化気泡関数有限要素流体解析装置は、三角形および四面体形状を使用した有限要素法による流体解析に有用であり、特に、流体−構造連成問題、逆解析などの定式化が非常に煩雑となる応用問題に適している。
【図面の簡単な説明】
【００９７】
【図１】この発明の本実施の形態にかかる安定化気泡関数有限要素流体解析装置のハードウェア構成を示すブロック図である。
【図２】この発明の本実施の形態にかかる安定化気泡関数有限要素流体解析装置の機能構成を示すブロック図である。
【図３】この発明の本実施の形態にかかる安定化気泡関数有限要素流体解析装置におけるデータ格納部に格納されている流体情報テーブルを示す説明図である。
【図４】この発明の本実施の形態にかかる安定化気泡関数有限要素流体解析装置における流体解析処理手順を示すフローチャートである。
【図５】三角形気泡関数要素の補間関数に関する説明図である。
【図６】三角形一次要素の補間関数に関する説明図である。
【図７】四面体気泡関数要素の補間関数に関する説明図である。
【図８】四面体一次要素の補間関数に関する説明図である。
【図９】三角形を使用した場合の気泡関数に関する説明図である。
【図１０】四面体を使用した場合の気泡関数に関する説明図である。
【図１１】２次元Ｃａｖｉｔｙ流れの解析モデルと境界条件を示す説明図である。
【図１２】３次元Ｃａｖｉｔｙ流れの解析モデルと境界条件を示す説明図である。
【図１３】図１０の解析モデルにおいて、三角形分割されたメッシュを示す説明図である。
【図１４】図１１の解析モデルにおいて、四面体分割されたメッシュを示す説明図である。
【図１５】図１０の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図１６】図１０の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図１７】図１０の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図１８】図１０の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図１９】図１０の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図２０】図１１の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図２１】図１１の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図２２】図１１の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図２３】図１１の解析モデルにおいて、Ｒｅ＝１を使用した場合の圧力図を示す説明図である。
【図２４】図１１の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図２５】図１１の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図２６】図１１の解析モデルにおいて、Ｒｅ＝４００を使用した場合の圧力図を示す説明図である。
【図２７】Ｓｔａｎｄｉｎｇ  Ｖｏｒｔｅｘ問題の解析モデルを示す説明図である。
【図２８】図２７の解析モデルにおいて、三角形分割されたメッシュを示す説明図である。
【図２９】図２７の解析モデルにおいて、初期流速ベクトル値を示す説明図である。
【図３０】図２７の解析モデルにおいて、各解析手法の運動エネルギー保存率を示す説明図である。
【図３１】流速と圧力において、三角形一次要素を使用した場合の、未知量の情報を示す説明図である。
【図３２】流速と圧力において、四面体一次要素を使用した場合の、未知量の情報を示す説明図である。
【図３３】３次元Ｃａｖｉｔｙ流れ解析における、圧力図を示す説明図である。
【図３４】３次元Ｃａｖｉｔｙ流れ解析における、圧力図を示す説明図である。
【図３５】流速と圧力において、三角形気泡関数要素、三角形一次要素を使用した場合の、未知量の情報を示す説明図である。
【図３６】流速と圧力において、四面体気泡関数要素、四面体一次要素を使用した場合の、未知量の情報を示す説明図である。
【図３７】３次元Ｃａｖｉｔｙ流れ解析における、圧力図を示す説明図である。
【符号の説明】
【００９８】
１  安定化気泡関数有限要素流体解析装置
２  データ格納部
３  データ抽出部
５  流体挙動値算出部
６  表示制御部
【図１】
2005115934.tif 000036
【図２】
2005115934.tif 000037
【図３】
2005115934.tif 000038
【図４】
2005115934.tif 000039
【図５】
2005115934.tif 000040
【図６】
2005115934.tif 000041
【図７】
2005115934.tif 000042
【図８】
2005115934.tif 000043
【図９】
2005115934.tif 000044
【図１０】
2005115934.tif 000045
【図１１】
2005115934.tif 000046
【図１２】
2005115934.tif 000047
【図１３】
2005115934.tif 000048
【図１４】
2005115934.tif 000049
【図１５】
2005115934.tif 000050
【図１６】
2005115934.tif 000051
【図１７】
2005115934.tif 000052
【図１８】
2005115934.tif 000053
【図１９】
2005115934.tif 000054
【図２０】
2005115934.tif 000055
【図２１】
2005115934.tif 000056
【図２２】
2005115934.tif 000057
【図２３】
2005115934.tif 000058
【図２４】
2005115934.tif 000059
【図２５】
2005115934.tif 000060
【図２６】
2005115934.tif 000061
【図２７】
2005115934.tif 000062
【図２８】
2005115934.tif 000063
【図２９】
2005115934.tif 000064
【図３０】
2005115934.tif 000065
【図３１】
2005115934.tif 000066
【図３２】
2005115934.tif 000067
【図３３】
2005115934.tif 000068
【図３４】
2005115934.tif 000069
【図３５】
2005115934.tif 000070
【図３６】
2005115934.tif 000071
【図３７】
2005115934.tif 000072