JP6900118B2 - シミュレーション方法、シミュレーション装置、及びプログラム - Google Patents

Description

本発明は、剛体の動きと複数の粒子の挙動とを連成解析するシミュレーション方法、シミュレーション装置、及びプログラムに関する。

分子動力学法（ＭＤ法）等の粒子法を用いたシミュレーションでは、剛体を複数の粒子で表して、各粒子の挙動を解析する手法が一般的である（例えば、特許文献１）。例えば、剛体の形状が複雑になると、剛体を表す粒子数が増大し、計算負荷が大きくなる。

また、粒子法において、一般的に、直交座標系の相互に直交する三方向に関して周期性を持つと仮定される系に対して、周期境界条件が適用される。剛体の形状が、相互に直交する三方向に関して周期性を持たない場合、このような周期境界条件を適用することができない。

特許第６０９８１９０号公報

本発明の目的は、剛体の動きと、その剛体と相互作用する複数の粒子の挙動とのシミュレーションにおいて、計算負荷を低減させることが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

本発明の一観点によると、
中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域内の複数の粒子の挙動をシミュレーションする方法であって、
前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義し、
前記剛体の形状に基づいて、前記剛体の形状を表す符号付距離関数を生成し、
前記部分解析領域の前記半平面に対応する面において、前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用し、前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションにより求めるシミュレーション方法が提供される。

本発明の他の観点によると、
中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域の内部に存在する複数の粒子の挙動をシミュレーションする処理装置を有し、
前記処理装置は、
入力装置から、前記剛体の形状を定義する形状データが入力されると、入力された前記形状データに基づいて、前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義し、
前記形状データに基づいて、前記剛体の形状を表す符号付距離関数を生成し、
前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用し、前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションにより求め、
シミュレーション結果を出力装置に出力するシミュレーション装置が提供される。

本発明のさらに他の観点によると、
中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域内の複数の粒子の挙動をシミュレーションする機能を、コンピュータに実現させるプログラムであって、
前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義する機能、
前記部分解析領域の前記半平面に対応する面において、前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用する機能、
前記剛体の形状に基づいて、前記剛体の形状を表す符号付距離関数を生成する機能、及び
前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションする機能
をコンピュータに実現させるプログラムが提供される。

剛体の動きと、その剛体と相互作用する複数の粒子の挙動とのシミュレーションにおいて、計算負荷を低減させることが可能である。

図１は、実施例によるシミュレーション方法のシミュレーション対象の一例を示す斜視図である。 図２Ａ及び図２Ｂは、シミュレーション対象の中心軸に対して垂直な断面図である。 図３は、剛体の形状を定義する符号付距離関数（ＳＤＦ）を説明するための模式図である。 図４は、実施例によるシミュレーション装置のブロック図である。 図５は、実施例によるシミュレーション方法のフローチャートである。 図６Ａは、シミュレーションにより求められた角速度の時間変化を示すグラフであり、図６Ｂは、計算に要した時間を示すグラフである。

図１を参照して、実施例によるシミュレーション方法のシミュレーション対象について説明する。

図１は、実施例によるシミュレーション方法のシミュレーション対象の一例を示す斜視図である。中心軸１０に関して回転対称な剛体１１、及び剛体が収容される空間（解析領域）内の流体に相当する複数の粒子１２がシミュレーション対象となる。すなわち、実施例では、剛体と流体との連成解析を行う。シミュレーションにおいて、剛体１１の回転運動、及び複数の粒子１２の挙動を求める。図１に示したシミュレーション対象の剛体１１は、円柱状の回転軸１１Ａと、この回転軸１１Ａの側面から半径方向に突出する複数の回転翼１１Ｂで構成される。回転翼１１Ｂの枚数をｎで表すと、剛体１１は中心軸１０の周りにｎ回回転対称である。ここでｎは２以上の整数である。

実施例によるシミュレーションでは、解析領域を、中心軸１０を境界とする２つの半平面１３、１４で分割して、一方の領域を部分解析領域１５として定義する。部分解析領域１５の、中心軸１０に垂直な断面は扇形である。シミュレーションの演算は部分解析領域１５内で行う。剛体１１の形状及び解析領域が中心軸１０に関してｎ回回転対称性（ｎは２以上の整数）を持つとき、２つの半平面１３、１４のなす角度が３６０°／ｎの整数倍になるように半平面１３、１４を定義する。

