JP5750091B2

JP5750091B2 - 流体シミュレーション方法

Info

Publication number: JP5750091B2
Application number: JP2012252498A
Authority: JP
Inventors: 角田　昌也; 昌也角田; ヤダフアルジュン
Original assignee: Sumitomo Rubber Industries Ltd
Current assignee: Sumitomo Rubber Industries Ltd
Priority date: 2012-11-16
Filing date: 2012-11-16
Publication date: 2015-07-15
Anticipated expiration: 2032-11-16
Also published as: US9449125B2; US20140142908A1; JP2014102574A

Description

本発明は、固体モデルと流体モデルとの間の境界面に作用する流体モデルのせん断応力を、短時間でかつ精度良く求めることができる流体シミュレーション方法に関する。

近年、様々な場面において、コンピュータを用いた流体シミュレーションが行われている。流体シミュレーションの一例としては、ゴルフボールの周囲の空気の流れを解析するものや、ゴム押出機の内部のゴムの流れを解析するもの等が知られている。このような流体シミュレーションは、飛行特性に優れたゴルフボールのディンプルの開発や、押出し抵抗の少ない押出機ないしゴム配合の開発等に役に立つ。

流体シミュレーションにおいては、例えば、図１２に示されるように、固体をモデル化した固体モデルａと、該固体の外側で流体が流れる空間を、有限個の要素ｃでモデル化した流体モデルｂとが定義される。また、流体モデルｂには、固体モデルａが配置される固体領域ｂ１と、流体モデルｂ側の流体領域ｂ２との境界面ｅが定義される。さらに、流体モデルｂの各要素ｃには、流体解析用の重心点ｄが設けられている。これらの重心点ｄでは、流体の速度、圧力、又は温度が計算される。そして、これらの物理量に基づいて、例えば、境界面ｅに作用する流体モデルｂのせん断応力等が計算される。関連する文献としては次のものがある。

特開２０１１−２１５８２３号公報

