粒子挙動の解析方法
【課題】大規模な設備における多数の粒子の挙動を計算するに際して、精度を落とすことなく、計算を高速化する計算方法及び粒子挙動解析方法を提供する。
【解決手段】多粒子系の粒子挙動解析において、粒子間の衝突過程に関する計算を簡略化するため、衝突した二粒子間の接触面法線方向の、時刻ｔpにおける粘性力を与える、二粒子間の相対速度において、一方の粒子の速度を時刻ｔp+1の予測速度に置き換える。粒子１に働く法線方向粘性力を、時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
とし、粒子２に働く法線方向粘性力を、時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表す。
【解決手段】多粒子系の粒子挙動解析において、粒子間の衝突過程に関する計算を簡略化するため、衝突した二粒子間の接触面法線方向の、時刻ｔpにおける粘性力を与える、二粒子間の相対速度において、一方の粒子の速度を時刻ｔp+1の予測速度に置き換える。粒子１に働く法線方向粘性力を、時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
とし、粒子２に働く法線方向粘性力を、時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表す。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は，離散要素法を用いた粒子挙動の計算方法に関するもの及びその製鉄分野への応用に関するものである。
【背景技術】
【０００２】
焼結原料は，図７の焼結原料分級設備概要に示すように，シャトルコンベアからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下させることで下部のパレット上に積載される。パレット上の積載層の底部には大きな粒径の焼結原料が，上部には小さな粒径の焼結原料が分布するようにこれらの装置配置は設計される。しかしながら、分級された焼結材料は必ずしも前述のような粒度分布で配置されない。そのような場合、焼結原料を払い出すにつれ、原料の集合体であったものが不定形に崩れやすく、効率的な原料の払い出しに支障をきたしたり、焼結原料が拡散し、割れたりして発塵源となったりして、深刻な生産効率の低下や環境への影響が懸念される場合があった。
【０００３】
また、シャトルコンベア，サージホッパー，ドラムフィーダ，フルイ等の要素からなり、焼結原料を分布させることを目的とした焼結原料装入装置を設計するにあたり、経験や試行錯誤で各機械要素の配置を決める場合も多かった。そのような装置の最適設計を行うに当たり、計算機を用いたシミュレーションを用いて精度の良い計算機実験を行うことは有効である。このようなシミュレーションにおける粒子群の挙動解析には、非特許文献１に示される離散要素法がよく用いられる。離散要素法とは、粒子群を構成する個々の粒子を要素とし、要素の集合体に対して個々の要素が運動方程式を満足することを条件として、集合体の動力学的挙動を数値解析する手法である。刻み時間をΔｔとし、運動方程式を差分方程式として数値計算を行う。
【０００４】
離散要素法においては、多数の粒子の挙動を解析する。粒子は、外場がなければ、衝突と等速度運動を繰り返して運動する。重力や電場等の外場があればそのような運動にさらに外場の力による運動が加わる。そこでは、多くの粒子の衝突は、二つの粒子の衝突が基本である。正確な結果を得るためには、粒子が衝突してから離れるまでの衝突過程は、細かく計算して精度を上げる必要があり、運動の計算のための刻み時間Δｔを細かくとる必要がある。
【０００５】
離散要素法においては、衝突している二つの粒子の取り扱いの概要は、次のようである。以下、時刻指標ｐを導入する。ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。時刻指標ｐにおける時刻をｔpで表す。同様に、時刻指標ｐにおける特定の物理量（例えばＶ）をＶpと表す。
【０００６】
時刻ｔpにおいて衝突している２個の粒子を粒子１及び粒子２であるとする。各々の粒子は、衝突している場合、粒子はお互いにめり込んでいると考える。そのときのめり込み量をδpとする。その様子は図１に示されている。めり込みによって形成される接触面の法線方向の力とその面内（接線）方向の成分に分けて、各々の粒子が受ける力を求める。どちらの方向の力も、図２に示されるように、バネ定数ｋのバネと粘性係数ηのダッシュポットの並列結合による力で表している。例えば、粒子１が受ける接触面法線方向の力は、めり込み量に比例するバネによる反発力と、衝突している二つの粒子の相対速度に比例する粘性力の和で表す。時刻ｔpにおける粒子１及び粒子２の法線方向速度を各々Ｖ1p、Ｖ2pであるとすれば、粒子１が、時刻ｔpにおいて衝突によって受ける法線方向の力Ｆ1pは、次のように表される。
Ｆ1p＝−ｋ×δp＋η×（Ｖ2p−Ｖ1p）
上式右辺第一項がバネによる反発力、第二項がダッシュポットによる粘性力である。重力などの外場があればそのような項が加わることになる。同様に、粒子２の受ける力は、上式における添え字１を添え字２に変えることによって表される。
【０００７】
法線方向の運動は、ニュートンの運動方程式に上式の力を代入して求める。計算は、刻み時間Δｔで離散的に計算が進められる。従って、粒子１の法線方向の運動、即ち、速度と位置の時間変化は、次に示す離散的運動方程式を解くことで得られる。
ｍ1（Ｖ1p+1−Ｖ1p）／Δｔ＝Ｆ1p
ｘ1p+1＝ｘ1p＋Ｖ1pΔｔ
ここで、ｍ1は粒子１の質量、ｘ1pは時刻ｔpにおける粒子１の中心の位置の接触面法線方向成分である。粒子２についても同様に求め、このようにして、刻み時間Δｔごとに速度と位置が変化する。二つの粒子の中心間距離が各々の半径の和であるｒ1＋ｒ2より大きくなった時点で衝突過程は終了し、等速運動あるいは、外場中の自由運動に移行することになる。
【０００８】
接触面面内方向（接線方向）の力も、同様なモデルで表され、その運動も同様に求められる。このほか、粒子の回転運動も加味されて、総合的な粒子運動として解析される。
【０００９】
このようにして精度よく粒子群の挙動が解析される。粒子が複数個の粒子と衝突している場合も、上記のような二個の粒子の衝突に分けて扱い、その合力を求めて運動方程式を解けばよい。更に、前記弾性係数ｋ、及び粘性係数ηについては、精度よく効率的に計算できるような表示式が非特許文献２及び非特許文献３等において開示されている。
【００１０】
また、特許文献１においては、高炉への原料装入装置用の旋回シュートを離散要素法を用いて設計することを開示している。
【００１１】
しかしながら、流動する粉粒体の巨視的挙動を離散要素法を用いて明らかにしようとすると、多数の粒子について、上記のような衝突過程を追跡する必要があり、すべての運動の計算を終了するまでに膨大な時間がかかるという問題がある。計算時間を短くするために、刻み時間Δｔを大きくして粗い時間スケールにして計算ステップ数を少なくすることは、効果的な計算時間短縮の方策の一つとして挙げられる。しかしながら、その場合、粒子が自由運動している場合は大きな問題はないが、粒子同士が衝突する際の運動解析については、上に示したように、細かい時間刻みで衝突過程を追跡しないと、大きな誤差が生じてしまうことが分かった。即ち、刻み時間Δｔを大きくすると、その衝突過程の計算で非常に大きい誤差が出てしまい、計算精度が大きく低下してしまうのである。
【００１２】
計算時間を短くするための技術の開示は、いくつか行われている。特許文献２においては、離散要素法を用いてコンベアベルトにより搬送される石炭や鉱石などの運搬物の挙動を解析しているが、時間短縮を図るため、ベルト搬送面等搬送物の挙動に直接的に影響する部分だけを記述して解析している。また、特許文献３においては、計算対象とする空間領域を分割し、分割領域内の代表粒子について計算を行うことを開示している。更に、特許文献４においては、複写機に用いられるトナーの挙動解析の例が開示されているが、移動量の大きい粒子については高精度の計算の対象とし、移動量の小さい粒子については低精度の計算の対象とするとしている。ここで、低精度の計算方法として、粒子間の接触判定を省略したり、計算間隔を長くとることで計算量を減らすことを開示している。
【００１３】
即ち、従来の技術では、解析対象となる粒子の位置や移動量を判定し、精度の異なる計算を適用することが必要となるため、煩雑な判定が必要となったり、低精度の計算の対象となった粒子の誤差が大きくなってしまい正しく粒子の挙動を解析できなくなる懸念があった。
【００１４】
このような懸念を回避する方策として、精度を保ったまま、計算量を減らすことのできる計算方法が求められている。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】特開２００７−２６２４５３号公報
【特許文献２】特開２００１−１３９１１６号公報
【特許文献３】特開２００７−２８６５１４号公報
【特許文献４】特開２０１０−７９４９３号公報
【非特許文献】
【００１６】
【非特許文献１】木山英郎，藤村尚，土木学会論文報告集，第３３３号 １３７−１４６ （１９８３）
【非特許文献２】S.P.Timoshenko and J.N.Goodier, Theory of Elasticity, p409, McGraw-Hill ｐ４０９（１９７０）
【非特許文献３】粉体工学会編，粉体シミュレーション入門，産業図書（株）ｐ３３（１９９８）
【発明の概要】
【発明が解決しようとする課題】
【００１７】
本発明は、大規模な設備における多数の粒子の挙動を計算するに際して、精度を落とすことなく、計算を高速化する計算方法及び、そのような計算方法を用いることで短時間での粒子の挙動の解析が可能となる粒子挙動解析方法を提供することを目的とするものである。
【課題を解決するための手段】
【００１８】
本発明者らは、本発明が対象とする多粒子系の粒子挙動解析の場合、粒子間の接触あるいは衝突が非常に多く、このような衝突過程に関する計算を省略できれば計算負荷を低下できると考えた。精度を保ったまま、時間刻みを大きくして衝突過程を省略できる計算方法について調査研究を行った。その結果、衝突した二粒子間の接触面法線方向の、時刻ｔpにおける粘性力を与える、二粒子間の相対速度において、一方の粒子の速度を時刻ｔp+1の予測速度に置き換えることにより、目的が達成されることが分かった。
【００１９】
本発明の要旨は、以下の通りである。
（１）粒子の運動の計算における任意の二つの粒子が衝突してから離れるまでの運動を計算する方法であって、時刻ｔpにおいて衝突している二粒子の接触面法線方向（以下、「法線方向」という。）に働く力を、弾性係数がｋであり、時刻ｔpにおける粒子のめり込み量δpに比例する弾性反発力
−ｋ×δp
と、粘性係数がηの粘性抵抗力の和で表し、
一方の粒子１の、時刻ｔpより刻み時間Δｔだけ進んだ時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
を、粒子１に働く法線方向粘性力とし、さらに、他方の粒子２に働く法線方向粘性力を、同様に時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表すことを特徴とする、二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
ここで、時刻指標ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。また、右上添え字ｐは、時刻指標ｐにおける当該物理指標を意味する。
（２）粒子１の質量をｍ1及び、粒子２の質量をｍ2として、時刻ｔp+1における、粒子１の予測法線方向速度Ｅ1p+1及び粒子２の予測法線方向速度Ｅ2p+1をそれぞれ
【数１】
及び
【数２】
で表すことを特徴とする、前記（１）記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
（３）時刻ｔpにおいて、粒子１及び粒子２が衝突から受ける法線方向の力Ｆ1p及びＦ2pを
【数３】
及び
【数４】
で表すことを特徴とする、前記（１）又は（２）に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
（４）多粒子系粒子の挙動を解析するための粒子挙動解析方法であって、個々の粒子の運動を運動方程式を用いて離散要素法で計算し、任意の２個の粒子が衝突してから離れるまでの運動の計算について上記（１）乃至（３）のいずれかに記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることを特徴とする粒子挙動解析方法。
（５）シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を求めるにあたり，上記（４）に記載の粒子挙動解析方法を用いることを特徴とする、焼結原料の粒子配置分布解析方法。
【図面の簡単な説明】
【００２０】
【図１】粒子が衝突している図である。
【図２】衝突している粒子に働く法線方向の力を説明する図である。
【図３】粒子の衝突過程の図である。
【図４】粒子の衝突過程の図である。
【図５】粒子の衝突過程の図である。
【図６】複数粒子の運動を計算するフローチャートである。
【図７】焼結原料分級設備概要の図である。
【図８】粒子挙動計算を説明する図である。
【発明を実施するための形態】
【００２１】
［第一の実施形態］
（粒子の衝突の判定とめりこみ量）
対象とする粒子の形状は特に限定しないが、ここでは球である場合を例にとって説明する。半径ｒ1の球である粒子１と、半径ｒ2の粒子２の二つの粒子が接触しているかどうかの判定は、種々の方法があるが、例えば粒子１の中心と粒子２の中心の距離がｒ1＋ｒ2以下であれば接触していると判定してもよい。粒子が球でなければ、粒子の重心をその中心とし、平均半径をその半径としてもよい。
【００２２】
めり込み量は、前記ｒ1＋ｒ2と接触している二つの粒子の中心間距離の差によって求められる。計算に用いる場合は、そのようにして得られた差を更に、前記ｒ1＋ｒ2で除した値を用いる場合もある。
【００２３】
粒子１と粒子２相互間のめり込み量をδとし、二つの粒子の中心の、その中心を結ぶ線上で計った位置をｘ1、ｘ2とすれば、二つの粒子の中心距離は、｜ｘ1−ｘ2｜と表される。これらは、時間的に変化する。各々の量の時刻ｔpにおける値であることを右肩添え字ｐで示せば、時刻ｔpにおけるめり込み量δpは次のように表すことができる。めり込みの例は図１に示されている。
δp＝｛（ｒ1＋ｒ2）−｜ｘ1p−ｘ2p｜｝／（ｒ1＋ｒ2） ［数式１］
【００２４】
ある時点で二つの粒子の接触判定を行った場合、接触と判定されたものについて、二つの粒子の中心距離は各々の半径の和よりもかなり小さく、即ち大きいめりこみ量である場合が多い。このようなめり込みが大きい場合、反発力も大きく、時間刻みを細かくとらないと、大きな反発力によって粒子が急速に反発してしまい、計算が不安定になってしまう場合がある。
【００２５】
（接触面法線方向に働く力）
接触した二つの粒子が離れるまでの衝突過程は、その接触面に対する法線方向に働く力（以下、接触法線力という。）が大きな影響を与える。即ち、該接触面から該法線に沿って大きく後退すれば衝突過程は早く終了する。従って、該法線方向に働く力を調整することで衝突過程を省略できるのである。また、ここで時間の推移は、刻み時間Δｔずつ推移するものとする。即ち、時刻ｔpにおいて、時刻指標ｐ＝０からｐΔｔだけ時間が経過している。
【００２６】
非特許文献１においては、時刻ｔpにおいて接触している二つの粒子に、その接触面の法線方向に働く力としては、粒子１に着目した場合、（１）めり込み量に比例する反発力−ｋδp 、および、（２）二つの粒子の相対速度に比例する粘性力、η（Ｖ2p−Ｖ1p）の和として与えられる。ここで右肩の添え字ｐは時刻指標を表す。δpは時刻ｔpにおける粒子１及び２のめり込み量である。また、Ｖ1pおよびＶ2pは各々、時刻ｔpにおける粒子１および２の前記接触面法線方向速度である。また、ｋおよびηは各々比例係数である。これらは図２に示すように、バネ定数ｋのバネと粘性係数ηのダッシュポットが並列につながったモデルで表される。即ち、非特許文献１にあるように、先行技術においては、粒子１に働く前記接触法線力Ｆ1pは
Ｆ1p＝−ｋδp＋η（Ｖ2p−Ｖ1p） ［数式２］
と表される。また、粒子２に働く該法線力は、［数式２］の下添え字の１と２を入れ替えればよい。
【００２７】
時間刻みΔｔを大きくすると、時刻ｔpで接触判定されても、粒子はすぐに離れ、衝突過程を省略し、計算が高速化できるが、上記［数式２］の反発力の項が大きくなり、後述するように粒子の速度が大きくなりすぎたりして、計算が不安定となったり、誤差が大きくなったりしてしまう。従って、先行技術においては、Δｔを細かくとって衝突過程を細かい時間間隔で追跡して計算する必要があった。本発明者等は、上記（２）の粘性力の表式に着目した。そこで、粒子１の時刻ｔpにおける該法線方向速度Ｖ1pを、粒子１の時刻ｔp+1における速度、即ちＶ1p+1で置き換えたのである。ここで、時刻ｔp+1とは、時刻ｔpから刻み時間Δｔだけ進んだ時刻である。しかしながら、時刻ｔpの時点で、時刻ｔp+1における正確な速度を把握することは困難なので、Ｖ1p+1は時刻ｔp+1における予測値Ｅ1p+1で表す。即ち、本発明においては、粒子１に働く接触法線力Ｆ1pは、前記［数式２］のＶ1pをＥ1p+1で置き換え、
Ｆ1p＝−ｋδp＋η（Ｖ2p−Ｅ1p+1） ［数式３］
と表される。粒子２に働く接触法線力Ｆ2pは、［数式３］の下添え字の１と２を入れ替えればよい。即ち
Ｆ2p＝−ｋδp＋η（Ｖ1p−Ｅ2p+1） ［数式４］
と表す。
【００２８】
予測法線方向速度の算出方法はいくつかある。例えば、本発明になる算出方法によれば、時間刻みΔｔを粗くしても、時刻ｔpにおいて衝突していても、時刻ｔp+1においては粒子は離れている場合が多く、かつ反発力も大きすぎず、精度を保ったまま計算量が減少し、［数式２］を用いた先行技術の場合と比較して計算時間が１／４以下になることが分かった。ここで、注意すべきなのは、粒子２に働く法線方向粘性力は、通常なら粒子１に働く粘性力の反力となるために同じ大きさで反対方向の力となり、それぞれ、作用、反作用の関係にあり、その和はゼロとなることである。本発明においては、［数式４］となり、これは、［数式３］と比較して必ずしも同じ大きさとはならない。従って、接触している２粒子において、これら、法線方向の粘性力の和はゼロではなく、有限な値を持った残留力となり、作用、反作用の関係が満たされず、厳密にいえば物理的に正しい描像を与えない。しかし、そのために生じる誤差は、無視できるほど小さいことが分かった。更に、多数の粒子の衝突を考えた場合、このような残留力は、ランダムな方向と大きさとなり、結果的にその総和はさらに小さくなり、精度に対する影響はほとんどないことが分かったのである。
【００２９】
（予測法線方向速度）
粒子１の予測法線速度Ｅ1p+1を、粒子１が衝突している粒子２のみから受ける力で運動した場合の時刻ｔp+1における粒子１の法線方向速度と考える。従って、次の離散的な運動方程式から求められる。
ｍ1（Ｅ1p+1−Ｖ1p）／Δｔ＝Ｆ1p ［数式５］
【００３０】
上式のＦ1pに［数式３］を代入し、整理すれば、
【数５】
が得られる。
【００３１】
［数式６］を［数式３］に再度代入すれば、時刻ｔpにおいて粒子１が粒子２より受ける法線方向力は、
【数６】
となる。同様に、粒子２について予測法線速度Ｅ2p+1 及び法線方向力Ｆ2pは以下の［数式８］、［数式９］のように表される。
【数７】
【数８】
【００３２】
以上より、予測法線方向速度Ｅ1p+1及びＥ2p+1が時刻ｔpの量で示され、法線方向力も時刻ｔpにおける値で示されたため、効率良いプログラム作成が可能となり、更に、前述のように精度を保ったまま、計算時間が短縮できるのである。
【００３３】
（接触判定された粒子の接触後の衝突運動の計算の省略）
時刻ｔpで二つの粒子が接触しているとした場合のその後の衝突過程の進展の概念図を図３〜図５に示す。従来技術である［数式２］で法線方向力を表した場合、時間刻みΔｔを細かくとり、図３に示すように、粒子が離れるまでの衝突過程を細かく計算する必要がある。例えば、同じ法線方向力のままΔｔを数倍に取った場合の衝突過程は、図４に示すように、接触後の粒子は次のステップで大きくめり込み、その次のステップでは反発力のために各々高速度で反発してしまい、計算が不安定になり、誤差が大きくなってしまう。次に、大きい時間刻みΔｔのままで［数式３］の本発明を用いると、図５に示すように、接触してから次のステップである程度めり込むが、その後、穏やかに離れ、精度を保ったまま衝突過程の計算を省略できる。これは、［数式３］を用いた場合、めり込み量が大きくなり、反発力がかなり大きくなっても、逆方向に働く粘性力も大きくなり、ほどよくバランスされるからである。即ち、［数式３］を用いれば、精度を保ったまま時間刻みを大きくとって衝突過程を省略し、計算速度を短くできるのである。
【００３４】
以上の計算ステップの大幅な省略のため、本発明になる計算方法を用いる場合、時間刻みΔｔは、従来技術の場合の時間刻みの４倍〜１０倍とっても、計算結果の精度はほとんど変わらないことが分かった。即ち、従来技術による計算の精度と同一の精度のまま計算時間が１／４〜１／１０になることが分かった。
【００３５】
（接触粒子の運動）
多数の粒子からなる粒子群においては、一つの粒子が同時に複数の粒子と衝突している場合が多い。その場合も、各々の二粒子の衝突について、前記の法線方向力をもとめその合力が粒子に働くと考える。
【００３６】
例えば、ある着目粒子ｉを固定し、他の対向粒子ｊについて粒子ｉと衝突しているかどうか判定する。衝突していると判定されたものについては、各々の接触面について法線方向力を［数式７］及び［数式９］において表される二粒子間の法線方向力を求める。これを繰り返し行い、粒子ｉについてすべての衝突粒子を求めたら、別の粒子に着目して同様な探索を行い、全ての粒子について衝突粒子と、各衝突における法線方向力を求め、これらの和を求める。但し、各々の衝突において接触面法線は異なる方向を向いているからベクトル的な和を求める必要がある。例えば着目粒子ｉにＮ個の対向粒子ｊが衝突しているとすれば、各々のｉとｊの衝突に対して異なる方向の法線速度を考え、前記［数式３］〜［数式９］を適用すればよい。即ち、着目粒子ｉの、各対向粒子ｊとの各接触面法線方向の粒子ｉの時刻ｔp+1での予測速度Ｅjip+1は、［数式６］と同様に考えればよく、
【数９】
となる。ここで、δijpは時刻ｔpにおける粒子ｉ，ｊの衝突によるめり込み量、Ｖijpは時刻ｔpにおける粒子ｊの、粒子ｉ、ｊの接触面法線方向の速度成分である。Ｖjipは、同様に時刻ｔpにおける粒子ｉの前記接触面法線方向の速度成分である。
【００３７】
従って、時刻ｔpにおいて粒子ｉと粒子ｊとの衝突から粒子ｉが受ける法線方向の力Ｆjipは、［数式７］と同様に考えて次のように表される。
【数１０】
【００３８】
この力は、前記のように、時刻ｔpにおける粒子ｉと粒子ｊの接触面に対する法線方向の量であり、そのような法線方向の単位ベクトルをｅjipであらわせば、上記［数式１１］の力の方向を示すベクトルとなる。
【００３９】
時刻ｔpにおいて粒子ｉがＮ個の衝突から受ける各々の法線方向力の合力は、［数式１１］で表される力を方向ｅjipを持つベクトルとして対向粒子ｊについて和を取ることで得る。
【数１１】
【００４０】
更に、各々の衝突における接触面の接線方向の力、及び各々の衝突による、回転モーメントも求め、その合力を求める。ここで、各接触面の接線方向の力及び回転モーメントは、非特許文献１にあるように、先行技術を用いて算出することができる。即ち、接触面の接線方向の力については、せん断抗力を与える弾性スプリングと粘性ダッシュポットの並列配置を仮定し、接触点近傍のせん断変形が主として要素間の摩擦力によって生ずるとして計算を行う。回転モーメントについては、要素毎に回転モーメントの運動方程式を立て、変位の運動方程式と同様に差分法による数値解析で求めることができる。これらの合力をもとめたら、先行技術と同様にして離散的な運動方程式を解き、時刻ｔp+1における速度と位置を求める。以上のステップのフローチャートを図６に示す。
【００４１】
（粒子挙動解析方法）
多粒子系粒子の挙動を解析するための粒子挙動解析方法において、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることができる。粒子挙動解析方法においては、個々の粒子の運動を運動方程式を用いて離散要素法で計算する。所定の空間領域内に存在するすべての粒子について、刻み時間をΔｔとし、各粒子の運動方程式を差分方程式として数値計算を行う。他の粒子と衝突していない粒子については、外場がなければ等速直線運動、重力や電場等の外場があればそのような運動にさらに外場の力による運動が加わる。ここにおいて、任意の２個の粒子が衝突してから離れるまでの運動の計算について、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いる。
【００４２】
系の初期条件と境界条件を定めた上で、以上の手順で、多粒子系粒子を構成するすべての粒子について挙動解析を行い、時間の経過と共に対象とする多粒子系粒子が全体としてどのように運動するのかを明らかにすることができる。
【００４３】
本発明になる計算方法、装置を用いて、大きさの異なる粒子によって構成される粒子群が、重力加速度のもとで種々の装置の中を通過し、最終的に運動を停止するまでを追跡し、最終的に粒子群が形成する塊における粒子の空間粒度分布をもとめることが可能である。この時、粒子群が通過する装置として、シャトルコンベヤや、ホッパー、ドラムフィーダ、フルイ等の機械要素を仮想的に配置し、粒子群を焼結原料の粒を模した粒径分布を持つように決めてやることによって実際の焼結原料装入装置をモデル化することができる。最終的な粒子の空間粒径分布としては、大きい粒子が下になり、小さい粒子が上にくるような分布が形成されることが作業効率上望ましい。しかしながら、前述のように、機械要素の空間的な位置関係、操業条件によってそのような望ましい分布から外れた分布ができてしまう場合が多い。本発明になる計算方法、解析方法を用いて、焼結原料装置をモデル化し、最終的に形成される空間粒径分布を求めることができる。さらに、機械要素の位置、操業条件を変えて計算を実施することによって理想的な空間粒径分布を求めることができる。即ち、焼結原料装入装置の設計、操業条件の決定を著しく高速で行うことができるのである。
（焼結原料の粒子配置分布解析方法）
即ち、上記本発明の粒子挙動解析方法を、焼結原料の粒子配置分布解析方法として用いることができる。図７に示すような焼結原料分級設備において、シャトルコンベヤ１からサージホッパー２に投入され，ドラムフィーダ３から切り出され，フルイ４の上を通って落下してパレット５上の堆積物６となるた焼結原料の粒子配置分布を、本発明の粒子挙動解析方法によって求める。まず、シャトルコンベヤ１、サージホッパー２、ドラムフィーダ３、フルイ４を空間に配置し、それぞれの装置に各粒子が接触した際の境界条件を定める。次いで、シャトルコンベア内に搭載された焼結原料の粒子配置分布を初期条件として与える。以上の準備を整えた上で、系内に存在するすべての粒子（要素）について、運動方程式を刻み時間Δｔの差分方程式として数値計算を行う。その際、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いる。
【００４４】
このような計算を行った結果として、シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を、時間経過と共に明らかにすることができる。
【００４５】
（接触面法線方向の力にかかる係数）
［数式３］に用いられる係数ｋとηはそれぞれ，弾性係数および粘性係数と考えることができるが，単純な定数とするよりも、非特許文献２及び非特許文献３等に開示されたように、次のように設定することで精度のよい計算を効率良く行うことができる。
【００４６】
衝突している２粒子を粒子ｉ及び粒子ｊとする。
【００４７】
【数１２】
【数１３】
ここで，各変数は以下のように表わされる。
１／Ｅ*＝（１−νi2）／Ｅi＋（１−νj2）／Ｅj ［数式１５］
１／ｒ*＝１／ｒi＋１／ｒj ［数式１６］
１／ｍ*＝１／ｍi＋１／ｍj ［数式１７］
【数１４】
上記［数式１３］〜［数式１８］の範囲内で，Ｅはヤング率、νはポアソン比，ｒは衝突粒子の半径，ｍは衝突粒子の質量、ｅは反発係数，πは円周率を表わす。
【００４８】
（画像処理用演算装置ＧＰＵを用いた解析）
粒子挙動解析の計算に用いる計算機として、画像処理を高速処理するための演算装置であるところのグラフィック・プロセシング・ユニット（ＧＰＵ）で計算を行うと、通常の中央演算装置（ＣＰＵ）を用いる汎用計算機を用いた同じ計算と比較して、計算時間が１／５０以下に短縮できることが分かった。即ち、本発明になる計算方法を用いて作成されたプログラムをＧＰＵに適用すると、従来技術の計算方法を用いたプログラムをＣＰＵに適用した場合に比較して計算時間を１／２００以下に短縮することができるのである。ただし、ＧＰＵで実施する場合は、計算機が実施するプログラム形態がＣＰＵに適用するプログラム形態とは異なるので、プログラム形態を変更しなければならない。しかしながら、計算のアルゴリズムは同一であり、本発明になる計算方法になるアルゴリズムを用いることは同じであるので、ＧＰＵを用いて計算しても、本発明になるものであることは明白である。また、本発明になる計算方法には、ＧＰＵによる計算の特徴である並列計算を阻害する要因は全くない。従ってＧＰＵによる計算の特徴を効率良く引き出すことが可能である。
【実施例】
【００４９】
［実施例１］（汎用計算機を用いた解析）
〔数式１３〕〜［数式１８］に示された係数を用いて、Δｔ＝１０-6秒の時間刻みで、５種類の大きさの球状粒子、２×１０4個からなる粒子群が、ホッパーからパレットに落下する計算を、汎用計算機を用いてシミュレーションした。従来技術で実施したものは、約５０時間でシミュレーションが終了した。同様なシミュレーションを、時間刻みを５０倍の５×１０-5秒にして実施したところ、従来技術では、正確な解が得られなかった。しかしながら、時間刻みが５×１０-5のままで本発明になる計算方法を用いると、約５時間でシミュレーションが終了し、従来技術で１０-6の時間刻みで行った計算とほぼ同様な粒子分布を得た。即ち、前記「予測速度」に起因する誤差は、ほとんど無視できることが分かった。さらに、本発明になるプログラムにより、精度を保ったまま、計算時間が１／１０に短縮されたことが分かった。
【００５０】
［実施例２］
（汎用計算機を用いた解析と画像処理用演算装置ＧＰＵを用いた解析との比較）
２×１０5個の粒子からなる粒子群の挙動の計算を、汎用計算機を用いて１×１０-6秒の時間刻みで従来の計算方法を用いてシミュレーションしたところ、計算終了まで約９８日かかった。同様な計算を、時間刻みを５×１０-5秒として、本発明の計算方法を用いて実施したところ、約８日かかった。時間刻み５×１０-5秒で、本発明の計算方法を適用したプログラムをＧＰＵを用いた装置で動作するように変換してＧＰＵを持つ装置で計算した。その結果１２時間で計算が終了し、従来法で１×１０-6秒の時間刻みの場合とほぼ同じ結果を得た。従って、汎用計算機を用いた従来技術に対して精度を維持したまま、計算時間が約１／２００に短縮された。
【００５１】
［実施例３］
図７に示す装置を用い、シャトルコンベヤ１からサージホッパー２に投入され，ドラムフィーダ３から切り出され，フルイ４の上を通って落下した焼結原料の分布をシミュレーションするにあたり，本発明による粒子運動の計算方法を、汎用計算機を用いてシミュレーションを行った。まず、シャトルコンベヤ、サージホッパー、ドラムフィーダ、フルイを空間に配置し、それぞれの装置に各粒子が接触した際の境界条件を定めた。次いで、シャトルコンベア内に搭載された焼結原料の粒子配置分布を初期条件として与えた。以上の準備を整えた上で、系内に存在するすべての粒子（要素）について、運動方程式を刻み時間Δｔの差分方程式として数値計算を行った。その際、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いた。
【００５２】
図８に，フルイ部がプレートと複数の条材からなる装置における粒子挙動の計算結果を示す。落下した粒子は下部のパレット上に積もる。下部のパレットは図中右から左に所定の速度で移動しており，堆積した粒子は図中左方向に運ばれて次の工程に進む。計算では，３種類の大きさの粒子を考慮した（大粒子（最も明るい、白 粒子径６ｍｍ），中粒子（最も暗い、黒 粒子径５ｍｍ），小粒子（中間の明るさ、灰色 粒子径３ｍｍ））。各々の粒子は、同じ重量密度、３．３ｇ／ｃｍ3とした。条材フルイの間隔は先端（右下）ほど広くなっており，大粒子ほど図右よりに落下する。落下した大粒子は堆積物斜面上を右方向に転がり，斜面右端から堆積物下部に巻き込まれ易くなる。一方小粒子は，フルイ部付け根あたりから落下し，下部バレットの堆積物の上面に降り積もる。これらにより，堆積物は，底から大粒子（最も明るい、白），中粒子（最も暗い、黒 ），小粒子（中間の明るさ、灰色）と分級されることが分かった。この結果を得るのに要した計算時間は１０時間であった。さらに、同様な装置構成で同様な粒子挙動を、従来の離散要素法を用いて汎用計算機を用いて計算したところ約１００時間かかった。本発明の手法により効率が向上したことが分かる。
【００５３】
［実施例４］
ＧＰＵを用いた画像処理専用計算機に、実施例３で用いた、本発明になるプログラムを移植し、実施例３と同じ要素配置での計算を行った。この時実施例３と同様な結果を得たが、計算に要した時間は２時間となり、汎用計算機を用いた実施例３に比較して１／５０以下となった。更に、画像処理専用の計算機を用いたために、途中結果や最終結果のグラフィック表示も効率良く行うことができ、設計効率が著しく向上した。
【００５４】
［実施例５］
ＧＰＵを用いた画像処理専用計算機を用いて、フルイ部がプレートと複数の条材からなる装置において，プレートの長さと角度，条材の開き角度を変えて計算を行った。その結果，堆積した粒子の分布は，プレート長さが０．５〜１．０ ｍ，角度が４０〜５０度，条材の開き角度が１．５〜４．０度の範囲が良いことが分かった。
【符号の説明】
【００５５】
１ シャトルコンベア
２ サージホッパー
３ ドラムフィーダ
４ フルイ
５ パレット
６ 堆積物
【技術分野】
【０００１】
本発明は，離散要素法を用いた粒子挙動の計算方法に関するもの及びその製鉄分野への応用に関するものである。
【背景技術】
【０００２】
焼結原料は，図７の焼結原料分級設備概要に示すように，シャトルコンベアからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下させることで下部のパレット上に積載される。パレット上の積載層の底部には大きな粒径の焼結原料が，上部には小さな粒径の焼結原料が分布するようにこれらの装置配置は設計される。しかしながら、分級された焼結材料は必ずしも前述のような粒度分布で配置されない。そのような場合、焼結原料を払い出すにつれ、原料の集合体であったものが不定形に崩れやすく、効率的な原料の払い出しに支障をきたしたり、焼結原料が拡散し、割れたりして発塵源となったりして、深刻な生産効率の低下や環境への影響が懸念される場合があった。
【０００３】
また、シャトルコンベア，サージホッパー，ドラムフィーダ，フルイ等の要素からなり、焼結原料を分布させることを目的とした焼結原料装入装置を設計するにあたり、経験や試行錯誤で各機械要素の配置を決める場合も多かった。そのような装置の最適設計を行うに当たり、計算機を用いたシミュレーションを用いて精度の良い計算機実験を行うことは有効である。このようなシミュレーションにおける粒子群の挙動解析には、非特許文献１に示される離散要素法がよく用いられる。離散要素法とは、粒子群を構成する個々の粒子を要素とし、要素の集合体に対して個々の要素が運動方程式を満足することを条件として、集合体の動力学的挙動を数値解析する手法である。刻み時間をΔｔとし、運動方程式を差分方程式として数値計算を行う。
【０００４】
離散要素法においては、多数の粒子の挙動を解析する。粒子は、外場がなければ、衝突と等速度運動を繰り返して運動する。重力や電場等の外場があればそのような運動にさらに外場の力による運動が加わる。そこでは、多くの粒子の衝突は、二つの粒子の衝突が基本である。正確な結果を得るためには、粒子が衝突してから離れるまでの衝突過程は、細かく計算して精度を上げる必要があり、運動の計算のための刻み時間Δｔを細かくとる必要がある。
【０００５】
離散要素法においては、衝突している二つの粒子の取り扱いの概要は、次のようである。以下、時刻指標ｐを導入する。ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。時刻指標ｐにおける時刻をｔpで表す。同様に、時刻指標ｐにおける特定の物理量（例えばＶ）をＶpと表す。
【０００６】
時刻ｔpにおいて衝突している２個の粒子を粒子１及び粒子２であるとする。各々の粒子は、衝突している場合、粒子はお互いにめり込んでいると考える。そのときのめり込み量をδpとする。その様子は図１に示されている。めり込みによって形成される接触面の法線方向の力とその面内（接線）方向の成分に分けて、各々の粒子が受ける力を求める。どちらの方向の力も、図２に示されるように、バネ定数ｋのバネと粘性係数ηのダッシュポットの並列結合による力で表している。例えば、粒子１が受ける接触面法線方向の力は、めり込み量に比例するバネによる反発力と、衝突している二つの粒子の相対速度に比例する粘性力の和で表す。時刻ｔpにおける粒子１及び粒子２の法線方向速度を各々Ｖ1p、Ｖ2pであるとすれば、粒子１が、時刻ｔpにおいて衝突によって受ける法線方向の力Ｆ1pは、次のように表される。
Ｆ1p＝−ｋ×δp＋η×（Ｖ2p−Ｖ1p）
上式右辺第一項がバネによる反発力、第二項がダッシュポットによる粘性力である。重力などの外場があればそのような項が加わることになる。同様に、粒子２の受ける力は、上式における添え字１を添え字２に変えることによって表される。
【０００７】
法線方向の運動は、ニュートンの運動方程式に上式の力を代入して求める。計算は、刻み時間Δｔで離散的に計算が進められる。従って、粒子１の法線方向の運動、即ち、速度と位置の時間変化は、次に示す離散的運動方程式を解くことで得られる。
ｍ1（Ｖ1p+1−Ｖ1p）／Δｔ＝Ｆ1p
ｘ1p+1＝ｘ1p＋Ｖ1pΔｔ
ここで、ｍ1は粒子１の質量、ｘ1pは時刻ｔpにおける粒子１の中心の位置の接触面法線方向成分である。粒子２についても同様に求め、このようにして、刻み時間Δｔごとに速度と位置が変化する。二つの粒子の中心間距離が各々の半径の和であるｒ1＋ｒ2より大きくなった時点で衝突過程は終了し、等速運動あるいは、外場中の自由運動に移行することになる。
【０００８】
接触面面内方向（接線方向）の力も、同様なモデルで表され、その運動も同様に求められる。このほか、粒子の回転運動も加味されて、総合的な粒子運動として解析される。
【０００９】
このようにして精度よく粒子群の挙動が解析される。粒子が複数個の粒子と衝突している場合も、上記のような二個の粒子の衝突に分けて扱い、その合力を求めて運動方程式を解けばよい。更に、前記弾性係数ｋ、及び粘性係数ηについては、精度よく効率的に計算できるような表示式が非特許文献２及び非特許文献３等において開示されている。
【００１０】
また、特許文献１においては、高炉への原料装入装置用の旋回シュートを離散要素法を用いて設計することを開示している。
【００１１】
しかしながら、流動する粉粒体の巨視的挙動を離散要素法を用いて明らかにしようとすると、多数の粒子について、上記のような衝突過程を追跡する必要があり、すべての運動の計算を終了するまでに膨大な時間がかかるという問題がある。計算時間を短くするために、刻み時間Δｔを大きくして粗い時間スケールにして計算ステップ数を少なくすることは、効果的な計算時間短縮の方策の一つとして挙げられる。しかしながら、その場合、粒子が自由運動している場合は大きな問題はないが、粒子同士が衝突する際の運動解析については、上に示したように、細かい時間刻みで衝突過程を追跡しないと、大きな誤差が生じてしまうことが分かった。即ち、刻み時間Δｔを大きくすると、その衝突過程の計算で非常に大きい誤差が出てしまい、計算精度が大きく低下してしまうのである。
【００１２】
計算時間を短くするための技術の開示は、いくつか行われている。特許文献２においては、離散要素法を用いてコンベアベルトにより搬送される石炭や鉱石などの運搬物の挙動を解析しているが、時間短縮を図るため、ベルト搬送面等搬送物の挙動に直接的に影響する部分だけを記述して解析している。また、特許文献３においては、計算対象とする空間領域を分割し、分割領域内の代表粒子について計算を行うことを開示している。更に、特許文献４においては、複写機に用いられるトナーの挙動解析の例が開示されているが、移動量の大きい粒子については高精度の計算の対象とし、移動量の小さい粒子については低精度の計算の対象とするとしている。ここで、低精度の計算方法として、粒子間の接触判定を省略したり、計算間隔を長くとることで計算量を減らすことを開示している。
【００１３】
即ち、従来の技術では、解析対象となる粒子の位置や移動量を判定し、精度の異なる計算を適用することが必要となるため、煩雑な判定が必要となったり、低精度の計算の対象となった粒子の誤差が大きくなってしまい正しく粒子の挙動を解析できなくなる懸念があった。
【００１４】
このような懸念を回避する方策として、精度を保ったまま、計算量を減らすことのできる計算方法が求められている。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】特開２００７−２６２４５３号公報
【特許文献２】特開２００１−１３９１１６号公報
【特許文献３】特開２００７−２８６５１４号公報
【特許文献４】特開２０１０−７９４９３号公報
【非特許文献】
【００１６】
【非特許文献１】木山英郎，藤村尚，土木学会論文報告集，第３３３号 １３７−１４６ （１９８３）
【非特許文献２】S.P.Timoshenko and J.N.Goodier, Theory of Elasticity, p409, McGraw-Hill ｐ４０９（１９７０）
【非特許文献３】粉体工学会編，粉体シミュレーション入門，産業図書（株）ｐ３３（１９９８）
【発明の概要】
【発明が解決しようとする課題】
【００１７】
本発明は、大規模な設備における多数の粒子の挙動を計算するに際して、精度を落とすことなく、計算を高速化する計算方法及び、そのような計算方法を用いることで短時間での粒子の挙動の解析が可能となる粒子挙動解析方法を提供することを目的とするものである。
【課題を解決するための手段】
【００１８】
本発明者らは、本発明が対象とする多粒子系の粒子挙動解析の場合、粒子間の接触あるいは衝突が非常に多く、このような衝突過程に関する計算を省略できれば計算負荷を低下できると考えた。精度を保ったまま、時間刻みを大きくして衝突過程を省略できる計算方法について調査研究を行った。その結果、衝突した二粒子間の接触面法線方向の、時刻ｔpにおける粘性力を与える、二粒子間の相対速度において、一方の粒子の速度を時刻ｔp+1の予測速度に置き換えることにより、目的が達成されることが分かった。
【００１９】
本発明の要旨は、以下の通りである。
（１）粒子の運動の計算における任意の二つの粒子が衝突してから離れるまでの運動を計算する方法であって、時刻ｔpにおいて衝突している二粒子の接触面法線方向（以下、「法線方向」という。）に働く力を、弾性係数がｋであり、時刻ｔpにおける粒子のめり込み量δpに比例する弾性反発力
−ｋ×δp
と、粘性係数がηの粘性抵抗力の和で表し、
一方の粒子１の、時刻ｔpより刻み時間Δｔだけ進んだ時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
を、粒子１に働く法線方向粘性力とし、さらに、他方の粒子２に働く法線方向粘性力を、同様に時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表すことを特徴とする、二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
ここで、時刻指標ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。また、右上添え字ｐは、時刻指標ｐにおける当該物理指標を意味する。
（２）粒子１の質量をｍ1及び、粒子２の質量をｍ2として、時刻ｔp+1における、粒子１の予測法線方向速度Ｅ1p+1及び粒子２の予測法線方向速度Ｅ2p+1をそれぞれ
【数１】
及び
【数２】
で表すことを特徴とする、前記（１）記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
（３）時刻ｔpにおいて、粒子１及び粒子２が衝突から受ける法線方向の力Ｆ1p及びＦ2pを
【数３】
及び
【数４】
で表すことを特徴とする、前記（１）又は（２）に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
（４）多粒子系粒子の挙動を解析するための粒子挙動解析方法であって、個々の粒子の運動を運動方程式を用いて離散要素法で計算し、任意の２個の粒子が衝突してから離れるまでの運動の計算について上記（１）乃至（３）のいずれかに記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることを特徴とする粒子挙動解析方法。
（５）シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を求めるにあたり，上記（４）に記載の粒子挙動解析方法を用いることを特徴とする、焼結原料の粒子配置分布解析方法。
【図面の簡単な説明】
【００２０】
【図１】粒子が衝突している図である。
【図２】衝突している粒子に働く法線方向の力を説明する図である。
【図３】粒子の衝突過程の図である。
【図４】粒子の衝突過程の図である。
【図５】粒子の衝突過程の図である。
【図６】複数粒子の運動を計算するフローチャートである。
【図７】焼結原料分級設備概要の図である。
【図８】粒子挙動計算を説明する図である。
【発明を実施するための形態】
【００２１】
［第一の実施形態］
（粒子の衝突の判定とめりこみ量）
対象とする粒子の形状は特に限定しないが、ここでは球である場合を例にとって説明する。半径ｒ1の球である粒子１と、半径ｒ2の粒子２の二つの粒子が接触しているかどうかの判定は、種々の方法があるが、例えば粒子１の中心と粒子２の中心の距離がｒ1＋ｒ2以下であれば接触していると判定してもよい。粒子が球でなければ、粒子の重心をその中心とし、平均半径をその半径としてもよい。
【００２２】
めり込み量は、前記ｒ1＋ｒ2と接触している二つの粒子の中心間距離の差によって求められる。計算に用いる場合は、そのようにして得られた差を更に、前記ｒ1＋ｒ2で除した値を用いる場合もある。
【００２３】
粒子１と粒子２相互間のめり込み量をδとし、二つの粒子の中心の、その中心を結ぶ線上で計った位置をｘ1、ｘ2とすれば、二つの粒子の中心距離は、｜ｘ1−ｘ2｜と表される。これらは、時間的に変化する。各々の量の時刻ｔpにおける値であることを右肩添え字ｐで示せば、時刻ｔpにおけるめり込み量δpは次のように表すことができる。めり込みの例は図１に示されている。
δp＝｛（ｒ1＋ｒ2）−｜ｘ1p−ｘ2p｜｝／（ｒ1＋ｒ2） ［数式１］
【００２４】
ある時点で二つの粒子の接触判定を行った場合、接触と判定されたものについて、二つの粒子の中心距離は各々の半径の和よりもかなり小さく、即ち大きいめりこみ量である場合が多い。このようなめり込みが大きい場合、反発力も大きく、時間刻みを細かくとらないと、大きな反発力によって粒子が急速に反発してしまい、計算が不安定になってしまう場合がある。
【００２５】
（接触面法線方向に働く力）
接触した二つの粒子が離れるまでの衝突過程は、その接触面に対する法線方向に働く力（以下、接触法線力という。）が大きな影響を与える。即ち、該接触面から該法線に沿って大きく後退すれば衝突過程は早く終了する。従って、該法線方向に働く力を調整することで衝突過程を省略できるのである。また、ここで時間の推移は、刻み時間Δｔずつ推移するものとする。即ち、時刻ｔpにおいて、時刻指標ｐ＝０からｐΔｔだけ時間が経過している。
【００２６】
非特許文献１においては、時刻ｔpにおいて接触している二つの粒子に、その接触面の法線方向に働く力としては、粒子１に着目した場合、（１）めり込み量に比例する反発力−ｋδp 、および、（２）二つの粒子の相対速度に比例する粘性力、η（Ｖ2p−Ｖ1p）の和として与えられる。ここで右肩の添え字ｐは時刻指標を表す。δpは時刻ｔpにおける粒子１及び２のめり込み量である。また、Ｖ1pおよびＶ2pは各々、時刻ｔpにおける粒子１および２の前記接触面法線方向速度である。また、ｋおよびηは各々比例係数である。これらは図２に示すように、バネ定数ｋのバネと粘性係数ηのダッシュポットが並列につながったモデルで表される。即ち、非特許文献１にあるように、先行技術においては、粒子１に働く前記接触法線力Ｆ1pは
Ｆ1p＝−ｋδp＋η（Ｖ2p−Ｖ1p） ［数式２］
と表される。また、粒子２に働く該法線力は、［数式２］の下添え字の１と２を入れ替えればよい。
【００２７】
時間刻みΔｔを大きくすると、時刻ｔpで接触判定されても、粒子はすぐに離れ、衝突過程を省略し、計算が高速化できるが、上記［数式２］の反発力の項が大きくなり、後述するように粒子の速度が大きくなりすぎたりして、計算が不安定となったり、誤差が大きくなったりしてしまう。従って、先行技術においては、Δｔを細かくとって衝突過程を細かい時間間隔で追跡して計算する必要があった。本発明者等は、上記（２）の粘性力の表式に着目した。そこで、粒子１の時刻ｔpにおける該法線方向速度Ｖ1pを、粒子１の時刻ｔp+1における速度、即ちＶ1p+1で置き換えたのである。ここで、時刻ｔp+1とは、時刻ｔpから刻み時間Δｔだけ進んだ時刻である。しかしながら、時刻ｔpの時点で、時刻ｔp+1における正確な速度を把握することは困難なので、Ｖ1p+1は時刻ｔp+1における予測値Ｅ1p+1で表す。即ち、本発明においては、粒子１に働く接触法線力Ｆ1pは、前記［数式２］のＶ1pをＥ1p+1で置き換え、
Ｆ1p＝−ｋδp＋η（Ｖ2p−Ｅ1p+1） ［数式３］
と表される。粒子２に働く接触法線力Ｆ2pは、［数式３］の下添え字の１と２を入れ替えればよい。即ち
Ｆ2p＝−ｋδp＋η（Ｖ1p−Ｅ2p+1） ［数式４］
と表す。
【００２８】
予測法線方向速度の算出方法はいくつかある。例えば、本発明になる算出方法によれば、時間刻みΔｔを粗くしても、時刻ｔpにおいて衝突していても、時刻ｔp+1においては粒子は離れている場合が多く、かつ反発力も大きすぎず、精度を保ったまま計算量が減少し、［数式２］を用いた先行技術の場合と比較して計算時間が１／４以下になることが分かった。ここで、注意すべきなのは、粒子２に働く法線方向粘性力は、通常なら粒子１に働く粘性力の反力となるために同じ大きさで反対方向の力となり、それぞれ、作用、反作用の関係にあり、その和はゼロとなることである。本発明においては、［数式４］となり、これは、［数式３］と比較して必ずしも同じ大きさとはならない。従って、接触している２粒子において、これら、法線方向の粘性力の和はゼロではなく、有限な値を持った残留力となり、作用、反作用の関係が満たされず、厳密にいえば物理的に正しい描像を与えない。しかし、そのために生じる誤差は、無視できるほど小さいことが分かった。更に、多数の粒子の衝突を考えた場合、このような残留力は、ランダムな方向と大きさとなり、結果的にその総和はさらに小さくなり、精度に対する影響はほとんどないことが分かったのである。
【００２９】
（予測法線方向速度）
粒子１の予測法線速度Ｅ1p+1を、粒子１が衝突している粒子２のみから受ける力で運動した場合の時刻ｔp+1における粒子１の法線方向速度と考える。従って、次の離散的な運動方程式から求められる。
ｍ1（Ｅ1p+1−Ｖ1p）／Δｔ＝Ｆ1p ［数式５］
【００３０】
上式のＦ1pに［数式３］を代入し、整理すれば、
【数５】
が得られる。
【００３１】
［数式６］を［数式３］に再度代入すれば、時刻ｔpにおいて粒子１が粒子２より受ける法線方向力は、
【数６】
となる。同様に、粒子２について予測法線速度Ｅ2p+1 及び法線方向力Ｆ2pは以下の［数式８］、［数式９］のように表される。
【数７】
【数８】
【００３２】
以上より、予測法線方向速度Ｅ1p+1及びＥ2p+1が時刻ｔpの量で示され、法線方向力も時刻ｔpにおける値で示されたため、効率良いプログラム作成が可能となり、更に、前述のように精度を保ったまま、計算時間が短縮できるのである。
【００３３】
（接触判定された粒子の接触後の衝突運動の計算の省略）
時刻ｔpで二つの粒子が接触しているとした場合のその後の衝突過程の進展の概念図を図３〜図５に示す。従来技術である［数式２］で法線方向力を表した場合、時間刻みΔｔを細かくとり、図３に示すように、粒子が離れるまでの衝突過程を細かく計算する必要がある。例えば、同じ法線方向力のままΔｔを数倍に取った場合の衝突過程は、図４に示すように、接触後の粒子は次のステップで大きくめり込み、その次のステップでは反発力のために各々高速度で反発してしまい、計算が不安定になり、誤差が大きくなってしまう。次に、大きい時間刻みΔｔのままで［数式３］の本発明を用いると、図５に示すように、接触してから次のステップである程度めり込むが、その後、穏やかに離れ、精度を保ったまま衝突過程の計算を省略できる。これは、［数式３］を用いた場合、めり込み量が大きくなり、反発力がかなり大きくなっても、逆方向に働く粘性力も大きくなり、ほどよくバランスされるからである。即ち、［数式３］を用いれば、精度を保ったまま時間刻みを大きくとって衝突過程を省略し、計算速度を短くできるのである。
【００３４】
以上の計算ステップの大幅な省略のため、本発明になる計算方法を用いる場合、時間刻みΔｔは、従来技術の場合の時間刻みの４倍〜１０倍とっても、計算結果の精度はほとんど変わらないことが分かった。即ち、従来技術による計算の精度と同一の精度のまま計算時間が１／４〜１／１０になることが分かった。
【００３５】
（接触粒子の運動）
多数の粒子からなる粒子群においては、一つの粒子が同時に複数の粒子と衝突している場合が多い。その場合も、各々の二粒子の衝突について、前記の法線方向力をもとめその合力が粒子に働くと考える。
【００３６】
例えば、ある着目粒子ｉを固定し、他の対向粒子ｊについて粒子ｉと衝突しているかどうか判定する。衝突していると判定されたものについては、各々の接触面について法線方向力を［数式７］及び［数式９］において表される二粒子間の法線方向力を求める。これを繰り返し行い、粒子ｉについてすべての衝突粒子を求めたら、別の粒子に着目して同様な探索を行い、全ての粒子について衝突粒子と、各衝突における法線方向力を求め、これらの和を求める。但し、各々の衝突において接触面法線は異なる方向を向いているからベクトル的な和を求める必要がある。例えば着目粒子ｉにＮ個の対向粒子ｊが衝突しているとすれば、各々のｉとｊの衝突に対して異なる方向の法線速度を考え、前記［数式３］〜［数式９］を適用すればよい。即ち、着目粒子ｉの、各対向粒子ｊとの各接触面法線方向の粒子ｉの時刻ｔp+1での予測速度Ｅjip+1は、［数式６］と同様に考えればよく、
【数９】
となる。ここで、δijpは時刻ｔpにおける粒子ｉ，ｊの衝突によるめり込み量、Ｖijpは時刻ｔpにおける粒子ｊの、粒子ｉ、ｊの接触面法線方向の速度成分である。Ｖjipは、同様に時刻ｔpにおける粒子ｉの前記接触面法線方向の速度成分である。
【００３７】
従って、時刻ｔpにおいて粒子ｉと粒子ｊとの衝突から粒子ｉが受ける法線方向の力Ｆjipは、［数式７］と同様に考えて次のように表される。
【数１０】
【００３８】
この力は、前記のように、時刻ｔpにおける粒子ｉと粒子ｊの接触面に対する法線方向の量であり、そのような法線方向の単位ベクトルをｅjipであらわせば、上記［数式１１］の力の方向を示すベクトルとなる。
【００３９】
時刻ｔpにおいて粒子ｉがＮ個の衝突から受ける各々の法線方向力の合力は、［数式１１］で表される力を方向ｅjipを持つベクトルとして対向粒子ｊについて和を取ることで得る。
【数１１】
【００４０】
更に、各々の衝突における接触面の接線方向の力、及び各々の衝突による、回転モーメントも求め、その合力を求める。ここで、各接触面の接線方向の力及び回転モーメントは、非特許文献１にあるように、先行技術を用いて算出することができる。即ち、接触面の接線方向の力については、せん断抗力を与える弾性スプリングと粘性ダッシュポットの並列配置を仮定し、接触点近傍のせん断変形が主として要素間の摩擦力によって生ずるとして計算を行う。回転モーメントについては、要素毎に回転モーメントの運動方程式を立て、変位の運動方程式と同様に差分法による数値解析で求めることができる。これらの合力をもとめたら、先行技術と同様にして離散的な運動方程式を解き、時刻ｔp+1における速度と位置を求める。以上のステップのフローチャートを図６に示す。
【００４１】
（粒子挙動解析方法）
多粒子系粒子の挙動を解析するための粒子挙動解析方法において、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることができる。粒子挙動解析方法においては、個々の粒子の運動を運動方程式を用いて離散要素法で計算する。所定の空間領域内に存在するすべての粒子について、刻み時間をΔｔとし、各粒子の運動方程式を差分方程式として数値計算を行う。他の粒子と衝突していない粒子については、外場がなければ等速直線運動、重力や電場等の外場があればそのような運動にさらに外場の力による運動が加わる。ここにおいて、任意の２個の粒子が衝突してから離れるまでの運動の計算について、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いる。
【００４２】
系の初期条件と境界条件を定めた上で、以上の手順で、多粒子系粒子を構成するすべての粒子について挙動解析を行い、時間の経過と共に対象とする多粒子系粒子が全体としてどのように運動するのかを明らかにすることができる。
【００４３】
本発明になる計算方法、装置を用いて、大きさの異なる粒子によって構成される粒子群が、重力加速度のもとで種々の装置の中を通過し、最終的に運動を停止するまでを追跡し、最終的に粒子群が形成する塊における粒子の空間粒度分布をもとめることが可能である。この時、粒子群が通過する装置として、シャトルコンベヤや、ホッパー、ドラムフィーダ、フルイ等の機械要素を仮想的に配置し、粒子群を焼結原料の粒を模した粒径分布を持つように決めてやることによって実際の焼結原料装入装置をモデル化することができる。最終的な粒子の空間粒径分布としては、大きい粒子が下になり、小さい粒子が上にくるような分布が形成されることが作業効率上望ましい。しかしながら、前述のように、機械要素の空間的な位置関係、操業条件によってそのような望ましい分布から外れた分布ができてしまう場合が多い。本発明になる計算方法、解析方法を用いて、焼結原料装置をモデル化し、最終的に形成される空間粒径分布を求めることができる。さらに、機械要素の位置、操業条件を変えて計算を実施することによって理想的な空間粒径分布を求めることができる。即ち、焼結原料装入装置の設計、操業条件の決定を著しく高速で行うことができるのである。
（焼結原料の粒子配置分布解析方法）
即ち、上記本発明の粒子挙動解析方法を、焼結原料の粒子配置分布解析方法として用いることができる。図７に示すような焼結原料分級設備において、シャトルコンベヤ１からサージホッパー２に投入され，ドラムフィーダ３から切り出され，フルイ４の上を通って落下してパレット５上の堆積物６となるた焼結原料の粒子配置分布を、本発明の粒子挙動解析方法によって求める。まず、シャトルコンベヤ１、サージホッパー２、ドラムフィーダ３、フルイ４を空間に配置し、それぞれの装置に各粒子が接触した際の境界条件を定める。次いで、シャトルコンベア内に搭載された焼結原料の粒子配置分布を初期条件として与える。以上の準備を整えた上で、系内に存在するすべての粒子（要素）について、運動方程式を刻み時間Δｔの差分方程式として数値計算を行う。その際、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いる。
【００４４】
このような計算を行った結果として、シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を、時間経過と共に明らかにすることができる。
【００４５】
（接触面法線方向の力にかかる係数）
［数式３］に用いられる係数ｋとηはそれぞれ，弾性係数および粘性係数と考えることができるが，単純な定数とするよりも、非特許文献２及び非特許文献３等に開示されたように、次のように設定することで精度のよい計算を効率良く行うことができる。
【００４６】
衝突している２粒子を粒子ｉ及び粒子ｊとする。
【００４７】
【数１２】
【数１３】
ここで，各変数は以下のように表わされる。
１／Ｅ*＝（１−νi2）／Ｅi＋（１−νj2）／Ｅj ［数式１５］
１／ｒ*＝１／ｒi＋１／ｒj ［数式１６］
１／ｍ*＝１／ｍi＋１／ｍj ［数式１７］
【数１４】
上記［数式１３］〜［数式１８］の範囲内で，Ｅはヤング率、νはポアソン比，ｒは衝突粒子の半径，ｍは衝突粒子の質量、ｅは反発係数，πは円周率を表わす。
【００４８】
（画像処理用演算装置ＧＰＵを用いた解析）
粒子挙動解析の計算に用いる計算機として、画像処理を高速処理するための演算装置であるところのグラフィック・プロセシング・ユニット（ＧＰＵ）で計算を行うと、通常の中央演算装置（ＣＰＵ）を用いる汎用計算機を用いた同じ計算と比較して、計算時間が１／５０以下に短縮できることが分かった。即ち、本発明になる計算方法を用いて作成されたプログラムをＧＰＵに適用すると、従来技術の計算方法を用いたプログラムをＣＰＵに適用した場合に比較して計算時間を１／２００以下に短縮することができるのである。ただし、ＧＰＵで実施する場合は、計算機が実施するプログラム形態がＣＰＵに適用するプログラム形態とは異なるので、プログラム形態を変更しなければならない。しかしながら、計算のアルゴリズムは同一であり、本発明になる計算方法になるアルゴリズムを用いることは同じであるので、ＧＰＵを用いて計算しても、本発明になるものであることは明白である。また、本発明になる計算方法には、ＧＰＵによる計算の特徴である並列計算を阻害する要因は全くない。従ってＧＰＵによる計算の特徴を効率良く引き出すことが可能である。
【実施例】
【００４９】
［実施例１］（汎用計算機を用いた解析）
〔数式１３〕〜［数式１８］に示された係数を用いて、Δｔ＝１０-6秒の時間刻みで、５種類の大きさの球状粒子、２×１０4個からなる粒子群が、ホッパーからパレットに落下する計算を、汎用計算機を用いてシミュレーションした。従来技術で実施したものは、約５０時間でシミュレーションが終了した。同様なシミュレーションを、時間刻みを５０倍の５×１０-5秒にして実施したところ、従来技術では、正確な解が得られなかった。しかしながら、時間刻みが５×１０-5のままで本発明になる計算方法を用いると、約５時間でシミュレーションが終了し、従来技術で１０-6の時間刻みで行った計算とほぼ同様な粒子分布を得た。即ち、前記「予測速度」に起因する誤差は、ほとんど無視できることが分かった。さらに、本発明になるプログラムにより、精度を保ったまま、計算時間が１／１０に短縮されたことが分かった。
【００５０】
［実施例２］
（汎用計算機を用いた解析と画像処理用演算装置ＧＰＵを用いた解析との比較）
２×１０5個の粒子からなる粒子群の挙動の計算を、汎用計算機を用いて１×１０-6秒の時間刻みで従来の計算方法を用いてシミュレーションしたところ、計算終了まで約９８日かかった。同様な計算を、時間刻みを５×１０-5秒として、本発明の計算方法を用いて実施したところ、約８日かかった。時間刻み５×１０-5秒で、本発明の計算方法を適用したプログラムをＧＰＵを用いた装置で動作するように変換してＧＰＵを持つ装置で計算した。その結果１２時間で計算が終了し、従来法で１×１０-6秒の時間刻みの場合とほぼ同じ結果を得た。従って、汎用計算機を用いた従来技術に対して精度を維持したまま、計算時間が約１／２００に短縮された。
【００５１】
［実施例３］
図７に示す装置を用い、シャトルコンベヤ１からサージホッパー２に投入され，ドラムフィーダ３から切り出され，フルイ４の上を通って落下した焼結原料の分布をシミュレーションするにあたり，本発明による粒子運動の計算方法を、汎用計算機を用いてシミュレーションを行った。まず、シャトルコンベヤ、サージホッパー、ドラムフィーダ、フルイを空間に配置し、それぞれの装置に各粒子が接触した際の境界条件を定めた。次いで、シャトルコンベア内に搭載された焼結原料の粒子配置分布を初期条件として与えた。以上の準備を整えた上で、系内に存在するすべての粒子（要素）について、運動方程式を刻み時間Δｔの差分方程式として数値計算を行った。その際、上記本発明の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いた。
【００５２】
図８に，フルイ部がプレートと複数の条材からなる装置における粒子挙動の計算結果を示す。落下した粒子は下部のパレット上に積もる。下部のパレットは図中右から左に所定の速度で移動しており，堆積した粒子は図中左方向に運ばれて次の工程に進む。計算では，３種類の大きさの粒子を考慮した（大粒子（最も明るい、白 粒子径６ｍｍ），中粒子（最も暗い、黒 粒子径５ｍｍ），小粒子（中間の明るさ、灰色 粒子径３ｍｍ））。各々の粒子は、同じ重量密度、３．３ｇ／ｃｍ3とした。条材フルイの間隔は先端（右下）ほど広くなっており，大粒子ほど図右よりに落下する。落下した大粒子は堆積物斜面上を右方向に転がり，斜面右端から堆積物下部に巻き込まれ易くなる。一方小粒子は，フルイ部付け根あたりから落下し，下部バレットの堆積物の上面に降り積もる。これらにより，堆積物は，底から大粒子（最も明るい、白），中粒子（最も暗い、黒 ），小粒子（中間の明るさ、灰色）と分級されることが分かった。この結果を得るのに要した計算時間は１０時間であった。さらに、同様な装置構成で同様な粒子挙動を、従来の離散要素法を用いて汎用計算機を用いて計算したところ約１００時間かかった。本発明の手法により効率が向上したことが分かる。
【００５３】
［実施例４］
ＧＰＵを用いた画像処理専用計算機に、実施例３で用いた、本発明になるプログラムを移植し、実施例３と同じ要素配置での計算を行った。この時実施例３と同様な結果を得たが、計算に要した時間は２時間となり、汎用計算機を用いた実施例３に比較して１／５０以下となった。更に、画像処理専用の計算機を用いたために、途中結果や最終結果のグラフィック表示も効率良く行うことができ、設計効率が著しく向上した。
【００５４】
［実施例５］
ＧＰＵを用いた画像処理専用計算機を用いて、フルイ部がプレートと複数の条材からなる装置において，プレートの長さと角度，条材の開き角度を変えて計算を行った。その結果，堆積した粒子の分布は，プレート長さが０．５〜１．０ ｍ，角度が４０〜５０度，条材の開き角度が１．５〜４．０度の範囲が良いことが分かった。
【符号の説明】
【００５５】
１ シャトルコンベア
２ サージホッパー
３ ドラムフィーダ
４ フルイ
５ パレット
６ 堆積物
【特許請求の範囲】
【請求項１】
粒子の運動の計算における任意の二つの粒子が衝突してから離れるまでの運動を計算する方法であって、時刻ｔpにおいて衝突している二粒子の接触面法線方向（以下、「法線方向」という。）に働く力を、弾性係数がｋであり、時刻ｔpにおける粒子のめり込み量δpに比例する弾性反発力
−ｋ×δp
と、粘性係数がηの粘性抵抗力の和で表し、
一方の粒子１の、時刻ｔpより刻み時間Δｔだけ進んだ時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
を、粒子１に働く法線方向粘性力とし、さらに、他方の粒子２に働く法線方向粘性力を、同様に時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表すことを特徴とする、二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
ここで、時刻指標ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。また、右上添え字ｐは、時刻指標ｐにおける当該物理指標を意味する。
【請求項２】
粒子１の質量をｍ1及び、粒子２の質量をｍ2として、時刻ｔp+1における、粒子１の予測法線方向速度Ｅ1p+1及び粒子２の予測法線方向速度Ｅ2p+1をそれぞれ
【数１】
及び
【数２】
で表すことを特徴とする、請求項１に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
【請求項３】
時刻ｔpにおいて、粒子１及び粒子２が衝突から受ける法線方向の力Ｆ1p及びＦ2pを
【数３】
及び
【数４】
で表すことを特徴とする、請求項１又は２に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
【請求項４】
多粒子系粒子の挙動を解析するための粒子挙動解析方法であって、個々の粒子の運動を運動方程式を用いて離散要素法で計算し、任意の２個の粒子が衝突してから離れるまでの運動の計算について請求項１乃至３のいずれかに記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることを特徴とする粒子挙動解析方法。
【請求項５】
シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を求めるにあたり，請求項４に記載の粒子挙動解析方法を用いることを特徴とする、焼結原料の粒子配置分布解析方法。
【請求項１】
粒子の運動の計算における任意の二つの粒子が衝突してから離れるまでの運動を計算する方法であって、時刻ｔpにおいて衝突している二粒子の接触面法線方向（以下、「法線方向」という。）に働く力を、弾性係数がｋであり、時刻ｔpにおける粒子のめり込み量δpに比例する弾性反発力
−ｋ×δp
と、粘性係数がηの粘性抵抗力の和で表し、
一方の粒子１の、時刻ｔpより刻み時間Δｔだけ進んだ時刻ｔp+1における法線方向速度Ｖ1p+1の予測値Ｅ1p+1及び他方の粒子２の時刻ｔpにおける法線方向速度Ｖ2pとの差に比例する力
η×（Ｖ2p−Ｅ1p+1）
を、粒子１に働く法線方向粘性力とし、さらに、他方の粒子２に働く法線方向粘性力を、同様に時刻ｔp+1における粒子２の法線方向速度の予測値Ｅ2p+1と時刻ｔpにおける法線方向速度Ｖ1pとの差に比例する力
η×（Ｖ1p−Ｅ2p+1）
で表すことを特徴とする、二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
ここで、時刻指標ｐは整数であり、時間がΔｔ進行するとｐの値が１増加する。また、右上添え字ｐは、時刻指標ｐにおける当該物理指標を意味する。
【請求項２】
粒子１の質量をｍ1及び、粒子２の質量をｍ2として、時刻ｔp+1における、粒子１の予測法線方向速度Ｅ1p+1及び粒子２の予測法線方向速度Ｅ2p+1をそれぞれ
【数１】
及び
【数２】
で表すことを特徴とする、請求項１に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
【請求項３】
時刻ｔpにおいて、粒子１及び粒子２が衝突から受ける法線方向の力Ｆ1p及びＦ2pを
【数３】
及び
【数４】
で表すことを特徴とする、請求項１又は２に記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法。
【請求項４】
多粒子系粒子の挙動を解析するための粒子挙動解析方法であって、個々の粒子の運動を運動方程式を用いて離散要素法で計算し、任意の２個の粒子が衝突してから離れるまでの運動の計算について請求項１乃至３のいずれかに記載の二粒子が衝突してから離れるまでの法線方向の運動の計算方法を用いることを特徴とする粒子挙動解析方法。
【請求項５】
シャトルコンベヤからサージホッパーに投入され，ドラムフィーダから切り出され，フルイの上を通って落下した焼結原料の粒子配置分布を求めるにあたり，請求項４に記載の粒子挙動解析方法を用いることを特徴とする、焼結原料の粒子配置分布解析方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２０１３−２３６９８（Ｐ２０１３−２３６９８Ａ）
【公開日】平成２５年２月４日（２０１３．２．４）
【国際特許分類】
【出願番号】特願２０１１−１５６４８２（Ｐ２０１１−１５６４８２）
【出願日】平成２３年７月１５日（２０１１．７．１５）
【出願人】（０００００６６５５）新日鐵住金株式会社 (6,474)
【Ｆターム（参考）】
【公開日】平成２５年２月４日（２０１３．２．４）
【国際特許分類】
【出願日】平成２３年７月１５日（２０１１．７．１５）
【出願人】（０００００６６５５）新日鐵住金株式会社 (6,474)
【Ｆターム（参考）】
[ Back to top ]