半平面１３、１４に対応する部分解析領域１５の境界において、中心軸１０を回転中心とした回転方向に関して周期性を持つように周期境界条件を定義する。周期境界条件については、後に図２Ａ及び図２Ｂを参照して説明する。部分解析領域１５の扇形断面の円弧部分に対応する境界には、例えば反射境界条件を適用する。

剛体１１の形状は、剛体１１の形状を表す符号付距離関数（ＳＤＦ）で表す。剛体１１と粒子１２との距離が、このＳＤＦと粒子１２の座標とから算出される。ＳＤＦについては、後に図３を参照して説明する。

次に、図２Ａを参照して、周期境界条件について説明する。
図２Ａは、中心軸１０に対して垂直な断面図である。図２Ａでは、回転翼１１Ｂが８枚の例を示している。また、２枚の半平面１３と１４とのなす角度（中心角）が１８０度の例を示している。図２Ａに示した例では、図の上半分を部分解析領域１５として定義し、部分解析領域１５内の剛体１１及び粒子１２の挙動についてシミュレーション演算を行う。シミュレーション演算には、例えば分子動力学法を適用する。

部分解析領域１５の境界、すなわち半平面１３、１４の近傍に位置する粒子１２Ａに関しては、粒子１２Ａの位置に対して回転対称位置１２ａの近傍に存在する粒子１２及び剛体１１との相互作用を考慮して演算を行う。粒子１２Ｂが周期境界を越えて部分解析領域１５の外側に移動した場合には、移動後の粒子の位置に対して回転対称となる位置に周期境界を越えて粒子１２ｂを進入させる。

図２Ｂに示すように、２枚の半平面１３と１４とのなす角度を９０度にしてもよい。この場合には、部分解析領域１５の周期境界が９０度の角度で交差することになる。すなわち、部分解析領域１５の中心角が９０度になる。半平面１３、１４に対応する２つの境界に対して周期境界条件が適用される。また、図２Ａ及び図２Ｂに示したように、剛体１１が８回回転対称性を有する場合には、部分解析領域１５の中心角を４５度にしてもよい。

次に、図３を参照して剛体１１（図１）の形状から生成されたＳＤＦについて説明する。
図３は、ＳＤＦを説明するための模式図である。剛体１１の周囲に三次元直交格子を配置する。直交格子の各格子点（ノード）２０から剛体までの最短距離Ｌを、格子点２０ごとに算出する。格子点２０ごとに算出された最短距離Ｌが、ＳＤＦである。格子点２０が剛体１１の外側に位置するとき、最短距離Ｌの符号を正と定義し、格子点２０が剛体１１の内側に位置するとき、最短距離Ｌの符号を負と定義する。

粒子１２から剛体１１までの距離は、粒子１２が存在する直交格子のセルの頂点のＳＤＦ値（頂点から剛体１１までの最短距離Ｌ）に基づいて補間演算を行うことにより算出することができる。

次に、図４を参照して実施例によるシミュレーション装置について説明する。
図４は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、処理装置３０、記憶装置５０、入力装置６１、及び出力装置６２を含む。処理装置３０は、ハードウェア的にはコンピュータの中央処理ユニット（ＣＰＵ）等で実現され、各ブロックの機能はＣＰＵがコンピュータプログラムを実行することにより実現される。

処理装置３０は、記憶装置５０、入力装置６１、及び出力装置６２に接続されている。入力装置６１は、処理装置３０が実行する処理に関連する種々のデータをユーザから受け付ける。入力装置６１として、例えばキーボード、ポインティングデバイス、リムーバブルメディアの読み取り装置、通信装置等を用いることができる。出力装置６２として、例えば画像表示装置、リムーバブルメディアにデータを書き込む書き込み装置、通信装置等を用いることができる。処理装置３０は、シミュレーション結果等の情報を出力装置６２に出力する。

処理装置３０は、部分解析領域取得部３１、粒子取得部３２、剛体取得部３３、出力部３４、及び数値演算部４０を含む。記憶装置５０に、部分解析領域情報保持部５１、剛体情報保持部５２、及び粒子情報保持部５３が確保されている。

部分解析領域取得部３１は、入力装置６１を介して、解析領域を分割して部分解析領域１５（図１）を定義するための情報を取得する。部分解析領域取得部３１は、取得した部分解析領域１５を仮想空間内で表現する情報を部分解析領域情報保持部５１に登録する。部分解析領域１５を仮想空間内で表現する情報は、平面、曲面、直線、曲線を表す関数を含み、特に中心軸１０、半平面１３、１４（図１）を表す情報を含む。

