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

Description

本発明は、シミュレーション装置、シミュレーションプログラム及びシミュレーション方法に関する。

近年の計算機パワーの向上に伴い、シミュレーションの手法も徐々に発展してきた。その結果として、様々な応用分野にシミュレーションが用いられるようになってきた。

流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。

しかしながら、有限要素法又は有限体積法等の格子を用いる手法では、自由表面などの界面の存在する問題や、流体・構造連成問題などの移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難である場合が多い。

これに対して格子を用いないＭＰＳ（Moving Particle Semi-implicit又はMoving Particle Simulation）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等の粒子法は、計算対象物である流体を粒子の集まりとして表し、移動境界の取り扱いに特別な処置を必要としない。それゆえ、粒子法は近年広く用いられるようになっている。これら、粒子分布を用いて解析対象を表現する手法においては、ある粒子から影響半径（カーネル半径）と呼ばれる距離ｈ以内にある粒子を近傍粒子してその粒子の情報を用いている。

例として、ＳＰＨ法における運動方程式の離散化の（式１）を示す。
ここで、下付添え字ａはある粒子を表し、下付添え字ｂは近傍粒子を表す。
は、それぞれ粒子ａの位置ベクトル、速度ベクトル、密度、圧力、質量である。
は、カーネル関数で、粒子の分布から連続場を構成するのに用い、以下の３次のスプライン関数（式２）等がよく使われる。

また、弾性体などにおいては、下記の（式３）であらわされる変形勾配テンソルを離散化して計算することで加速度を求める。
ここで、Ｘ_jは位置間の距離のｊ成分、ｕ_iは変位のｉ成分、Ｉは恒等テンソルである。

粒子法においては、ある粒子ａの変形勾配テンソルは、近傍粒子ｂからの影響を重ね合わせて以下の（式４）で表すことができる。
ここで、ｄｕは変位を表す。

しかしながら、ここでカーネル関数を基準点のもののまま用いると、粒子の相対位置が影響半径以上に離れた場合相互作用がなくなってしまう。そのため、変形量に応じて影響半径の形状を変形させ、基準となる粒子との相互作用を持ち続ける必要がある。

そこで、弾性体などの解析において、最初の近傍の計算が進み粒子位置が変更した後でも、引き続き近傍として取り扱う手法が知られている（例えば、非特許文献１参照。）。

また、各粒子に対して相互作用を及ぼす粒子を識別する技術が開示されている（特許文献１参照。）。

また、ある粒子が近傍に分布する粒子を探索するためには、他の全ての粒子との距離を求めて近傍かどうかの判定を繰り返す必要があるが、空間中に分布する粒子を、近傍粒子を探索して新しい属性情報を計算する粒子と、近傍粒子探索を行わず粒子自身の属性情報から新しい属性情報を計算する粒子の二種類に分け、それぞれに対応する計算処理を実行することによって、近傍粒子探索の回数を減らす技術が開示されている（特許文献２参照。）。

また、粒子間の距離に対応した力またはポテンシャルの関数計算において。関数の微分係数を用い、補間法で新たなポテンシャルを求める技術が開示されている（特許文献３参照。）。

特開平４−１８１３５０号公報 特許第２９２９８４０号明細書 特許第３２７７７９９号明細書

皆木浩幸、野口裕久、「高精度弾塑性大変形解析のためのＳＰＨの拡張に関する検討」、日本機械学会第１７回計算力学講演会講演論文集、２００４年、ｐ．７６７-７６８

基準の近傍粒子を用いて計算を行う場合には、初期位置における距離で求めた値と変形後の値は一致するが、その微分値は初期位置のものとは一致せず、相互作用を求めるときに基準位置に直して計算する必要がある。そのため、たとえば応力を基準位置における応力値であるピオラキルヒホッフ（Piola-Kirchhoff）応力として計算する必要がある。

