JP5255714B2 - 三次元の流体シミュレーション方法 - Google Patents

Description

本発明は、三次元の流体シミュレーション方法に関する。

従来より、物体の回りを流れる流体の状態（速度、圧力等）をコンピュータを用いて解析する流体シミュレーションが種々行われている。このような流体シミュレーションは、物体の表面形状、例えばゴルフボールのディンプルについて、空気抵抗の低減や飛行性能の向上のための開発に役立っている。

ところで、流体シミュレーションでは、流体が流れる空間が、格子状にメッシュ分割される。このメッシュ分割に関して、物体から比較的離れた領域では、流体の流れが単純になりやすいため、大きなメッシュ（即ち、格子密度が小）で構成されるのが計算効率の面で望ましい。他方、物体の近くの領域では、流体の流れが複雑になりやすいため、小さなメッシュ（即ち、格子密度を大）を用いて空間の解像度を高め、詳細な流体の流れを調べることが望ましいものとなる。しかしながら、三次元の３つの自由度（３軸）それぞれについて、空間解像度に差を設けて不均一なメッシュで分割して計算を行うと、ＦＦＴで二次元化し、ブロックサイクリックリダクション法といった直接法を用いて解くことができず、Algebraic Multigrid Method（ＡＭＧ法）やFull Multigrid Method（ＦＭＧ法）等の反復法で計算しなければならず、計算時間が著しく増大するという問題がある。他方、三次元の各自由度を均一な構造格子でメッシュ分割すると、十分な解析精度を確保するためには格子数が巨大化し、やはり計算コストが増大するという問題がある。

本発明は、以上のような問題点に鑑み案出なされたもので、流体が流れる空間が３自由度全て不均一にメッシュ分割された第１メッシュを定義して前記コンピュータに入力するステップと、この第１メッシュの３自由度のうち１自由度についてのみ均一にメッシュ分割され他の２自由度については第１メッシュと同一にメッシュ分割された第２メッシュとを定義し、これらの間で計算結果等のデータを受け渡し（カップリング）させることを基本として、詳細な解析が必要な領域は細かく、かつ、詳細な解析が不要な領域は粗く不均一にメッシュ分割しつつも計算コストの増加を抑制しうる三次元の流体シミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、物体の回りの流体の流れを、コンピュータを用いて解析するための三次元の流体シミュレーション方法であって、前記コンピュータに、前記流体が流れる空間が３自由度全て不均一にメッシュ分割された第１メッシュを定義して入力するステップａ、前記第１メッシュの３自由度のうち１自由度についてのみ均一にメッシュ分割され他の２自由度については第１メッシュと同一にメッシュ分割された第２メッシュを定義して前記コンピュータに入力するステップｂ、前記物体が有限個の要素で分割された物体モデルを前記コンピュータに入力するステップｃ、前記第１メッシュ内に前記物体モデルを配置して予め定めた境界条件を与えて前記物体モデルの回りの流体に関する運動方程式を設定して前記流体の速度を計算するステップｄ、前記ステップｄで得られた流体の速度に基づいて各メッシュ内の流量アンバランスを計算するステップｅ、前記流量アンバランスに基づいて前記流体に関する圧力補正方程式を設定するステップｆ、前記得られた流量アンバランスを前記第２メッシュに写像し、前記流体の圧力補正量を計算するステップｇ、及び前記ステップｇで得られた圧力補正量を前記第１メッシュに写像するステップｈを含み、前記ステップｅで得られる流量アンバランスが、収束するまで上記ステップｄ乃至ｈを繰り返すとともに、前記ステップｂにおいて、均一にメッシュ分割される１自由度は、前記流体の流れの方向と同一でないことを特徴とする。

また請求項２記載の発明は、前記写像は、数値補間法を用いて行われることを特徴とする。

また請求項３記載の発明は、前記物体が球体であり、前記流体が空気であることを特徴とする。

また請求項４記載の発明は、前記第１メッシュ及び前記第２メッシュは、ともに直交格子で構成されることを特徴とする。