粒子取得部３２は、ユーザから入力装置６１を介して、シミュレーション対象の粒子系を定義する情報を取得する。粒子系を定義する情報には、粒子１２（図１）の個数、各粒子の位置及び速度を定義するための情報が含まれる。粒子取得部３２は、粒子１２の各々に、その位置及び速度を対応付けて粒子情報保持部５３に登録する。

剛体取得部３３は、ユーザから入力装置６１を介して、剛体１１（図１）の形状を定義する形状データを取得する。さらに、この形状データに基づいて、部分解析領域１５（図１）内に位置する部分についてＳＤＦを生成し、剛体情報保持部５２に登録する。

出力部３４は、シミュレーション結果を出力装置６２に出力する。その他に、ユーザに通知すべき情報、例えばデータの入力を促す情報等を出力装置６２に出力する。

次に、数値演算部４０の各ブロックの機能について説明する。
粒子剛体距離演算部４１が、粒子１２（図１）の現在位置、及び剛体１１（図１）の形状を定義するＳＤＦに基づいて、粒子１２から剛体１１までの距離を算出する。

粒子座標変換部４２が、剛体１１の回転に対応して粒子１２の座標を変換する。具体的には、剛体１１及び粒子１２が配置されている仮想空間内で剛体１１を回転させる代わりに、粒子１２の座標変換を行う。これにより、剛体１１と粒子１２との相対位置関係は、仮想空間内で剛体１１を回転させたときの相対位置関係と等価になる。

力演算部４３が、粒子１２の間に作用する力、及び粒子１２と剛体１１との間に作用する力を算出する。この力の算出には、粒子１２の間の相互作用ポテンシャル、剛体１１が周囲の粒子１２に与えるポテンシャルを用いる。このポテンシャルとして、例えばレナードジョーンズ型ポテンシャルを用いることができる。

粒子状態演算部４４が、複数の粒子１２（図１）の各々について、運動方程式に基づいて、現在の状態から１タイムステップ後の状態を求める。粒子１２の状態には、位置及び速度が含まれる。

粒子状態更新部４５が、粒子１２の各々の状態を更新する。具体的には、粒子状態演算部４４で求められた粒子１２の状態を、粒子１２の現在の状態とし、粒子情報保持部５３に登録する。更新後の粒子１２の位置が、部分解析領域１５（図１）の周期境界を越えて外に移動してしまった場合には、周期境界条件を適用して部分解析領域１５の内部に戻す。具体的には、部分解析領域１５の外に移動した粒子１２の位置ベクトル及び速度ベクトルを、中心軸１０（図１）を回転中心として部分解析領域１５の中心角と等しい角度だけ回転させる。

剛体状態演算部４６が、運動方程式に基づいて、剛体１１（図１）の現在の状態から１タイムステップ後の剛体１１の状態を求める。剛体１１の状態には、中心軸１０（図１）を回転中心とした剛体１１の回転速度、及び基準位置からの回転角を含む。

剛体状態更新部４７が、剛体１１の状態を更新する。具体的には、剛体状態演算部４６で求められた剛体１１の状態を、剛体１１の現在の状態とし、剛体１１の状態を剛体情報保持部５２に登録する。

終了条件判定部４８が、シミュレーション処理を終了するか否かを判定する。例えば、予め設定されているタイムステップ数の演算を行うとシミュレーション処理を終了させる。

次に、図５を参照して実施例によるシミュレーション方法について説明する。
図５は、実施例によるシミュレーション方法のフローチャートである。まず、部分解析領域取得部３１が部分解析領域１５（図１）の情報を取得し、この情報を部分解析領域情報保持部５１に登録する（ステップＳ０１）。剛体取得部３３が、剛体１１（図１）の形状を定義する情報を取得する（ステップＳ０２）。取得した情報に基づいて、剛体取得部３３がＳＤＦを生成し、剛体情報保持部５２に登録する（図テップＳ０３）。

粒子取得部３２が粒子１２の初期状態を取得し、粒子情報保持部５３に登録する（ステップＳ０４）。粒子１２の初期状態には、粒子１２の初期位置及び初期速度が含まれる。なお、粒子１２の質量も同時に取得される。力演算部４３が、粒子間距離に基づいて粒子１２の各々に働く力を求める（ステップＳ０５）。さらに、粒子剛体距離演算部４１が、粒子１２の各々について粒子剛体間距離を求め、力演算部４３が、粒子剛体間距離に基づいて粒子１２及び剛体１１に働く力を求める（ステップＳ０６）。

