JP2018147315A - 流体シミュレーション方法及び流体シミュレーションのプログラム - Google Patents

Abstract

【課題】解析精度を確保しつつも計算負荷を抑えることが可能な流体シミュレーション方法を提供する。【解決手段】粒子法における固体物と流体とが隣接する状況下での流体シミュレーション方法である。そして、固体物の境界をポリゴンと固体粒子でモデル化するステップと、流体粒子の影響半径内にあるポリゴンの数を数えるステップＳ１２と、ポリゴンの数が１枚の場合は固体物モデルの積分をポリゴンによる積分領域で計算し（ステップＳ１５）、ポリゴンの数が複数の場合は固体物モデルの積分を固体粒子により計算する（ステップＳ１６）ことを特徴とする。【選択図】図１

Description

本発明は、粒子法における固体物と流体とが隣接する状況下での流体シミュレーション方法及び流体シミュレーションのプログラムに関するものである。

特許文献１に開示されているように、水や雪や泥などのタイヤによって飛び跳ねる介在物を、複数の粒子モデルに置き換えて、粒子法によって解析することで飛び跳ね特性を評価する方法が知られている。

粒子法による数値解析では、モデル化する粒子の数が多くなると、メモリの使用量が増大するため計算機に大きな負荷がかかる。このため、特許文献２に開示されているように、解析領域を限定することで、メモリ使用量を低減するプログラムが開発されている。

また、特許文献３に開示されているように、流体に接する壁境界をポリゴン距離関数によって作成することで、効率的な解析が可能となるようにする方法も開発されている。

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

しかしながらポリゴンを使用して境界表現を行うと、計算負荷は小さくなるが、形状が複雑な箇所では計算のアルゴリズムが複雑になり、簡略化を行えば解析精度が低下するという問題がある。一方、固体壁などもすべて粒子でモデル化すれば、計算負荷が増大して効率が低下する。

そこで、本発明は、解析精度を確保しつつも計算負荷を抑えることが可能な流体シミュレーション方法及び流体シミュレーションのプログラムを提供することを目的としている。

前記目的を達成するために、本発明の流体シミュレーション方法は、粒子法における固体物と流体とが隣接する状況下での流体シミュレーション方法であって、固体物の境界をポリゴンと固体粒子でモデル化するステップと、流体粒子の影響半径内にある前記ポリゴンの数を数えるステップと、前記ポリゴンの数が１枚の場合は固体物モデルの積分を前記ポリゴンによる積分領域で計算し、前記ポリゴンの数が複数の場合は固体物モデルの積分を前記固体粒子により計算するステップとを備えたことを特徴とする。

また、流体シミュレーションのプログラムの発明は、粒子法における固体物と流体とが隣接する状況下で流体シミュレーションを実行させる演算過程においてコンピュータを、固体物の境界をポリゴンと固体粒子でモデル化した解析データを読み込ませる手段と、流体粒子の影響半径内にある前記ポリゴンの数をカウントさせる手段と、前記ポリゴンの数が１枚の場合は固体物モデルの積分を前記ポリゴンによる積分領域で計算させ、前記ポリゴンの数が複数の場合は固体物モデルの積分を前記固体粒子により計算させる手段として機能させる。

このように構成された本発明の流体シミュレーション方法及び流体シミュレーションのプログラムは、固体物の境界をポリゴンでモデル化したうえで、複数のポリゴンが交差するような箇所では、固体粒子によるモデル化を行う。

このため、ポリゴンが交差するような複雑な形状となる箇所でも簡略化が行われず解析精度を確保することができる。さらに、単純な形状の箇所ではポリゴンによって関数表現されるので、計算負荷を抑えることができる。