流体シミュレーションにおいて、非圧縮性流体のNavier-Stokes方程式を、圧力ベースの分離型解法で解く場合、運動方程式と圧力方程式とを計算する必要がある。ここで、運動方程式に関して、メッシュ分割が均一であるか又は不均一であるかは計算時間に殆ど影響がない。一方、圧力方程式に関しては、メッシュ分割を不均一にすると、ＡＭＧ法やＦＭＧ法等で解く必要性が生じ、計算時間が増大する。即ち、圧力方程式の計算に際しては、例えばＦＦＴをかけて次元数を減らし、ブロックサイクリックリダクション法等を用いて解く必要がある。

本発明では、流体に関する運動方程式が、３自由度全てが不均一にメッシュ分割された第１メッシュを使用して構築され、かつ、そこから流体の速度が計算される。また、この運動方程式の計算によって得られた流体の速度に基づく流量アンバランスは、第２メッシュに写像される。第２メッシュは、第１メッシュの１自由度に関して等間隔で分割されている。そして、第２メッシュでは、前記第１メッシュでの計算結果を利用した圧力補正方程式が構築され、これを計算して流体の圧力補正量が求められる。このように、圧力補正方程式が１自由度を均一に分割した第２メッシュで、また運動方程式は３自由度全てが不均一にメッシュ分割された第１メッシュでそれぞれ計算され、かつ、流量アンバランスが収束するまで上記計算が繰り返される。このような手順を用いた本発明の流体シミュレーション方法では、シミュレーションにおけるトータルでの計算時間及び計算コストを低減させることができる。

本発明のシミュレーション方法を実行するためのコンピュータ装置の斜視図である。 本発明のシミュレーション方法の全体フローチャートである。 第１メッシュを視覚化して示す斜視図である。 第２メッシュを視覚化して示す斜視図である。 物体モデルを視覚化して示す図である。 第１メッシュに物体モデルを配置した状態を示す側面図である。 本実施形態の計算処理のさらに詳細なフローチャートである。 （ａ）は第１メッシュのｚ値とそのソース項との関係を示すグラフ、（ｂ）は写像を説明するグラフ、（ｃ）は、第２メッシュに写像されたソース項とｚ値との関係を示すグラフである。 （ａ）は第１メッシュのｚ値とそのソース項との関係を示すグラフ、（ｂ）は他の形態の写像を説明するグラフ、（ｃ）は、第２メッシュに写像されたソース項とｚ値との関係を示すグラフである。 実施例のシミュレーションにおいて、物体モデルと第１メッシュのｚ軸との関係を示す。 実施例のシミュレーションにおいて、第１メッシュのｚ軸方向に関する格子サイズの分布を示すグラフである。 実施例のシミュレーションにおいて、第１メッシュのｚ軸のソース項の値を示すグラフである。 実施例のシミュレーションにおいて、第２メッシュのｚ軸に写像されたソース項の値を示すグラフである。 実施例及び比較例について、ソース項を比較したグラフである。 格子間隔とシミュレーション精度との関係を示すグラフである。 流体の速度を色彩の変化で表してシミュレーションの結果を視覚化して示す側面図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本発明は、物体の回りの流体の流れを、コンピュータを用いて解析するための三次元の流体シミュレーション方法である。

前記流体には、非圧縮性流体、即ち、圧力や温度による密度の変化が無視できるほど小さな流体が用いられる。我々の身の回りにある、空気等の流れに関する現象の多くは、流体を非圧縮性と仮定できるケースが多く、流れ場を支配する方程式は非圧縮性を前提に単純化することができる。一般にマッハ数が０．３以下の条件でシミュレーションが行われる場合、空気は非圧縮性と見なすことができる。また、前記物体には、その外面形状を特定できるものであれば、種々のものが対象となる。本実施形態では、流体に空気が、物体にゴルフボールがそれぞれ適用されている。これにより、本実施形態では、空気中を飛行するゴルフボールの周囲の空気の流れを解析するシミュレーションが行われる。

