JP6703743B2 - 粒子法を用いた解析装置及びプログラム - Google Patents

Description

本発明は、粒子法を用いた解析装置及びプログラムに関する。

特許文献１は、粒子法等を用いた数値解析で使用する解析モデルの境界に位置する粒子を探索する境界粒子探索システムであって、解析モデルを構成する全粒子の位置座標より境界に位置する粒子をひとつ、または複数個抽出して探索基準粒子とする探索基準粒子設定部と、前記探索基準粒子設定部が抽出した探索基準粒子を始点とした探索する方向を示す探索方向ベクトルを設定するベクトル設定部と、探索基準粒子の位置座標を元に定義される領域内に存在する他の粒子を抽出する領域内粒子抽出部と、探索基準粒子に対する前記領域内粒子抽出部が抽出した領域内の各粒子の各位置ベクトルを作成し、その位置ベクトルの大きさと、位置ベクトルと前記ベクトル設定部が設定した探索方向ベクトルの角度を計算するベクトル計算部と、前記ベクトル計算部が計算したベクトルの大きさと角度により境界粒子を決定する境界粒子判断部と、を備える構成を開示し、特定の表面粒子を起点に近接する表面粒子を境界に沿ったベクトルに基づいて順次探索する方法が提案されている。

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

本発明は、表面方向を検出しないで判定する場合と比較して、表面粒子を高精度に判定することができる粒子法を用いた解析装置及びプログラムを提供することを目的とする。

請求項１に係る本発明は、対象粒子から解析対象の表面方向を検出する表面方向検出手段と、前記対象粒子の前記表面方向に存在する所定領域内の他粒子の数が閾値以上である場合には前記対象粒子を内部粒子であると判定し、閾値未満である場合には前記対象粒子を表面粒子であると判定する判定手段と、を有する粒子法を用いた解析装置である。

請求項２に係る本発明は、前記表面方向検出手段は、前記対象粒子の一定範囲内に存在する他粒子の座標平均値から重心を算出する重心算出手段と、前記対象粒子と前記重心を結んだベクトルを設定するベクトル設定手段と、を有する請求項１記載の粒子法を用いた解析装置である。

請求項３に係る本発明は、前記判定手段は、前記対象粒子を頂点とし、前記ベクトルに平行な直線を軸として任意の中心角を有する円錐内の他粒子の数が閾値未満である場合に前記対象粒子を表面粒子であると判定する請求項２に記載の粒子法を用いた解析装置である。

請求項４に係る本発明は、前記判定手段は、前記対象粒子を頂点とし、前記ベクトルに平行な直線を軸として対称な任意の角度を有する三角形内の他粒子の数が閾値未満である場合に前記対象粒子を表面粒子であると判定する請求項２に記載の粒子法を用いた解析装置である。

請求項５に係る本発明は、対象粒子から解析対象の表面方向を検出するステップと、前記対象粒子の前記表面方向に存在する所定領域内の他粒子の数が閾値以上である場合には前記対象粒子を内部粒子であると判定し、閾値未満である場合には前記対象粒子を表面粒子であると判定するステップと、をコンピュータに実行させるためのプログラムである。

請求項１に係る本発明によれば、表面方向を検出しないで判定する場合と比較して、表面粒子を高精度に判定することができる粒子法を用いた解析装置を提供することができる。

請求項２に係る発明によれば、複雑なアルゴリズムなしで表面方向を検出することができる。

請求項３に係る本発明によれば、三次元の場合であっても、ミクロなバラつきに強く、誤判定を抑制することができる。

請求項４に係る本発明によれば、二次元の場合であっても、ミクロなバラつきに強く、誤判定を抑制することができる。

請求項５に係る本発明によれば、表面方向を検出しないで判定する場合と比較して、表面粒子を高精度に判定することができるプログラムを提供することができる。