本実施の形態の流体シミュレーション方法の境界処理の各ステップを説明するフロー図である。 本実施の形態の流体シミュレーション方法の全体の処理の流れを説明するフロー図である。 固体壁を固体壁粒子で表現したモデルを模式的に示した説明図である。 固体壁を固体壁ポリゴンで表現したモデルを模式的に示した説明図である。 流体粒子の影響半径内に存在する固体壁粒子を積分点とすることを模式的に示した説明図である。 流体粒子の影響半径内に存在する固体壁ポリゴンを積分領域とすることを模式的に示した説明図である。 流体粒子の影響半径内に交差する固体壁ポリゴンによって凸形状の積分領域が形成されたケースを示した説明図である。 凸形状の積分領域に対して簡略的に設定される積分領域を示した説明図である。 流体粒子の影響半径内に交差する固体壁ポリゴンによって凹形状の積分領域が形成されたケースを示した説明図である。 凹形状の積分領域に対して簡略的に設定される積分領域を示した説明図である。 流体粒子の影響半径内に固体壁ポリゴンが１枚存在する場合の積分領域を示した説明図である。 流体粒子の影響半径内に固体壁ポリゴンが２枚存在する場合の積分点を示した説明図である。 本実施の形態の流体シミュレーション方法の境界条件の設定を説明する模式図である。 本実施の形態の流体シミュレーション方法の境界処理を他の処理方法と比較して一覧表にした図である。 車輪とレールとの間に介在される水膜の挙動を解析する例の説明図である。 車輪とレールと水膜を粒子で表現した場合の解析モデルを比較例として例示した説明図である。 固体粒子とポリゴンとを併用した併用モデルを例示した説明図である。 車輪の併用モデルを拡大して示した説明図である。 レールの併用モデルを拡大して示した説明図である。

以下、本発明の実施の形態について図面を参照して説明する。図１は、本実施の形態の流体シミュレーション方法の要部となる境界処理の流れを各ステップにして説明するフロー図である。一方、図２は、流体シミュレーション方法の全体の処理の流れを各ステップにして説明するフロー図である。

本実施の形態の流体シミュレーションは、水や雪や泥などを粒子としてモデル化する粒子法を用いて実行される。粒子法は、計算対象物を複数の粒子モデルの集まりとし、分離及び移動可能にして行われる数値解析の離散化手法の一つである。粒子法では、移動規定条件の下に粒子モデルを移動させ、粒子モデルの持つ速度、密度又は全エネルギー量などによって、粒子モデルの状態を評価する。

粒子法には、MPS（Moving Particle Simulation）法、SPH（Smoothed Particle Hydrodynamics）法などがあり、図２は非圧縮性流れを解くための粒子法の概略的な処理の流れを示している。

粒子法では、図３Ａに示すように、水や雪や泥などの流体として扱えるものを流体粒子２，・・・としてモデル化する。ここで流体粒子２は、２次元解析では円形、３次元解析では球体の一様な大きさに設定される。

一方、流体粒子２，・・・が接する壁などの固体物は、境界条件として設定する必要がある。その場合に２種類の境界表現の方法があり、一つ目の表現方法としては、固体物モデルである固体壁１Ａを、固体粒子である固体壁粒子１１，・・・で表現する方法が挙げられる。

２つ目の表現方法としては、図３Ｂに示すように、固体物モデルである固体壁１Ｂを、ポリゴンである固体壁ポリゴン１２，・・・で表現する方法である。ポリゴンは、流体から外力を受けても移動しない壁などの固体物の境界を、ポリゴンによって表したものである。ここでは、固体壁ポリゴン１２を直線（３次元解析では平面）で示している。

流体粒子２は、図４Ａに示すように、流体粒子２の影響半径２１内に入る粒子の物理量を積分する必要がある。このため、固体壁１Ａ付近に存在する流体粒子２は、影響半径２１内の固体壁１Ａ内部の物理量を積分することになる。

一方、固体物内部の物理量は、流体部分の物理量から決まるため、必ずしも固体粒子を配置して物理量を与える必要はない。このため、図４Ｂに示すように、流体粒子２の物理量から固体壁１Ｂ内部の積分をすることが可能となり、境界条件を固体壁ポリゴン１２で表現することができる。

そこで、図４Ａに示すように、固体壁１Ａを固体壁粒子１１，・・・でモデル化した場合は、流体粒子２の影響半径２１内に入るものは積分点３，・・・として、積分計算に組み込む。

一方、図４Ｂに示すように、固体壁１Ｂを固体壁ポリゴン１２でモデル化した場合は、流体粒子２の影響半径２１内に入る範囲を積分領域４として、積分計算に組み込む。この図のケースでは、影響半径２１内を直線状（３次元解析では平面状）の固体壁ポリゴン１２で仕切るだけなので、弓形（２次元解析の場合）の積分領域４の積分は、単純なアルゴリズムで精度よく計算することができる。

これに対して、図５Ａに示すように、２枚の固体壁ポリゴン１２Ａ，１２Ａが交差して凸形状の積分領域４Ａが形成される場合がある。この場合は、積分領域４Ａの計算のアルゴリズムが複雑になり、負荷のかかる計算が必要になる。