本実施形態のシミュレーションは、図１に示されるようなコンピュータ装置１を用いて行われる。このコンピュータ装置１は、大規模計算が可能な計算サーバーＳと、該計算サーバーＳに接続されかつ該計算サーバーＳに所定の計算を行わせるクライアントコンピュータ１ａとを含む。クライアントコンピュータ１ａは、キーボード１ｂ、マウス１ｃ、及びディスプレイ装置１ｄを含んで構成されている。クライアントコンピュータ１ａの内部には、ＣＰＵ、ＲＯＭ、作業用メモリー及び磁気ディスク等のなどの大容量記憶装置が設けられる。また、クライアントコンピュータ１ａには、ＣＤ−ＲＯＭやフレキシブルディスクのドライブ装置１ａ１、１ａ２が設けられる。そして、前記大容量記憶装置には後述する本実施形態のシミュレーション方法を実行するためのプログラムが記憶されている。

図２には、本実施形態のシミュレーション方法を実施するためのフローチャートが記載されており、以下、具体的な各ステップについて説明する。

［ステップａ（図２のＳ１）］
ステップａでは、前記コンピュータ装置１に、第１メッシュが定義されかつ入力される。図３には、第１メッシュＭ１が視覚化されて表示される。第１メッシュＭ１は、前記流体が流れる空間の大きさとその物理量の計算位置とを定める。このため、第１メッシュＭ１は、有限空間として定義されるとともに、内部が多数の格子としてメッシュ分割される。流体の物理量は、例えば、前記メッシュの各格子点で得られる。従って、コンピュータ装置１には、このような空間及びメッシュ分割を特定しうるパラメータ（この実施形態では、第１メッシュＭ１のｘ、ｙ及びｚの大きさや、各格子点の座標値など）が入力される。なお、第１メッシュＭ１の大きさ等は、解析対象や解析目的等に応じて種々定めることができる。

本実施形態の第１メッシュＭ１は、直交座標系の中で直方体状に定義されるが、このような態様に限定されるものではなく、三次元の空間を規定しうるものであれば、円柱座標系や極座標系等が採用されても良い。

また、第１メッシュＭ１は、三次元空間の３自由度（この例ではｘ、ｙ及びｚの各軸）全てにおいて不均一にメッシュ分割されている。本実施形態の第１メッシュＭ１は、ｘ軸方向において、格子密度が小さくなるように分割された両側の粗分割領域ｘａと、これよりも格子密度が大きくなるように分割された中央の細分割領域ｘｂとを有する。同様に、第１メッシュＭ１は、ｙ軸方向及びｚ軸後方においても、それぞれ粗分割領域ｙａ、ｚａと、細分割領域ｙｂ、ｚｂとを有する。このような分割により、本実施形態の第１メッシュＭ１は、空間の中央部分の格子密度が大きく構成される。従って、第１メッシュＭ１の中央部分に物体モデル（後述）を配置することにより、物体周辺の流体の流れがより高精度に解析され得る。

［ステップｂ（図２のＳ２）］
ステップｂでは、前記コンピュータ装置１に、第２メッシュが定義されかつ入力される。図４には、第２メッシュＭ２が視覚化されて表示される。第２メッシュＭ２も、流体が流れる空間を定めるとともに、該空間が多数の格子でメッシュ分割されている。従って、コンピュータ装置１には、このような空間及びメッシュ分割を特定しうるパラメータ（この実施形態では、第２メッシュＭ２のｘ、ｙ及びｚの大きさや、各格子点の座標値など）が入力される。

また、第２メッシュＭ２は、第１メッシュＭ１と同一の空間（外形形状）を定めているが、第１メッシュＭ１の３自由度のうち１自由度についてのみ均一にメッシュ分割され他の２自由度については第１メッシュと同一にメッシュ分割されている。本実施形態の第２メッシュＭ２は、ｘ軸方向及びｙ軸方向は、第１メッシュＭ１と同様に不均一な分割がなされている。しかしながら、本実施形態の第２メッシュＭ２は、ｚ軸方向に関しては等間隔でかつ細分割領域ｚｂのピッチで分割されている。このため、第２メッシュＭ２の格子数は、第１メッシュの格子数よりも多い。

［ステップｃ（図２のＳ３）］
ステップｃでは、前記コンピュータ装置１に、物体モデルが定義されかつ入力される。図５には、物体モデルＭ３が視覚化して示される。本実施形態の物体モデルＭ３は、ゴルフボールに近似させて球形で作られている。また、ディンプルが気流に与える影響を調べるためには、該物体モデルＭ３の表面にディンプルに従った凹凸（図示省略）を形成しておくことが望ましい。