本発明の一実施形態における解析装置の構成を示すブロック図である。 本発明の一実施形態における解析装置の動作を示すフローチャートである。 ステップ１０２の表面粒子判定処理の詳細を示すフローチャートである。 本発明の一実施形態における表面粒子の判定方法を説明するための模式図である。 本発明の一実施形態における表面粒子の判定方法を説明するための模式図である。 本発明の実施形態に係る解析装置による表面粒子の判定結果を示す図である。 本発明の実施形態に係る解析装置を用いた数値解析例を示す図である。 比較例に係る解析装置による表面粒子の判定結果を示す図である。

粒子法は、粒子を用いた連続体の数値解析手法であって、連続体を有限個の粒子によって表し、連続体の挙動を粒子の運動によって計算する方法である。

代表的な粒子法として、非圧縮性流れの数値解析手法であるＭＰＳ（Ｍｏｖｉｎｇ Ｐａｒｔｉｃｌｅ Ｓｅｍｉ−ｉｍｐｌｉｃｉｔ）法が用いられている。ＭＰＳ法では、非圧縮条件を粒子数密度一定の条件とし、ここから圧力のポアソン方程式を導出している。自由液面は、粒子数密度の低下によって判定するので液面形状を描く必要はなく、流体の分裂や合体も容易に計算でき、界面の大変形を扱うことが容易、煩雑な格子生成作業が必要ない等の利点がある。ＭＰＳ法を用いた数値解析では、計算粒子が流体の表面に存在するか内部に存在するかを判定し、圧力勾配や表面張力といった相互作用力の計算対象を決めている。このため、表面粒子の判定を高精度に行うことが重要である。

従来から該当粒子の周辺に他粒子がどのくらい存在するかを表す数密度を用いて、数密度が閾値以下であれば連続体の表面に存在する表面粒子と判定する方法が用いられている。
しかしながら、この方法では、表面に存在しない粒子であっても表面粒子として判定される場合があり、粒子座標や液滴形状により判定精度がばらついて、３〜３０パーセント程度の誤判定が発生していた。図８は、液滴の濡れ広がりについて、数密度を用いて表面粒子か否か判定した結果を示しており、表面粒子であると判定された粒子と、表面に存在しない内部粒子であると判定された粒子が示されている。このように、上述した判定方法では、内部粒子であっても表面粒子であると誤判定される場合がある。

本発明は、上述した実情に鑑みてなされたものであり、本発明の一実施形態について図面を参照して詳細に説明する。
図１は、本発明の一実施形態における解析装置１０の全体構成を示すブロック図である。

解析装置１０は、図１に示されるように、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）１２、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）１４、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）１６、記憶部１８、操作部２０及び表示部２２を備えている。これらは、バス２４に接続されている。

ＣＰＵ１２は、ＲＯＭ１４または記憶部１８に格納されたプログラムに基づいて所定の処理を実行して、解析装置１０の各部の動作を制御する。

なお、本実施形態では、ＣＰＵ１２は、ＲＯＭ１４または記憶部１８内に格納されたプログラムを読み出して実行するものとして説明するが、プログラムをＣＤ−ＲＯＭ等の記憶媒体に格納してＣＰＵ１２に提供し、又はプログラムを不図示の通信装置を介して提供することも可能である。

ＲＯＭ１４は、半導体素子等で構成された読み出し専用の不揮発性記憶装置である。

ＲＡＭ１６は、半導体素子等で構成された揮発性の記憶手段であり、ＣＰＵ１２がプログラムを実行する際のワークエリアとして利用される。

記憶部１８は、ハードディスクドライブなどの大容量の記憶手段であり、ＣＰＵ１２に読み込まれるプログラムを記憶する。

操作部２０は、タッチセンサやキーボードやマウスなどを備えており、ユーザによる操作を受け付けてその操作内容に応じた信号をＣＰＵ１２に供給する。

表示部２２は、ＣＰＵ１２の制御の下、種々の情報を表示するモニタ等の出力手段である。

