生命科学関連特許情報
| タイトル： | 公開特許公報(A)_粒子挙動シミュレーション装置、粒子挙動シミュレーション方法、制御プログラム、および記録媒体 |
| 出願番号： | 2010264857 |
| 年次： | 2012 |
| IPC分類： | G06F 19/00,G01N 11/00 |
特許情報キャッシュ
児林 智成 島田 直樹 田中 敏嗣 JP 2012118579 公開特許公報(A) 20120621 2010264857 20101129 粒子挙動シミュレーション装置、粒子挙動シミュレーション方法、制御プログラム、および記録媒体 住友化学株式会社 000002093 国立大学法人大阪大学 504176911 特許業務法人原謙三国際特許事務所 110000338 児林 智成 島田 直樹 田中 敏嗣 G06F 19/00 20110101AFI20120525BHJP G01N 11/00 20060101ALI20120525BHJP JPG06F19/00 110G01N11/00 Z 9 4 ＯＬ 23 本発明は、多数の粒子の挙動を予測する粒子挙動シミュレーションに関する。 化学合成プロセスにおいて、気体から反応生成物を得る触媒反応装置として、流動層装置が広く用いられている。流動層装置は、処理容器内に粉体等を入れておき、粉体の下部から処理容器内に気体を吹き込むことで、粉体を流動化させる装置である。例えば、流動層装置の処理容器内に触媒粒子を入れておき、処理容器の下部から原料ガスを供給し、流動化状態の触媒粒子と原料ガスとを効率よく接触させることによって、ガス状の反応生成物を得ることができる。 効率よく反応を起こさせる流動層装置を設計するため、または触媒粒子の量および気体の流量等のよりよい条件を見つけるために、流動層における粒子の挙動をシミュレーションすることが重要である。 このような流動層触媒反応に用いられる触媒粒子としては、平均粒径が３０〜１５０μｍで、密度が７００〜２５００ｋｇ／ｍ３程度の粒子、いわゆるGeldart線図においてＡ粒子に分類される粒子が広く用いられている。Geldart線図における粒子の分類は、流動層の粒子の流動化状態を予測する指標となる。Geldart線図においてＡ粒子に分類される粒子を触媒粒子として用いれば、流動層が小さい気泡を形成するので、反応効率を高くすることができる。 流動層におけるこのような粒子の挙動をシミュレーションする方法の１つに、離散要素法（個別要素法：ＤＥＭ：Discrete Element Method）がある。離散要素法を用いた数値解析では、分割された時間刻み毎に、各粒子に作用する力に基づいて各粒子に対して運動方程式を計算し、時刻を時間ステップだけ進めた次の時刻における粒子の状態を数値計算によって求める。そして、時間刻み毎の計算を繰り返し、粒子の挙動を求める。ある粒子に作用する力としては、他の粒子または処理容器の内壁との衝突による反発を表す接触力、粒子の周囲の気体（流体）から受ける流体力、他の粒子または処理容器の内壁との間に働く付着力、および重力等がある。 離散要素法において、シミュレーションの時間刻みの幅（時間ステップ）が小さければ、精度のよいシミュレーション結果を得ることができる。しかしながら、時間ステップを小さくすると、所定の時間（例えば１秒間）の粒子の挙動を求めるための計算量が多くなる。そのため、シミュレーションにかかる計算時間が長くなる。時間ステップを大きくすれば計算時間を短くすることができるが、一方で、計算を安定に行う（妥当なシミュレーション結果を得る）ための時間ステップには上限がある。 非特許文献１には、粒子の衝突の接触力を表す線形ばねモデルを仮定した場合、妥当なシミュレーション結果を得るための時間ステップの経験的な上限が、ばねの固有振動の周期Ｔの約１／１０であることが記載されている。この場合、半周期Ｔ／２が粒子同士の接触時間になるので、粒子同士の接触時間の約１／５が時間ステップの経験的な上限になる。そして、固有振動の周期Ｔはばね定数ｋに依存するので、ばね定数が大きな粒子、すなわち硬い粒子の場合、時間ステップを小さくする必要がある。そのため、ばね定数が小さい軟らかい粒子に比べて、硬い粒子のシミュレーションは、計算時間が長くなる。 また、非特許文献１には、付着性を持たない粒子からなる系について、本来のばね定数よりも非常に小さい修正ばね定数を用いて、離散要素法によって粒子の挙動を計算しても、修正ばね定数がある程度大きければ、結果に大きな影響がないことが記載されている。すなわち、付着性を持たない粒子からなる系について、本来のばね定数よりも小さい修正ばね定数を用いれば、シミュレーションの計算時間を短くして、妥当なシミュレーション結果を得ることができる。社団法人日本粉体工業技術協会、「流動層ハンドブック」、株式会社培風館、１９９９年３月２５日、p.125-144 しかしながら、Geldart線図においてＡ粒子に分類される粒子のように、粒子径が小さい場合、ファンデルワールス力および液架橋力等に起因する粒子同士の付着力が大きくなり、付着力の影響が無視できなくなる。付着力とは、粒子が他の粒子または処理容器の内壁等と近接した場合に受ける引力である。付着力が大きく、付着力の影響が無視できない系において、本来のばね定数よりもばね定数を小さく設定して時間ステップを長くすると、シミュレーションの精度が悪くなり、妥当な結果が得られないという問題がある。そのため、付着力の影響を無視できない系では、シミュレーションの計算に長い時間がかかっていた。 本発明は、上記の問題点に鑑みてなされたものであり、その目的は、付着力の影響を考慮しながら、離散要素法を用いて複数の粒子の挙動を計算する時間を短縮することができる粒子挙動シミュレーションを実現することにある。 本発明に係る粒子挙動シミュレーション装置は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、上記の課題を解決するために、粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定部と、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴としている。 本発明に係る粒子挙動シミュレーション方法は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション方法であって、粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、シミュレーション装置の接触力特定部が、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定ステップと、上記シミュレーション装置の挙動計算部が、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算ステップとを含むことを特徴としている。 上記の構成によれば、ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて修正接触力を特定し、修正接触力を用いて離散要素法によって粒子の挙動を計算する。よって、粒子同士の接触時間が長くなるので、時間ステップを長くし、挙動計算に要する時間を短縮することができる。また、挙動計算において、修正ばね定数ｋｄを用いると同時に、ばね定数ｋｒと修正ばね定数ｋｄとの比に応じて付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いる。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。 なお、粒子に働く接触力を表すばねモデルは、線形ばねモデルでも非線形ばねモデルでもよい。また、上記挙動計算部は、粒子同士の付着力を、粒子同士が所定の距離より近接しているときのみに粒子に働く力として扱ってもよい。 また、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根とほぼ同じになるようにしてもよい。 上記の構成によれば、小さい修正ばね定数ｋｄを用いて挙動を計算しても、粒子同士の衝突後に、２つの粒子が付着したままになるか、離れていくかの閾値となる速度ｖｃが変化しない。よって、小さい修正ばね定数ｋｄを用いて挙動を計算しても、粒子の挙動が変化しない。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。 また、上記修正付着力Ｆａｄは、条件として、後述の式（１５）を満たすようにしてもよい。 上記の構成によれば、修正付着力Ｆａｄの値は、修正ばね定数ｋｄをばね定数ｋｒで除したものの平方根に付着力Ｆａｒを乗じた値の、０．７倍から１．３倍までの範囲になる。 このように、修正付着力Ｆａｄが、修正ばね定数ｋｄをばね定数ｋｒで除したものの平方根に付着力Ｆａｒを乗じた値から少しずれてもよい。 また、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根と同じオーダーになるようにしてもよい。 また、上記修正ばね定数ｋｄは、上記ばね定数ｋｒの０．１倍以下としてもよい。 上記の構成によれば、修正ばね定数ｋｄに応じて時間ステップを長くして、挙動計算に要する時間を短縮することができる。 また、上記付着力Ｆａｒと上記修正ばね定数ｋｄと上記ばね定数ｋｒとに基づき、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根と同じになるような上記修正付着力Ｆａｄを特定する修正付着力特定部を備える構成であってもよい。 上記の構成によれば、修正付着力特定部が、小さい修正ばね定数ｋｄを用いて挙動を計算しても粒子の挙動が変化しない修正付着力Ｆａｄを特定する。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。 なお、上記粒子挙動シミュレーション装置は、一部をコンピュータによって実現してもよく、この場合には、コンピュータを上記各部として動作させる制御プログラム、および上記制御プログラムを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に入る。 本発明に係る粒子挙動シミュレーション装置は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定部と、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴としている。 そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。シミュレーションの対象とする流動層装置の概略構成を示す断面図である。２つの粒子が衝突する様子を時刻毎に示した図である。２つの粒子が衝突した後の２つの状態を示す図である。本発明の一実施形態の粒子挙動シミュレーション装置の機能的構成を示すブロック図である。上記粒子挙動シミュレーション装置の処理を示すフロー図である。本発明の別の実施形態の粒子挙動シミュレーション装置の機能的構成を示すブロック図である。実験例における粒子層の状態を撮影した画像である。従来技術による比較例のシミュレーション結果を示す画像である。従来技術による比較例のシミュレーション結果を示す画像である。本発明の実施例のシミュレーション結果を示す画像である。 以下の実施の形態では、主に、流動層装置における粒子の挙動をシミュレーションする場合について説明するが、本発明はこれに限定されるものではない。 ［実施の形態１］ 以下、本実施の形態について、図１〜図５を参照して詳細に説明する。本実施の形態では、図１に示す流動層装置における粒子をシミュレーションの対象とする。 ＜シミュレーション系＞ 図１は、流動層装置１の概略構成を示す断面図である。流動層装置１は、処理容器２、分散板３、気体供給口４、および排気口５を備える。分散板３は、処理容器２内の下部に容器内を分割するように水平に配置されており、微細な貫通孔を有する。処理用に２内の分散板３の上には、触媒粒子６が堆積されている。分散板３は、気体を通過させるが触媒粒子６は通過させない。気体供給口４は、処理容器２の分散板３の下部に設けられている。原料ガスは、気体供給口４から処理容器２内に供給され、分散板３の下から堆積した粒子の隙間を通って噴き上げる。原料ガスおよびその反応物とは処理容器２の上部に設けられた排気口５から排出される。 ＜シミュレーション方法＞ 離散要素法では、着目する粒子に働く力から、次時刻における粒子の位置および運動状態を計算する。着目する粒子に働く力Ｆは、次式で与えられる。ここで、Ｆｄは流体力、Ｆｇは重力、Ｆｃは接触力、Ｆａは付着力をそれぞれ表す。 流体力Ｆｄは、粒子が吹き上げる気体（原料ガス）から受ける、主に上向きの力である。重力Ｆｇは自重によって粒子に働く力である。接触力Ｆｃは、粒子同士の接触、または粒子と処理容器の内壁（または分散板）との接触による反発力を表す。付着力Ｆａは、粒子が他の粒子または処理容器の内壁と近接または接触した際に受ける、ファンデルワールス力および液架橋力等に起因する引力を表す。流体力Ｆｄおよび重力Ｆｇは、各粒子が常に受ける力である。ただし、流体力Ｆｄは、一定ではなく、位置および周囲の粒子の分布（空間密度）等によって変化する。一方、接触力Ｆｃおよび付着力Ｆａは、例えば粒子が他の粒子または内壁と接触したときにのみ働く。 粒子が他の粒子と接触した場合の接触力Ｆｃおよび付着力Ｆａについて、説明する。ここでは、粒子が接触した場合、粒子が弾性変形をして、反発力を受けるモデルを用いて接触力を表す。半径ｒの粒子同士が衝突した際に粒子が受ける接触力について、接触面の法線方向の力を、弾性反発力を表す線形ばねと、粘性減衰力を表すダッシュポットによって表現する。粘性減衰係数は、粒子の反発係数から決定される。なお、接触面の接線方向の力は、粒子の変形量が小さい場合は、法線方向と同様の線形ばねとダッシュポットの組み合わせによって表現され、粒子の変形量が大きい場合は、滑りを表すための摩擦スライダを用いて表現される。 法線方向の接触力Ｆｃは、次式で表される。ここで、ｋは粒子の弾性反発力を表す法線方向のばね定数、ηはダッシュポットによる減衰を表す法線方向の粘性減衰係数、ｘは接触時を基準とした粒子の中心位置の変位、ｖは粒子の相対速度ベクトル、ｎは法線方向の単位ベクトルを表す。 図２は、２つの粒子が衝突する様子を時刻毎に示した図である。ここでは２つの粒子が相対速度の方向に正面衝突する場合を例にとって考える。粒子ｊに対する粒子ｉの相対速度をｖｉｊとする。時刻ｔ０において、相対速度をｖｉｊにて粒子ｉと粒子ｊとが接触する。時刻ｔ０における粒子ｉおよび粒子ｊの中心（質量中心）間の距離は粒子の半径の２倍である。それから粒子ｉと粒子ｊとが衝突により弾性変形して粒子の中心がさらに近づく（時刻ｔ１）。そして、反発力によって互いの粒子ｉｊが相対速度ｖｉｊ’で離れていく（時刻ｔ２）。接触力Ｆｃは、２つの粒子ｉｊが接触して（時刻ｔ０）から２つの粒子ｉｊが離れ始める（時刻ｔ２）まで粒子に働く。そして、この線形ばね系の振動周期Ｔは次式で表される。ここで、ｍは粒子ｉ（粒子ｊ）の質量である。すなわち、半周期のＴ／２が、粒子同士が接触している時間（接触時間Ｔｃ）になる。 なお、本実施形態では、付着力Ｆａは、粒子同士が接触している間、一定の力で粒子ｉｊが引き合うように作用すると仮定するモデルを採用する。なお、付着力Ｆａが、粒子ｉと粒子ｊとの接触面積に応じて変化するモデルを用いてもよい。 流動層の離散要素法シミュレーションにおいて、接触時間Ｔｃを５分割以上に分割した時間刻みでシミュレーションを行えば、安定かつ妥当な計算結果が得られることが経験的に知られている（非特許文献１参照）。よって、妥当な時間ステップΔｔの範囲は、次式で表される。上式を見れば分かるように、時間ステップΔｔの上限は、粒子の質量ｍとばね定数ｋとに依存し、粒子の質量ｍが小さいまたは粒子が硬い（ばね定数ｋが大きい）場合、小さくなり、計算時間が長くなる。 上述したように、付着力を考慮しない流動層の離散要素法シミュレーションにおいて、実際の粒子の反発を表すばね定数より小さい、修正ばね定数を用いて計算を行っても、問題ないことが知られている。例えば、硬い粒子の衝突において、粒子の弾性変形量が粒子径の数％以内（例えば１％以内）に収まるように修正ばね定数を設定して計算しても、妥当なシミュレーション結果が得られることが経験的に知られている。衝突の際の相対速度には幅があるが、衝突の平均相対速度の場合に、上記条件が満たされるようにしてもよいし、例えば、９５％の衝突において、上記条件が満たされるようにしてもよい。また、気体の吹き上げがないときに、堆積した粒子の重さによって一番下の粒子の弾性変形量が粒子径の数％以内に収まるように、修正ばね定数を設定してもよい。このような条件で修正ばね定数を設定する場合、硬い粒子であれば、本来のばね定数ｋｒの０．１倍以下の値に修正ばね定数ｋｄを設定できる。例えば、修正ばね定数ｋｄが本来のばね定数の０．１倍であれば、時間ステップを約３．１６倍（１０1/2倍）に設定することができる。 ただし付着力を考慮する場合、単に本来のばね定数より小さい修正ばね定数を用いて計算を行うと、妥当な結果が得られない。以下に、付着力の影響について説明する。 図３は、２つの粒子が衝突した後の２つの状態を示す図である。付着力と接触力とが働く２つの粒子ｉｊが衝突する場合、衝突時の相対速度ｖｉｊがある速度ｖｃ以下の場合、付着力が勝り、衝突後は２つの粒子ｉｊが付着した状態が維持される（状態Ａ）。相対速度ｖｉｊがある速度ｖｃより大きい場合、接触力による反発力が勝り、衝突後は２つの粒子ｉｊが速度Ｖｉｊ’で離れていく（状態Ｂ）。この閾値となる速度ｖｃを、求めることができる。 粒子同士が接触した場合の、接触力と付着力とを考慮した粒子ｉの運動方程式は、次式で表される。なお、流体力および重力は省略する。ここで、ｋ＞０、η＞０、Ｆａ＞０である。ｘは法線方向の粒子ｉの相対位置を表す。初期条件（ｔ＝０）で、ｘ＝０、衝突時の粒子ｉｊの相対速度（初期速度）をｖ0として上記運動方程式を解くと、次式が得られる。ここで、ｅpは粒子同士の衝突における反発係数である。 衝突した後のｘがｘ＜０（粒子が離れた状態）となりうる、初期速度ｖ0の下限ｖｃは、次式で表される。すなわち、初期速度ｖ0が閾値となる速度ｖｃ以下であれば、衝突後の粒子は付着したままとなり（状態Ａ）、初期速度ｖ0が閾値となる速度ｖｃより小さければ、衝突後の粒子は離れていく（状態Ｂ）。 ここで、速度ｖｃは、ばね定数ｋの−１／２乗に比例していると同時に、付着力Ｆａの１乗に比例している。そのため、付着力を考慮する場合、計算負荷を軽減するためにばね定数のみを小さく設定して計算を行うと、速度ｖｃの値が変化してしまう。速度ｖｃは、衝突した粒子同士が付着するか否かを分ける閾値となるので、粒子挙動が変化してしまう。そのため、従来のシミュレーションでは、妥当なシミュレーション結果が得られなかった。 本実施の形態では、この問題を回避するために、ばね定数として本来のばね定数ｋｒよりも小さい修正ばね定数ｋｄを用いると同時に、付着力として本来の付着力Ｆａｒよりも小さい修正付着力Ｆａｄを用いて、シミュレーションを行う。修正付着力Ｆａｄは、粒子の挙動を変えないために、閾値となる速度ｖｃが変化しないように選ぶのが好ましい。すなわち、本来のばね定数ｋｒと修正ばね定数ｋｄとの比の平方根と、本来の付着力Ｆａｒと修正付着力Ｆａｄとの比とが、ほぼ同じになるように修正付着力Ｆａｄを設定する。この関係を式に表せば以下のようになる。上式の関係を満たすような修正ばね定数ｋｄおよび修正付着力Ｆａｄを用いれば、閾値となる速度ｖｃが変化しない、すなわち、粒子の付着に関する挙動が変化しないので、妥当なシミュレーション結果を得ることができる。よって、計算に用いるばね定数を小さくすると共に、シミュレーションにおける時間ステップΔｔを大きくすることができるので、計算時間を短縮することができる。 なお、修正ばね定数ｋｄと修正付着力Ｆａｄを用いて粒子挙動を計算する本実施の形態のシミュレーション方法は、本来の付着力Ｆａｒが重力Ｆｇの０．１倍以上である場合に、好適に用いることができる。 ＜粒子挙動シミュレーション装置の構成＞ 図４は、本実施の形態の粒子挙動シミュレーション装置１０の機能的構成を示すブロック図である。粒子挙動シミュレーション装置１０は、条件入力部１１、挙動計算部１２、および接触力計算部（接触力特定部）１３を備える。 条件入力部１１は、修正ばね定数取得部１４、修正付着力取得部１５、および時間ステップ取得部１６を備える。条件入力部１１は、シミュレーションに必要な各種の条件が記載された設定ファイルを読み込み、シミュレーションの条件の入力を受け付ける。シミュレーションの条件としては、容器の形状、粒子数、粒子径、粒子質量、粒子同士の衝突における修正ばね定数、粒子同士の衝突における修正付着力、粒子同士の衝突における反発係数、分散板の下側から流入する気体の流速（流量）、気体の密度、気体の粘性係数、および時間ステップ等がある。シミュレーションの条件は他にもあるが省略する。なお、条件入力部１１は、キーボード等の入力装置を備えて、ユーザから直接シミュレーションの条件の入力を受け付けてもよい。 修正ばね定数取得部１４は、入力されたシミュレーション条件のうち、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、修正ばね定数ｋｄおよび粘性減衰係数ηを接触力計算部１３に出力する。修正ばね定数ｋｄは、本来のばね定数ｋｒより小さい値である。 修正付着力取得部１５は、入力されたシミュレーション条件のうち、修正付着力Ｆａｄを取得し、修正付着力Ｆａｄを挙動計算部１２に出力する。修正付着力Ｆａｄは、本来のばね定数ｋｒと修正ばね定数ｋｄとの比（の平方根）に応じて本来の付着力Ｆａｒより値を小さくしたものである。 時間ステップ取得部１６は、入力されたシミュレーション条件のうち、時間ステップｔｄを取得し、時間ステップｔｄを挙動計算部１２に出力する。時間ステップｔｄは、小さくした修正ばね定数ｋｄに応じて値を大きくした時間ステップであり、修正時間ステップと呼ぶことができるものである。修正時間ステップｔｄは、本来のばね定数ｋｒと修正ばね定数ｋｄとの比の平方根に応じて値を大きく設定できる。安定かつ妥当な結果を得るためには、衝突における粒子の接触時間Ｔｃの１／５以下の時間ステップで計算を行うという条件を満たせばよい。小さい修正ばね定数ｋｄを用いる場合、接触時間Ｔｃは長くなる。よって、計算時間を短縮し、安定かつ妥当な結果を得るための最大限の修正時間ステップｔｄは、次式で表される。 条件入力部１１は、取得したシミュレーションの他の条件についても挙動計算部１２に出力する。 挙動計算部１２は、各シミュレーション条件を用いて、修正時間ステップｔｄ毎に、各時刻における粒子の挙動（位置および運動）の計算を行う。挙動計算部１２は、ある時刻ｔにおいて、各粒子に働く力を求める。ここで、着目する１つの粒子に働く力Ｆは、次式で表される。挙動計算部１２は、着目する粒子が他の粒子（または壁）と接触している（着目する粒子と他の粒子との距離が所定の距離以内である）場合、着目する粒子に働く接触力を求めるために、着目する粒子と他の粒子（または壁）との距離および相対速度を接触力計算部１３に出力する。 接触力計算部１３は、着目する粒子と他の粒子との距離および相対速度と、修正ばね定数ｋｄと、粘性減衰係数ηとを用いて、着目する粒子に働く接触力（修正接触力）Ｆｃｄを求める。法線方向の修正接触力Ｆｃｄは次式で表される。接触力計算部１３は、算出した着目する粒子に働く修正接触力Ｆｃｄを、挙動計算部１２に出力する。 挙動計算部１２は、接触力計算部１３から受け取った着目する粒子に働く修正接触力Ｆｃｄ、および、修正付着力Ｆａｄ、着目する粒子に働く流体力Ｆｄ、重力Ｆｇを用いて、着目する粒子の修正時間ステップｔｄ後の位置および速度を計算する。挙動計算部１２は、同様に、全ての粒子について、修正時間ステップｔｄ後の位置および速度を計算する。なお、挙動計算部１２は、粒子の運動状態として粒子の自転角速度をさらに求めてもよい。 挙動計算部１２は、このようにして、修正時間ステップｔｄ毎の各時刻における各粒子の位置および運動を求める。挙動計算部１２は、各時刻の各粒子の位置および運動の情報を、記憶装置（図示せず）に出力する。例えば、粒子挙動シミュレーション装置１０は、各粒子の位置および運動を、時刻毎に順次表示装置（図示せず）に点や丸で可視化して表示させ、粒子の挙動をユーザに提示してもよい。 ＜粒子挙動シミュレーションのフロー＞ 次に、粒子挙動シミュレーションのパラメータの決定方法および粒子挙動シミュレーション装置１０における処理のフローについて説明する。 （ばね定数ｋｒの特定） 本実施の形態の粒子挙動シミュレーションを行うために、修正ばね定数ｋｄおよび修正付着力Ｆａｄが必要になる。そのために、本来のばね定数ｋｒおよび本来の付着力Ｆａｒが必要になる。 粒子同士が接触しているときに粒子に働く接触力を、線形ばねによる弾性反発力モデルで表した場合のばね定数ｋｒを、以下の手順で求める。 ヘルツの接触理論によれば、材質と大きさが等しい２つの球（粒子）が衝突するときの接触時間Ｔｈは、次式によって表される。ここで、ρは粒子密度、νは粒子のポアソン比、Ｅは粒子のヤング率、Ｒは粒子半径である。またｖ0は、粒子が衝突するときの相対速度である。ｖ0としては、粒子同士が衝突するときの相対速度の代表的な値を用いる。例えば、仮のばね定数を用いて付着力を考慮しない線形ばねモデルによって粒子挙動シミュレーションを行い、粒子同士の衝突時の平均相対速度を求め、その平均相対速度をｖ0として用いることができる。 一方、線形ばねとエネルギー減衰とによる反発力モデルにおける接触時間Ｔｌは、次式で表される。本実施の形態では、ヘルツの接触理論から求められる接触時間Ｔｈと、線形ばねとエネルギー減衰とによる反発力モデルにおける接触時間Ｔｌとが同じになるような、線形ばねのばね定数ｋｒを特定する。なお、ばね定数ｋｒを、実際の測定値に基づいて特定（仮定）してもよい。 なお、粘性減衰係数ηは、反発係数の実測値から特定する。 （付着力Ｆａｒの特定） 付着力Ｆａｒは、実測によって特定することができる。例えば、２つの粒子を押しつけて付着させてから、それを引きはがす際に必要な力を付着力とすることができる。なお、付着力Ｆａｒを、Hamaker理論等に基づいて特定（仮定）してもよい。 （修正ばね定数ｋｄの特定） 次に、ばね定数ｋｒより小さい修正ばね定数ｋｄを特定（設定）する。上述したように、例えば、硬い粒子の衝突において、粒子の弾性変形量が粒子径の数％以内（例えば３％以内または１％以内）に収まるように修正ばね定数を設定して計算しても、妥当なシミュレーション結果が得られることが経験的に知られている。本実施の形態では、主な衝突による粒子の弾性変形量（図２における時刻ｔ１の距離ｘ（粒子の重なりの距離））が粒子径の１％以内になるように、修正ばね定数ｋｄを設定する。また、例えば、付着力を考慮しない線形ばねモデルによって粒子挙動シミュレーションを行ったときの衝突時の平均相対速度から、弾性変形量が所定の範囲内に収まるような、修正ばね定数ｋｄを設定してもよい。 （修正付着力Ｆａｄの特定） 次に、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比（の平方根）に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを特定する。修正付着力Ｆａｄは、式（９）によって求めることができる。 （修正時間ステップｔｄの特定） 修正ばね定数ｋｄを用いて、計算に用いる修正時間ステップｔｄを求める。本実施の形態では、衝突における粒子の接触時間Ｔｃの１／５以下を、修正時間ステップｔｄとする。よって、修正時間ステップｔｄは、式（１０）によって決定される。 （粒子挙動の計算） 図５は、粒子挙動シミュレーション装置１０の処理を示すフロー図である。 粒子挙動シミュレーション装置１０の条件入力部１１は、上記のようにして特定した修正ばね定数ｋｄ、修正付着力Ｆａｄおよび修正時間ステップｔｄを含むシミュレーションの条件を、取得する（Ｓ１）。 修正ばね定数取得部１４は、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、接触力計算部１３に出力する。 修正付着力取得部１５は、修正付着力Ｆａｄを取得し、挙動計算部１２に出力する。 時間ステップ取得部１６は、時間ステップｔｄを取得し、挙動計算部１２に出力する。 着目する粒子が他の粒子と接触している場合（Ｓ２でＹｅｓ）、挙動計算部１２は、着目する粒子と他の粒子との距離および相対速度を接触力計算部１３に出力する（Ｓ３）。 接触力計算部１３は、着目する粒子と他の粒子との距離および相対速度と、修正ばね定数ｋｄと、粘性減衰係数ηとを用いて、式（１２）から、着目する粒子に働く接触力（修正接触力）Ｆｃｄを求める（Ｓ４）。接触力計算部１３は、着目する粒子に働く修正接触力Ｆｃｄを挙動計算部１２に出力する。 着目する粒子が他の粒子と接触していない場合（Ｓ２でＮｏ）、接触力および付着力は働かない（Ｆｃｄ＝Ｆａｄ＝０）。 Ｓ２でＮｏ、またはＳ４の後、挙動計算部１２は、修正接触力Ｆｃｄ、着目する粒子に働く修正付着力Ｆａｄ、流体力Ｆｄ、および重力Ｆｇを用いて、修正時間ステップｔｄ後における着目する粒子の位置および速度を求める（Ｓ５）。 修正時間ステップｔｄ後における位置および速度を求めていない粒子が未だある場合（Ｓ６でＮｏ）、残りの粒子に着目してＳ２〜Ｓ５の処理を繰り返す。 ある時刻について全ての粒子について修正時間ステップｔｄ後における位置および速度を求めた場合（Ｓ６でＹｅｓ）、所定の終了時刻（時間刻み）まで計算を行ったか否かを判定する（Ｓ７）。ここで、１秒間の粒子の挙動のシミュレーションを行う場合、終了時刻は最初の時刻から１秒後である。 所定の終了時刻まで計算していない場合（Ｓ７でＮｏ）、修正時間ステップｔｄだけ時刻を進め、次の時刻における計算（Ｓ２〜Ｓ６）を繰り返す。 所定の終了時刻まで計算が完了した場合（Ｓ７でＹｅｓ）、処理を終了する。 （実施の形態１のまとめ） 本実施の形態によれば、本来のばね定数ｋｒより小さい修正ばね定数ｋｄと、それに応じて本来の付着力Ｆａｒより小さくした修正付着力Ｆａｄとを用いて、離散要素法によって粒子の挙動を計算する。よって、小さい修正ばね定数ｋｄに応じて時間ステップを大きくすることができる。そのため、付着力を考慮する場合においても、シミュレーションの計算時間を短くし、かつ、妥当なシミュレーション結果を得ることができる。 なお、平均相対速度ｖ0、ばね定数ｋｒ、および修正ばね定数ｋｄを特定する際に、付着力を考慮しない線形ばねモデルによる粒子挙動シミュレーションを利用しているが、最初に一度利用するだけなので長い計算時間は問題にならない。上記従来の粒子挙動シミュレーションによって一度、平均相対速度ｖ0、ばね定数ｋｒ、および修正ばね定数ｋｄ等を決定した後は、類似の条件のシミュレーションにおいて、それらの値を何度も利用することができる。例えば、流動層装置を設計する場合、最適な条件を模索するために、処理容器の形状や流速等を変化させて何度も粒子挙動のシミュレーションを行う必要がある。一度修正ばね定数ｋｄ等を求めた後は、本実施の形態のシミュレーションを用いて、短時間で多くの条件のシミュレーションを行うことができる。 また、修正付着力Ｆａｄは、修正ばね定数ｋｄに応じて小さくなっていればよく、厳密に式（９）の関係を満たさなくても、ある程度の幅以内であれば、妥当な結果を得ることができる。例えば、修正付着力Ｆａｄを、次式の範囲内で設定してもよい。または、修正付着力Ｆａｄと付着力Ｆａｒとの比が、修正ばね定数ｋｄとばね定数ｋｒとの比の平方根と同じオーダーとなるように、修正付着力Ｆａｄを設定してもよい。 なお、本実施の形態では、修正接触力のモデルとして線形ばねモデルを仮定し、その線形ばね定数を特定しているが、これに限らず、ヘルツの接触理論による非線形ばねモデルを仮定し、その反発力を表す非線形ばね定数を用いてもよい。 また、上記は粒子同士の衝突について説明したが、粒子と壁との衝突について、別途修正ばね定数、反発係数、修正付着力を求めて、計算に用いてもよい。 ［実施の形態２］ 以下、本実施の形態について、図６を参照して詳細に説明する。なお、説明の便宜上、実施の形態１にて説明した図面と同じ機能を有する部材・構成については、同じ符号を付記し、その詳細な説明を省略する。本実施の形態では、図１に示す流動層装置における粒子をシミュレーションの対象とする。 ＜粒子挙動シミュレーション装置の構成＞ 図６は、本実施の形態の粒子挙動シミュレーション装置２０の機能的構成を示すブロック図である。粒子挙動シミュレーション装置２０は、条件入力部２１、挙動計算部１２、接触力計算部（接触力特定部）１３、および修正付着力特定部２２を備える。 条件入力部２１は、ばね定数取得部２３、付着力取得部２４、修正ばね定数取得部２５、および時間ステップ取得部１６を備える。条件入力部２１は、シミュレーションに必要な各種の条件が記載された設定ファイルを読み込み、シミュレーションの条件の入力を受け付ける。 ばね定数取得部２３は、入力されたシミュレーション条件のうち、本来のばね定数ｋｒを取得し、ばね定数ｋｒを修正付着力特定部２２に出力する。 付着力取得部２４は、入力されたシミュレーション条件のうち、本来の付着力Ｆａｒを取得し、付着力Ｆａｒを修正付着力特定部２２に出力する。 修正ばね定数取得部２５は、入力されたシミュレーション条件のうち、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、修正ばね定数ｋｄおよび粘性減衰係数ηを接触力計算部１３に出力する。また、修正ばね定数取得部２５は、修正ばね定数ｋｄを修正付着力特定部２２に出力する。 修正付着力特定部２２は、修正付着力Ｆａｄと付着力Ｆａｒとの比が、修正ばね定数ｋｄとばね定数ｋｒとの比の平方根と同じになるような、修正付着力Ｆａｄを特定する。すなわち、修正付着力特定部２２は、式（９）に従って修正付着力Ｆａｄを特定する。修正付着力特定部２２は、特定した修正付着力Ｆａｄを、挙動計算部１２に出力する。 条件入力部２１は、取得したシミュレーションの他の条件についても挙動計算部１２に出力する。 挙動計算部１２および接触力計算部１３の構成は、実施の形態１と同様であるので、その説明を省略する。 本実施の形態では、本来のばね定数ｋｒ、本来の付着力Ｆａｒ、および修正ばね定数ｋｄに基づき、修正付着力特定部２２が適切な修正付着力Ｆａｄを求める。よって、付着力を考慮する場合において、短い計算時間で、妥当なシミュレーション結果を得ることができる。 本発明の粒子挙動シミュレーションに関する実施例について、以下に説明する。 ＜実験例＞ まず、対比に用いる実験結果について説明する。流動層装置として、図１に示す構成のもので、処理容器内側の幅は１８ｍｍ、奥行きは３．６ｍｍのものを用いた。処理容器は透明のガラス製である。粒子として、平均径６０μｍの球形ガラスビーズ（日本粉体技術協会ＪＩＳ試験用粉体２（ＧＢＬ−６０））を用いる。初期の（気流のない状態での）粒子層の堆積高さは２０ｍｍである。シリンジポンプ（有限会社ミナトコンセプト製）を用いて、流動層内での流速が０．０２ｍ／ｓになるように、帆布製の分散板の下から空気を供給した。 図７は、処理容器の背面から照明を当て、流動層状態になったときの粒子層の状態を撮影した画像である。粒子層の中に気泡ができ、流動層状態になっていることが分かる。 ＜シミュレーション例＞ 本実施例のシミュレーションの条件について、説明する。想定する流動層装置は、実験例と同じ、処理容器内側の幅は１８ｍｍのものである。ただし、奥行きが粒子１層からなる２次元流動層として計算を行っている。粒子径は６０μｍで、その物性としてガラスの物性を仮定し、粒子のヤング率Ｅ＝７１．３ＧＰａ、ポアソン比ν＝０．２２とした。また、粒子の密度を２２５０ｋｇ／ｍ3、反発係数を０．９４、摩擦係数を０．１８とした。粒子数は１０万個である。初期の粒子層の堆積高さは２０ｍｍである。分散板の下から流入する気体の流速を、０．０２ｍ／ｓ、気体の密度を１．２０５ｋｇ／ｍ3、気体の粘性率を１．８１×１０-5Ｐａ・ｓとした。 上記条件から、ヘルツの接触理論による接触時間と線形ばねモデルによる接触時間が同じになるばね定数ｋｒを求めた。上記条件と式（１３）および式（１４）とから、ばね定数ｋｒを約７０００Ｎ／ｍと特定した。 付着力Ｆａｒは、微小力測定装置（ＮＳ−Ａ１００：株式会社ナノシーズ製）を用いて測定した。その結果、粒子同士の間には、平均して１０７ｎＮの付着力が作用することが分かった。これは自重Ｆｇの約４０倍である。なお付着力の測定値の変動幅は２７ｎＮ以内であり、測定された付着力は、付着時に粒子にかけた押し込み荷重には依存していなかった。 （比較例１） 比較例１では、ばね定数ｋｒ＝７０００Ｎ／ｍと、付着力Ｆａｒ＝４０Ｆｇ（重力の４０倍）とを用いて、従来の離散要素法によって計算を行った。図８は、比較例１のシミュレーション結果を示す画像である。比較例１は、気泡が形成される流動層状態を再現しており、図７に示す実験結果を良好に再現している。ただし、比較例１において、用いた時間ステップは１．０×１０-7秒であり、１秒間のシミュレーションの計算に要した時間は、４２０時間であった（ＣＰＵ：Intel Xeon（登録商標） 3.4GHz、メモリ：4GB のＰＣを用いて計算）。 （比較例２） 比較例２では、修正ばね定数ｋｄ＝１０Ｎ／ｍと、付着力Ｆａｒ＝４０Ｆｇとを用いて、従来の離散要素法によって計算を行った。すなわち、本来のばね定数ｋｒの１．４×１０-3倍の値を修正ばね定数として用いている。図９は、比較例２のシミュレーション結果を示す画像である。比較例２では、気泡の代わりに亀裂状の空隙の形成が見られ、粒子の挙動が大きく変化してしまっている。 （実施例１） 実施例１では、修正ばね定数ｋｄ＝１０Ｎ／ｍと、修正付着力Ｆａｄ＝１．５１Ｆｇとを用いて、本発明の離散要素法によって計算を行った。修正付着力Ｆａｄは、式（９）に従って求めた。図１０は、実施例１のシミュレーション結果を示す画像である。実施例１は、気泡が形成される流動層状態を再現しており、図７に示す実験結果を良好に再現している。さらに、実施例１において、用いた時間ステップは２．０×１０-6秒であり、１秒間のシミュレーションの計算に要した時間は、６０時間であった（計算機は比較例１と同条件）。すなわち、修正ばね定数ｋｄと修正付着力Ｆａｄとを用いることで、妥当な結果を得るための計算時間を約１／７に短縮することができた。なお、離散要素法においては、粒子に関する計算と流体に関する計算とを行っている。そのため、時間ステップと全体の計算時間とは、必ずしも比例しない。 最後に、粒子挙動シミュレーション装置１０・２０の各ブロック、特に条件入力部１１、挙動計算部１２、接触力計算部１３、修正ばね定数取得部１４、修正付着力取得部１５、時間ステップ取得部１６、条件入力部２１、修正付着力特定部２２、ばね定数取得部２３、付着力取得部２４、および修正ばね定数取得部２５は、ハードウェアロジックによって構成してもよいし、次のようにＣＰＵ（central processing unit）を用いてソフトウェアによって実現してもよい。 すなわち、粒子挙動シミュレーション装置１０・２０は、各機能を実現する制御プログラムの命令を実行するＣＰＵ、上記プログラムを格納したＲＯＭ（read only memory）、上記プログラムを展開するＲＡＭ（random access memory）、上記プログラムおよび各種データを格納するメモリ等の記憶装置（記録媒体）などを備えている。そして、本発明の目的は、上述した機能を実現するソフトウェアである粒子挙動シミュレーション装置１０・２０の制御プログラムのプログラムコード（実行形式プログラム、中間コードプログラム、ソースプログラム）をコンピュータで読み取り可能に記録した記録媒体を、粒子挙動シミュレーション装置１０・２０に供給し、そのコンピュータ（またはＣＰＵやＭＰＵ（microprocessor unit））が記録媒体に記録されているプログラムコードを読み出し実行することによっても、達成可能である。 上記記録媒体としては、例えば、磁気テープやカセットテープ等のテープ系、フロッピー（登録商標）ディスク／ハードディスク等の磁気ディスクやＣＤ−ＲＯＭ（compact disc read-only memory）／ＭＯ（magneto-optical）／ＭＤ（Mini Disc）／ＤＶＤ（digital versatile disk）／ＣＤ−Ｒ（CD Recordable）等の光ディスクを含むディスク系、ＩＣカード（メモリカードを含む）／光カード等のカード系、あるいはマスクＲＯＭ／ＥＰＲＯＭ（erasable programmable read-only memory）／ＥＥＰＲＯＭ（electrically erasable and programmable read-only memory）／フラッシュＲＯＭ等の半導体メモリ系などを用いることができる。 また、粒子挙動シミュレーション装置１０・２０を通信ネットワークと接続可能に構成し、上記プログラムコードを通信ネットワークを介して供給してもよい。この通信ネットワークとしては、特に限定されず、例えば、インターネット、イントラネット、エキストラネット、ＬＡＮ（local area network）、ＩＳＤＮ（integrated services digital network）、ＶＡＮ（value-added network）、ＣＡＴＶ（community antenna television）通信網、仮想専用網（virtual private network）、電話回線網、移動体通信網、衛星通信網等が利用可能である。また、通信ネットワークを構成する伝送媒体としては、特に限定されず、例えば、ＩＥＥＥ（institute of electrical and electronic engineers）１３９４、ＵＳＢ、電力線搬送、ケーブルＴＶ回線、電話線、ＡＤＳＬ（asynchronous digital subscriber loop）回線等の有線でも、ＩｒＤＡ（infrared data association）やリモコンのような赤外線、Ｂｌｕｅｔｏｏｔｈ（登録商標）、８０２．１１無線、ＨＤＲ（high data rate）、携帯電話網、衛星回線、地上波デジタル網等の無線でも利用可能である。 本発明は上述した各実施形態に限定されるものではなく、請求項に示した範囲で種々の変更が可能であり、異なる実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても本発明の技術的範囲に含まれる。 本発明は、粒子の挙動を予測する粒子挙動シミュレーションに利用することができる。 １ 流動層装置 ２ 処理容器 ３ 分散板 ４ 気体供給口 ５ 排気口 ６ 触媒粒子１０、２０ 粒子挙動シミュレーション装置１１、２１ 条件入力部１２ 挙動計算部１３ 接触力計算部（接触力特定部）１４、２５ 修正ばね定数取得部１５ 修正付着力取得部１６ 時間ステップ取得部２２ 修正付着力特定部２３ ばね定数取得部２４ 付着力取得部 離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、 粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、 上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定部と、 接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴とする粒子挙動シミュレーション装置。 上記修正付着力Ｆａｄは、条件としてを満たすことを特徴とする請求項１に記載の粒子挙動シミュレーション装置。 上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比は、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根と同じオーダーであることを特徴とする請求項１に記載の粒子挙動シミュレーション装置。 上記修正ばね定数ｋｄは、上記ばね定数ｋｒの０．１倍以下であることを特徴とする請求項１から３のいずれか一項に記載の粒子挙動シミュレーション装置。 上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比は、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根とほぼ同じであることを特徴とする請求項１に記載の粒子挙動シミュレーション装置。 上記付着力Ｆａｒと上記修正ばね定数ｋｄと上記ばね定数ｋｒとに基づき、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の平方根と同じになるような上記修正付着力Ｆａｄを特定する修正付着力特定部を備えることを特徴とする請求項１に記載の粒子挙動シミュレーション装置。 離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション方法であって、 粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、 シミュレーション装置の接触力特定部が、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定ステップと、 上記シミュレーション装置の挙動計算部が、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算ステップとを含むことを特徴とする粒子挙動シミュレーション方法。 請求項１から６のいずれか一項に記載の粒子挙動シミュレーション装置としてコンピュータを機能させるための制御プログラムであって、コンピュータを上記の各部として機能させるための制御プログラム。 請求項８に記載の制御プログラムを記録したコンピュータ読み取り可能な記録媒体。 【課題】付着力の影響を考慮しながら、離散要素法を用いて複数の粒子の挙動を計算する時間を短縮することができる粒子挙動シミュレーションを実現する。【解決手段】本発明に係る粒子挙動シミュレーション装置１０は、粒子同士が接触しているときに上記粒子に働く接触力をばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記ばねモデルに基づいた修正接触力を求める接触力特定部１３と、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって複数の粒子の挙動を計算する挙動計算部１２とを備える。【選択図】図４