しかしながら、上述のような従来の技術においては、基準位置に近傍でなかった粒子が新たに接触してくる影響を考慮できないという接触問題など、弾性体の変形後の座標で作用する応力値を考慮した計算を行うのが困難である、という問題点がある。

１つの側面では、本発明は、弾性体の変形後の座標で作用する応力値を考慮した計算を容易に行うことが可能なシミュレーション装置、シミュレーションプログラム又はシミュレーション方法を提供することを目的とする。

本発明は、上記課題を解決するため、下記のような構成を採用した。
すなわち、本発明の一態様によれば、シミュレーションプログラムにおいて、コンピュータに、連続体を粒子の集まりで表現した連続体モデルの基準粒子から影響半径内の１つ以上の第１の近傍粒子について、前記連続体モデルの変位後の前記第１の近傍粒子の位置とカーネル関数に基づいて、前記基準粒子と前記第１の近傍粒子との第１の相互作用を算出させ、前記変位後に前記基準粒子から影響半径内に入った１つ以上の第２の近傍粒子について、前記連続体モデルの変位後の位置とカーネル関数に基づいて、前記基準粒子と前記第２の近傍粒子との第２の相互作用を算出させ、前記第１の相互作用及び前記第２の相互作用に基づいて、前記基準粒子と前記基準粒子の現時点の近傍粒子との相互作用を算出させ、前記第１の相互作用を算出させる場合、前記カーネル関数の変形勾配テンソルを用いることを特徴とするシミュレーションプログラムが提供される。

開示されたシミュレーション装置、シミュレーションプログラム及びシミュレーション方法によれば、変形勾配テンソルを用いて粒子が移動した後の変形に応じて影響半径が変化することを想定した相互作用を計算し、変形後の座標での応力値であるコーシー応力を算出することができ、変形後に新たに近傍に来た粒子との相互作用の計算を組み合わせることにより自己接触問題等を解くことを可能とする。特に、弾性体を取り扱う場合においては、弾性体の挙動を示すために変形勾配テンソルを用いることが多くあるため、別途特別な計算を行うことなく算出が可能となり、特別計算量を増やすことなく計算が可能である。

本発明を適用したシミュレーション装置の構成例を示す図である。 第１の実施の形態におけるシミュレーション処理の流れを示すフローチャートである。 弾性体の変形に伴う影響半径の変化を示す図である。 弾性体が大変形した場合の粒子の状態を示す図である。 ２種類の近傍を用いた場合の概念を示す図である。 自己接触計算を行う場合の概念を示す図である。 第２の実施の形態におけるシミュレーション処理の流れを示すフローチャートである。 情報処理装置の構成図である。

以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
本発明は、シミュレーション装置、コンピュータに実行させるシミュレーションプログラム及びシミュレーション方法であって、シミュレーションの対象である連続体を粒子の集まりとして捕らえる。

本発明のシミュレーション装置、シミュレーションプログラム及びシミュレーション方法は、粒子法を用いた数値計算プログラムにおける、近傍粒子の取り扱い手法である。特に、弾性体などの連続体の運動を粒子法を用いて数値的に計算する手法であって、さらに詳しくは近傍粒子の取り扱いにおいて基準位置における近傍粒子を考慮したまま変形後のコーシー応力を求め、座標変形後の近傍と組み合わせて計算を行う手法である。

図１は、本発明を適用したシミュレーション装置の構成例を示す図である。
図１において、シミュレーション装置１００は、処理部１０１、格納部１０２及び出力部１０３を備え、初期条件１１１に基づいて粒子法による数値計算を行って、シミュレーション結果１１２を出力する。

格納部１０２は、本発明を適用したシミュレーションプログラムを実行するための各計算式の情報を格納する。