本実施形態では、解析装置１０によって、粒子法（ＭＰＳ法）による流体の解析を実行する。ＭＰＳ法では、粒子の密度、動粘性係数等の物性に関するパラメータ、重力加速度等の外力に関するパラメータ、初期配置における粒子速度、粒子座標等のパラメータ、後述する表面粒子判定に用いられる円錐Ｃの形状を決めるパラメータ等のパラメータが初期値として解析装置１０に入力される。
そして、解析装置１０は、記憶部１８に格納されたプログラムをＣＰＵ１２が実行することによって入力されたパラメータにより所定の計算が行われ、粒子の流速及び座標を所定の時間毎に更新していく。

次に、本発明の一実施形態における解析装置１０の動作について詳述する。
図２は、本発明の一実施形態における解析装置１０の動作を示すフローチャートである。

ステップ１００（Ｓ１００）において、操作部２０における操作者の入力により、計算条件が入力される。具体的には、粒子の密度、動粘性係数等の物性に関するパラメータ、重力加速度等の外力に関するパラメータ、境界条件等のパラメータを入力して解析を実行する。

ステップ１０１（Ｓ１０１）において、操作部２０における操作者の入力により、判定対象である対象粒子を所定の形状に初期配置し、粒子速度、粒子座標等を定義してパラメータを入力する。

ステップ１０２（Ｓ１０２）において、対象粒子が表面粒子か否かを判定する表面粒子判定処理を実行する。この表面粒子判定処理について図３乃至図５を用いて詳述する。図３は、ステップ１０２の表面粒子判定処理の詳細を示すフローチャートである。図４及び図５は、表面粒子の判定方法を説明するための模式図である。

ステップ２００（Ｓ２００）において、判定対象である対象粒子ｉを決定する。

ステップ２０１（Ｓ２０１）において、図４（ａ）及び図４（ｂ）に示すように、対象粒子ｉの有効半径ｒ内に存在する他粒子の座標平均値から重心ｇを算出する。ここで、有効半径ｒは任意であって、一定範囲内であればよく、一定範囲内の他粒子の座標平均値から重心ｇを算出する。

ステップ２０２（Ｓ２０２）において、図４（ｃ）に示すように、対象粒子ｉから重心ｇに向かうベクトルνを定義する。このベクトルνの方向が対象粒子ｉに対する表面方向を示しており、ベクトルνは、表面方向の法線ベクトルとなっている。

そして、図４（ｄ）に示すように、このベクトルνに平行な線を回転軸ｌとして、頂点を対象粒子ｉとして任意の中心角θで形成される円錐Ｃを形成する。そして、円錐Ｃ内に存在する他粒子の数を算出する。具体的には、後述するステップ２０３からステップ２０７に基づいて算出する。

ステップ２０３（Ｓ２０３）において、対象粒子ｉの近傍に存在する相手粒子ｊを決定する。

ステップ２０４（Ｓ２０４）において、対象粒子ｉと相手粒子ｊとの距離ｄ１を求める。距離ｄ１は、対象粒子ｉの中心の座標（ｘ_i，ｙ_i，ｚ_i）と相手粒子ｊの中心の座標（ｘ_j，ｙ_j，ｚ_j）とを用いて次のように表される。

次に、相手粒子ｊとベクトルνとの距離ｄ２を求める。すなわち、対象粒子ｉの中心の座標（ｘ_i，ｙ_i，ｚ_i）を通り、ベクトルν（ν_x，ν_y，ν_z）で方向が表されている軸ｌと、相手粒子ｊの中心の座標（ｘ_j，ｙ_j，ｚ_j）との距離ｄ２を求める。
ここで、ベクトルν（ν_x，ν_y，ν_z）は方向ベクトルであり、対象粒子ｉから重心ｇに向かうベクトルであるから、ベクトルνは、対象粒子ｉの中心の座標（ｘ_i，ｙ_i，ｚ_i）と重心ｇの座標（ｘ_g，ｙ_g，ｚ_g）、を用いて次のように表される。