［ステップｄ（図２のＳ６及び図７のＳ６１〜ＳＳ６２）］
ステップｄでは、図６に視覚化して示されるように、第１メッシュＭ１内に物体モデルＭ３が配置され、両者の相対位置関係が特定される。即ち、シミュレーションを行うに際して、第１メッシュＭ１の内部に物体モデルＭ３が位置するよう、両者の相対位置関係が定義される。本実施形態では、第１メッシュＭ１と物体モデルＭ３との中心が、揃えられて配置される。

次に、ステップｄでは、予め定めた計算条件が与えられる（図２のＳ５）。計算条件としては、シミュレーションの計算を行う時間ステップΔｔや、流体の密度などの他、メッシュの各壁面での境界条件が含まれる。本実施形態の境界条件として、物体モデル３の表面はNo Slip Wallとした。また第１メッシュＭ１において、流体の入口Ｉには流体の速度ｖが、流体の出口Ｏには圧力０がそれぞれ定義されるとともに、入口Ｉと出口Ｏとの間のメッシュ外面ＷにはそれぞれFull Slip Wallの条件が与えられている。また、計算条件として、例えば、圧力補正方程式を用いて速度場の補正を行う際の繰り返し計算（サブイタレーション）等の収束条件や反復回数なども与えられる。これらの条件は、コンピュータ装置１に入力され、シミュレーションの際に参照される。

さらに、ステップｄでは、図２のステップＳ６の詳細を示す図７から明らかなように、前記コンピュータ装置１において、シミュレーションの時間ステップが更新され（図７のＳ６１）、次に、相対位置関係が定義づけられた第１メッシュＭ１と物体モデルＭ３とを用いて流体に関する運動方程式（下記式（１）参照）が設定され、該流体の速度が計算される（同Ｓ６２）。これにより、第１メッシュＭ１内に、前記条件で流体を流したときに、第１メッシュＭ１内での流体の各速度が計算される。

ここで、ｖは速度ベクトル、ρは流体の密度、ｐは流体の圧力、ｔは時刻、ｇは重力加速度である。

［ステップｅ（図７のＳ６３）］
ステップｅでは、上記ステップｄで得られた流体の速度に基づいて、各メッシュ内の流量アンバランスがコンピュータ装置１により計算される。即ち、各格子（セル）毎に、質量流束（Mass Flux）が計算され、単位時間当たりの入ってくる流量と出て行く流量との差である流量アンバランスΔｍ’が計算される。

［ステップｆ（図７のＳ６４）］
ステップｆでは、コンピュータ装置１が、前記流量アンバランスΔｍ’に基づいて、流体に関する圧力補正方程式を設定する。即ち、コンピュータ装置１は、前記流量アンバランスΔｍ’を各セルの体積で除し、有限差分法で離散化した形式に変換し、下記式（２）の圧力補正方程式を立てる。

ここで、符号は次の通りである。
ｋ：構造格子によって定まる定数
ｐ’：式１のｐの圧力補正量（Pressure Correction）
Δｍ’：流量アンバランス

［ステップｇ（図７のＳ６５）］
ステップｇでは、前記コンピュータ装置１が、各格子（セル）毎に得られた流量アンバランスΔｍ’を前記第２メッシュＭ２の対応する格子（セル）に写像し、解くべきマトリックスを構築して、流体の圧力補正量を計算する。

［写像］
第１メッシュＭ１は、３軸とも不均一に分割されている。一方、本実施形態の第２メッシュＭ２では、２軸（ｘ軸及びｙ軸）のみが第１メッシュＭ１と同様に分割されているものの、残りの１軸（ｚ軸）については、均等に分割されている。従って、第１メッシュＭ１の各格子点の流量アンバランスΔｍ’を第２メッシュＭ２で用いる場合には、第２メッシュＭ２で均等に分割されている軸に関しては、写像（補間）する必要が生じる。

写像は、種々の方法が存在するが、好ましくは数値補間法を用いて行われるのが良い。数値補間法としては、その中でも代表的な線形補間法（linear interpolation）や加重距離法（weighted inverse distance）が好適である。