処理部１０１は、本発明を適用したシミュレーション処理を実行する。
処理部１０１が実行するシミュレーション処理は、連続体を粒子の集まりで表現した連続体モデルの基準粒子から影響半径内の１つ以上の第１の近傍粒子について、この連続体モデルの変位後の第１の近傍粒子の位置とカーネル関数に基づいて、基準粒子と第１の近傍粒子との第１の相互作用を算出する。また、シミュレーション処理は、変位後に基準粒子から影響半径内に入った１つ以上の第２の近傍粒子について、連続体モデルの変位後の位置とカーネル関数に基づいて、基準粒子と第２の近傍粒子との第２の相互作用を算出する。そして、これら第１の相互作用及び第２の相互作用に基づいて、基準粒子と基準粒子の現時点の近傍粒子との相互作用を算出する。なお、第１の相互作用を算出させる場合、カーネル関数の変形勾配テンソルを用いる。

このシミュレーション処理における基準粒子と現時点の近傍粒子との相互作用の算出は、第１の相互作用と第２の相互作用を任意の割合で算出させることができる。また、第１の相互作用と第２の相互作用の算出の割合は、連続体モデルの変位の大きさに応じて定めることができる。

出力部１０３は、処理部１０１が実行したシミュレーション結果１１２を出力する。
本発明を適用したシミュレーション処理は、基準位置での近傍粒子を取り扱う場合に、「弾性体の変形によって、影響半径は弾性体の粒子の変形にしたがって変化する」と仮定することにより、変形後のカーネルに変形の大きさを示すパラメータである変形勾配テンソルを用いた補正項を導入し、下記（式５）として計算する。

なお、添え字iniがついた変数は基準位置での値であり、
は、ａ粒子の（ｎ−１）ステップでの変形勾配テンソルである。

これにより基準位置でのカーネルと影響半径を用いて、弾性体の変形後のカーネルの微分の正確な値を求めることができ、変形後の座標系における応力であるコーシー応力の算出が可能となり、現在の座標系における他の応力との計算を容易にすることができる。

また、弾性体が大変形をした場合、その基準位置の近傍による近傍粒子の寄与と、変形後の近傍による近傍粒子の寄与を一定の割合で変化させながら使うことにより、基準位置に近傍でなかった粒子が新たに接触してくる影響等を考慮し、弾性体大変形時の挙動を求める。このとき、この割合を決める方法としては変形勾配テンソルを用いることができる。

図２は、第１の実施の形態におけるシミュレーション処理の流れを示すフローチャートである。

まず、ステップＳ２０１において、入力データの取得を行う。入力データは、弾性体等をモデル化した連続体粒子モデルを表現するために必要なデータである。例えば、各粒子の中心座標、移動速度、影響半径（カーネル半径）、密度、質量、変形勾配テンソル、材料物性値、材料温度等である。

ステップＳ２０２において、各粒子について、ステップＳ２０１で入力された中心座標に基づき、各粒子から所定の距離内に存在する近傍粒子の情報を取得しメモリに格納する。

図３は、弾性体の変形に伴う影響半径の変化を示す図である。
図３において、弾性体を構成する任意の粒子である基準位置（以下、基準粒子３０１ともいう。）から影響半径３０２内には、複数の近傍粒子３０３が定義されている。弾性体が変形すると、その変形に伴い影響半径３０２も変形する。

図３に示すように、影響半径３０２が弾性体の変形に応じて変化していると仮定すると、例えばある粒子間の距離が特定の方向のＮ倍離れた場合、影響半径はその方向に同じくＮ倍伸びると想定することができる。このとき、これらの２粒子間の体積の寄与は、基準粒子３０１で計ったものと同じ値を持つことになる。そこで、図２のステップＳ２０３において、この想定によりカーネル関数を変形させることなく基準粒子３０１の相対位置と基準粒子３０１のカーネル関数にて変形後の影響を計算する。ただし、微分の値はその限りではないため、上記（式４）の補正項を用いて補正を行い、あるステップ数ｎにおける変形勾配テンソルを下記（式６）とする。