剛体状態演算部４６が、剛体１１の現在の回転速度及び回転角を初期値とし、剛体１１に働く力及び剛体１１の慣性モーメントに基づいて運動方程式を解くことにより、剛体１１の１タイムステップ後の回転速度及び回転角を求める（ステップＳ０７）。

粒子状態演算部４４が、粒子１２の各々について、現在の位置及び速度を初期値とし、粒子１２に働く力及び粒子１２の質量に基づいて運動方程式を解くことにより、粒子１２の１タイムステップ後の位置及び速度を求める。さらに、粒子座標変換部４２が、剛体１１の回転角に基づいて粒子１２の各々の位置及び速度の座標変換を行う。粒子１２の位置及び速度を更新する（ステップＳ０８）。具体的には、座標変換後の粒子１２の位置及び速度を、粒子１２の現在の状態とする。

剛体状態更新部４７が、剛体１１の回転速度及び回転角を更新する（ステップＳ０９）。具体的には、ステップＳ０７で求められた剛体１１の回転速度及び回転角を、剛体１１の現在の状態とする。

粒子状態更新部４５が、部分解析領域１５の外に粒子１２が存在するか否かを判定する（ステップＳ１０）。部分解析領域１５の外に粒子１２が存在する場合には、部分解析領域１５の外の粒子１２に周期境界条件を適用して、位置及び速度を変換する（ステップＳ１１）。部分解析領域１５の外に存在するすべての粒子１２に対してステップＳ１１を実行する。

終了条件判定部４８が、終了条件が満たされるか否かを判定する（ステップＳ１２）。終了条件が満たされない場合には、ステップＳ０５からの処理を繰り返す。終了条件が満たされる場合には、出力部３４がシミュレーション結果を出力装置６２に出力する（ステップＳ１３）。

次に、実施例によるシミュレーション方法の持つ優れた効果について説明する。
本実施例によると、剛体１１（図１）を多数の粒子の集合体として表す代わりに、ＳＤＦを用いてその形状を表している。このため、剛体１１を構成する多数の粒子について１つずつ演算を行う場合に比べて、演算負荷が軽減される。

また、本実施例では剛体１１の回転対称性に着目して解析領域を分割し、部分解析領域１５内について演算を行っている。このため、解析領域の全域について演算を行う場合と比べて、演算負荷が軽減される。

さらに、本実施例では、剛体１１及び複数の粒子１２を配置する仮想空間内で剛体１１が回転した場合に、剛体１１の形状を定義するＳＤＦを再生成する代わりに、粒子１２の位置及び速度に対して座標変換を行っている。剛体１１の回転に応じてＳＤＦを再生成する必要がないため、ＳＤＦの生成に要する演算負荷が軽減される。

実施例によるシミュレーション方法の効果を確認するために、実際にシミュレーションを行った。次に、図６Ａ及び図６Ｂを参照して、シミュレーション結果について説明する。シミュレーションでは、風車形状の回転剛体を流れ場に設置し、流体からの作用による剛体の角速度の時間変化求めた。シミュレーションのタイムステップ数を約１００万とし、解析領域全域内の粒子数を約３２０００個とした。相互作用ポテンシャルとしてレナードジョーンズ型ポテンシャルを用いた。

図６Ａは、シミュレーションにより求められた角速度の時間変化を示すグラフである。横軸は経過時間（タイムステップ数）を無次元で表し、縦軸は角速度を単位「ラジアン／タイムステップ」で表す。

図６Ａのグラフ中の太い実線、細い実線は、それぞれ部分解析領域１５（図１）の中心角を１８０度及び９０度にした場合のシミュレーション結果を示す。参考のために、解析領域全域についてシミュレーションを行った結果を破線で示す。

図６Ａに示すように、解析領域全体でシミュレーションした場合と、解析領域全体から切り出した部分解析領域１５内でシミュレーションした場合とで、有意な差は見られなかった。

図６Ｂは、計算に要した時間を示すグラフである。横軸は、部分解析領域１５（図１）の中心角を単位「度」で表し、縦軸は計算時間を単位「秒」で表す。中心角を３６０度にしたときのシミュレーションは、解析領域全域をシミュレーション対象にした場合に相当する。部分解析領域１５の中心角を小さくするに従って計算時間が短くなっていることがわかる。

図６Ａに示したように、解析領域全体から切り出した部分解析領域１５（図１）についてシミュレーションを行っても、十分なシミュレーションの精度が得られることが確認された。さらに、図６Ｂに示したように、部分解析領域１５についてシミュレーションを行うことにより、計算時間を短縮できることが確認された。