そこで凸形状の積分領域４Ａを簡略化するために、図５Ｂに示すように、最近接の１枚の固体壁ポリゴン１２０Ａ面のみを考慮し、それで仕切られた範囲を積分領域４０Ａとして、計算を簡単にすることができる。しかしながらこの方法では、積分領域４Ａが過大に評価された積分領域４０Ａとなって、簡略化した分だけ精度は低下する。

同様に、図６Ａに示すように、２枚の固体壁ポリゴン１２Ｂ，１２Ｂが交差して凹形状の積分領域４Ｂが形成される場合がある。この場合も、積分領域４Ｂの計算のアルゴリズムが複雑になり、負荷のかかる計算が必要になる。

そこで凹形状の積分領域４Ｂを簡略化するために、図６Ｂに示すように、最近接の１枚の固体壁ポリゴン１２０Ｂ面のみを考慮し、それで仕切られた範囲を積分領域４０Ｂとして、計算を簡単にすることができる。しかしながらこの場合は、積分領域４Ｂが過小に評価された積分領域４０Ｂとなって、簡略化した分だけ精度は低下する。

そこで、本実施の形態の流体シミュレーション方法の境界処理では、固体壁粒子１１と固体壁ポリゴン１２を併用した固体物モデルの境界表現を行う。すなわち、図７に示すように、固体壁ポリゴン１２が１枚だけ流体粒子２の影響半径２１内に存在する場合は、固体壁ポリゴン１２の境界表現をそのままにして、直線状の固体壁ポリゴン１２によって仕切られた円弧と弦に囲まれた領域（３次元解析の場合は球面と割平面とに囲まれた領域）を積分領域４とする。

一方、図８に示すように、２枚の固体壁ポリゴン１２，１２が流体粒子２の影響半径２１内に存在する場合は、固体壁ポリゴン１２，１２と影響半径２１の円弧（３次元解析の場合は球面）で仕切られた領域を固体壁粒子１１，・・・で置き換え、固体壁粒子１１を積分点３とする。

要するに、図９に示すように、形状が複雑となる固体壁１の隅角部には予め固体壁粒子１１，・・・を配置しておき、直線状（３次元解析の場合は平面状）の隅角部間には固体壁ポリゴン１２を配置する。こうすることによって、固体壁粒子１１を配置した箇所では精度の高い計算が行われ、固体壁ポリゴン１２のみを配置した箇所では効率的な計算が行われる。

図１０は、ここまで説明してきた計算方法のメリットとデメリットを対比した一覧表である。この一覧表に示したように、固体壁１Ａを固体壁粒子１１のみで表現して計算した場合（図３Ａ参照）、計算コスト（計算負荷）が大きくなる。

また、固体壁１Ｂを固体壁ポリゴン１２のみで表現して正確に計算した場合（図３Ｂ参照）、計算のアルゴリズムが複雑になる。また、簡略に計算すると（図５Ｂ，図６Ｂ参照）、計算コストは小さくできるが、精度が低くなる。
そこで、計算コストが小さく、精度が良いうえに、計算のアルゴリズムを単純にできる固体壁粒子１１と固体壁ポリゴン１２とを併用した固体壁１で計算することとする。

次に、図１，２を参照しながら、本実施の形態の流体シミュレーション方法の各ステップを説明する。
まず、図２のフローチャートに示すように、粒子法を用いた流体シミュレーションを行うにあたって、ステップＳ０で解析データを読み込ませる。この解析データは、流体を流体粒子２でモデル化し、固体壁１を固体壁粒子１１と固体壁ポリゴン１２とによってモデル化させたものである。

例えば、解析データを単純化して説明すると、図９に示すように、容器の中に充填された流体粒子２，・・・と、固体壁１の側面及び底面に配置される固体壁ポリゴン１２，１２，１２と、固体壁１の隅角部周辺に配置される固体壁粒子１１，・・・とがモデルデータとなる。
その他にも、移動規定条件となる流体の物性（質量、粘性など）、初期条件としての速度や配置位置などの解析に必要となるデータが読み込まれる。