これによりカーネル関数の微分を用いる変形勾配テンソルを変形後の位置を用いずに算出することができ、コーシー応力を算出ができる。

図４は、弾性体が大変形した場合の粒子の状態を示す図であり、図５は、２種類の近傍を用いた場合の概念を示す図である。

上述の手法を用いた場合、弾性体の変形が進むに連れて、図４に示したように、当初は遠方にあった粒子（遠方粒子４０１）が基準粒子３０１のごく近くまで来た場合や、たとえ粒子同士が重なり合った場合にも影響を及ぼしあわないことがあり、これは現実の状態と即していないと考えられる。

そこで、ステップＳ２０４において、図５に示すように、当初から基準粒子３０１と近傍粒子３０３とが近接していた場合と、弾性体の大変形後に基準粒子３０１と遠方粒子４０１とが近接した場合の２種類を用いて、両者を一定の割合で重ね合わせる。なお、最初の１ステップ目においては、両影響は同量となる。また、変形途中の近傍粒子３０３を新たな影響を及ぼしあう基準のものとして３つ以上の近傍状態を重ね合わせても良い。重ね合わせる指標としては、例えば以下の通りである。

基準粒子３０１からの変形量が小さい場合には、最初遠方にあった遠方粒子４０１からの影響力は小さいと考えられるので、変形量の指標となる１つである粒子間の相対位置座標を使ってもよい。例えば、基準粒子３０１の近傍粒子３０３の相対位置座標のうち、最も近いものｘ_nearと影響半径ｈを用いて、基準の近傍の割合Ｐを以下のように決める。
Ｐ＝１ （ｘ_near／２ｈ＜０．５）
Ｐ＝ｘ_near／２ｈ （０．５≦ｘ_near／２ｈ＜１）
Ｐ＝０．５ （ｘ_near／２ｈ＞１）

次に、ステップＳ２０５において、弾性体粒子に設定された物理モデルに従って、弾性体粒子の運動を解き応力を求める。

具体的には、弾性体においては、上記（式４）より求められる変形勾配テンソルから、
上記（式７）（式８）（式９）を解くことにより、コーシー応力σを算出する。Ｃは、右コーシーグリーン変形テンソル（right Cauchy-Green deformation tensor）、Ｅは、グリーンラグランジュ歪テンソル（Green-Lagrange strain tensor）に値し、ｊは、体積変化率（＝｜Ｆ｜）、Ｑは、弾性ポテンシャルである。

弾性ポテンシャルとしては、例えば（式１０）で表されるMooney-Rivlinポテンシャルが良く知られている。

ここで、ｃ₁、ｃ₂は、Mooney-Rivlinパラメータで、Ｉ_c、ＩＩ_cは、右コーシーグリーン変形テンソルＣの主不変量であり、

そして、ステップＳ２０６において、（式１２）までで求まったコーシー応力より加速度を算出する。
ここで、σはコーシー応力テンソルである。

ステップＳ２０７において、最後に求められた加速度を含む時間微分項を用いて連続体粒子の物理量を時間積分し、時刻を時間刻みの分だけ進めると１ステップ分の計算が終了する。

ステップＳ２０８において、更新された近傍粒子３０３の情報を取得する。
そして、ステップＳ２０９において、一定ステップ数の経過、又は一定のシミュレーション内時間の経過など、必要に応じて計算結果を外部ファイルとして出力する。

ステップＳ２１０において、指定したステップ数が終了したかを判断し、終了していなければ（ステップＳ２１０：Ｎｏ）、ステップＳ２０３以降を指定したステップ数繰り返すことで第１の実施の形態におけるシミュレーション処理が終了する。