次に、実施例の変形例について説明する。
上記実施例では、シミュレーション対象の剛体１１（図１）として、回転軸１１Ａと回転翼１１Ｂとからなる剛体を採用したが、その他に、ｎ回回転対称性（ｎは２以上の整数）を持つ種々の形状の剛体をシミュレーション対象とすることができる。また、上記実施例では、剛体１１と粒子１２との挙動のシミュレーションに分子動力学法を適用したが、その他の方法、例えば繰り込み群分子動力学法（ＲＭＤ法）、離散要素法（ＤＥＭ法）、ＳＰＨ法、ＭＰＳ法等の粒子法を用いてもよい。

上述の各実施例は例示であり、異なる実施例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例の同様の構成による同様の作用効果については実施例ごとには逐次言及しない。さらに、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

１０ 中心軸
１１ 剛体
１１Ａ 回転軸
１１Ｂ 回転翼
１２ 粒子
１２Ａ 周期境界の近傍の粒子
１２ａ 周期境界の近傍の粒子の位置の回転対称位置
１２Ｂ 周期境界を越えて部分解析領域の外に移動する粒子
１２ｂ 粒子１２Ｂの移動後の位置に対する回転対称位置に侵入する粒子
１３、１４ 半平面
１５ 部分解析領域
２０ ＳＤＦを定義するための直交格子のノード
３０ 処理装置
３１ 部分解析領域取得部
３２ 粒子取得部
３３ 剛体取得部
３４ 出力部
４０ 数値演算部
４１ 粒子剛体距離演算部
４２ 粒子座標変換部
４３ 力演算部
４４ 粒子状態演算部
４５ 粒子状態更新部
４６ 剛体状態演算部
４７ 剛体状態更新部
４８ 終了条件判定部
５０ 記憶装置
５１ 部分解析領域情報保持部
５２ 剛体情報保持部
５３ 粒子情報保持
６１ 入力装置
６２ 出力装置

Claims (7)

1. 中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域内の複数の粒子の挙動をシミュレーションする方法であって、
   前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義し、
   前記剛体の形状に基づいて、前記剛体の形状を表す符号付距離関数を生成し、
   前記部分解析領域の前記半平面に対応する面において、前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用し、前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションにより求めるシミュレーション方法。
2. 前記剛体の形状、前記剛体の挙動、及び前記解析領域が前記中心軸に関してｎ回回転対称性（ｎは２以上の整数）を持ち、２つの前記半平面のなす角度が３６０°／ｎの整数倍である請求項１に記載のシミュレーション方法。
3. 前記粒子の挙動のシミュレーションに分子動力学法を適用する請求項１または２に記載のシミュレーション方法。
4. 中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域の内部に存在する複数の粒子の挙動をシミュレーションする処理装置を有し、
   前記処理装置は、
   入力装置から、前記剛体の形状を定義する形状データが入力されると、入力された前記形状データに基づいて、前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義し、
   前記形状データに基づいて、前記剛体の形状を表す符号付距離関数を生成し、
   前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用し、前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションにより求め、
   シミュレーション結果を出力装置に出力するシミュレーション装置。
5. 前記剛体の形状、前記剛体の挙動、及び前記解析領域が前記中心軸に関してｎ回回転対称性（ｎは２以上の整数）を持ち、前記処理装置は、２つの前記半平面のなす角度が３６０°／ｎの整数倍になるように前記部分解析領域を定義する請求項４に記載のシミュレーション装置。
6. 前記処理装置は、前記粒子の挙動のシミュレーションに分子動力学法を適用する請求項４または５に記載のシミュレーション装置。
7. 中心軸に関して回転対称な形状を持つ剛体の前記中心軸のまわりの回転運動、及び前記剛体が収容される解析領域内の複数の粒子の挙動をシミュレーションする機能を、コンピュータに実現させるプログラムであって、
   前記解析領域を、前記中心軸を境界とする２つの半平面で分割して、一方の領域を部分解析領域として定義する機能、
   前記部分解析領域の前記半平面に対応する面において、前記中心軸を回転中心とした回転方向に関して周期性を持つように周期境界条件を適用する機能、
   前記剛体の形状に基づいて、前記剛体の形状を表す符号付距離関数を生成する機能、及び
   前記部分解析領域内において、前記符号付距離関数を用いて前記部分解析領域の内部の前記剛体及び複数の粒子の挙動をシミュレーションする機能
   をコンピュータに実現させるプログラム。

Images