続いて流体シミュレーションの演算を開始すると、重力と粘性による力の計算が行われる（ステップＳ１）。粘性による力の計算を行う際に、速度の空間２階微分を計算する必要があるため、流体粒子２の影響半径２１内の物理量の積分を行う必要がある。そして、このステップＳ１の計算の中で、流体粒子２の影響半径２１内の物理量の積分を行うに際して、境界処理（ステップＳ１０）が行われる。

図１に境界処理のフローチャートを示した。流体粒子２の影響半径２１内の物理量の積分計算（ステップＳ１１）は、影響半径２１内に入る固体壁１内部の物理量を重み付け加算することで行われる。

まずステップＳ１２で、流体粒子２の影響半径２１内にある固体壁ポリゴン１２の数Ｎを数える。例えば、図７に示した流体粒子２の影響半径２１内にある固体壁ポリゴン１２の数Ｎは、１枚とカウントされる。一方、図８に示した流体粒子２の影響半径２１内にある固体壁ポリゴン１２，１２の数Ｎは、２枚（複数）とカウントされる。

そこで、ステップＳ１３ではカウントされた数Ｎを判定し、流体粒子２の影響半径２１内に固体壁ポリゴン１２が存在しなかった場合（Ｎ＝０）は、固体壁１の影響がその流体粒子２には及ばないとして固体壁処理を停止する（ステップＳ１４）。

一方、図７に示すように流体粒子２の影響半径２１内に固体壁ポリゴン１２が１枚だけ存在した場合は、積分計算を固体壁ポリゴン１２によって仕切られた積分領域４で行う（ステップＳ１５）。

さらに、図８に示すように流体粒子２の影響半径２１内で固体壁ポリゴン１２，１２が交差して２枚存在した場合は、２枚の固体壁ポリゴン１２，１２と影響半径２１の球面で囲まれた領域内の固体壁粒子１１を積分点３として、積分計算を行う（ステップＳ１６）。

そして、すべての流体粒子２，・・・の積分計算を終えた後に、ステップＳ１に戻る（図１参照）。このようにして行われた重力と粘性による力の計算結果に基づいて、ステップＳ２では、初期条件で与えられた流体粒子２，・・・の速度と位置の更新を行う。

続いてステップＳ３では、圧力計算を行う。圧力計算は、速度の空間１階微分や圧力の空間２階微分や密度を計算する必要があるため、流体粒子２の影響半径２１内の物理量の積分を行う必要がある。この圧力計算の際にも、境界処理（ステップＳ１０）が行われる。

またステップＳ４では、圧力計算の結果に基づいて圧力勾配による力の計算が行われる。圧力勾配による力の計算は、圧力の空間１階微分を計算する必要があるため、流体粒子２の影響半径２１内の物理量の積分を行う必要がある。この圧力勾配による力の計算の際にも、境界処理（ステップＳ１０）が行われる。

そして、ステップＳ５では、ここまでの計算結果に基づいて、流体粒子２，・・・の速度と位置の修正を行う。さらにステップＳ６で、演算を続けるか否かの終了判定をし、続ける場合はステップＳ１からの演算処理を繰り返す。

次に、具体的に解析モデルの一例を示して説明する。
図１１に例示するように、走行する車両の車輪５とレール６との間に、水膜７が介在した場合の影響を流体シミュレーションによって解明する水膜７の挙動解析について説明する。ここで、車輪５とレール６は固体物に該当し、水膜７は流体に該当する。

解析は、車両が20m/s(72km/hr)で等速に移動している状態を想定し、車輪５とレール６は新幹線用に準拠してモデル化した。また、車輪５とレール６の接触部分の要素サイズは、レール６側が3mm、車輪５側が6.23mmとした。そして、車輪５とレール６の接触部分は、有限要素法の弾塑性解析の結果により動的に変形するものとした。

まず、比較例として、図１２Ａに示すように、車輪５とレール６と水膜７のすべてを粒子で表現した解析モデルを考える。この解析モデルでは、流体粒子７１の直径を0.05mmとすると、流体粒子７１は2560万粒子となる。また、車輪５の固体粒子５１は、直径を0.05mmとすると6890万粒子となり、レール６の固体粒子６１は、直径を0.05mmとすると2390万粒子となった。要するに、固体粒子５１，６１の合計は、9280万粒子となる。

これに対して、図１２Ｂに示すように、車輪５とレール６とを固体粒子５１，６１とポリゴン５２，６２との併用モデル５０，６０とすると、車輪５の固体粒子５１は553万粒子となり、レール６の固体粒子６１は377万粒子となった。固体粒子５１，６１の合計は、930万粒子となり、すべてを粒子で表現した場合と比較して固体粒子５１，６１の数を１／１０に低減することができた。

