粒子モデルの離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置
【課題】転がり抵抗を導入した多角形型粒子モデル、即ち多角形型粒子の回転運動を導入した離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置を提供する。
【解決手段】多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、仮想粒子に物性値を与えるステップと、前記仮想粒子を配置する仮想粒子配置ステップと、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、回転運動を基に立てた運動方程式を解いて個々の粒子の挙動を求める。
【解決手段】多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、仮想粒子に物性値を与えるステップと、前記仮想粒子を配置する仮想粒子配置ステップと、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、回転運動を基に立てた運動方程式を解いて個々の粒子の挙動を求める。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子の転がり抵抗を導入した多角形型粒子の挙動の離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置に関する。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、粒子の挙動を把握することが重要な課題である。従来は、このような粒子の挙動を試行錯誤の実験により把握することが多かった。このような場合、現実に用いる粒子の組成や、実験装置や実験環境を用意しなければ、所望の粒子の挙動解析を行うことが出来ない。また、これらの仕様を変更するたびに、実験しなければならず、その都度時間やコストの面を考慮しなくてはならなかった。
【０００３】
近年では、粒子の挙動解析を実験的にではなく、数式化若しくはモデル化し、ほぼ同じ法則に支配されるシステムの挙動をコンピュータ上で模擬すること、すなわちシミュレーションが広く利用されている。シミュレーション技術を利用することにより、粒子の挙動を予測することができる。また、粒子（粉体）を投入する条件やパラメータを変えて、同じシミュレーション計算を繰り返し行うことにより、様々な粒子組成や装置設計・制御体型の性能を評価することができる。
【０００４】
数式化した場合の粉体の挙動解析シミュレーション手段として、粉体毎の挙動を解析する離散要素（若しくは個別要素）法（Ｄｉｓｃｒｅｔｅ Ｅｌｅｍｅｎｔｓ Ｍｅｔｈｏｄ、若しくはＤｉｓｔｉｎｃｔ Ｅｌｅｍｅｎｔｓ Ｍｅｔｈｏｄ）解析シミュレーション手法（以下本明細書ではＤＥＭシミュレーションと表記する。）と、解析対象の粒子の形状をモデル化して挙動を解析する手法などが挙げられる。
【０００５】
ＤＥＭシミュレーションの原理は、球状粒子に見立てた２つのボール間に作用する法線方向（図１（ａ）参照）、剪断方向（図１（ｂ）参照）のボール間接触力について、図１に示すような、弾性的性質を表すバネ２，４と、非弾性的性質を表すダッシュポット３，５と、ボール１間の摩擦的性質を表すスライダ６（剪断方向）とで構成されるフォークトモデルによって説明される。前記フォークトモデルにより、円筒型の容器の壁とボールの接触、ボール同士の接触を判定し、接触している場合のボールと壁、またはボール同士との法線、剪断方向の接触力と摩擦力を図１（ａ）、図１（ｂ）のフォークトモデルに基づいて演算し、バネ２，４の弾性係数はヘルツの弾性接触理論より求められ、ダッシュポット３，５の粘性減衰係数（法線方向）、摩擦係数（剪断方向）は実験により求められる。そして、その接触力に基づいて全ボールの加速度、速度、変位を各離散化時間で演算し、この演算を繰り返すことによって任意の時間のボール群全体の運動挙動がシミュレーションされる。
【０００６】
上記ＤＥＭシミュレーションによれば、粒子と粒子との間の衝突判定が容易且つ、計算負荷を小さくできる。ＤＥＭシミュレーションを使用した粒子の挙動解析の例として、特開２００６−２５９９１１号公報（特許文献１参照）に開示されているようなトナーに使用する粒子の挙動解析等が挙げられる。
【０００７】
しかしながら、粒子の形状は現実的には球形ではないことが多く、球形を仮定すると、接触角度が精度良く形成されないこと、粉体の流動が実際よりも大きく見積もられてしまうといったことなどが課題である。また、特許文献１に記載の方法は、あくまで粒子を球形で仮定していること、粒子間の点若しくは面接触を考慮していないため、例えば金属粒子の最密充填（高濃度充填）には適していない。
【０００８】
一方、従来のＤＥＭシミュレーションによる粒子挙動解析の改善策として、例えばマツシマら（非特許文献１参照）が提案した非球形粒子に複数の球形粒子をつなぎ合わせたクラスターモデルとＤＥＭシミュレーションをベースにした非球形粒子挙動解析方法がある。非特許文献１に記載の方法では、特許文献１のような従来のＤＥＭシミュレーションに比べると、モデル粉体が実際の粉体により近い流動を見積もることが可能になったが、クラスターモデルを使用するため計算体系が小さくなる、即ち一領域あたりのモデル粉体の数が少なく見積もられるといったことが問題である。
【０００９】
近年では、ＤＥＭシミュレーションと、多角形型粒子モデルとを組み合わせたジョン−ポール レイサムら（Ｊｏｈｎ−Ｐａｕｌ Ｌａｔｈａｍ， ｅｔ ａｌ）の提案モデル（非特許文献２参照）がある。非特許文献２に記載の方法は、従来のＤＥＭシミュレーションに比べ、粒子と壁面との接触（粒子姿勢）角度が精度良く形成されないといった問題は解消されたが、該シミュレーションで必要な運動方程式については球の回転運動を考慮しなくてはならないため、計算負荷がかかってしまうなどの問題点がある。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特開２００６−２５９９１１号公報
【非特許文献】
【００１１】
【非特許文献１】マツシマ他（Ｔ． Ｍａｔｓｕｓｈｉｍａ， ｅｔ ａｌ）， ジャーナル オブ エーロスペース エンジニアリング（Ｊｏｕｒｎａｌ ｏｆ Ａｅｒｏｓｐａｃｅ Ｅｎｇｉｎｅｅｒｉｎｇ）， ２００９年， １５−２３頁．
【非特許文献２】ジョン−ポール レイサム他（Ｊｏｈｎ−Ｐａｕｌ Ｌａｔｈａｍ， ｅｔ ａｌ）， ミネラルズ エンジニアリング（Ｍｉｎｅｒａｌｓ Ｅｎｇｉｎｅｅｒｉｎｇ）， ２１巻（２００８）， ７９７−８０５頁．
【発明の概要】
【発明が解決しようとする課題】
【００１２】
上記に述べた問題点に対して望まれる解決案は、計算上は球形粒子を使用して、粒子の転がり抵抗、即ち多角形型粒子の回転運動を導入することである。
【００１３】
以上の事情を鑑み、本発明は、転がり抵抗を導入した多角形型粒子モデル、即ち多角形型粒子の回転運動を導入した多角形型粒子モデルの離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置を提供する。
【課題を解決するための手段】
【００１４】
本発明の上記目的は、多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、仮想粒子に物性値を与えるステップと、前記仮想粒子を配置する仮想粒子配置ステップと、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、並進運動を基に立てた運動方程式を解いて個々の粒子の挙動を求めることを特徴とすることにより効果的に達成される。
【００１５】
更に、前記仮想粒子に与える物性値は、ヤング率、粘性減衰係数、摩擦係数のうち少なくとも１つであることにより、或いは前記解析対象粒子が多角形型粒子であることにより、或いは前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が０のときに、前記粒子同士が面接触と判定されることにより、或いは前記粒子同士が面接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触面を考慮して回転角度を求め、個々の粒子の挙動を求めることにより、或いは前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が測定されるときに、前記粒子同士が点接触と判定されることにより、或いは前記粒子同士が点接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触点及び接触力を考慮して回転角度を求め、個々の粒子の挙動を求めることにより、より効果的に達成される。
【００１６】
また、本発明は、前記離散要素法解析シミュレーション方法を実行するための離散要素法解析シミュレーションプログラム、或いは前記離散要素法解析シミュレーションプログラムを記録した記録媒体により、或いは前記記録媒体を搭載した離散要素法解析シミュレーション装置であって、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備することを特徴とすることにより、効果的に達成される。
【発明の効果】
【００１７】
本発明の離散要素法解析シミュレーション方法は、粒子姿勢角度を測定することで粒子間の面若しくは点接触を判定することが容易であり、精度良く粉体の流動を模擬することができる。
【００１８】
また、従来のＤＥＭシミュレーションと同程度の計算負荷で粉体の挙動をシミュレートすることが可能である。
【図面の簡単な説明】
【００１９】
【図１】従来技術におけるＤＥＭシミュレーションを力学的に模したフォークトモデルである。
【図２】本発明の粒子挙動解析を行う前の前処理を示すフローチャートである。
【図３】本発明における、ＤＥＭシミュレーション方法のフローチャートである。
【図４】本発明の実施形態における、多角形（正六角形）型仮想粒子の概略図である。
【発明を実施するための形態】
【００２０】
以下、本発明を実施するための形態について、図面を参照しながら説明する。
【００２１】
図２は、本発明における粒子挙動解析を行う前の前処理の流れを示すフローチャートである。
【００２２】
先ず、挙動解析を行う前に、本発明に必要なシミュレーション用粒子、即ち仮想粒子に物性値を与える（ステップＳ１）。具体的には、バネ定数、粘性減衰係数、摩擦係数を仮想粒子に与えるようにする。バネ定数については、線形バネと非線形バネがある。非線形バネの場合、ヤング率などを入力する。
【００２３】
次に、仮想粒子を生成するための処理を行う（ステップＳ２）。処理手順は、従来のＤＥＭシミュレーションと同様に、仮想粒子を配置する範囲などを入力し、仮想粒子平均径や仮想粒子平均間隔を設定する。
【００２４】
次いで、仮想粒子数を設定（ステップＳ３）し、各仮想粒子の直径の設定、各仮想粒子を配置する座標位置の設定、並びに各仮想粒子の直径及び座標位置を出力（ステップＳ４）し、前処理を終了する。
【００２５】
図３は、本発明における粒子挙動解析の処理手順をフローチャートの形式で示している。また、図４には本発明の実施形態における、２次元体系の多角形（正六角形）型仮想粒子の概略図を示す。なお、図４において、正六角形型仮想粒子を示したが、実際に本発明のシミュレーション方法を使用する際は、特に粒子の形状は限定されない。
【００２６】
まず、解析対象となる各粒子に初期条件を与えて、初期化する（ステップＳ１１）。ここで言う初期条件とは、前述の仮想粒子の処理で予め計算した物性値等のことである。次に、粒子同士の接触（衝突）を判定する（ステップＳ１２）。ここで、前記判定は、球形粒子を仮定して行う（図４のパターンＡ及びＢの「計算」参照）。
【００２７】
ここで言う粒子には、解析対象粒子と、図４で表す仮想粒子の双方を含む。なお、前記解析対象粒子及び仮想粒子は、金属、薬剤、トナーなど、粉体ならば特に制限無くシミュレーション可能である。
【００２８】
次に、粒子同士の接触が検出（判定）された場合には、更に粒子姿勢角度を測定すること（ステップＳ１３）により、粒子姿勢角度θ（θについては図４参照）が０の場合は、面接触（ステップＳ１３１）と判定され、θが観測される場合は点接触（ステップＳ１３２）と判定される。回転運動の計算にあたり、面接触の場合と、点接触の場合とで、回転モーメントの計算の仕方を変更する。
【００２９】
ここで言う粒子姿勢角度θとは、図４に示すように、粒子を平面に見立てた場合、正六角形ＡＢＣＤＥＦと、正六角形Ａ´Ｂ´Ｃ´Ｄ´Ｅ´Ｆ´において、辺ＡＢと辺Ａ´Ｂ´とが重なる場合（図４パターンＡ「想定」参照）をθ＝０、即ち面接触と定義され、辺Ａ´Ｂ´に点Ａが接触（図４パターンＢ「想定」参照）した場合を角度ＢＡＢ´＝θとして角度が測定され、点接触と定義される角度のことをいう。
【００３０】
ステップＳ１３、ステップＳ１３１及びステップＳ１３２で粒子姿勢角度を測定した後、粒子間の距離を計算し（ステップＳ１４）、接触による弾性反発力や粘性力などの作用力を計算する（ステップＳ１５）。さらに回転運動を計算し、粒子の角速度および粒子姿勢角度を計算する。そして、接触粒子における計算が終了するまで（ステップＳ１６）、これらの計算処理（ステップＳ１３〜Ｓ１５）を繰り返し行う。
【００３１】
なお、ステップＳ１２にて接触が判定されない場合は、後述するステップＳ１７に進む。
【００３２】
次いで、粒子に作用する外力の計算を行う（ステップＳ１７）。ここで言う外力には、重力、ファンデルワールス力などが挙げられる。
【００３３】
そして、すべての粒子について、ステップＳ１３〜Ｓ１７の処理を繰り返し実行する（ステップＳ１８）。
【００３４】
次いで、解析対象である各粒子についての運動方程式を立てて、それぞれの加速度、速度、並びに変位の計算を行う（ステップＳ１９）。なお、本発明で用いる運動方程式は下記の通りである（下記数１参照）。
【００３５】
【数１】
上記数１にて、並進運動は、ニュートンの第二法則より、式（Ａ）のように表される。ここで、式（Ａ）におけるＦＣ及びＦｇは、接触力及び重力である。
【００３６】
次に、接触力ＦＣについてはそれぞれ、法線方向成分ＦＣｎ及び接線方向成分ＦＣｔについて、上記数１の式（１）及び（２）のように表される。式（１）及び（２）において、ｋ、δ、η、ｖ、ｒ及びωは、それぞれ、バネ定数、変位、粘性減衰係数、速度、粒子半径及び角速度である。添字のｎ及びｔは、法線方向成分及び接線方向成分を意味する。
【００３７】
一方、上記数１にて、回転運動は、上記数１の式（３）のように表される。ここで、Ｉは慣性モーメントであり、ｒ及びＦＣｔは、上述のように式（１）及び式（２）で定義したものである。
【００３８】
ステップＳ１９の運動方程式を解くにあたって、粒子同士が面接触と判定された場合（ステップＳ１３１）は、面接触を考慮して回転運動を計算する。回転運動は、面接触を考慮して式（３）を用いて計算する。面接触において角加速度を計算するにあたり、接触している面とそれに作用する接触力を考慮する。
【００３９】
ステップＳ１９の運動方程式を解くにあたって、粒子同士が点接触と判定された場合（ステップＳ１３２）は、該点接触において角加速度を計算する場合、接触点および接触力を考慮する。さらに、回転運動の計算において、粒子姿勢角度を考慮する。これは、たとえば、２体粒子の点接触の場合、接触頂点が粒子の重心を結ぶ直線上にない場合、回転することを意味する。球形を仮定した場合これができない。
【００４０】
そして、時間ステップが終了するまで、ステップＳ１２〜Ｓ１９の処理を繰り返し実行する（ステップＳ２０）。
【００４１】
本実施形態によれば、解析対象粒子に現実の諸物理特性値を与えるとともに、仮想粒子にも諸物理特性値を与えているので、ステップＳ１９において運動方程式を計算して粒子間の機械的相互作用を求める際には、粒子同士の受ける効果を再現することができ、解析対象粒子の挙動を高精度に解析することができる。
【００４２】
本実施形態では、２次元体系の正六角形型粒子の挙動解析について記載したが、この形態に限られるものではなく、粒子姿勢角度さえ特定できれば、転がり抵抗の小さい粒子、即ち正十二角形などの多角形型粒子のシミュレーションも可能である。また、３次元体系モデルへの応用も可能である。
【００４３】
本発明に係る離散要素法解析シミュレーション方法を実行させるためのプログラム（以下「離散要素法解析シミュレーションプログラム」と記す）は、任意の装置外部の記録媒体（ＣＤ−ＲＯＭやフラッシュメモリなど）、又は装置内部の記録媒体（ハードディスクやメモリなど）に格納されており、当該装置により読み出されて実行される。ここで言う装置は、パーソナルコンピュータなど、前記離散要素法解析シミュレーションプログラムの実行制御が可能な情報処理装置であれば良く、ハードウェア構成は汎用的なものを適用できる。なお、本発明に係る離散要素法解析シミュレーションプログラム、前記プログラムを記録するための記録媒体及び装置は種類や数に限定されるものではない。また、前記プログラムの搭載（格納）については装置の形態によって適宜変更可能である。
【００４４】
また、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備するが、前記離散要素法解析シミュレーションプログラムを記録した記録媒体の搭載方法は適宜変更可能である。
【産業上の利用可能性】
【００４５】
本発明のシミュレーション方法を使用することにより、金属粒子の最密充填、トナー粒子、製薬分野などで取り扱われる粉体の挙動解析が期待される。
【符号の説明】
【００４６】
１ ボール
２ 法線方向バネ
３ 法線方向ダッシュポット
４ 剪断方向バネ
５ 剪断方向ダッシュポット
６ スライダ
【技術分野】
【０００１】
本発明は、粒子の転がり抵抗を導入した多角形型粒子の挙動の離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置に関する。
【背景技術】
【０００２】
粉体や粒体などの粒子を取り扱う分野では、粒子の挙動を把握することが重要な課題である。従来は、このような粒子の挙動を試行錯誤の実験により把握することが多かった。このような場合、現実に用いる粒子の組成や、実験装置や実験環境を用意しなければ、所望の粒子の挙動解析を行うことが出来ない。また、これらの仕様を変更するたびに、実験しなければならず、その都度時間やコストの面を考慮しなくてはならなかった。
【０００３】
近年では、粒子の挙動解析を実験的にではなく、数式化若しくはモデル化し、ほぼ同じ法則に支配されるシステムの挙動をコンピュータ上で模擬すること、すなわちシミュレーションが広く利用されている。シミュレーション技術を利用することにより、粒子の挙動を予測することができる。また、粒子（粉体）を投入する条件やパラメータを変えて、同じシミュレーション計算を繰り返し行うことにより、様々な粒子組成や装置設計・制御体型の性能を評価することができる。
【０００４】
数式化した場合の粉体の挙動解析シミュレーション手段として、粉体毎の挙動を解析する離散要素（若しくは個別要素）法（Ｄｉｓｃｒｅｔｅ Ｅｌｅｍｅｎｔｓ Ｍｅｔｈｏｄ、若しくはＤｉｓｔｉｎｃｔ Ｅｌｅｍｅｎｔｓ Ｍｅｔｈｏｄ）解析シミュレーション手法（以下本明細書ではＤＥＭシミュレーションと表記する。）と、解析対象の粒子の形状をモデル化して挙動を解析する手法などが挙げられる。
【０００５】
ＤＥＭシミュレーションの原理は、球状粒子に見立てた２つのボール間に作用する法線方向（図１（ａ）参照）、剪断方向（図１（ｂ）参照）のボール間接触力について、図１に示すような、弾性的性質を表すバネ２，４と、非弾性的性質を表すダッシュポット３，５と、ボール１間の摩擦的性質を表すスライダ６（剪断方向）とで構成されるフォークトモデルによって説明される。前記フォークトモデルにより、円筒型の容器の壁とボールの接触、ボール同士の接触を判定し、接触している場合のボールと壁、またはボール同士との法線、剪断方向の接触力と摩擦力を図１（ａ）、図１（ｂ）のフォークトモデルに基づいて演算し、バネ２，４の弾性係数はヘルツの弾性接触理論より求められ、ダッシュポット３，５の粘性減衰係数（法線方向）、摩擦係数（剪断方向）は実験により求められる。そして、その接触力に基づいて全ボールの加速度、速度、変位を各離散化時間で演算し、この演算を繰り返すことによって任意の時間のボール群全体の運動挙動がシミュレーションされる。
【０００６】
上記ＤＥＭシミュレーションによれば、粒子と粒子との間の衝突判定が容易且つ、計算負荷を小さくできる。ＤＥＭシミュレーションを使用した粒子の挙動解析の例として、特開２００６−２５９９１１号公報（特許文献１参照）に開示されているようなトナーに使用する粒子の挙動解析等が挙げられる。
【０００７】
しかしながら、粒子の形状は現実的には球形ではないことが多く、球形を仮定すると、接触角度が精度良く形成されないこと、粉体の流動が実際よりも大きく見積もられてしまうといったことなどが課題である。また、特許文献１に記載の方法は、あくまで粒子を球形で仮定していること、粒子間の点若しくは面接触を考慮していないため、例えば金属粒子の最密充填（高濃度充填）には適していない。
【０００８】
一方、従来のＤＥＭシミュレーションによる粒子挙動解析の改善策として、例えばマツシマら（非特許文献１参照）が提案した非球形粒子に複数の球形粒子をつなぎ合わせたクラスターモデルとＤＥＭシミュレーションをベースにした非球形粒子挙動解析方法がある。非特許文献１に記載の方法では、特許文献１のような従来のＤＥＭシミュレーションに比べると、モデル粉体が実際の粉体により近い流動を見積もることが可能になったが、クラスターモデルを使用するため計算体系が小さくなる、即ち一領域あたりのモデル粉体の数が少なく見積もられるといったことが問題である。
【０００９】
近年では、ＤＥＭシミュレーションと、多角形型粒子モデルとを組み合わせたジョン−ポール レイサムら（Ｊｏｈｎ−Ｐａｕｌ Ｌａｔｈａｍ， ｅｔ ａｌ）の提案モデル（非特許文献２参照）がある。非特許文献２に記載の方法は、従来のＤＥＭシミュレーションに比べ、粒子と壁面との接触（粒子姿勢）角度が精度良く形成されないといった問題は解消されたが、該シミュレーションで必要な運動方程式については球の回転運動を考慮しなくてはならないため、計算負荷がかかってしまうなどの問題点がある。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特開２００６−２５９９１１号公報
【非特許文献】
【００１１】
【非特許文献１】マツシマ他（Ｔ． Ｍａｔｓｕｓｈｉｍａ， ｅｔ ａｌ）， ジャーナル オブ エーロスペース エンジニアリング（Ｊｏｕｒｎａｌ ｏｆ Ａｅｒｏｓｐａｃｅ Ｅｎｇｉｎｅｅｒｉｎｇ）， ２００９年， １５−２３頁．
【非特許文献２】ジョン−ポール レイサム他（Ｊｏｈｎ−Ｐａｕｌ Ｌａｔｈａｍ， ｅｔ ａｌ）， ミネラルズ エンジニアリング（Ｍｉｎｅｒａｌｓ Ｅｎｇｉｎｅｅｒｉｎｇ）， ２１巻（２００８）， ７９７−８０５頁．
【発明の概要】
【発明が解決しようとする課題】
【００１２】
上記に述べた問題点に対して望まれる解決案は、計算上は球形粒子を使用して、粒子の転がり抵抗、即ち多角形型粒子の回転運動を導入することである。
【００１３】
以上の事情を鑑み、本発明は、転がり抵抗を導入した多角形型粒子モデル、即ち多角形型粒子の回転運動を導入した多角形型粒子モデルの離散要素法解析シミュレーション方法、離散要素法解析シミュレーションプログラム、及び離散要素法解析シミュレーション装置を提供する。
【課題を解決するための手段】
【００１４】
本発明の上記目的は、多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、仮想粒子に物性値を与えるステップと、前記仮想粒子を配置する仮想粒子配置ステップと、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、並進運動を基に立てた運動方程式を解いて個々の粒子の挙動を求めることを特徴とすることにより効果的に達成される。
【００１５】
更に、前記仮想粒子に与える物性値は、ヤング率、粘性減衰係数、摩擦係数のうち少なくとも１つであることにより、或いは前記解析対象粒子が多角形型粒子であることにより、或いは前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が０のときに、前記粒子同士が面接触と判定されることにより、或いは前記粒子同士が面接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触面を考慮して回転角度を求め、個々の粒子の挙動を求めることにより、或いは前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が測定されるときに、前記粒子同士が点接触と判定されることにより、或いは前記粒子同士が点接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触点及び接触力を考慮して回転角度を求め、個々の粒子の挙動を求めることにより、より効果的に達成される。
【００１６】
また、本発明は、前記離散要素法解析シミュレーション方法を実行するための離散要素法解析シミュレーションプログラム、或いは前記離散要素法解析シミュレーションプログラムを記録した記録媒体により、或いは前記記録媒体を搭載した離散要素法解析シミュレーション装置であって、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備することを特徴とすることにより、効果的に達成される。
【発明の効果】
【００１７】
本発明の離散要素法解析シミュレーション方法は、粒子姿勢角度を測定することで粒子間の面若しくは点接触を判定することが容易であり、精度良く粉体の流動を模擬することができる。
【００１８】
また、従来のＤＥＭシミュレーションと同程度の計算負荷で粉体の挙動をシミュレートすることが可能である。
【図面の簡単な説明】
【００１９】
【図１】従来技術におけるＤＥＭシミュレーションを力学的に模したフォークトモデルである。
【図２】本発明の粒子挙動解析を行う前の前処理を示すフローチャートである。
【図３】本発明における、ＤＥＭシミュレーション方法のフローチャートである。
【図４】本発明の実施形態における、多角形（正六角形）型仮想粒子の概略図である。
【発明を実施するための形態】
【００２０】
以下、本発明を実施するための形態について、図面を参照しながら説明する。
【００２１】
図２は、本発明における粒子挙動解析を行う前の前処理の流れを示すフローチャートである。
【００２２】
先ず、挙動解析を行う前に、本発明に必要なシミュレーション用粒子、即ち仮想粒子に物性値を与える（ステップＳ１）。具体的には、バネ定数、粘性減衰係数、摩擦係数を仮想粒子に与えるようにする。バネ定数については、線形バネと非線形バネがある。非線形バネの場合、ヤング率などを入力する。
【００２３】
次に、仮想粒子を生成するための処理を行う（ステップＳ２）。処理手順は、従来のＤＥＭシミュレーションと同様に、仮想粒子を配置する範囲などを入力し、仮想粒子平均径や仮想粒子平均間隔を設定する。
【００２４】
次いで、仮想粒子数を設定（ステップＳ３）し、各仮想粒子の直径の設定、各仮想粒子を配置する座標位置の設定、並びに各仮想粒子の直径及び座標位置を出力（ステップＳ４）し、前処理を終了する。
【００２５】
図３は、本発明における粒子挙動解析の処理手順をフローチャートの形式で示している。また、図４には本発明の実施形態における、２次元体系の多角形（正六角形）型仮想粒子の概略図を示す。なお、図４において、正六角形型仮想粒子を示したが、実際に本発明のシミュレーション方法を使用する際は、特に粒子の形状は限定されない。
【００２６】
まず、解析対象となる各粒子に初期条件を与えて、初期化する（ステップＳ１１）。ここで言う初期条件とは、前述の仮想粒子の処理で予め計算した物性値等のことである。次に、粒子同士の接触（衝突）を判定する（ステップＳ１２）。ここで、前記判定は、球形粒子を仮定して行う（図４のパターンＡ及びＢの「計算」参照）。
【００２７】
ここで言う粒子には、解析対象粒子と、図４で表す仮想粒子の双方を含む。なお、前記解析対象粒子及び仮想粒子は、金属、薬剤、トナーなど、粉体ならば特に制限無くシミュレーション可能である。
【００２８】
次に、粒子同士の接触が検出（判定）された場合には、更に粒子姿勢角度を測定すること（ステップＳ１３）により、粒子姿勢角度θ（θについては図４参照）が０の場合は、面接触（ステップＳ１３１）と判定され、θが観測される場合は点接触（ステップＳ１３２）と判定される。回転運動の計算にあたり、面接触の場合と、点接触の場合とで、回転モーメントの計算の仕方を変更する。
【００２９】
ここで言う粒子姿勢角度θとは、図４に示すように、粒子を平面に見立てた場合、正六角形ＡＢＣＤＥＦと、正六角形Ａ´Ｂ´Ｃ´Ｄ´Ｅ´Ｆ´において、辺ＡＢと辺Ａ´Ｂ´とが重なる場合（図４パターンＡ「想定」参照）をθ＝０、即ち面接触と定義され、辺Ａ´Ｂ´に点Ａが接触（図４パターンＢ「想定」参照）した場合を角度ＢＡＢ´＝θとして角度が測定され、点接触と定義される角度のことをいう。
【００３０】
ステップＳ１３、ステップＳ１３１及びステップＳ１３２で粒子姿勢角度を測定した後、粒子間の距離を計算し（ステップＳ１４）、接触による弾性反発力や粘性力などの作用力を計算する（ステップＳ１５）。さらに回転運動を計算し、粒子の角速度および粒子姿勢角度を計算する。そして、接触粒子における計算が終了するまで（ステップＳ１６）、これらの計算処理（ステップＳ１３〜Ｓ１５）を繰り返し行う。
【００３１】
なお、ステップＳ１２にて接触が判定されない場合は、後述するステップＳ１７に進む。
【００３２】
次いで、粒子に作用する外力の計算を行う（ステップＳ１７）。ここで言う外力には、重力、ファンデルワールス力などが挙げられる。
【００３３】
そして、すべての粒子について、ステップＳ１３〜Ｓ１７の処理を繰り返し実行する（ステップＳ１８）。
【００３４】
次いで、解析対象である各粒子についての運動方程式を立てて、それぞれの加速度、速度、並びに変位の計算を行う（ステップＳ１９）。なお、本発明で用いる運動方程式は下記の通りである（下記数１参照）。
【００３５】
【数１】
上記数１にて、並進運動は、ニュートンの第二法則より、式（Ａ）のように表される。ここで、式（Ａ）におけるＦＣ及びＦｇは、接触力及び重力である。
【００３６】
次に、接触力ＦＣについてはそれぞれ、法線方向成分ＦＣｎ及び接線方向成分ＦＣｔについて、上記数１の式（１）及び（２）のように表される。式（１）及び（２）において、ｋ、δ、η、ｖ、ｒ及びωは、それぞれ、バネ定数、変位、粘性減衰係数、速度、粒子半径及び角速度である。添字のｎ及びｔは、法線方向成分及び接線方向成分を意味する。
【００３７】
一方、上記数１にて、回転運動は、上記数１の式（３）のように表される。ここで、Ｉは慣性モーメントであり、ｒ及びＦＣｔは、上述のように式（１）及び式（２）で定義したものである。
【００３８】
ステップＳ１９の運動方程式を解くにあたって、粒子同士が面接触と判定された場合（ステップＳ１３１）は、面接触を考慮して回転運動を計算する。回転運動は、面接触を考慮して式（３）を用いて計算する。面接触において角加速度を計算するにあたり、接触している面とそれに作用する接触力を考慮する。
【００３９】
ステップＳ１９の運動方程式を解くにあたって、粒子同士が点接触と判定された場合（ステップＳ１３２）は、該点接触において角加速度を計算する場合、接触点および接触力を考慮する。さらに、回転運動の計算において、粒子姿勢角度を考慮する。これは、たとえば、２体粒子の点接触の場合、接触頂点が粒子の重心を結ぶ直線上にない場合、回転することを意味する。球形を仮定した場合これができない。
【００４０】
そして、時間ステップが終了するまで、ステップＳ１２〜Ｓ１９の処理を繰り返し実行する（ステップＳ２０）。
【００４１】
本実施形態によれば、解析対象粒子に現実の諸物理特性値を与えるとともに、仮想粒子にも諸物理特性値を与えているので、ステップＳ１９において運動方程式を計算して粒子間の機械的相互作用を求める際には、粒子同士の受ける効果を再現することができ、解析対象粒子の挙動を高精度に解析することができる。
【００４２】
本実施形態では、２次元体系の正六角形型粒子の挙動解析について記載したが、この形態に限られるものではなく、粒子姿勢角度さえ特定できれば、転がり抵抗の小さい粒子、即ち正十二角形などの多角形型粒子のシミュレーションも可能である。また、３次元体系モデルへの応用も可能である。
【００４３】
本発明に係る離散要素法解析シミュレーション方法を実行させるためのプログラム（以下「離散要素法解析シミュレーションプログラム」と記す）は、任意の装置外部の記録媒体（ＣＤ−ＲＯＭやフラッシュメモリなど）、又は装置内部の記録媒体（ハードディスクやメモリなど）に格納されており、当該装置により読み出されて実行される。ここで言う装置は、パーソナルコンピュータなど、前記離散要素法解析シミュレーションプログラムの実行制御が可能な情報処理装置であれば良く、ハードウェア構成は汎用的なものを適用できる。なお、本発明に係る離散要素法解析シミュレーションプログラム、前記プログラムを記録するための記録媒体及び装置は種類や数に限定されるものではない。また、前記プログラムの搭載（格納）については装置の形態によって適宜変更可能である。
【００４４】
また、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備するが、前記離散要素法解析シミュレーションプログラムを記録した記録媒体の搭載方法は適宜変更可能である。
【産業上の利用可能性】
【００４５】
本発明のシミュレーション方法を使用することにより、金属粒子の最密充填、トナー粒子、製薬分野などで取り扱われる粉体の挙動解析が期待される。
【符号の説明】
【００４６】
１ ボール
２ 法線方向バネ
３ 法線方向ダッシュポット
４ 剪断方向バネ
５ 剪断方向ダッシュポット
６ スライダ
【特許請求の範囲】
【請求項１】
多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、
仮想粒子に物性値を与えるステップと、
前記仮想粒子を配置する仮想粒子配置ステップと、
前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、
前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、回転運動を基に立てた運動方程式を解いて個々の粒子の挙動を求めることを特徴とする離散要素法解析シミュレーション方法。
【請求項２】
前記仮想粒子に与える物性値は、バネ定数、粘性減衰係数、摩擦係数のうち少なくとも１つである請求項１に記載の離散要素法解析シミュレーション方法。
【請求項３】
前記解析対象粒子が多角形型粒子である請求項１又は２に記載の離散要素法解析シミュレーション方法。
【請求項４】
前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が０のときに、前記粒子同士が面接触と判定される請求項１乃至３のいずれか１項に記載の離散要素法解析シミュレーション方法。
【請求項５】
前記粒子同士が面接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触面を考慮して回転角度を求め、個々の粒子の挙動を求める請求項４に記載の離散要素法解析シミュレーション方法。
【請求項６】
前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が測定されるときに、前記粒子同士が点接触と判定される請求項１乃至３のいずれか１項に記載の離散要素法解析シミュレーション方法。
【請求項７】
前記粒子同士が点接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触点及び接触力を考慮して回転角度を求め、個々の粒子の挙動を求める請求項６に記載の離散要素法解析シミュレーション方法。
【請求項８】
請求項１乃至７のいずれか１項に記載の離散要素法解析シミュレーション方法を実行するための離散要素法解析シミュレーションプログラム。
【請求項９】
請求項８に記載の離散要素法解析シミュレーションプログラムを記録した記録媒体。
【請求項１０】
請求項９に記載の記録媒体を搭載した離散要素法解析シミュレーション装置であって、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備することを特徴とする離散要素法解析シミュレーション装置。
【請求項１】
多数の解析対象粒子の挙動を解析する離散要素法解析シミュレーション方法であって、
仮想粒子に物性値を与えるステップと、
前記仮想粒子を配置する仮想粒子配置ステップと、
前記解析対象粒子の挙動をシミュレーションする挙動シミュレーションステップとを具備し、
前記挙動シミュレーションステップは、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップを具備し、前記解析対象粒子及び／又は仮想粒子について、回転運動を基に立てた運動方程式を解いて個々の粒子の挙動を求めることを特徴とする離散要素法解析シミュレーション方法。
【請求項２】
前記仮想粒子に与える物性値は、バネ定数、粘性減衰係数、摩擦係数のうち少なくとも１つである請求項１に記載の離散要素法解析シミュレーション方法。
【請求項３】
前記解析対象粒子が多角形型粒子である請求項１又は２に記載の離散要素法解析シミュレーション方法。
【請求項４】
前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が０のときに、前記粒子同士が面接触と判定される請求項１乃至３のいずれか１項に記載の離散要素法解析シミュレーション方法。
【請求項５】
前記粒子同士が面接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触面を考慮して回転角度を求め、個々の粒子の挙動を求める請求項４に記載の離散要素法解析シミュレーション方法。
【請求項６】
前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定するステップにおいて、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度が測定されるときに、前記粒子同士が点接触と判定される請求項１乃至３のいずれか１項に記載の離散要素法解析シミュレーション方法。
【請求項７】
前記粒子同士が点接触と判定されるときは、回転運動を基に運動方程式を立て、粒子姿勢角度の情報に基づき回転モーメント及び角加速度の計算に、粒子同士の接触点及び接触力を考慮して回転角度を求め、個々の粒子の挙動を求める請求項６に記載の離散要素法解析シミュレーション方法。
【請求項８】
請求項１乃至７のいずれか１項に記載の離散要素法解析シミュレーション方法を実行するための離散要素法解析シミュレーションプログラム。
【請求項９】
請求項８に記載の離散要素法解析シミュレーションプログラムを記録した記録媒体。
【請求項１０】
請求項９に記載の記録媒体を搭載した離散要素法解析シミュレーション装置であって、前記装置は、仮想粒子に物性値を与える手段と、前記仮想粒子を配置する仮想粒子配置手段と、前記解析対象粒子の挙動をシミュレーションする挙動シミュレーション手段とを具備し、前記挙動シミュレーション手段は、前記解析対象粒子及び／又は仮想粒子間の粒子姿勢角度を測定する手段を具備することを特徴とする離散要素法解析シミュレーション装置。
【図１】
【図２】
【図３】
【図４】
【図２】
【図３】
【図４】
【公開番号】特開２０１１−８１５３０（Ｐ２０１１−８１５３０Ａ）
【公開日】平成２３年４月２１日（２０１１．４．２１）
【国際特許分類】
【出願番号】特願２００９−２３２１５９（Ｐ２００９−２３２１５９）
【出願日】平成２１年１０月６日（２００９．１０．６）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【出願人】（５０５０１１５３４）プロメテック・ソフトウェア株式会社 (5)
【出願人】（５９１１６０５１２）金属技研株式会社 (14)
【公開日】平成２３年４月２１日（２０１１．４．２１）
【国際特許分類】
【出願日】平成２１年１０月６日（２００９．１０．６）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【出願人】（５０５０１１５３４）プロメテック・ソフトウェア株式会社 (5)
【出願人】（５９１１６０５１２）金属技研株式会社 (14)
[ Back to top ]