図６は、自己接触計算を行う場合の概念を示す図である。
自己接触問題を解く場合、図６に示すように、当初の近傍粒子３０３においては上述のように連続体粒子の運動を解き、弾性体の変形後に新たに近接した近傍粒子６０１に関しては、基準粒子３０１と近傍粒子６０１の接触が起こったとして反発力を求める。この場合、当初及び変形後の両方とも、若しくは当初のみ近接している近傍粒子３０３に関しては、上述したように何らかの割合により異なる影響力を与えてもよい。そして、反発力の算出にはペナルティ法などを用いてもよい。この場合、それぞれの影響力は弾性体の運動と接触という異なるものであるため、その影響力の割合は変化させないでもよい。

図７は、第２の実施の形態におけるシミュレーション処理の流れを示すフローチャートである。

第２の実施の形態におけるシミュレーション処理は、第１の実施の形態におけるシミュレーション処理において用いた（式５）を用いず、基準座標での影響半径を求めるための座標変換を行う処理である。

ステップＳ２０１乃至Ｓ２０３は、第１の実施の形態におけるシミュレーション処理と同様である。

次に、基準座標を用いた場合、その座標系は初期位置を基準としたものとなる。そこで、ステップＳ７０１において、影響半径３０２を変換して基準座標での影響半径３０２に変換を行う。このとき、変形勾配テンソルをＦとして、影響半径３０２は、
となる（
は、基準座標系を用いた場合の影響半径）。

次に、ステップＳ２０５において、第１の実施の形態と同様に、当初及び変形後の近傍粒子３０３の両方を用いて応力を計算する。計算式については基本的には第１の実施の形態と同様であるが、変形勾配テンソルについて、補正を入れていない変形勾配テンソルを用いる。そのため、応力は基準座標系における応力となり、それは第１ピオラキルヒホッフ応力となる。

他方、当初及び変形後の二つの影響半径３０２を重ね合わせる場合に、現在の座標系で同様の式を解いて求めた応力や、接触を考えた場合の応力は現在の座標系の応力値であるコーシー応力となっており、その二つを単純に重ね合わせることができない。

そこで、ステップＳ７０２において、
（σ：Cauchy応力、
：第１ピオラキルヒホッフ応力）という（式１５）を用いて、基準座標系を用いた応力に変換して計算を行う。

その後のステップＳ２０６以降については、第１の実施の形態と同様である。
このように、本発明の実施の形態は、粒子法を用いた数値計算において広く利用可能である。例えば、鋳造過程を解析することで製品設計に用いることができる。弾性体に対して適用することで、製品設計の際に封止ゲルの形状などを適切に決定することができる。

図１のシミュレーション装置は、例えば、図８に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図８の情報処理装置は、ＣＰＵ（Central Processing Unit）８０１、メモリ８０２、入力装置８０３、出力装置８０４、外部記録装置８０５、媒体駆動装置８０６及びネットワーク接続装置８０７を備える。これらはバス８０８により互いに接続されている。

メモリ８０２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ８０１は、メモリ８０２を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ８０２は、図１の格納部１０２としても使用できる。

入力装置８０３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置８０４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置８０４は、図１の出力部１０３としても使用できる。

外部記録装置８０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記録装置８０５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記録装置８０５にプログラム及びデータを格納しておき、それらをメモリ８０２にロードして使用することができる。

媒体駆動装置８０６は、可搬型記録媒体８０９を駆動し、その記録内容にアクセスする。可搬型記録媒体８０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体８０９には、Compact Disk Read Only Memory （ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体８０９にプログラム及びデータを格納しておき、それらをメモリ８０２にロードして使用することができる。

このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ８０２、外部記録装置８０５、及び可搬型記録媒体８０９のような、物理的な（非一時的な）記録媒体が含まれる。

ネットワーク接続装置８０７は、通信ネットワーク８１０に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置８０７を介して受け取り、それらをメモリ８０２にロードして使用することができる。ネットワーク接続装置８０７は、図１の出力部１０３としても使用できる。

開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。