また、距離ｄ２は、相手粒子ｊからベクトルνへ下した垂線の長さであるから、相手粒子ｊの中心の座標（ｘ_j，ｙ_j，ｚ_j）からベクトルνへ下した垂線のベクトルν上の座標を点Ｔ（ｘ_t，ｙ_t，ｚ_t）とすると、対象粒子ｉの中心の座標（ｘ_i，ｙ_i，ｚ_i）を通り、方向ベクトルが（ν_x，ν_y，ν_z）であるような直線の方程式は媒介変数ｔを用いて次のように表される。

したがって、相手粒子ｊの中心の座標（ｘ_j，ｙ_j，ｚ_j）と点Ｔ（ｘ_t，ｙ_t，ｚ_t）から相手粒子ｊから点Ｔへ向かうベクトルｊＴが次のように算出される。

また、このベクトルｊＴとベクトルνは直交するので、これらのベクトルの内積は０となるから、

が算出される。
また、相手粒子ｊとベクトルνとの距離ｄ２は、相手粒子ｊの中心の座標（ｘ_j，ｙ_j，ｚ_j）と点Ｔ（ｘ_t，ｙ_t，ｚ_t）から次のように表される。

以上から対象粒子ｉと相手粒子ｊと重心ｇの中心の座標をパラメータとして入力することにより距離ｄ１と距離ｄ２は算出される。

ステップ２０５（Ｓ２０５）において、相手粒子ｊがｄ１・ｓｉｎθ＞ｄ２を満たすか否か判定する。満たす場合には、相手粒子ｊが円錐Ｃ内に存在するとカウントされて、次のステップ２０６へ進む。満たさない場合には、相手粒子ｊは円錐Ｃ内に存在しないのでカウントされずに次のステップ２０６へ進む。

ステップ２０６（Ｓ２０６）において、全ての相手粒子ｊについて実行されたか否か判定する。実行された場合には、次のステップ２０７へ進み、実行されていない場合には、ステップ２０３へ戻る。

ステップ２０７（Ｓ２０７）において、円錐Ｃ内に存在するとカウントされた相手粒子ｊの個数Ｎが閾値Ｎｏ以上であるか否か判定する。個数Ｎが閾値Ｎｏ以上の場合には、対象粒子ｉは内部粒子であると判定され（ステップ２０８（Ｓ２０８））、閾値Ｎｏ未満の場合には、対象粒子ｉは表面粒子であると判定される（ステップ２０９（Ｓ２０９））。

ステップ１０３（Ｓ１０３）において、判定結果に応じて圧力勾配や表面張力といった相互作用力を計算する。具体的には、対象粒子ｉが内部粒子であると判定された場合には圧力として所定のパラメータを入力し、対象粒子ｉが表面粒子であると判定された場合には圧力を０として相互作用力を計算する。

ステップ１０４（Ｓ１０４）において、圧力勾配による力をもとに粒子座標を更新する。

ステップ１０５（Ｓ１０５）において、指定した計算時間を経過したか等の終了条件を満たすか否かが判定され、満たす場合には計算を終了し、満たさない場合にはステップ１０２へ戻る。

図６は、上述した本発明の実施形態により液滴の表面粒子判定を行った結果を示しており、表面粒子であると判定された粒子と、表面に存在しない内部粒子であると判定された粒子が示されている。このように、上述した判定方法では、液滴形状に対してロバストな判定が可能であって、誤判定率が０〜２％程度に大幅に低減された。

図７は、本発明の実施形態に係る解析装置１０を用いて粒子法シミュレーション解析を行った例を示している。図７は、液滴の濡れ広がりについて、液滴形状の接触角θを（ａ）は７５度、（ｂ）は４３度、（ｃ）は１７度、（ｄ）は７度に設定して表示されたシミュレーション解析結果である。
親水条件では、三相界面での表面張力の釣り合いをより厳密に扱う必要があるため、高精度に表面粒子を判定しなければならず、従来では、接触角３０度程度までしかシミュレーション解析で再現できなかったが、本発明により接触角７度の親水条件までの形状を再現することができた。