線形補間法を用いた写像は、図８（ｂ）に示されるように、写像先の第２メッシュＭ２が、写像元の第１のメッシュＭ１よりも細かく分割されている場合に有効である。線形補間法では、図８（ａ）に示されるように、写像元の前記１自由度（ここではｚ軸）における第１メッシュＭ１の隣り合う２つのソース項φ（Δｍ’）の間は、線形に変化すると仮定して補間が行われ、第２メッシュＭ２の前記１自由度（ここではｚ’軸とする）の格子点にその値が割り当てられる。図８（ｃ）には、第２メッシュＭ２に写像された流量アンバランスを白丸のプロットで示している（なおグレーで着色されたプロットは、第１メッシュＭ１のソース項である。）。

また、加重距離法を用いた写像は、図９（ｂ）に示されるように、写像先の第２メッシュＭ２が、写像元の第１メッシュＭ１よりも粗く分割されている場合、かつ、図９（ａ）に示されるように、第１メッシュＭ１の持つ情報（ソース項）がスムースでは無く細かく振動している場合に有効である。このような加重距離法は、第２メッシュＭ２と第１メッシュＭ１との距離の逆数で重み付けをし、関係する情報全てを反映させて補間するものである。例えば下式で補間値を求めることができる。
φz'＝Σ（ωi・φz） （ｉ＝１〜Ｎ）
ωi＝［（１／ｄi）／Σ（１／ｄi）］
ここで、ｄi は、ｚとｚ'との距離である。従って、ｄi＝０の場合、その値がそのまま写像され、他の値はほとんど無視されるという特徴がある。また、図９（ｃ）には、第２メッシュＭ２に写像された流量アンバランスを白丸のプロットで示している（なおグレーで着色されたプロットは、第１メッシュＭ１のソース項である。）。

［解くべきマトリックス］
第２メッシュＭ２の各格子に流量アンバランスΔｍ’が写像されると、コンピュータ装置１は、第２メッシュＭ２で解くべきマトリックスが設定される。解くべきマトリックスは、例えば下式（３）のように設定される。

ここで、符号は次の通りである。
Ｔ_ｉ（ｉ＝１〜ｎ）：各スライスの持つマトリックス
ａ_ｉ、ｂ_ｉ（ｉ＝１〜ｎ）：カップリング項
Φ_ｉ（ｉ＝１〜ｎ）：求めたい解ｐ’
ｒ_ｉ（ｉ＝１〜ｎ）：−Δｍ’
ｎ：図４のＺ軸の格子分割数
ここで、上記スライスとは、計算領域を、並列計算する個数分だけ分割するときの一つを指す（なお、全計算領域で解く必要がある方程式は式（３）であるが、分割された領域一つ一つで個々に計算する方程式は後述の式（５）である。）。

［流体の圧力の計算（図７のＳ６６）］
第２メッシュＭ２での流体の圧力の計算は、コンピュータ装置１が上記マトリックスをフーリエ変換し、複数の対角行列の集まりに分割する（下式（４）参照）。式（４）と式（３）の比較から明らかなように、式（４）では、フーリエ変換によって式（３）で存在していたカップリング項ａ_ｉ、ｂ_ｉを無くすことができる。この結果、式（４）では、各スライスにおいて、他のスライスを参照することなく独立して計算することができる。そして、コンピュータ装置１は、分割された各行列式（下式（５）参照）を、例えばブロックサイクリングリダクション法等の直接法を用いて計算することができる。そして、コンピュータ装置１は、計算された解Φを逆フーリエ変換し、元の未知数である圧力ｐ’へと戻す。

ここで、符号は次の通りである。
Ｔ〜i（ｉ＝１〜ｎ）：各スライスの持つマトリックス
Φ〜i（ｉ＝１〜ｎ）：求めたい解（ｐ’）にＦＦＴをかけたもの
ｒ〜i（ｉ＝１〜ｎ）：−Δｍ’にＦＦＴをかけたもの
ｎ：図４のｚ軸の格子分割数