そこで、本実施の形態の併用モデル５０，６０の詳細についてさらに説明する。図１３Ａには、車輪５を固体粒子５１，・・・とポリゴン５２，・・・とを併用してモデル化した併用モデル５０を拡大して示した。この図に示すように、複数の部材面が交差して辺となる箇所や隅角部や頂点となる箇所にのみ固体粒子５１を配置し、それ以外の部材が一面で形成されている内部の箇所はポリゴン５２でモデル化した。

一方図１３Ｂには、レール６を固体粒子６１，・・・とポリゴン６２，・・・とを併用してモデル化した併用モデル６０を拡大して示した。この図に示すように、複数の部材面が交差して辺となる箇所や隅角部や頂点となる箇所にのみ固体粒子６１を配置し、それ以外の部材が一面で形成されている内部の箇所はポリゴン６２でモデル化した。

次に、本実施の形態の流体シミュレーション方法及び流体シミュレーションのプログラムの作用について説明する。
このように構成された本実施の形態の流体シミュレーション方法及び流体シミュレーションのプログラムは、固体壁１の境界を固体壁ポリゴン１２でモデル化したうえで、固体壁ポリゴン１２が交差するような隅角部などの複雑な形状となる箇所では、固体壁粒子１１によるモデル化を行う。

このため、固体壁ポリゴン１２が凸形状や凹形状のような複雑な形状となる箇所では固体壁粒子１１を積分点３とする計算が行われ、簡略化が行われないので解析精度を確保することができる。さらに、直線や平面などの単純な形状の箇所ではポリゴンによって関数表現されるので、計算負荷を抑えることができる。

以上、図面を参照して、本発明の実施の形態を詳述してきたが、具体的な構成は、この実施の形態に限らず、本発明の要旨を逸脱しない程度の設計的変更は、本発明に含まれる。

例えば前記実施の形態で説明した流体シミュレーション方法及び流体シミュレーションのプログラムは、２次元解析であっても３次元解析であっても適用することができる。

また、前記実施の形態では、車輪５とレール６との間に流体である水膜７が介在される場合を例に説明したが、これに限定されるものではなく、粒子法を用いるあらゆる流体シミュレーション方法及び流体シミュレーションのプログラムに本発明を適用することができる。

１ 固体壁（固体物モデル）
１１ 固体壁粒子
１２ 固体壁ポリゴン
２ 流体粒子
２１ 影響半径
３ 積分点
４ 積分領域
５ 車輪（固体物）
５０ 併用モデル（固体物モデル）
５１ 固体粒子
５２ ポリゴン
６ レール（固体物）
６０ 併用モデル（固体物モデル）
６１ 固体粒子
６２ ポリゴン
７ 水膜（流体）
７１ 流体粒子
Ｎ ポリゴンの数

Claims (3)

1. 粒子法における固体物と流体とが隣接する状況下での流体シミュレーション方法であって、
   固体物の境界をポリゴンと固体粒子でモデル化するステップと、
   流体粒子の影響半径内にある前記ポリゴンの数を数えるステップと、
   前記ポリゴンの数が１枚の場合は固体物モデルの積分を前記ポリゴンによる積分領域で計算し、前記ポリゴンの数が複数の場合は固体物モデルの積分を前記固体粒子により計算するステップとを備えたことを特徴とする流体シミュレーション方法。
2. 前記ポリゴンの数が１枚の場合は、前記ポリゴンによって仕切られた前記影響半径内を積分領域とし、前記ポリゴンの数が複数の場合は、前記影響半径内の前記固体粒子を積分点とすることを特徴とする請求項１に記載の流体シミュレーション方法。
3. 粒子法における固体物と流体とが隣接する状況下で流体シミュレーションを実行させる演算過程においてコンピュータを、
   固体物の境界をポリゴンと固体粒子でモデル化した解析データを読み込ませる手段と、
   流体粒子の影響半径内にある前記ポリゴンの数をカウントさせる手段と、
   前記ポリゴンの数が１枚の場合は固体物モデルの積分を前記ポリゴンによる積分領域で計算させ、前記ポリゴンの数が複数の場合は固体物モデルの積分を前記固体粒子により計算させる手段として機能させるための流体シミュレーションのプログラム。