なお、表面粒子と判定された粒子のうち、壁面から一定距離内にある粒子を三相界面粒子として判定し、濡れ挙動の計算時の界面張力の算出に活用することができる。その際、三相界面近傍の流体高さが低い場合に三相界面粒子の誤判定が発生しやすくなる。その場合は、さらに壁面に平行な方向で表面粒子の判定アルゴリズムを適用して、三相界面粒子か否かを判定する。

なお、本実施形態においては、表面粒子判定処理において、円錐Ｃ内に相手粒子ｊの中心点が入る個数を閾値にして表面粒子か否かを判定したが、これに限らず、円錐Ｃ内に相手粒子の全部又は一部が入る個数を閾値にしてもよく、また、円錐Ｃ内に入る相手粒子の体積を閾値にしてもよい。

また、本実施形態においては、円錐Ｃの中心角と高さを任意としたが、計算目的によって使い分けることが可能である。すなわち、中心角を小さくすればするほど厳密な判定が可能となる。また、高さを高くするほど判定精度は向上される。

また、本実施形態においては、三次元の流体について、円錐Ｃ内に存在するとカウントされた相手粒子ｊの個数を閾値にして表面粒子か否かを判定したが、二次元の場合には、対象粒子ｉを頂点として、対象粒子ｉから重心ｇに向かうベクトルνに平行な直線を軸ｌとして対称な任意の角度を有する三角形を用いて、この三角形内に存在するとカウントされた相手粒子の個数を閾値にして表面粒子か否かを判定する。

以上、本発明は、粒子法シミュレーション等に用いられるが、上記の実施形態に限定されるものではなく、種々の変更が可能であり、数密度法と組み合わせて用いることもできる。

なお、本実施形態においては、流体を粒子によって表し、流体の挙動を粒子の運動によって計算する方法について説明したが、これに限らず、粉体や固体等の連続体に適用される。

１０ 解析装置
１２ ＣＰＵ
１４ ＲＯＭ
１６ ＲＡＭ
１８ 記憶部
２０ 操作部
２２ 表示部

Claims (5)

1. 対象粒子から解析対象の表面方向を検出する表面方向検出手段と、
   前記対象粒子の前記表面方向に存在する所定領域内の他粒子の数が閾値以上である場合には前記対象粒子を内部粒子であると判定し、閾値未満である場合には前記対象粒子を表面粒子であると判定する判定手段と、
   を有する粒子法を用いた解析装置。
2. 前記表面方向検出手段は、
   前記対象粒子の一定範囲内に存在する他粒子の座標平均値から重心を算出する重心算出手段と、
   前記対象粒子と前記重心を結んだベクトルを設定するベクトル設定手段と、
   を有する請求項１記載の粒子法を用いた解析装置。
3. 前記判定手段は、前記対象粒子を頂点とし、前記ベクトルに平行な直線を軸として任意の中心角を有する円錐内の他粒子の数が閾値未満である場合に前記対象粒子を表面粒子であると判定する請求項２に記載の粒子法を用いた解析装置。
4. 前記判定手段は、前記対象粒子を頂点とし、前記ベクトルに平行な直線を軸として対称な任意の角度を有する三角形内の他粒子の数が閾値未満である場合に前記対象粒子を表面粒子であると判定する請求項２に記載の粒子法を用いた解析装置。
5. 対象粒子から解析対象の表面方向を検出するステップと、
   前記対象粒子の前記表面方向に存在する所定領域内の他粒子の数が閾値以上である場合には前記対象粒子を内部粒子であると判定し、閾値未満である場合には前記対象粒子を表面粒子であると判定するステップと、
   をコンピュータに実行させるためのプログラム。

Images