１００ シミュレーション装置
１０１ 処理部
１０２ 格納部
１０３ 出力部
１１１ 初期条件
１１２ シミュレーション結果
３０１ 基準粒子
３０２ 影響半径
３０３ 近傍粒子
４０１ 遠方粒子
６０１ 近傍粒子
８０１ ＣＰＵ（Central Processing Unit）
８０２ メモリ
８０３ 入力装置
８０４ 出力装置
８０５ 外部記録装置
８０６ 媒体駆動装置
８０７ ネットワーク接続装置
８０８ バス
８０９ 可搬型記録媒体
８１０ 通信ネットワーク

Claims (5)

1. コンピュータに、
   連続体を粒子の集まりで表現した連続体モデルの基準粒子から影響半径内の１つ以上の第１の近傍粒子について、前記連続体モデルの変位後の前記第１の近傍粒子の位置とカーネル関数に基づいて、前記基準粒子と前記第１の近傍粒子との第１の相互作用を算出させ、
   前記変位後に前記基準粒子から影響半径内に入った１つ以上の第２の近傍粒子について、前記連続体モデルの変位後の位置とカーネル関数に基づいて、前記基準粒子と前記第２の近傍粒子との第２の相互作用を算出させ、
   前記第１の相互作用及び前記第２の相互作用に基づいて、前記基準粒子と前記基準粒子の現時点の近傍粒子との相互作用を算出させ、
   前記第１の相互作用を算出させる場合、前記カーネル関数の変形勾配テンソルを用いることを特徴とするシミュレーションプログラム。
2. 前記基準粒子と前記基準粒子の現時点の近傍粒子との相互作用の算出は、前記第１の相互作用の算出結果と前記第２の相互作用の算出結果を任意の割合で算出させることを特徴とする請求項１記載のシミュレーションプログラム。
3. 前記第１の相互作用の算出結果と前記第２の相互作用の算出結果の算出の割合は、前記連続体モデルの変位の大きさに応じて定めることを特徴とする請求項２記載のシミュレーションプログラム。
4. コンピュータが、
   連続体を粒子の集まりで表現した連続体モデルの基準粒子から影響半径内の１つ以上の第１の近傍粒子について、前記連続体モデルの変位後の前記第１の近傍粒子の位置とカーネル関数に基づいて、前記基準粒子と前記第１の近傍粒子との第１の相互作用を算出し、
   前記変位後に前記基準粒子から影響半径内に入った１つ以上の第２の近傍粒子について、前記連続体モデルの変位後の位置とカーネル関数に基づいて、前記基準粒子と前記第２の近傍粒子との第２の相互作用を算出し、
   前記第１の相互作用及び前記第２の相互作用に基づいて、前記基準粒子と前記基準粒子の現時点の近傍粒子との相互作用を算出し、
   前記第１の相互作用を算出する場合、前記カーネル関数の変形勾配テンソルを用いることを特徴とするシミュレーション方法。
5. 連続体を粒子の集まりで表現した連続体モデルの基準粒子から影響半径内の１つ以上の第１の近傍粒子について、前記連続体モデルの変位後の前記第１の近傍粒子の位置とカーネル関数に基づいて、前記基準粒子と前記第１の近傍粒子との第１の相互作用を算出する第１相互作用算出手段と、
   前記変位後に前記基準粒子から影響半径内に入った１つ以上の第２の近傍粒子について、前記連続体モデルの変位後の位置とカーネル関数に基づいて、前記基準粒子と前記第２の近傍粒子との第２の相互作用を算出する第２相互作用算出手段と、
   前記第１相互作用算出手段によって算出した第１の相互作用及び前記第２相互作用算出手段によって算出した前記第２の相互作用に基づいて、前記基準粒子と前記基準粒子の現時点の近傍粒子との相互作用を算出する相互作用算出手段と、
   を備え、
   前記第１相互作用算出手段は、前記カーネル関数の変形勾配テンソルを用いて前記第１の相互作用を算出することを特徴とするシミュレーション装置。