DARTZI PAN and TZUNG-TZA SHEN 著、「 A Ghost Cell Method for the Computation of Incompressible Flows with Immersed Bodies 」、 6th IASME/WSEAS International Conference on FLUID MECHANICS and AERODYNAMICS (FMA'08)Rhodes, Greece, August 20-22, 2008

しかしながら、従来の流体シミュレーション方法では、境界面ｅに作用する流体モデルｂのせん断応力の計算過程が、複雑であったため、計算時間の増大や、計算精度の低下を招きやすいという問題があった。従って、流体モデルｂのせん断応力を、短時間、かつ、精度良く求めることができる流体シミュレーション方法が強く求められていた。

本発明は、以上のような実状に鑑み案出されたもので、固体モデルと流体モデルとの間の境界面に作用する流体モデルのせん断応力を、短時間、かつ、精度良く求めることができる流体シミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、固体をモデル化した固体モデルと、前記固体の外側で流体が流れる空間を有限個の要素でモデル化した流体モデルとを用いて、前記固体の周囲での前記流体の流れを解析する流体シミュレーション方法であって、前記流体モデルの各要素には、流体解析用の重心点が設けられ、前記コンピュータが、前記流体モデルにおいて、前記固体モデルが配置される固体領域と、前記固体領域の外側の流体領域との境界面を定義するステップ、前記コンピュータが、前記境界面を挟んで隣り合う前記固体領域の前記要素と、前記流体領域の前記要素との間において、前記各重心点を連結する直線を定義するステップ、前記コンピュータが、前記境界面と前記直線との交点を定義するステップ、前記コンピュータが、前記流体領域の前記要素の前記重心点と、前記交点との間の距離を計算するステップ、及び、前記コンピュータが、前記境界面に作用する前記流体モデルのせん断応力を、下記式（１）を用いて計算するステップを含むことを特徴とする。

ここで、
τ_correct：境界面に作用する流体モデルのせん断応力
Ｕ_live：流体領域の要素の重心点での速度
Ｕ_boundary：交点での速度
dist _{live cell to boundary}：流体領域の要素の重心点と交点との間の距離
μ：流体のせん断粘度
Ａ：境界面の面積

また、請求項２記載の発明は、上記式（１）の dist _{live cell to boundary}には、下記式（２）で修正されるdist’ _{live cell to boundary}が用いられる請求項１に記載の流体シミュレーション方法である。

ここで、
Ｌ１：直線の長さ
α：直線の長さＬ１に対するdist _{live cell to boundary}の比

また、請求項３記載の発明は、前記境界面には、圧力境界条件としてノイマン境界条件が設定される請求項１又は２に記載の流体シミュレーション方法である。

本発明の流体シミュレーション方法は、コンピュータが、境界面を挟んで隣り合う固体領域の要素と、流体領域の要素との間において、各重心点を結ぶ直線を定義するステップ、及び、境界面と直線との交点を定義するステップを含んでいる。

さらに、流体シミュレーション方法は、コンピュータが、流体領域の要素の重心点と、交点との間の距離を計算するステップ、及び、境界面に作用する流体モデルのせん断応力を、下記式（１）を用いて計算するステップを含んでいる。

ここで、
τ_correct：境界面に作用する流体モデルのせん断応力
Ｕ_live：流体領域の要素の重心点での速度
Ｕ_boundary：交点での速度
dist _{live cell to boundary}：流体領域の要素の重心点と交点との間の距離
μ：流体のせん断粘度
Ａ：境界面の面積

本発明の流体シミュレーション方法では、境界面を挟んで隣り合う固体領域の要素の重心点と、流体領域の要素の重心点との間で連結される直線を基準にして、境界面に作用する流体モデルのせん断応力を求めている。従って、本発明の流体シミュレーション方法では、計算が複雑になるのを防ぐことができるため、流体モデルのせん断応力を短時間で、かつ、精度良く求めることができる。

本実施形態の流体シミュレーション方法を実行するコンピュータの斜視図である。本実施形態の流体シミュレーション方法を示すフローチャートである。固体モデルを視覚化して示す正面図である。流体モデルを視覚化して示す斜視図である。図４の側面図である。本実施形態の流体シミュレーションステップを示すフローチャートである。（ａ）、（ｂ）は、従来のせん断応力の計算方法を説明する側面図である。（ａ）、（ｂ）は、本実施形態のせん断応力の計算方法を説明する側面図である。（ａ）は、αとαの逆数（１／α）との関係を示すグラフ、（ｂ）は、αとα_newの逆数との関係を示すグラフである。 αとα_newとの関係を示すグラフである。（ａ）は、実施例のシミュレーション結果を視覚化した側面図、（ｂ）は、比較例のシミュレーション結果を視覚化した側面図である。従来の流体シミュレーション方法を説明する側面図である。

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

図１に示されるように、コンピュータ１は、例えば、大規模計算が可能な計算サーバー１Ｓと、該計算サーバー１Ｓに接続されかつ該計算サーバー１Ｓに所定の計算を行わせるクライアントコンピュータ１Ｃとを含んで構成されている。

また、クライアントコンピュータ１Ｃは、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含む。この本体１ａには、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２などが設けられる。記憶装置には、本実施形態の流体シミュレーション方法を実行するための処理手順（プログラム）が予め記憶されている。

図２には、本発明の流体シミュレーションの処理手順の一例が示される。以下、順に説明する。

流体シミュレーション方法では、先ず、コンピュータ１に、固体モデルが入力される（ステップＳ１）。図３に示されるように、本実施形態の固体モデル２は、固体（非流体物）としてのゴルフボールを、数値解析法により取り扱い可能な有限個の要素５でモデル化（離散化）されたゴルフボールモデル２Ａからなる。数値解析法としては、例えば、有限要素法等が採用される。また、ゴルフボールモデル２Ａには、ゴルフボールの表面に形成されるディンプルが形成されるのが望ましい。

本実施形態の固体モデル２の要素５としては、変位を未知数とするラグランジュ（Lagrange）要素が用いられている。また、各要素５には、要素番号、節点番号、全体座標系Ｘ−Ｙ−Ｚの節点座標値及び材料特性（例えば密度、ヤング率及び／又は減衰係数等）などの数値データが定義される。これらの数値は、コンピュータ１に記憶される。

次に、コンピュータ１に、流体モデルが入力される（ステップＳ２）。図４に示されるように、流体モデル３は、固体の外側で流体（空気）が流れる空間を、数値解析法により取り扱い可能な有限個の要素６でモデル化（離散化）することによって定義されている。

本実施形態の流体モデル３は、固体モデル２を完全に覆う大きさを有した直方体状に設定されている。また、流体モデル３は、ｘ軸方向の一端に配置される前面３ａ、ｘ軸方向の他端側に配置される後面３ｂ、及び、前面３ａと後面３ｂとの間をのびる側面３ｓが定義されている。

流体モデル３の各要素６は、直方体状に設定された流体モデル３において、ｘ軸方向にｓ分割、ｙ軸方向にｍ分割、及び、ｚ軸方向にｎ分割した直交格子状に定義されている。ｓ、ｍ及びｎは、いずれも自然数である。これらの自然数ｓ、ｍ及びｎは、モデルの大きさ、解析の精度、又は、流速等に応じて適宜定められる。例えば、ゴルフボールモデル２Ａの外径が４２．７mmの場合には、各軸ｘ、ｙ、ｚにおいて、要素６の間隔Ｗ１が２．５×１０^−２〜１mm程度に設定されるのが望ましい。

図５に示されるように、本実施形態の流体モデル３の各要素６としては、流体解析用の重心点８を有するオイラー（Euler）要素からなる。この要素６には、空気の初期の比重や、粘性、圧力及び温度といったパラメータが割り当てられる。また、各重心点８では、流体シミュレーションにおいて、空気の速度、圧力又は温度等が計算される。

次に、固体モデル２の計算条件が、コンピュータ１に入力される（ステップＳ３）。この計算条件には、固体モデル２の流体モデル３に対する位置情報、固体モデル２の回転を定義する角速度Ｒ１等が含まれる。

また、固体モデル２の位置情報は、固体モデル２の中心２ｃの座標と、流体モデル３の中心３ｃの座標とが一致するように設定させている。これにより、固体モデル２が流体モデル３の中央に配置される状態を定義することができる。これらの計算条件は、コンピュータ１に記憶される。

次に、流体モデル３の計算条件が、コンピュータ１に入力される（ステップＳ４）。この計算条件には、流体モデル３の境界条件、圧力補正方程式を用いた速度場計算の収束条件、及び、反復回数等が含まれる。

また、流体モデル３の境界条件の設定では、流体モデル３の前面３ａに、流体モデル３の要素６の流入（流速Ｖ１）が定義される。また、流体モデル３の後面３ｂには、前面３ａでの要素６の流入と同量の流出量になるように繰り返し収束計算される。さらに、後面３ｂでの要素６は、全ての変数（圧力、速度等）が拡散しないものとして定義される。また、流体モデル３の側面３ｓには、流体モデル３の内部への要素６の流入出が不能に定義されている。これらの計算条件は、コンピュータ１に記憶される。

次に、コンピュータ１が、固体モデル２の周囲での流体モデル３の流れを解析する流体シミュレーションステップＳ５が行われる。図６には、本実施形態の流体シミュレーションステップＳ５の具体的な処理手順が示されている。

本実施形態の流体シミュレーションステップＳ５では、先ず、図５に示されるように、流体モデル３において、固体モデル２が配置される固体領域Ｔ１と、固体領域Ｔ１の外側の流体領域Ｔ２との境界面７が、コンピュータ１によって定義される（ステップＳ５１）。

本実施形態の境界面７は、流体モデル３において、固体モデル２の表面を特定する座標値である。このような座標値は、固体モデル２の表面を構成する各要素５（図３に示す）の座標値を追跡して設定される。これらの座標値は、コンピュータ１に記憶される。

次に、流体モデル３の物理量が計算される（ステップＳ５２）。このステップＳ５２では、固体モデル２が配置される固体領域Ｔ１を除く流体領域Ｔ２の全ての要素６を対象に、流体モデル３の物理量が計算される。

本実施形態では、流体モデル３の運動が、Immersed Boundary法に基づいて、ナビエ・ストークスの式によって表される。このナビエ・ストークスの式は、例えば、コンピュータ１で計算可能な近似式に変換して計算されることにより、空気の運動状態を表すパラメータ、即ち、流体モデル３の各要素６の重心点８での圧力や速度等が計算される。なお、シミュレーションの単位時間Ｔｘの増加により、固体領域Ｔ１から流体領域Ｔ２に変化した要素６は、該要素６に隣接する流体領域Ｔ２の要素６の圧力を、固体モデル２の移動速度の初期値として計算される。

このような流体モデル３の運動状態の計算は、例えば、特開２０１２−８３９５８号公報に記載されている「構造格子を用いたシミュレーション方法」に準じて計算することができるが、ANSYS社のFLUNETなどの市販の流体解析用のアプリケーションソフトを用いて計算することもできる。

なお、Immersed Boundary法に基づく流体シミュレーションでは、固体領域Ｔ１にも圧力場が存在するものとして、固体領域Ｔ１の要素６ａの全範囲を対象に、圧力補正方程式が適用される。このため、固体領域Ｔ１の圧力場が、流体領域Ｔ２の圧力場に影響し、シミュレーション精度が低下しやすい。このため、境界面７には、固体領域Ｔ１内の圧力勾配がゼロになるように、ノイマン境界条件が設定されるのが望ましい。本実施形態では、先ず、固体領域Ｔ１の要素６ａにノイマン条件を設定する。さらに、境界面７が横切る要素６ｃでは、流体領域Ｔ２側の境界面７に、ノイマン境界条件が定義される。

ノイマン境界条件は、隣り合う要素６、６の微分値に対して、拘束条件を与える（例えば、要素６、６の微分値の差をゼロにする）ものである。これにより、本実施形態では、固体領域Ｔ１を計算対象から除外して、流体領域Ｔ２のみを計算対象とすることができる。従って、本実施形態では、固体領域Ｔ１の圧力場が、流体領域Ｔ２の圧力場に影響するのを抑制でき、シミュレーション精度を高めることができる。

ところで、流体モデル３の物理量のうち、境界面７に作用する流体モデル３のせん断応力を求める計算方法としては、種々提案されている。上記非特許文献１に基づく従来の方法では、例えば、下記式（３）が用いられる。なお、固体モデル２及び流体モデル３が三次元モデルからなる場合には、ｘ軸、ｙ軸及びｚ軸の各方向において、せん断応力τが求められる。

ここで、
μ：せん断粘度
Ａ：境界面の面積
ｖ：流体の速度
ｎ：法線ベクトル

上記式（３）を用いるには、境界面７から垂直にのびる法線ベクトルにおいて、境界面７での速度をゼロとした流体モデル３の速度勾配（ｄｖ／ｄｎ）を求める必要がある。この速度勾配（ｄｖ／ｄｎ）を求めるには、先ず、図７（ａ）に示されるように、境界面７に隣接する流体モデル３の要素６、６をのび、かつ、境界面７に直交する垂線１２が定義される。

この垂線１２の一端には、流体領域Ｔ２の要素６に配置される流体側点１３が定義される。さらに、垂線１２の他端には、固体領域Ｔ１の要素６の重心点８によって定義された固体側点１４が設定される。また、境界面７には、垂線１２との交点１５が定義されている。

次に、図７（ｂ）に示されるように、交点１５から流体側点１３までの速度勾配が、下記式（４）によって求められる。なお、流体側点１３及び交点１５は、流体領域Ｔ２の要素６の重心点８によって定義されていない。このため、流体側点での流体の速度Ｖ_i、及び、交点での流体の速度Ｖ_pは、それぞれの流体側点１３及び交点１５において、隣接する複数の要素６の重心点８で計算された速度から推定される。

ここで、
Ｖ_i：流体側点での流体の速度
Ｖ_p：交点での流体の速度
δ：流体側点から交点までの距離

次に、下記式（５）を用いて、固体側点１４での速度Ｖ_gが求められる。下記式（５）では、速度Ｖ_gが、固体側点１４の位置ベクトルｒ_gと、交点１５の位置ベクトルｒ_pとの差（ｒ_g−ｒ_p）から、上記式（４）で求めた速度勾配を用いて計算されている。そして、速度Ｖ_gの値を変更して、速度Ｖ_pを相対的にゼロして、境界面７での速度をゼロとした流体の速度勾配（ｄｖ／ｄｎ）が求められる。

ここで、
Ｖ_g：固体側点での流体の速度
Ｖ_p：交点での流体の速度
ｒ_g：固体側点の位置ベクトル
ｒ_p：交点の位置ベクトル

このように、従来の方法では、速度Ｖ_i及び速度Ｖ_pを、複数の要素６の重心点８の速度から推定したり、速度Ｖ_gを変更して速度Ｖ_pを相対的にゼロにしたりする必要がある。このため、従来の方法では、計算過程が複雑であったため、計算時間の増大や、計算精度の低下を招きやすいという問題があった。

本実施形態では、流体モデル３のせん断応力を計算するのに先立ち、図８（ａ）に示されるように、境界面７を挟んで隣り合う固体領域Ｔ１の要素６ａと、流体領域Ｔ２の要素６ｂとの間において、各重心点８、８を結ぶ直線９が定義される（ステップＳ５３）。

このステップＳ５３では、先ず、固体領域Ｔ１側の要素６ａと、流体領域Ｔ２の要素６ｂとの間において、境界面７を挟んで最も隣接する要素６ａ、６ｂの重心点８、８が探索される。そして、それらの重心点８、８間に、境界面７と交わる直線９が定義される。この直線９は、各重心点８、８の座標を特定する数値データであり、コンピュータ１に記憶される。

また、最も隣接する要素６ａ、６ｂの重心点８、８を探索する方法としては、特に限定されるわけではないが、例えば、特開２０１１−１０８１１３号公報に開示されている「任意の実数群から探索点の最近傍値を探索する方法」が好適に採用される。これにより、最も隣接する要素６ａ、６ｂの重心点８、８を効率的に探索することができる。

次に、境界面７と直線９との交点１０が、コンピュータ１によって定義される（ステップＳ５４）。この交点１０は、流体モデル３において特定される座標値である。この座標値は、コンピュータ１に記憶される。

次に、図８（ｂ）に示されるように、境界面７を挟んで最も隣接する要素６ａ、６ｂにおいて、流体領域Ｔ２の要素６ｂの重心点８と、交点１０との間の距離Ｌ２が、コンピュータ１によって計算される（ステップＳ５５）。この距離Ｌ２は、流体領域Ｔ２の要素６ｂの重心点８と、交点１０との間において、直線９上で特定される距離である。この距離Ｌ２は、コンピュータ１に記憶される。

次に、境界面７に作用する流体モデル３のせん断応力τ_correctが、下記式（１）を用いて計算される（ステップＳ５６）。本実施形態のように、固体モデル２及び流体モデル３が三次元モデルからなる場合には、ｘ軸、ｙ軸及びｚ軸の各方向において、せん断応力τ_correctが求められる。

流体領域Ｔ２の重心点８での速度Ｕ_liveは、ステップＳ５２で実施される流体モデル３の物理量の計算結果に基づいて設定される。また、交点１０での速度Ｕ_boundaryは、交点１０に隣接する流体領域Ｔ２の要素６ｂにおいて、少なくとも一つの重心点８で計算された速度から推定される。さらに、距離dist _{live cell to boundary}は、ステップＳ５５で計算された重心点８と交点１０との間の距離Ｌ２である。

また、流体のせん断粘度μは、解析対象の流体が持つ予め定められたせん断粘度が設定される。本実施形態では、空気のせん断粘度が用いられる。境界面７の面積Ａは、図８（ａ）に示されるように、境界面７が要素６に交わる断面１７で特定される。

このように、本実施形態のステップＳ５６では、境界面７を挟んで隣り合う固体領域Ｔ１の要素６ａの重心点８と、流体領域Ｔ２の要素６ｂの重心点８との間で連結される直線９を基準にして、境界面７に作用する流体モデル３のせん断応力τ_correctが求められている。このため、本発明では、図７（ａ）に示される従来の方法のように、境界面７に交わる垂線１２や、流体側点１３を定義する必要がない。

さらに、本発明では、従来のように、隣接する流体モデル３の重心点８の速度から流体側点１３の速度Ｖ_iを推定したり、固体側点１４での速度Ｖ_gを変更して、交点１０での速度Ｖ_pを相対的にゼロにしたりする必要もない。従って、本発明では、計算が複雑になるのを防ぐことができるため、流体モデル３のせん断応力を短時間で、かつ、精度良く求めることができる。

なお、上記式（１）では、分母である距離dist _{live cell to boundary}の値が非常に小さくなると、せん断応力τ_correctが過度に大きくなり、計算が不安定になるおそれがある。最悪の場合には、ゼロ割りが発生し、計算が発散するおそれがある。このため、上記式（１）の距離dist _{live cell to boundary}は、下記式（２）で修正される距離dist’ _{live cell to boundary}が用いられるのが望ましい。

ここで、
Ｌ１：直線の長さ
α：直線の長さＬ１に対するdist _{live cell to boundary}の比（０≦α≦１）

図９（ａ）に示されるように、距離dist _{live cell to boundary}の比であるαがゼロに収束すると、αの逆数１／αが∞に発散する。これは、上記式（１）において、距離dist _{live cell to boundary} が非常に小さくなると、せん断応力τ_correctが過度に大きくなることを示している。

一方、図９（ｂ）に示されるように、上記式（２）においてαの値が修正されたα_newでは、その逆数１／α_newが１２付近で収束している。これは、距離dist _{live cell to boundary}が非常に小さくても、直線９の長さＬ１にα_newを乗じた距離dist’ _{live cell to boundary}が、上記式（１）に代入されることにより、せん断応力τ_correctが過大になるのを防止できることを示している。従って、上記式（２）で修正される距離dist’ _{live cell to boundary}が、上記式（１）の距離dist _{live cell to boundary}に代えて用いられることにより、計算安定性を高めることができる。

なお、上記式（２）のα_newについては、適宜変更することができるが、図１０に示されるように、α_newとαとの関係が、αの増加とともにα_newが増加する単調増加の連続関数１８となるように設定されるのが望ましい。さらに、連続関数１８は、α_new＝αを示す直線１９に近似するのが望ましい。これにより、α_newとαとの誤差を最小限に抑えることができるため、計算制度を維持することができる。

次に、コンピュータ１が、予め定められた終了時間が経過したか判断する（ステップＳ５７）。このステップＳ５７では、終了時間が経過したと判断した場合、流体シミュレーションステップＳ５を終了させ、流体モデル３の物理量が出力される（ステップＳ６）。

一方、コンピュータ１が、終了時間が経過していないと判断した場合には、単位時間Ｔｘを一つ進めて（ステップＳ５８）、前記ステップＳ５１〜Ｓ５６が行われる。これにより、流体シミュレーションの開始から終了までの流体モデル３の物理量を、単位時間Ｔｘ毎に取得することができる。

次に、流体モデル３の物理量が許容範囲内であるかが判断される（ステップＳ７）。このステップＳ７では、流体モデル３の物理量が許容範囲内であると判断された場合、上記固体モデル２（ゴルフボールモデル２Ａ）に基づいて固体（ゴルフボール）が設計される（ステップＳ８）。一方、物理量が許容範囲内でないと判断された場合は、固体モデル２が設計変更され（ステップＳ９）、再度流体シミュレーションが行われる（工程Ｓ１〜Ｓ７）。

このように、本実施形態の流体シュミュレーション方法では、流体モデル３の物理量が許容範囲内になるまで、固体モデル２が設計変更されるため、性能の優れた固体（ゴルフボール）を効率良く設計することができる。

以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。

図２に示した手順に従って、図３に示した固体モデル（ゴルフボールモデル）、及び、図４に示した流体モデルがコンピュータに設定された。そして、図６に示した手順に従って、境界面に作用する流体モデルのせん断応力が計算された（実施例）。図１１（ａ）には、実施例のシミュレーション結果を視覚化した側面図が示されている。

また、比較のために、上記式（３）〜（５）を用いた従来の方法（ゴーストセル法）に基づいて、境界面に作用する流体モデルのせん断応力が計算された（比較例）。図１１（ｂ）には、比較例のシミュレーション結果を視覚化した側面図が示されている。なお、流体モデルの運動状態の計算（ソルバー）は、特開２０１２−８３９５８号公報に記載されている「構造格子を用いたシミュレーション方法」に準じて実施された。また、パラメータ等の共通仕様は、次のとおりである。
コンピュータ：ＳＧＩ社のワークステーションＸＥ３４０
ＣＰＵのコア数：１６０コア
搭載メモリ：１．４４ＴＢ（７２GB／ノード）
処理方法：ＭＰＩによる並列処理
固体モデル（ゴルフボールモデル）：
外径：４２．７mm
角速度Ｒ１：２５００rpm
流体モデル
要素の間隔Ｗ１：０．０８mm（ゴルフボール周辺）
流速Ｖ１：５８ｍ／Ｓ

テストの結果、実施例の計算時間が、比較例の計算時間の約８０％であった。従って、本発明の流体シミュレーション方法では、従来の流体シミュレーション方法に比べて、せん断応力の計算を短時間で計算できることが確認できた。

また、図１１（ａ）に示されるように、実施例の流体シミュレーション方法では、バックスピンがかかったゴルフボールに作用する揚力を忠実に再現できることが確認できた。一方、図１１（ｂ）に示されるように、比較例の流体シミュレーション方法では、バックスピンがかかったゴルフボールに作用する揚力を十分に再現できないことが確認できた。従って、本発明の流体シミュレーション方法では、従来の流体シミュレーション方法に比べて、せん断応力を精度良く求めうることを確認できた。

２固体モデル
３流体モデル
５要素
６要素
７境界面
８重心点
９直線
１０交点
Ｔ１固体領域
Ｔ２流体領域

Claims

固体をモデル化した固体モデルと、前記固体の外側で流体が流れる空間を有限個の要素でモデル化した流体モデルとを用いて、前記固体の周囲での前記流体の流れを解析する流体シミュレーション方法であって、
前記流体モデルの各要素には、流体解析用の重心点が設けられ、
前記コンピュータが、前記流体モデルにおいて、前記固体モデルが配置される固体領域と、前記固体領域の外側の流体領域との境界面を定義するステップ、
前記コンピュータが、前記境界面を挟んで隣り合う前記固体領域の前記要素と、前記流体領域の前記要素との間において、前記各重心点を連結する直線を定義するステップ、
前記コンピュータが、前記境界面と前記直線との交点を定義するステップ、
前記コンピュータが、前記流体領域の前記要素の前記重心点と、前記交点との間の距離を計算するステップ、及び、
前記コンピュータが、前記境界面に作用する前記流体モデルのせん断応力を、下記式（１）を用いて計算するステップを含むことを特徴とする流体シミュレーション方法。

ここで、
τ_correct：境界面に作用する流体モデルのせん断応力
Ｕ_live：流体領域の要素の重心点での速度
Ｕ_boundary：交点での速度
dist _{live cell to boundary}：流体領域の要素の重心点と交点との間の距離
μ：流体のせん断粘度
Ａ：境界面の面積
上記式（１）の dist _{live cell to boundary}には、下記式（２）で修正されるdist’ _{live cell to boundary}が用いられる請求項１に記載の流体シミュレーション方法。

ここで、
Ｌ１：直線の長さ
α：直線の長さＬ１に対するdist _{live cell to boundary}の比
前記境界面には、圧力境界条件としてノイマン境界条件が設定される請求項１又は２に記載の流体シミュレーション方法。