［ステップｈ（図７のＳ６７、６８）］
ステップｈでは、前記コンピュータ装置１が、上記ステップｇで得られた第２メッシュＭ２での圧力補正量ｐ’を、第１メッシュＭ１に写像する。ここでの写像も、先に述べた方法が好適に使用される。これにより、第１メッシュＭ１での流体の圧力が、第２メッシュＭ２の圧力補正方程式によって計算された前記圧力補正量ｐ’（ｐの修正値）によって更新される。

［ループステップ］
次に、コンピュータ装置１は、ステップｅで得られる流量アンバランスΔｍ’が収束するまで、上記ステップｄ乃至ｈを繰り返して計算する。即ち、コンピュータ装置１は、流量アンバランスΔｍ’が予め定めた収束条件（例えば前回計算値との差が閾値以下）を満たしたと判断した場合（Ｓ７０でＹ）、図２のＳ７へと戻る。他方、コンピュータ装置１は、流量アンバランスΔｍ’が収束条件を満たしていないと判断した場合（Ｓ７０でＮ）、収束条件を満たすまで再度、Ｓ６２〜Ｓ６８までを繰り返し実行する。なお、このループステップの収束条件は、種々設定することができ、例えばループの繰り返し回数等によって規定されても良い。ステップｄ乃至ｈの繰り返しは、例え、流量アンバランスが完全に取り除かれたとしても、一度、圧力が修正されれば、新たに補正された速度は、ステップｄの運動方程式を確認していないという事実から必要になる。そして、ステップｄ乃至ｈの繰り返しによって、圧力と速度の解の集合として、運動方程式および質量保存方程式(連続方程式)の両方が満たされるものが得られる。

［時間経過］
次に、コンピュータ装置１は、予め定めた所定の時間又は予め定めた計算ステップ数が経過しているか否かを判断し、真の場合（Ｓ７でＹ）には、処理を終えるが、偽（Ｓ７でＮ）の場合には、図７の処理を繰り返す。

以下、本発明のより具体的な実施例について述べるが、本発明はこのような実施例に限定して解釈されるものではない。

この実施例では、図３及び図４に示した第１メッシュ及び第２メッシュを使用してゴルフボールの空力シミュレーションが行われた。

図１０には、第１メッシュのｚ軸と、物体モデルＭ３との位置関係を示す。第１メッシュのｚ方向の巾は、０．２ｍである。また、物体モデルＭ３の直径は４２．７mmである。物体モデルＭ３は、第１メッシュのｚ軸方向の中央部に配置される。

図１１には、第１メッシュの格子サイズの分布を示し、横軸は、第１メッシュのｚ軸、縦軸は、格子サイズをそれぞれ示している。図１１から明らかなように、物体モデルＭ３の周辺で格子サイズが最も小さくなっており、そこから外側に向かって格子サイズは漸増している。ｚ軸の０．１〜−０．１ｍの領域の格子点は、１０２４個である。

これらのメッシュＭ１及び物体モデルＭ３を使用して、下記式（６）のポアソン方程式が計算された。

ただし、ソース項（ｒ）は、図１２に示されるような値をとっているものとする。このソース項は、本実施形態の様に、物体モデルＭ３が配置された場合を想定し、ゴルフボール（物体モデルＭ３）が存在する周辺に、ランダムに値を発生させた。

次に第２メッシュが設定された。第２メッシュのｚ軸の領域は、−０．１〜０．１ｍであり、この領域に７００個の格子点が均等に設けられた。従って、格子間隔は、約０．０００２８６ｍ（＝（０．１−（−０．１））／７００）となる。この実施形態では、加重距離法を用いて、第１メッシュから第２メッシュへと写像が行われた。その結果は図１３に示される。

次に、ｚ軸が等分割された上記第２メッシュについて、圧力補正方程式（ポアソン方程式）が計算され、求まった解が、第１メッシュに写像されて戻された。なお、圧力補正方程式の計算には、上述の通り、ＦＦＴを使用して次元数を落とし、ブロックサイクリックリダクション法が好適に用いられた。また、得られた解を第１メッシュに戻す写像には、加重距離法が採用された。

図１４には、本発明によって計算（本発明方法）された解が、写像することなく第１メッシュのみを使用して計算（以下、「直接法」という。）された解と、比較して表示されている。この様に、本発明の計算結果と、従来の直接法で求めた計算結果とは、ほぼ一致する事が確認できる。

