(57)【要約】【課題】付着力の影響を考慮しながら、計算時間を短縮することができる粒子挙動シミュレーションを実現する。【解決手段】本発明に係る粒子挙動シミュレーション装置１０は、接触力を非線形ばねモデルで表した場合のばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、修正接触力を求める接触力特定部１３と、接触力として上記修正接触力を用い、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、粒子の挙動を計算する挙動計算部１２とを備える。【選択図】図４
2013210918.tif 000002
【特許請求の範囲】
【請求項１】
離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、
粒子同士が接触しているときに上記粒子に働く接触力を非線形ばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、
上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記非線形ばねモデルに基づいた修正接触力を求める接触力特定部と、
接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴とする粒子挙動シミュレーション装置。
【請求項２】
上記非線形ばねモデルは、ヘルツ接触モデルであることを特徴とする請求項１に記載の粒子挙動シミュレーション装置。
【請求項３】
上記修正付着力Ｆａｄは、条件として、
【数１】
2013210918.tif 000034
を満たすことを特徴とする請求項１または２に記載の粒子挙動シミュレーション装置。
【請求項４】
上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比は、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗と同じオーダーであることを特徴とする請求項１または２に記載の粒子挙動シミュレーション装置。
【請求項５】
上記修正ばね定数ｋｄは、上記ばね定数ｋｒの０．１倍以下であることを特徴とする請求項１から４のいずれか一項に記載の粒子挙動シミュレーション装置。
【請求項６】
上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比は、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗とほぼ同じであることを特徴とする請求項１または２に記載の粒子挙動シミュレーション装置。
【請求項７】
上記付着力Ｆａｒと上記修正ばね定数ｋｄと上記ばね定数ｋｒとに基づき、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗と同じになるような上記修正付着力Ｆａｄを特定する修正付着力特定部を備えることを特徴とする請求項１または２に記載の粒子挙動シミュレーション装置。
【請求項８】
離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション方法であって、
粒子同士が接触しているときに上記粒子に働く接触力を非線形ばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、
シミュレーション装置の接触力特定部が、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記非線形ばねモデルに基づいた修正接触力を求める接触力特定ステップと、
上記シミュレーション装置の挙動計算部が、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を
計算する挙動計算ステップとを含むことを特徴とする粒子挙動シミュレーション方法。
【請求項９】
請求項１から７のいずれか一項に記載の粒子挙動シミュレーション装置としてコンピュータを機能させるための制御プログラムであって、コンピュータを上記の各部として機能させるための制御プログラム。
【請求項１０】
請求項９に記載の制御プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、多数の粒子の挙動を予測する粒子挙動シミュレーションに関する。
【背景技術】
【０００２】
化学合成プロセスにおいて、気体から反応生成物を得る触媒反応装置として、流動層装置が広く用いられている。流動層装置は、処理容器内に粉体等を入れておき、粉体の下部から処理容器内に気体を吹き込むことで、粉体を流動化させる装置である。例えば、流動層装置の処理容器内に触媒粒子を入れておき、処理容器の下部から原料ガスを供給し、流動化状態の触媒粒子と原料ガスとを効率よく接触させることによって、ガス状の反応生成物を得ることができる。
【０００３】
効率よく反応を起こさせる流動層装置を設計するため、または触媒粒子の量および気体の流量等のよりよい条件を見つけるために、流動層における粒子の挙動をシミュレーションすることが重要である。
【０００４】
このような流動層触媒反応に用いられる触媒粒子としては、平均粒径が３０〜１５０μｍで、密度が７００〜２５００ｋｇ／ｍ3程度の粒子、いわゆるGeldart線図においてＡ粒子に分類される粒子が広く用いられている。Geldart線図における粒子の分類は、流動層の粒子の流動化状態を予測する指標となる。Geldart線図においてＡ粒子に分類される粒子を触媒粒子として用いれば、流動層が小さい気泡を形成するので、反応効率を高くすることができる。
【０００５】
流動層におけるこのような粒子の挙動をシミュレーションする方法の１つに、離散要素法（個別要素法：ＤＥＭ：Discrete Element Method）がある。離散要素法を用いた数値解析では、分割された時間刻み毎に、各粒子に作用する力に基づいて各粒子に対して運動方程式を計算し、時刻を時間ステップだけ進めた次の時刻における粒子の状態を数値計算によって求める。そして、時間刻み毎の計算を繰り返し、粒子の挙動を求める。ある粒子に作用する力としては、他の粒子または処理容器の内壁との衝突による反発を表す接触力、粒子の周囲の気体（流体）から受ける流体力、他の粒子または処理容器の内壁との間に働く付着力、および重力等がある。
【０００６】
離散要素法において、シミュレーションの時間刻みの幅（時間ステップ）が小さければ、精度のよいシミュレーション結果を得ることができる。しかしながら、時間ステップを小さくすると、所定の時間（例えば１秒間）の粒子の挙動を求めるための計算量が多くなる。そのため、シミュレーションにかかる計算時間が長くなる。時間ステップを大きくすれば計算時間を短くすることができるが、一方で、計算を安定に行う（妥当なシミュレーション結果を得る）ための時間ステップには上限がある。
【０００７】
非特許文献１には、粒子の衝突の接触力を表す線形ばねモデルを仮定した場合、妥当なシミュレーション結果を得るための時間ステップの経験的な上限が、ばねの固有振動の周期Ｔの約１／１０であることが記載されている。この場合、半周期Ｔ／２が粒子同士の接触時間になるので、粒子同士の接触時間の約１／５が時間ステップの経験的な上限になる。そして、固有振動の周期Ｔはばね定数ｋに依存するので、ばね定数が大きな粒子、すなわち硬い粒子の場合、時間ステップを小さくする必要がある。そのため、ばね定数が小さい軟らかい粒子に比べて、硬い粒子のシミュレーションは、計算時間が長くなる。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】社団法人日本粉体工業技術協会、「流動層ハンドブック」、株式会社培風館、１９９９年３月２５日、p.125-144
【非特許文献２】田中敏嗣、石田俊哉、辻裕、「水平管内粒状体プラグ流の直接数値シミュレーション（付着力がない場合）」、日本機械学会論文集（Ｂ編）、５７巻、５３４号、１９９１年、p.456-463
【発明の概要】
【発明が解決しようとする課題】
【０００９】
線形ばねの代わりにヘルツ接触モデルを適用してシミュレーションを行う場合でも同様に、ばね定数が大きな硬い粒子の場合、粒子同士の接触時間が短くなるため、時間ステップを小さくする必要がある。
【００１０】
しかしながら、Geldart線図においてＡ粒子に分類される粒子のように、粒子径が小さい場合、ファンデルワールス力および液架橋力等に起因する粒子同士の付着力が大きくなり、付着力の影響が無視できなくなる。付着力とは、粒子が他の粒子または処理容器の内壁等と近接した場合に受ける引力である。付着力が大きく、付着力の影響が無視できない系において、本来のばね定数よりもばね定数を小さく設定して時間ステップを長くすると、シミュレーションの精度が悪くなり、妥当な結果が得られないという問題がある。そのため、付着力の影響を無視できない系では、シミュレーションの計算に長い時間がかかっていた。
【００１１】
本発明は、上記の問題点に鑑みてなされたものであり、その目的は、付着力の影響を考慮しながら、離散要素法を用いて複数の粒子の挙動を計算する時間を短縮することができる粒子挙動シミュレーションを実現することにある。
【課題を解決するための手段】
【００１２】
本発明に係る粒子挙動シミュレーション装置は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、上記の課題を解決するために、粒子同士が接触しているときに上記粒子に働く接触力を非線形ばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記非線形ばねモデルに基づいた修正接触力を求める接触力特定部と、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴としている。
【００１３】
本発明に係る粒子挙動シミュレーション方法は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション方法であって、粒子同士が接触しているときに上記粒子に働く接触力を非線形ばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、シミュレーション装置の接触力特定部が、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記非線形ばねモデルに基づいた修正接触力を求める接触力特定ステップと、上記シミュレーション装置の挙動計算部が、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算ステップとを含むことを特徴としている。
【００１４】
上記の構成によれば、ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて修正接触力を特定し、修正接触力を用いて離散要素法によって粒子の挙動を計算する。よって、粒子同
士の接触時間が長くなるので、時間ステップを長くし、挙動計算に要する時間を短縮することができる。また、挙動計算において、修正ばね定数ｋｄを用いると同時に、ばね定数ｋｒと修正ばね定数ｋｄとの比に応じて付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いる。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。
【００１５】
なお、粒子に働く接触力を表す非線形ばねモデルは、ヘルツ接触モデルでもよい。また、上記挙動計算部は、粒子同士の付着力を、粒子同士が所定の距離より近接しているときのみに粒子に働く力として扱ってもよい。
【００１６】
また、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗とほぼ同じになるようにしてもよい。
【００１７】
上記の構成によれば、小さい修正ばね定数ｋｄを用いて挙動を計算しても、粒子同士の衝突後に、２つの粒子が付着したままになるか、離れていくかの閾値となる速度ｖｃが変化しない。よって、小さい修正ばね定数ｋｄを用いて挙動を計算しても、粒子の挙動が変化しない。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。
【００１８】
また、上記修正付着力Ｆａｄは、条件として、後述の式（２３）を満たすようにしてもよい。
【００１９】
上記の構成によれば、修正付着力Ｆａｄの値は、修正ばね定数ｋｄをばね定数ｋｒで除したものの２／５乗に付着力Ｆａｒを乗じた値の、０．７倍から１．３倍までの範囲になる。
【００２０】
このように、修正付着力Ｆａｄが、修正ばね定数ｋｄをばね定数ｋｒで除したものの２／５乗に付着力Ｆａｒを乗じた値から少しずれてもよい。
【００２１】
また、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗と同じオーダーになるようにしてもよい。
【００２２】
また、上記修正ばね定数ｋｄは、上記ばね定数ｋｒの０．１倍以下としてもよい。
【００２３】
上記の構成によれば、修正ばね定数ｋｄに応じて時間ステップを長くして、挙動計算に要する時間を短縮することができる。
【００２４】
また、上記付着力Ｆａｒと上記修正ばね定数ｋｄと上記ばね定数ｋｒとに基づき、上記修正付着力Ｆａｄと上記付着力Ｆａｒとの比が、上記修正ばね定数ｋｄと上記ばね定数ｋｒとの比の２／５乗と同じになるような上記修正付着力Ｆａｄを特定する修正付着力特定部を備える構成であってもよい。
【００２５】
上記の構成によれば、修正付着力特定部が、小さい修正ばね定数ｋｄを用いて挙動を計算しても粒子の挙動が変化しない修正付着力Ｆａｄを特定する。そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。
【００２６】
なお、上記粒子挙動シミュレーション装置は、一部をコンピュータによって実現してもよく、この場合には、コンピュータを上記各部として動作させる制御プログラム、および上記制御プログラムを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に
入る。
【発明の効果】
【００２７】
本発明に係る粒子挙動シミュレーション装置は、離散要素法を用いて複数の粒子の挙動を予測する粒子挙動シミュレーション装置であって、粒子同士が接触しているときに上記粒子に働く接触力を非線形ばねモデルで表した場合のばね定数をｋｒとし、粒子同士の付着力をＦａｒとしたとき、上記ばね定数ｋｒより小さい修正ばね定数ｋｄを用いて、上記非線形ばねモデルに基づいた修正接触力を求める接触力特定部と、接触力として、上記修正接触力を用い、かつ、付着力として、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを用いて、離散要素法によって粒子の挙動を計算する挙動計算部とを備えることを特徴としている。
【００２８】
そのため、付着力を考慮すべき粒子系について、挙動計算に要する時間を短縮しながら、妥当なシミュレーション結果を得ることができる。
【図面の簡単な説明】
【００２９】
【図１】シミュレーションの対象とする流動層装置の概略構成を示す断面図である。
【図２】２つの粒子が衝突する様子を時刻毎に示した図である。
【図３】粒子が壁面に衝突した後の２つの状態を示す図である。
【図４】本発明の一実施形態の粒子挙動シミュレーション装置の機能的構成を示すブロック図である。
【図５】上記粒子挙動シミュレーション装置の処理を示すフロー図である。
【図６】本発明の別の実施形態の粒子挙動シミュレーション装置の機能的構成を示すブロック図である。
【図７】粘性減衰および付着力を考慮した運動方程式に基づいて一粒子の動きを数値計算して、付着した状態を維持する限界の初期相対速度ｖｃがある値になるときの、ばね定数ｋと付着力Ｆａとをプロットしたグラフである。
【図８】参考例として付着力を一定にしてばね定数のみを変化させた場合の、閾値となる速度ｖｃの値をプロットしたグラフである。
【図９】ばね定数と付着力とを変化させた場合の、閾値となる速度ｖｃの値をプロットしたグラフである。
【発明を実施するための形態】
【００３０】
以下の実施の形態では、主に、流動層装置における粒子の挙動をシミュレーションする場合について説明するが、本発明はこれに限定されるものではない。
【００３１】
［実施の形態１］
以下、本実施の形態について、図１〜図５、図７〜図９を参照して詳細に説明する。本実施の形態では、図１に示す流動層装置における粒子をシミュレーションの対象とする。
【００３２】
＜シミュレーション系＞
図１は、流動層装置１の概略構成を示す断面図である。流動層装置１は、処理容器２、分散板３、気体供給口４、および排気口５を備える。分散板３は、処理容器２内の下部に容器内を分割するように水平に配置されており、微細な貫通孔を有する。処理用に２内の分散板３の上には、触媒粒子６が堆積されている。分散板３は、気体を通過させるが触媒粒子６は通過させない。気体供給口４は、処理容器２の分散板３の下部に設けられている。原料ガスは、気体供給口４から処理容器２内に供給され、分散板３の下から堆積した粒子の隙間を通って噴き上げる。原料ガスおよびその反応物とは処理容器２の上部に設けられた排気口５から排出される。
【００３３】
＜シミュレーション方法＞
離散要素法では、着目する粒子に働く力から、次時刻における粒子の位置および運動状態を計算する。着目する粒子に働く力Ｆは、次式で与えられる。
【００３４】
【数１】
2013210918.tif 000003
【００３５】
ここで、Ｆｄは流体力、Ｆｇは重力、Ｆｃは接触力、Ｆａは付着力をそれぞれ表す。
【００３６】
流体力Ｆｄは、粒子が吹き上げる気体（原料ガス）から受ける、主に上向きの力である。重力Ｆｇは自重によって粒子に働く力である。接触力Ｆｃは、粒子同士の接触、または粒子と処理容器の内壁（または分散板）との接触による反発力を表す。付着力Ｆａは、粒子が他の粒子または処理容器の内壁と近接または接触した際に受ける、ファンデルワールス力および液架橋力等に起因する引力を表す。流体力Ｆｄおよび重力Ｆｇは、各粒子が常に受ける力である。ただし、流体力Ｆｄは、一定ではなく、位置および周囲の粒子の分布（空間密度）等によって変化する。一方、接触力Ｆｃおよび付着力Ｆａは、例えば粒子が他の粒子または内壁と接触したときにのみ働く。
【００３７】
粒子が他の粒子と接触した場合の接触力Ｆｃおよび付着力Ｆａについて、説明する。ここでは、粒子が接触した場合、粒子が弾性変形をして、反発力を受けるモデルを用いて接触力を表す。半径ｒの粒子同士が衝突した際に粒子が受ける接触力について、接触面の法線方向の力を、弾性反発力を表すばねと、粘性減衰力を表すダッシュポットによって表現する。粘性減衰係数は、粒子の反発係数から決定される。なお、接触面の接線方向の力は、粒子の変形量が小さい場合は、法線方向と同様のばねとダッシュポットの組み合わせによって表現され、粒子の変形量が大きい場合は、滑りを表すための摩擦スライダを用いて表現される。
【００３８】
ヘルツの接触理論を用いた非線形ばねモデル（ヘルツ接触モデル）における法線方向の接触力Ｆｃは、次式で表される。
【００３９】
【数２】
2013210918.tif 000004
【００４０】
ここで、ｋは粒子の弾性反発力を表す法線方向のヘルツ接触モデルの（非線形ばねの）ばね定数、ηはダッシュポットによる減衰を表す法線方向の粘性減衰係数、ｘは接触時を基準とした粒子の中心位置の変位、ｖは粒子の相対速度ベクトル、ｎは法線方向の単位ベクトルを表す。
【００４１】
図２は、２つの粒子が衝突する様子を時刻毎に示した図である。ここでは２つの粒子が相対速度の方向に正面衝突する場合を例にとって考える。粒子ｊに対する粒子ｉの相対速度をｖｉｊとする。時刻ｔ０において、相対速度をｖｉｊにて粒子ｉと粒子ｊとが接触する。時刻ｔ０における粒子ｉおよび粒子ｊの中心（質量中心）間の距離は粒子の半径の２倍である。それから粒子ｉと粒子ｊとが衝突により弾性変形して粒子の中心がさらに近づく（時刻ｔ１）。そして、反発力によって互いの粒子ｉｊが相対速度ｖｉｊ’で離れていく（時刻ｔ２）。接触力Ｆｃは、２つの粒子ｉｊが接触して（時刻ｔ０）から２つの粒子ｉｊが離れ始める（時刻ｔ２）まで粒子に働く。
【００４２】
なお、本実施形態では、付着力Ｆａは、粒子同士が接触している間、一定の力で粒子ｉ
ｊが引き合うように作用すると仮定するモデルを採用する。
【００４３】
付着力を考慮しない場合、本来のばね定数よりある程度小さい修正ばね定数を用いて計算を行っても、妥当なシミュレーション結果を得ることができる。流動層の計算結果に関して、ばね定数の影響は認められないことが確認されている（非特許文献１参照）。
【００４４】
しかしながら、付着力を考慮する場合、単に本来のばね定数より小さい修正ばね定数を用いて計算を行うと、妥当な結果が得られない。接触時間が長くなる分、付着力の力積が大きくなってしまい、粒子の挙動が変化してしまうためである。以下に、付着力の影響について説明する。
【００４５】
（付着力の影響）
図３は、粒子が壁面に衝突した後の２つの状態を示す図である。付着力と接触力とが働く粒子ｉが壁面に衝突する場合、衝突時の相対速度ｖｉｊがある速度ｖｃ以下の場合、付着力が勝り、衝突後は粒子ｉが壁面に付着した状態が維持される（状態Ａ）。相対速度ｖｉｊがある速度ｖｃより大きい場合、接触力による反発力が勝り、衝突後は粒子ｉが壁面から速度Ｖｉｊ’で離れていく（状態Ｂ）。速度ｖｃは、衝突した後の粒子の状態（挙動）を分ける閾値となる速度である。
【００４６】
このように、接触力および付着力が、粒子の反発・付着挙動にとって重要な因子である。そこで、接触している粒子と壁面との間に作用する接触力および付着力に着目すると、粒子が壁面と衝突するときの運動量の変化と力積との関係は次式で表される。
【００４７】
【数３】
2013210918.tif 000005
【００４８】
ここで、ｍは粒子の質量、ｖｉｊは衝突前の粒子の相対速度、ｖｉｊ’は衝突後の粒子の相対速度、Ｔｃは接触時間、Ｆｃは接触力、Ｆａは付着力を表す。左辺は衝突における粒子の運動量の変化を表し、右辺は衝突において粒子が受けた力積を表す。なお、接触力および付着力の力が働く方向は互いに逆向きである。
【００４９】
衝突後に粒子が付着した状態を維持する限界の初期相対速度をｖｃとすると、ｖｉｊ＝ｖｃのとき、ｖｉｊ’＝０となる。このとき、運動量と力積との関係は次式で表される。
【００５０】
【数４】
2013210918.tif 000006
【００５１】
さらに、粒子と壁面とが接触している間に作用する付着力Ｆａが一定であると仮定すると、次式が得られる。
【００５２】
【数５】
2013210918.tif 000007
【００５３】
一方、ヘルツ接触モデルを適用する場合の、接触力Ｆｃおよび付着力Ｆａを考慮した衝
突における粒子ｉの運動方程式は、以下で表される。
【００５４】
【数６】
2013210918.tif 000008
【００５５】
ここで、ｋは粒子の弾性反発力を表すヘルツ接触モデルの（非線形ばねの）ばね定数、ηはダッシュポットによる減衰を表す法線方向の粘性減衰係数、ｘは接触時を基準とした粒子の中心位置の変位である。ここで、ｋ＞０、η＞０、Ｆａ＞０である。付着力Ｆａまで考慮した上式から、例えば接触時間Ｔｃを解析的に求めるのは困難である。
【００５６】
そこでまず、粘性減衰および付着力を考慮しない場合の粒子の運動を考える。粘性減衰および付着力を考慮しない場合の衝突における粒子ｉの運動方程式は、次式で表される。
【００５７】
【数７】
2013210918.tif 000009
【００５８】
上式の力が働くとき、ヘルツの接触理論によれば、材質および大きさが等しい２つの球（粒子）が相対速度ｖ0で衝突するときの接触時間Ｔｃは、次式によって表される。
【００５９】
【数８】
2013210918.tif 000010
【００６０】
ここで、ρは粒子密度、νは粒子のポアソン比、Ｅは粒子のヤング率、Ｒは粒子半径である。また、衝突の初期条件（ｔ＝０）における、粒子ｉの位置ｘを０、粒子ｉｊの相対速度をｖ0とした。上式の定数を省略すると、接触時間Ｔｃを次式のように表すことができる。
【００６１】
【数９】
2013210918.tif 000011
【００６２】
一方、同じ材質の球と壁面とが接触する際の法線方向のばね定数ｋは、次式で表される。
【００６３】
【数１０】
2013210918.tif 000012
【００６４】
また、参考までに材質および大きさが等しい２つの球（粒子）が接触する際の法線方向のばね定数ｋは、次式で表される。
【００６５】
【数１１】
2013210918.tif 000013
【００６６】
式（１０）または式（１１）の定数を省略すると、いずれの場合でも、ヘルツ接触モデルによるばね定数ｋを次式のように表すことができる。
【００６７】
【数１２】
2013210918.tif 000014
【００６８】
上式と式（９）とにより、接触時間Ｔｃとばね定数ｋとの関係は次式で表される。
【００６９】
【数１３】
2013210918.tif 000015
【００７０】
すなわち、接触時間Ｔｃは、ヘルツ接触モデルのばね定数ｋの−２／５乗に比例していることが分かる。
【００７１】
一方、ヘルツの接触理論によれば、材質および大きさが等しい２つの球（粒子）が接触する際の最大オーバーラップ量ｘmaxは、次式によって表される。
【００７２】
【数１４】
2013210918.tif 000016
【００７３】
最大オーバーラップ量ｘmaxは、粒子同士が最も接近したときの粒子同士の半径が重なっている長さを表す。最大オーバーラップ量ｘmaxは、衝突時をｘ＝０として、粒子同士が最も接近したときの粒子ｉの変位と言い換えることもできる。
【００７４】
式（７）より、粒子同士が最も接近したときの接触力Ｆｃmaxは、最大オーバーラップ量ｘmaxを用いて以下のように表すことができる。
【００７５】
【数１５】
2013210918.tif 000017
【００７６】
式（１４）および式（１５）より、粒子同士が最も接近したときの（最大の）接触力Ｆｃmaxとばね定数ｋとの関係は、次式のようになる。
【００７７】
【数１６】
2013210918.tif 000018
【００７８】
また、式（１３）より、接触時間Ｔｃとばね定数ｋとの関係は、次式のようになる。
【００７９】
【数１７】
2013210918.tif 000019
【００８０】
式（１６）および式（１７）より、最大の接触力Ｆｃmaxは、ばね定数ｋの２／５乗に比例し、接触時間Ｔｃは、ばね定数ｋの−２／５乗に比例することが分かる。すなわち、接触力の力積（式（５）の右辺のＦｃの積分項）はばね定数ｋによらず一定とみなすことができる。
【００８１】
ここで、シミュレーションの計算時間の短縮のために本来のばね定数ｋよりも小さい修正ばね定数を使用する場合、妥当なシミュレーション結果を得るためには、修正ばね定数を使用しても衝突後の粒子の挙動が変化しないことが重要である。衝突における粒子の挙動（付着するか、離れるか）が変化しなければ、本来のばね定数ｋより小さい修正ばね定数を使用してシミュレーションを行っても、妥当なシミュレーション結果を得ることができる。すなわち、粒子系の状態を正確にシミュレーションすることができる。そして、小さい修正ばね定数を使用すれば、衝突の接触時間が長くなるので、シミュレーションの計算において時間ステップ（時間刻み）を大きく設定しても、妥当なシミュレーション結果が得られる。時間ステップを大きく設定することで、シミュレーションの計算負荷を小さくすることができる。
【００８２】
ばね定数を変えても衝突後の粒子の挙動が変化しないということは、式（３）の左辺が表す衝突における運動量の変化が変わらないということである。言い換えれば、ばね定数を変えても粒子の挙動が変化する閾値となる相対速度ｖｃ（衝突後に粒子が付着した状態を維持する限界の初期相対速度ｖｃ）が一定になるようにすればよい。すなわち、式（５）において左辺が一定であるとすると、右辺の接触力の力積（Ｆｃの積分項）はばね定数によらずに一定なので、右辺の付着力の力積もばね定数の変化に関わらず一定にならなければならない。
【００８３】
すなわち、本来のばね定数ｋよりも小さい修正ばね定数を使用する場合、次式が成り立つようにしなければならない。
【００８４】
【数１８】
2013210918.tif 000020
【００８５】
上式および式（１３）より、次式が得られる。
【００８６】
【数１９】
2013210918.tif 000021
【００８７】
ここで、ρ、Ｒ、ｖｃは、ばね定数に無関係な定数なので、付着力Ｆａとばね定数ｋとの関係を表す次式が得られる。
【００８８】
【数２０】
2013210918.tif 000022
【００８９】
すなわち、本来のばね定数ｋよりも小さい修正ばね定数を用いてシミュレーションを行う場合、上式に従って付着力Ｆａも小さくしなければならない。さもないと、閾値となる初期相対速度ｖｃが変化し、粒子の付着／反発の挙動が変化してしまう。
【００９０】
なお、式（６）のように粘性減衰および付着力を考慮した場合も、粒子の挙動を変化させないようなばね定数と付着力の関係式（式（２０））は同様に成り立つ。以下にこれについて説明する。
【００９１】
１つの粒子が同じ材質の壁面に衝突した時に粒子に働く力を式（６）とし、時間刻みを細かくして１つの粒子について数値計算を行うことで、ある相対速度ｖ0で衝突した粒子が衝突後に壁面に付着するか、離れるかを計算することができる。それゆえ、付着した状態を維持する限界の初期相対速度ｖｃを求めることができる。
【００９２】
図７は、粘性減衰および付着力を考慮した運動方程式（式（６））に基づいて一粒子の動きを数値計算して、付着した状態を維持する限界の初期相対速度ｖｃがある値になるときの、ばね定数ｋと付着力Ｆａとをプロットしたグラフである。プロットした点を曲線でフィッティングすると、Ｆａ＝４×１０-11×ｋ0.4 の関係が得られた。すなわち、付着した状態を維持する限界の初期相対速度ｖｃが一定になるような、付着力Ｆａとばね定数ｋとの関係は次式で表される。
【００９３】
【数２１】
2013210918.tif 000023
【００９４】
これは、粘性減衰および付着力を考慮しない場合の関係式（式（２０））と一致する。また、参考までに壁面と衝突する際の初期相対速度ｖ0と接触時間Ｔｃ2とについて、次式の関係が得られた。
【００９５】
【数２２】
2013210918.tif 000024
【００９６】
これは、粘性減衰および付着力を考慮しない場合の関係式（式（１３））と一致する。
【００９７】
本実施の形態では、ヘルツ接触モデルのばね定数として本来のばね定数ｋｒよりも小さい修正ばね定数ｋｄを用いると同時に、付着力として本来の付着力Ｆａｒよりも小さい修正付着力Ｆａｄを用いて、シミュレーションを行う。修正付着力Ｆａｄは、粒子の挙動を変えないために、閾値となる速度ｖｃが変化しないように選ぶのが好ましい。すなわち、本来のばね定数ｋｒと修正ばね定数ｋｄとの比の２／５乗と、本来の付着力Ｆａｒと修正付着力Ｆａｄとの比とが、ほぼ同じになるように修正付着力Ｆａｄを設定する。この関係を式に表せば以下のようになる。
【００９８】
【数２３】
2013210918.tif 000025
【００９９】
上式の関係を満たすような修正ばね定数ｋｄおよび修正付着力Ｆａｄを用いれば、閾値となる速度ｖｃが変化しない、すなわち、粒子の付着に関する挙動が変化しないので、妥当なシミュレーション結果を得ることができる。よって、計算に用いるばね定数を小さくすると共に、シミュレーションにおける時間ステップΔｔを大きくすることができるので、計算時間を短縮することができる。
【０１００】
図８は、参考例として付着力を一定にしてばね定数のみを変化させた場合の、閾値となる速度ｖｃの値をプロットしたグラフである。ここでは、粘性減衰および付着力を考慮した運動方程式（式（６））に基づいて一粒子の衝突時の動きを数値計算している。なお、粒子として直径６０μｍのガラスビーズを仮定している。付着力としては、当該ガラスビーズの物性値を用いている。図８によると、付着力を変化させずにばね定数ｋを変化させた場合（修正ばね定数を使用した場合）、閾値となる速度ｖｃが変わってしまうことが分かる。すなわち、図８に示す結果は、付着力を一定にして修正ばね定数を使用して計算すると、粒子の挙動が変わってしまい、妥当なシミュレーション結果が得られないことを示している。
【０１０１】
図９は、式（２３）の関係にしたがってばね定数と付着力とを変化させた場合の、閾値となる速度ｖｃの値をプロットしたグラフである。ここでは、粘性減衰および付着力を考慮した運動方程式（式（６））に基づいて一粒子の衝突時の動きを数値計算している。物性値等は図８の場合と同じである。図９によると、式（２３）の関係にしたがってばね定数と付着力とを変化させた場合（修正ばね定数と修正付着力とを使用した場合）、閾値となる速度ｖｃは変化しないことが分かる。閾値となる速度ｖｃが変わらなければ、衝突後の粒子の挙動（付着するか、離れるか）は変わらない。すなわち、図９に示す結果は、式（２３）の関係にしたがった修正ばね定数および修正付着力を使用して計算すると、粒子の挙動が本来の挙動と変わらず、妥当なシミュレーション結果を得ることができることを示している。
【０１０２】
なお、修正ばね定数ｋｄと修正付着力Ｆａｄを用いて粒子挙動を計算する本実施の形態のシミュレーション方法は、本来の付着力Ｆａｒが重力Ｆｇの０．１倍以上である場合に、好適に用いることができる。
【０１０３】
（時間ステップ）
妥当なシミュレーション結果を得るための、適切な時間ステップについて、以下に説明する。
【０１０４】
接触力のモデルとして線形ばねモデルを適用した場合、線形ばね定数をｋｌとすると、粒子の接触時間Ｔｃ3は次式で表される。
【０１０５】
【数２４】
2013210918.tif 000026
【０１０６】
ここで、ｍは粒子の質量である。
【０１０７】
線形ばねモデルを適用した流動層の離散要素法シミュレーションにおいて、接触時間Ｔｃ3を５分割以上に分割した時間刻みでシミュレーションを行えば、安定かつ妥当な計算結果が得られることが経験的に知られている（非特許文献１参照）。よって、線形ばねモデルの場合の妥当な時間ステップΔｔの範囲は、次式で表される。
【０１０８】
【数２５】
2013210918.tif 000027
【０１０９】
上式を見れば分かるように、時間ステップΔｔの上限は、粒子の質量ｍと線形ばね定数ｋｌとに依存し、粒子の質量ｍが小さいまたは粒子が硬い（線形ばね定数ｋｌが大きい）場合、小さくなり、計算時間が長くなる。
【０１１０】
一方、接触モデルとしてヘルツ接触モデルを適用した場合、粒子の接触時間は、接触時の初期相対速度により変化し、解析的に求めることができない。そこで、次式で示す一自由度振動系を考える。
【０１１１】
【数２６】
2013210918.tif 000028
【０１１２】
この非線形の振動系は、振幅が大きくなるほど振動数が増大する。そのため、時間を刻んだシミュレーション（数値積分）の安定性は、初期条件と時間ステップの長さとに依存する。そこで、実際にシミュレーションを行う粒子系の予想される最大の接触時の初期相対速度ｖ0を仮定して、この振動系（式（２６））のシミュレーションを行い、計算が安定的に行える限界の時間ステップを特定することで、妥当な時間ステップを決定することができる（参考：非特許文献２）。なお、限界を求めなくとも、上記振動系のシミュレーションが安定的に行える時間ステップを、実際にシミュレーションを行う粒子系に使用する時間ステップとして採用すればよい。
【０１１３】
また、接触時間Ｔｃを５分割以上に分割した時間ステップを採用してもよい。上記振動系における粒子の接触時間Ｔｃは、次式で表される。
【０１１４】
【数２７】
2013210918.tif 000029
【０１１５】
ここで、ρは粒子密度、νは粒子のポアソン比、Ｅは粒子のヤング率、Ｒは粒子半径である。接触時間Ｔｃは、付着力の影響により長くなる。よって、この接触時間Ｔｃを５分割以上に分割した時間ステップでシミュレーションを行えば、妥当なシミュレーション結果
を得ることができる。よって、ヘルツ接触モデルの場合の妥当な時間ステップΔｔの範囲は、次式で表される。
【０１１６】
【数２８】
2013210918.tif 000030
【０１１７】
ヘルツ接触モデルを適用する場合のこの非線形の振動系は、衝突の初期相対速度ｖ0が大きいほど、接触時間Ｔｃが短くなる。そのため、予想される最大の接触時の初期相対速度ｖ0を仮定して、そこから得られる予想される最小の接触時間Ｔｃを５分割したものを、時間ステップΔｔとしてもよい。
【０１１８】
上述したように、付着力を考慮しない線形ばねモデルの流動層の離散要素法シミュレーションにおいて、実際の粒子の反発を表すばね定数より小さい、修正ばね定数を用いて計算を行っても、問題ないことが知られている。例えば、硬い粒子の衝突において、粒子の弾性変形量が粒子径の数％以内（例えば１％以内）に収まるように線形ばねモデルの修正ばね定数を設定して計算しても、妥当なシミュレーション結果が得られることが経験的に知られている。ヘルツ接触モデルの場合でも、硬い粒子の衝突において、粒子の弾性変形量が粒子径の数％以内（例えば１％以内）に収まるようにヘルツ接触モデルの修正ばね定数を設定して計算すればよい。衝突の際の相対速度には幅があるが、衝突の平均相対速度の場合に、上記条件が満たされるようにしてもよいし、例えば、９５％の衝突において、上記条件が満たされるようにしてもよい。また、気体の吹き上げがないときに、堆積した粒子の重さによって一番下の粒子の弾性変形量が粒子径の数％以内に収まるように、ヘルツ接触モデルの修正ばね定数を設定してもよい。このような条件でヘルツ接触モデルの修正ばね定数を設定する場合、硬い粒子であれば、本来のばね定数ｋｒの０．１倍以下の値に修正ばね定数ｋｄを設定できる。例えば、ヘルツ接触モデルの修正ばね定数ｋｄが本来のばね定数の０．１倍であれば、時間ステップを約２．５倍（１０2/5倍）に設定することができる。
【０１１９】
＜粒子挙動シミュレーション装置の構成＞
図４は、本実施の形態の粒子挙動シミュレーション装置１０の機能的構成を示すブロック図である。粒子挙動シミュレーション装置１０は、条件入力部１１、挙動計算部１２、および接触力計算部（接触力特定部）１３を備える。
【０１２０】
条件入力部１１は、修正ばね定数取得部１４、修正付着力取得部１５、および時間ステップ取得部１６を備える。条件入力部１１は、シミュレーションに必要な各種の条件が記載された設定ファイルを読み込み、シミュレーションの条件の入力を受け付ける。シミュレーションの条件としては、容器の形状、粒子数、粒子径、粒子質量、粒子同士の衝突における修正ばね定数、粒子同士の衝突における修正付着力、粒子同士の衝突における反発係数、分散板の下側から流入する気体の流速（流量）、気体の密度、気体の粘性係数、および時間ステップ等がある。シミュレーションの条件は他にもあるが省略する。なお、条件入力部１１は、キーボード等の入力装置を備えて、ユーザから直接シミュレーションの条件の入力を受け付けてもよい。
【０１２１】
修正ばね定数取得部１４は、入力されたシミュレーション条件のうち、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、修正ばね定数ｋｄおよび粘性減衰係数ηを接触力計算部１３に出力する。修正ばね定数ｋｄは、本来のばね定数ｋｒより小さい値である。
【０１２２】
修正付着力取得部１５は、入力されたシミュレーション条件のうち、修正付着力Ｆａｄを取得し、修正付着力Ｆａｄを挙動計算部１２に出力する。修正付着力Ｆａｄは、本来のばね定数ｋｒと修正ばね定数ｋｄとの比（の２／５乗）に応じて本来の付着力Ｆａｒより値を小さくしたものである。
【０１２３】
時間ステップ取得部１６は、入力されたシミュレーション条件のうち、時間ステップｔｄを取得し、時間ステップｔｄを挙動計算部１２に出力する。時間ステップｔｄは、小さくした修正ばね定数ｋｄに応じて値を大きくした時間ステップであり、修正時間ステップと呼ぶことができるものである。修正時間ステップｔｄは、本来のばね定数ｋｒと修正ばね定数ｋｄとの比の平方根に応じて値を大きく設定できる。安定かつ妥当な結果を得るためには、例えば、衝突における粒子の接触時間Ｔｃの１／５以下の時間ステップで計算を行うという条件を満たせばよい。小さい修正ばね定数ｋｄを用いる場合、接触時間Ｔｃは長くなる。
【０１２４】
条件入力部１１は、取得したシミュレーションの他の条件についても挙動計算部１２に出力する。
【０１２５】
挙動計算部１２は、各シミュレーション条件を用いて、修正時間ステップｔｄ毎に、各時刻における粒子の挙動（位置および運動）の計算を行う。挙動計算部１２は、ある時刻ｔにおいて、各粒子に働く力を求める。ここで、着目する１つの粒子に働く力Ｆは、次式で表される。
【０１２６】
【数２９】
2013210918.tif 000031
【０１２７】
挙動計算部１２は、着目する粒子が他の粒子（または壁）と接触している（着目する粒子と他の粒子との距離が所定の距離以内である）場合、着目する粒子に働く接触力を求めるために、着目する粒子と他の粒子（または壁）との距離および相対速度を接触力計算部１３に出力する。
【０１２８】
接触力計算部１３は、着目する粒子と他の粒子との距離および相対速度と、修正ばね定数ｋｄと、粘性減衰係数ηとを用いて、着目する粒子に働く接触力（修正接触力）Ｆｃｄを求める。法線方向の修正接触力Ｆｃｄは次式で表される。
【０１２９】
【数３０】
2013210918.tif 000032
【０１３０】
接触力計算部１３は、算出した着目する粒子に働く修正接触力Ｆｃｄを、挙動計算部１２に出力する。
【０１３１】
挙動計算部１２は、接触力計算部１３から受け取った着目する粒子に働く修正接触力Ｆｃｄ、および、修正付着力Ｆａｄ、着目する粒子に働く流体力Ｆｄ、重力Ｆｇを用いて、着目する粒子の修正時間ステップｔｄ後の位置および速度を計算する。挙動計算部１２は、同様に、全ての粒子について、修正時間ステップｔｄ後の位置および速度を計算する。なお、挙動計算部１２は、粒子の運動状態として粒子の自転角速度をさらに求めてもよい。
【０１３２】
挙動計算部１２は、このようにして、修正時間ステップｔｄ毎の各時刻における各粒子の位置および運動を求める。挙動計算部１２は、各時刻の各粒子の位置および運動の情報を、記憶装置（図示せず）に出力する。例えば、粒子挙動シミュレーション装置１０は、各粒子の位置および運動を、時刻毎に順次表示装置（図示せず）に点や丸で可視化して表示させ、粒子の挙動をユーザに提示してもよい。
【０１３３】
＜粒子挙動シミュレーションのフロー＞
次に、粒子挙動シミュレーションのパラメータの決定方法および粒子挙動シミュレーション装置１０における処理のフローについて説明する。
【０１３４】
（ばね定数ｋｒの特定）
本実施の形態の粒子挙動シミュレーションを行うために、修正ばね定数ｋｄおよび修正付着力Ｆａｄが必要になる。そのために、本来のばね定数ｋｒおよび本来の付着力Ｆａｒが必要になる。
【０１３５】
粒子同士が接触しているときに粒子に働く接触力を、非線形ばねによる弾性反発力モデルで表した場合のばね定数ｋｒを、以下の手順で求める。
【０１３６】
ヘルツの接触理論によれば、球（粒子）が同じ材質の壁面と衝突するときの法線方向のばね定数は、式（１０）で表される。また、材質と大きさが等しい２つの球（粒子）が衝突するときの法線方向のばね定数は、式（１１）で表される。すなわち、ヘルツ接触モデルの本来のばね定数ｋｒは、粒子のヤング率Ｅ、粒子のポアソン比ν、および粒子の半径Ｒから決定することができる。
【０１３７】
なお、粘性減衰係数ηは、反発係数の実測値から特定する。
【０１３８】
（付着力Ｆａｒの特定）
付着力Ｆａｒは、実測によって特定することができる。例えば、２つの粒子を押しつけて付着させてから、それを引きはがす際に必要な力を付着力とすることができる。なお、付着力Ｆａｒを、Hamaker理論等に基づいて特定（仮定）してもよい。
【０１３９】
（修正ばね定数ｋｄの特定）
次に、ばね定数ｋｒより小さい修正ばね定数ｋｄを特定（設定）する。上述したように、例えば、硬い粒子の衝突において、粒子の弾性変形量が粒子径の数％以内（例えば３％以内または１％以内）に収まるように修正ばね定数を設定して計算しても、妥当なシミュレーション結果が得られることが経験的に知られている。本実施の形態では、主な衝突による粒子の弾性変形量（図２における時刻ｔ１の距離ｘ（粒子の重なりの距離））が粒子径の１％以内になるように、修正ばね定数ｋｄを設定する。
【０１４０】
（修正付着力Ｆａｄの特定）
次に、上記ばね定数ｋｒと上記修正ばね定数ｋｄとの比（の２／５乗）に応じて上記付着力Ｆａｒより値を小さくした修正付着力Ｆａｄを特定する。修正付着力Ｆａｄは、式（２３）によって求めることができる。
【０１４１】
（修正時間ステップｔｄの特定）
想定される衝突時の最大の初期相対速度ｖ0を仮定し、計算に用いる修正時間ステップｔｄを求める。本実施の形態では、衝突における粒子の接触時間Ｔｃの１／５以下を、修正時間ステップｔｄとする。よって、修正時間ステップｔｄは、式（２８）によって決定される。
【０１４２】
（粒子挙動の計算）
図５は、粒子挙動シミュレーション装置１０の処理を示すフロー図である。
【０１４３】
粒子挙動シミュレーション装置１０の条件入力部１１は、上記のようにして特定した修正ばね定数ｋｄ、修正付着力Ｆａｄおよび修正時間ステップｔｄを含むシミュレーションの条件を、取得する（Ｓ１）。
【０１４４】
修正ばね定数取得部１４は、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、接触力計算部１３に出力する。
【０１４５】
修正付着力取得部１５は、修正付着力Ｆａｄを取得し、挙動計算部１２に出力する。
【０１４６】
時間ステップ取得部１６は、時間ステップｔｄを取得し、挙動計算部１２に出力する。
【０１４７】
着目する粒子が他の粒子と接触している場合（Ｓ２でＹｅｓ）、挙動計算部１２は、着目する粒子と他の粒子との距離および相対速度を接触力計算部１３に出力する（Ｓ３）。
【０１４８】
接触力計算部１３は、着目する粒子と他の粒子との距離および相対速度と、修正ばね定数ｋｄと、粘性減衰係数ηとを用いて、式（１２）から、着目する粒子に働く接触力（修正接触力）Ｆｃｄを求める（Ｓ４）。接触力計算部１３は、着目する粒子に働く修正接触力Ｆｃｄを挙動計算部１２に出力する。
【０１４９】
着目する粒子が他の粒子と接触していない場合（Ｓ２でＮｏ）、接触力および付着力は働かない（Ｆｃｄ＝Ｆａｄ＝０）。
【０１５０】
Ｓ２でＮｏ、またはＳ４の後、挙動計算部１２は、修正接触力Ｆｃｄ、着目する粒子に働く修正付着力Ｆａｄ、流体力Ｆｄ、および重力Ｆｇを用いて、修正時間ステップｔｄ後における着目する粒子の位置および速度を求める（Ｓ５）。
【０１５１】
修正時間ステップｔｄ後における位置および速度を求めていない粒子が未だある場合（Ｓ６でＮｏ）、残りの粒子に着目してＳ２〜Ｓ５の処理を繰り返す。
【０１５２】
ある時刻について全ての粒子について修正時間ステップｔｄ後における位置および速度を求めた場合（Ｓ６でＹｅｓ）、所定の終了時刻（時間刻み）まで計算を行ったか否かを判定する（Ｓ７）。ここで、１秒間の粒子の挙動のシミュレーションを行う場合、終了時刻は最初の時刻から１秒後である。
【０１５３】
所定の終了時刻まで計算していない場合（Ｓ７でＮｏ）、修正時間ステップｔｄだけ時刻を進め、次の時刻における計算（Ｓ２〜Ｓ６）を繰り返す。
【０１５４】
所定の終了時刻まで計算が完了した場合（Ｓ７でＹｅｓ）、処理を終了する。
【０１５５】
（実施の形態１のまとめ）
本実施の形態によれば、本来のばね定数ｋｒより小さい修正ばね定数ｋｄと、それに応じて本来の付着力Ｆａｒより小さくした修正付着力Ｆａｄとを用いて、離散要素法によって粒子の挙動を計算する。よって、小さい修正ばね定数ｋｄに応じて時間ステップを大きくすることができる。そのため、付着力を考慮する場合においても、シミュレーションの計算時間を短くし、かつ、妥当なシミュレーション結果を得ることができる。
【０１５６】
また、修正付着力Ｆａｄは、修正ばね定数ｋｄに応じて小さくなっていればよく、厳密
に式（２３）の関係を満たさなくても、ある程度の幅以内であれば、妥当な結果を得ることができる。例えば、修正付着力Ｆａｄを、次式の範囲内で設定してもよい。
【０１５７】
【数３１】
2013210918.tif 000033
【０１５８】
または、修正付着力Ｆａｄと付着力Ｆａｒとの比が、修正ばね定数ｋｄとばね定数ｋｒとの比の２／５乗と同じオーダーとなるように、修正付着力Ｆａｄを設定してもよい。
【０１５９】
なお、本実施の形態では、修正接触力のモデルとしてヘルツ接触モデルを仮定している。さらに、接触力として、複数のヘルツ接触モデルのばねを組み合わせた非線形のばねを考慮してもよい。
【０１６０】
また、上記は粒子同士の衝突について説明したが、粒子と壁との衝突について、別途修正ばね定数、および修正付着力を求めて、計算に用いてもよい。
【０１６１】
［実施の形態２］
以下、本実施の形態について、図６を参照して詳細に説明する。なお、説明の便宜上、実施の形態１にて説明した図面と同じ機能を有する部材・構成については、同じ符号を付記し、その詳細な説明を省略する。本実施の形態では、図１に示す流動層装置における粒子をシミュレーションの対象とする。
【０１６２】
＜粒子挙動シミュレーション装置の構成＞
図６は、本実施の形態の粒子挙動シミュレーション装置２０の機能的構成を示すブロック図である。粒子挙動シミュレーション装置２０は、条件入力部２１、挙動計算部１２、接触力計算部（接触力特定部）１３、および修正付着力特定部２２を備える。
【０１６３】
条件入力部２１は、ばね定数取得部２３、付着力取得部２４、修正ばね定数取得部２５、および時間ステップ取得部１６を備える。条件入力部２１は、シミュレーションに必要な各種の条件が記載された設定ファイルを読み込み、シミュレーションの条件の入力を受け付ける。
【０１６４】
ばね定数取得部２３は、入力されたシミュレーション条件のうち、本来のばね定数ｋｒを取得し、ばね定数ｋｒを修正付着力特定部２２に出力する。
【０１６５】
付着力取得部２４は、入力されたシミュレーション条件のうち、本来の付着力Ｆａｒを取得し、付着力Ｆａｒを修正付着力特定部２２に出力する。
【０１６６】
修正ばね定数取得部２５は、入力されたシミュレーション条件のうち、修正ばね定数ｋｄおよび粘性減衰係数ηを取得し、修正ばね定数ｋｄおよび粘性減衰係数ηを接触力計算部１３に出力する。また、修正ばね定数取得部２５は、修正ばね定数ｋｄを修正付着力特定部２２に出力する。
【０１６７】
修正付着力特定部２２は、修正付着力Ｆａｄと付着力Ｆａｒとの比が、修正ばね定数ｋｄとばね定数ｋｒとの比の２／５乗と同じになるような、修正付着力Ｆａｄを特定する。すなわち、修正付着力特定部２２は、式（２３）に従って修正付着力Ｆａｄを特定する。修正付着力特定部２２は、特定した修正付着力Ｆａｄを、挙動計算部１２に出力する。
【０１６８】
条件入力部２１は、取得したシミュレーションの他の条件についても挙動計算部１２に出力する。
【０１６９】
挙動計算部１２および接触力計算部１３の構成は、実施の形態１と同様であるので、その説明を省略する。
【０１７０】
本実施の形態では、本来のばね定数ｋｒ、本来の付着力Ｆａｒ、および修正ばね定数ｋｄに基づき、修正付着力特定部２２が適切な修正付着力Ｆａｄを求める。よって、付着力を考慮する場合において、短い計算時間で、妥当なシミュレーション結果を得ることができる。
【０１７１】
最後に、粒子挙動シミュレーション装置１０・２０の各ブロック、特に条件入力部１１、挙動計算部１２、接触力計算部１３、修正ばね定数取得部１４、修正付着力取得部１５、時間ステップ取得部１６、条件入力部２１、修正付着力特定部２２、ばね定数取得部２３、付着力取得部２４、および修正ばね定数取得部２５は、ハードウェアロジックによって構成してもよいし、次のようにＣＰＵ（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）、携帯電話網、衛星回線、地上波デジタル網等の無線でも利用可能である。
【０１７５】
本発明は上述した各実施形態に限定されるものではなく、請求項に示した範囲で種々の変更が可能であり、異なる実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても本発明の技術的範囲に含まれる。
【産業上の利用可能性】
【０１７６】
本発明は、粒子の挙動を予測する粒子挙動シミュレーションに利用することができる。
【符号の説明】
【０１７７】
１    流動層装置
２    処理容器
３    分散板
４    気体供給口
５    排気口
６    触媒粒子
１０、２０    粒子挙動シミュレーション装置
１１、２１    条件入力部
１２    挙動計算部
１３    接触力計算部（接触力特定部）
１４、２５    修正ばね定数取得部
１５    修正付着力取得部
１６    時間ステップ取得部
２２    修正付着力特定部
２３    ばね定数取得部
２４    付着力取得部
【図１】
2013210918.tif 000035
【図２】
2013210918.tif 000036
【図３】
2013210918.tif 000037
【図４】
2013210918.tif 000038
【図５】
2013210918.tif 000039
【図６】
2013210918.tif 000040
【図７】
2013210918.tif 000041
【図８】
2013210918.tif 000042
【図９】
2013210918.tif 000043