次に、第２メッシュの格子密度と精度の関係について述べる。本実施例の場合、格子サイズ（間隔）（単位：ｍ）は、格子数をＮとすると、下式で表示することができる。
格子間隔＝０．２／Ｎ
ここでは、Ｎ＝１５０、３００、５００、１０００及び２０００ の場合、精度がどのように変化するかが確認された。結果は、図１５に示す通りである。

図１５から明らかなように、バックグランドメッシュの格子密度が高い（格子間隔が小さい）ほど、直接法の結果に近づいており、解の精度が向上する事が分かる。

図１６には、上記シミュレーションを実行した結果に基づき、流体の流速を求め、その大きさを色彩の変化を用いて視覚化した結果が示されている。

以上説明したように、本発明によれば、計算効率の向上効果が期待できる。具体的な例を挙げると、ｚ軸方向に等間隔に格子を配置しなければならない場合、直径０．０４２７ｍのゴルフボールをｚ軸の中心に配置すると、これまでの経験則上、解析精度を維持するためには、格子間隔は７．５０Ｅ^−５（ｍ）程度にする必要がある。一方、解析領域を−０．２５〜＋０．２５ｍと仮定すると、Ｚ方向の分割数は、６６７０個になる。ｘ，ｙ軸方向は、不等分割で格子を配置できるため、通常１０００個程度の分割数で十分である。よって、このようなメッシュの格子の総数は、６６７０×１０２４×１０２３＝６．９８億個になる。つまり、従来の方法では、約７億個の格子（セル）に対して運動方程式を解く必要があり、コンピュータ装置のメモリー容量や計算時間も膨大になる。しかしながら、本発明の第１メッシュは、ｚ方向も不等間隔で配置できるため、不必要な解析場を粗く分割することによって、ｚ方向も１００個程度の分割数に抑えることができる。このようなメッシュでは、格子の総数は、１０００×１０２４×１０２３＝１．０５億個と大幅に低減させる事ができ、計算規模を１／７程度まで圧縮することが可能である。

Ｍ１ 第１メッシュ
Ｍ２ 第２メッシュ
Ｍ３ 物体モデル

Claims (4)

1. 物体の回りの流体の流れを、コンピュータを用いて解析するための三次元の流体シミュレーション方法であって、
   前記コンピュータに、前記流体が流れる空間が３自由度全て不均一にメッシュ分割された第１メッシュを定義して入力するステップａ、
   前記第１メッシュの３自由度のうち１自由度についてのみ均一にメッシュ分割され他の２自由度については第１メッシュと同一にメッシュ分割された第２メッシュを定義して前記コンピュータに入力するステップｂ、
   前記物体が有限個の要素で分割された物体モデルを前記コンピュータに入力するステップｃ、
   前記第１メッシュ内に前記物体モデルを配置して予め定めた境界条件を与えて前記物体モデルの回りの流体に関する運動方程式を設定して前記流体の速度を計算するステップｄ、
   前記ステップｄで得られた流体の速度に基づいて各メッシュ内の流量アンバランスを計算するステップｅ、
   前記流量アンバランスに基づいて前記流体に関する圧力補正方程式を設定するステップｆ、
   前記得られた流量アンバランスを前記第２メッシュに写像し、前記流体の圧力補正量を計算するステップｇ、及び
   前記ステップｇで得られた圧力補正量を前記第１メッシュに写像するステップｈを含み、
   前記ステップｅで得られる流量アンバランスが、収束するまで上記ステップｄ乃至ｈを繰り返すとともに、
   前記ステップｂにおいて、均一にメッシュ分割される１自由度は、前記流体の流れの方向と同一でないことを特徴とする三次元の流体シミュレーション方法。
   
2. 前記写像は、数値補間法を用いて行われる請求項１記載の三次元の流体シミュレーション方法。
3. 前記物体が球体であり、前記流体が空気である請求項１又は２記載の三次元の流体シミュレーション方法。
4. 前記第１メッシュ及び前記第２メッシュは、ともに直交格子で構成される請求項１乃至３のいずれかに記載の三次元の流体シミュレーション方法。