JP7101387B2

JP7101387B2 - シミュレーション装置、シミュレーション方法、プログラム

Info

Publication number: JP7101387B2
Application number: JP2020162719A
Authority: JP
Inventors: 元彬猿渡; 英也仲村
Original assignee: Sumitomo Metal Mining Co Ltd; University Public Corporation Osaka
Current assignee: Sumitomo Metal Mining Co Ltd; University Public Corporation Osaka
Priority date: 2020-06-01
Filing date: 2020-09-28
Publication date: 2022-07-15
Anticipated expiration: 2040-09-28
Also published as: JP2021190060A

Description

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

特許文献１には、シミュレーション条件の入力を行う入力装置と、
シミュレーション結果を出力する出力装置と、
前記入力装置から入力されたシミュレーション条件に基づいて、大きさが異なる複数の
粒子を含む粉粒体の挙動を解析する処理装置とを有し、
前記処理装置は、
前記入力装置から入力されたシミュレーション対象の粉粒体の粒径分布を規定するパラメータの値、及び粒子を粗視化する基準となる粗視化係数の値に基づいて、粗視化された粉粒体の挙動をシミュレーションにより求め、
シミュレーションによって求められた粒子の挙動と、入力された粗視化係数の値とを関連付けて前記出力装置に出力するシミュレーション装置等に開示されているように、粉粒体についてのシミュレーション装置が開示されている。

特開２０２０－５７１３５号公報

工場等におけるプロセスの改善や、製造工程を検討する際の試験工数の削減等を目的として、複数の粒子を含む粉体（粉粒体）の挙動の解析を離散要素法（ＤＥＭ：ＤｉｓｃｒｅｔｅＥｌｅｍｅｎｔＭｅｔｈｏｄ）計算等により行うことが従来からなされてきた。

離散要素法計算は、個々の粒子について運動方程式を解くことで、粉体全体の運動を記述するシミュレーション技術である。

しかしながら、離散要素法計算では、取り扱う粒子数が多くなるほど計算負荷が高くなる。このため、例えば工場で用いるプラント等のような大きな規模での粉体の挙動について解析を行う場合、計算量が膨大になり、現実的には計算を行うことが困難となる。

そこで、複数個の粒子から構成される粒子群を１個の粗視化粒子とする粗視化の手法を用いたシミュレーション装置について従来から検討されていた（特許文献１を参照）。粗視化の手法を用いたシミュレーション装置では、正確な解析結果を得るため、計算で用いる粗視化粒子に関するパラメータを適切に選択することが求められる。このため、粗視化粒子に関するパラメータを新たな手法により選択、設定し、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置が求められていた。

上記従来技術が有する問題に鑑み、本発明の一側面では、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置を提供することを目的とする。

上記課題を解決するため本発明の一態様によれば、
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置を提供する。

本発明の一態様によれば、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置を提供することができる。

複数個の粒子から構成される粒子群と、粒子群と壁面との衝突時の説明図である。複数個の粒子から構成される粒子群を粗視化した粗視化粒子と、粗視化粒子と壁面との衝突時の説明図である。本発明の実施形態に係るシミュレーション装置のハードウェア構成図である。本発明の実施形態に係るシミュレーション装置の機能を示すブロック図である。本発明の実施形態に係るシミュレーション方法を説明するフローチャートである。実験例１における粉体の平均温度の変化を示すグラフである。実施例２－１におけるキルン内での粉体の混合の状態を示す図である。比較例２－１におけるキルン内での粉体の混合の状態を示す図である。比較例２－２におけるキルン内での粉体の混合の状態を示す図である。実施例２－１におけるキルン内の粉体の温度分布を示す図である。比較例２－１におけるキルン内の粉体の温度分布を示す図である。比較例２－２におけるキルン内の粉体の温度分布を示す図である。実験例２における粉体の平均温度の変化を示すグラフである。

本開示の一実施形態（以下「本実施形態」と記す）に係るシミュレーション装置、シミュレーション方法、プログラムの具体例を、以下に図面を参照しながら説明する。なお、本発明はこれらの例示に限定されるものではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内での全ての変更が含まれることが意図される。
１．第１の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
（１－１）粒子の粗視化について
本実施形態のシミュレーション装置の詳細について説明する前に、本実施形態のシミュレーション装置で用いることができる、複数個の粒子から構成される粒子群の粗視化、および粗視化した粒子である粗視化粒子に関連するパラメータの算出方法について以下に説明する。

既述のように、離散要素法計算では、取り扱う粒子数が多くなるほど計算負荷が高くなる。このため、例えば工場で用いるプラント等のような大きな規模での粉体の挙動について解析を行う場合、計算量が膨大になり、現実的には計算を行うことが困難となる。

そこで、多量の粒子を含む粉体について挙動を解析する場合、計算量を抑制するため、例えば図１（Ａ）に示した複数個の粒子から構成される粒子群１１を、図２（Ａ）に示すような１個の大きな粒子である粗視化粒子２１として取り扱う粗視化の技術が必要となる。

ただし、粗視化前の個別の粒子と、粗視化粒子とでは比表面積等が異なるため、計算に要する一部のパラメータは変化することになる。このため、粗視化粒子のパラメータを適切に決定する必要がある。
（１－２）粗視化粒子の粒子挙動の計算に用いるパラメータについて
そこで、本発明の発明者は、粗視化粒子に関するパラメータを決定する方法について検討を行った。計算にあたって、図１（Ａ）に示した粗視化する前の複数個の粒子から構成される粒子群１１が壁面１２に衝突する場合と、図２（Ａ）に示した粗視化粒子２１が壁面１２に衝突する場合とをモデルに用いた。以下の説明では壁面と粒子とが衝突する場合を例に粗視化粒子に関するパラメータを決定する方法を記載するが、粒子同士が衝突する場合でも同様の議論となるため説明を省略する。

図１（Ａ）に示すように、複数個の粒子から構成される粒子群１１が、立方体状に、縦方向、横方向、高さ方向に２個ずつ、合計で２^３個並んでいるとする。後述するように、係る８個の粒子をまとめて１個の粗視化粒子とする場合、一辺の方向に並べた粒子の数、すなわち２を粗視化倍率とする。

図１（Ａ）に示した複数個の粒子１１Ａ、１１Ｂを含む粒子群１１が壁面１２に衝突する際に、粒子群１１のうち壁面１２側に位置する粒子１１Ａが壁面あるいは外部粒子から受ける力を図１（Ｂ）に示すようにＦ_ｗとする。また、図１（Ｂ）に示すように粒子１１Ａの壁面１２あるいは外部粒子とのオーバーラップ量をδ_ｗ、粒子１１Ｂの隣接する粒子１１Ａとのオーバーラップ量をδ_ｐとする。なお、図１（Ｂ）は、粒子群１１が壁面１２に衝突する際の様子を側面側から見た図になる。

この場合、粒子群１１に加わる力の大きさは、以下の式（１）で表すことができる。

なお、式（１）におけるαは粗視化倍率であり、既述のように粒子群１１を１個の粗視化粒子とした場合に一辺の方向に並べた粒子の数を意味する。図１（Ａ）に示した粒子群１１を図２（Ａ）に示した１個の粗視化粒子２１とする場合、α＝２となる。

また、ｍは各粒子１１Ａ、１１Ｂの質量を、ａ_Ｇは粒子群１１の重心の加速度を、η_ｗは壁面１２あるいは外部粒子と粒子１１Ａとの反発係数から算出される粘性係数をそれぞれ意味している。作用反作用の法則により粒子間の接触力が打ち消されるため、壁面１２とは直接接しない粒子１１Ｂが、隣接する粒子１１Ａから受ける力Ｆ_ｐは式（１）には表れないことになる。

なお、上述のη_ｗ等を反発係数から算出する際に用いる、反発係数ｅと粘性係数ηとの関係は、以下の式（Ａ）で表すことができる。式（Ａ）中のｍ_＊は換算質量を、Ｋはバネ係数を意味している。

次に、図２（Ａ）に示すように、図１（Ａ）に示した８個の粒子からなる粒子群１１を１つの粗視化粒子２１としたとする。この場合において、粗視化粒子２１が、壁面１２に衝突した場合に、粗視化粒子２１が受ける力は以下の式（２）で表すことができる。

式（２）中のＦ_ｃｗは、図２（Ｂ）に示すように粗視化粒子２１が壁面１２あるいは外部粒子から受ける力を、δ_ｃｗは粗視化粒子２１の壁面１２あるいは外部粒子とのオーバーラップ量を、η_ｃｗは粗視化粒子２１と壁面１２あるいは外部粒子との反発係数から算出される粘性係数をそれぞれ意味している。なお、図２（Ｂ）は、粗視化粒子２１が壁面１２に衝突する際の様子を側面側から見た図になる。

既述のように、粗視化は、離散要素法計算において、計算量を抑制するために行うものである。このため、粗視化粒子２１についての計算結果と、粗視化粒子２１とする前の粒子群１１についての計算結果とは一致することになる。

従って、粒子群１１について計算を行った既述の式（１）と、粒子群を粗視化した粗視化粒子２１について計算を行った既述の式（２）とから、対応するパラメータが一致することを示す以下の式（３）、式（４）が導出される。

また、Ｈｅｒｔｚ－Ｍｉｎｄｌｉｎの接触モデルにより、粒子のオーバーラップ量δ_ｗ、δ_ｐ、δ_ｃｗを用いて、各粒子に加わる力を以下の式（５）～式（７）のように表すことができる。式（５）～式（７）におけるＫ_ｗは粒子１１Ａと壁面１２あるいは外部粒子とのバネ係数を、Ｋ_ｐは粒子群１１の内部粒子のバネ係数を、Ｋ_ｃｗは粗視化粒子２１と壁面１２あるいは外部粒子とのバネ係数をそれぞれ意味している。

そうすると、式（３）、式（５）、式（７）から以下の式（８）の関係が導かれる。

ここで、以下の式（９）のようにＫ_ｒを定義すると、上記式（８）は以下の式（１０）のように表すことができる。

なお、粗視化前の粒子群１１と、粗視化粒子２１とで重心が一致しているとすると、以下の式（１１）の関係を満たすことになる。

従って、Ｋ_ｒを適切に設定することで、粗視化粒子２１の壁面１２あるいは外部粒子とのオーバーラップ量δ_ＣＷから、粗視化前の粒子群１１を構成する粒子のオーバーラップ量を算出できることもわかる。

そして、Ｋ_ｒは、衝突中の粗視化前の粒子群１１の弾性エネルギーと、粗視化粒子２１の弾性エネルギーとの関係を用いた特性方程式により算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、特性方程式を作成してＫ_ｒを算出できる。

壁面１２と衝突中の、粒子群１１の弾性エネルギーと、粗視化粒子の弾性エネルギーは、既述の粒子群１１を構成する粒子１１Ａ、１１Ｂに加わる力や、粗視化粒子に加わる力を表した式（５）～式（７）をオーバーラップ量の距離で積分することで算出できる。

従って、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーとを用いて、以下の式（１２）が得られる

上記式（１２）は、既述の式（８）～式（１１）を用いて以下の式（１３）に変形できる。

式（１３）は垂直方向のＫ_ｒの特性方程式である。そして、式（９）の定義式から明らかなように、Ｋ_ｒは粗視化粒子、および粗視化前の粒子群１１を構成する粒子のオーバーラップ量に関するパラメータであり、粗視化粒子の挙動を支配するパラメータである。このため、特性方程式によりＫ_ｒを事前に求めておくことによって、粗視化後粒子のオーバーラップ量から粗視化前粒子群のオーバーラップ量を算出するなど、粗視化粒子についてのパラメータを算出したり、粗視化粒子の挙動を計算したりすることが可能になる。

ここまで、壁面１２に対する垂直方向の運動方程式を用いて説明したが、接線方向の運動方程式や、回転の運動方程式についても同様のことが言える。

具体的には、接線方向の運動方程式は、式（１４）で表すことができる。

この場合も、式（１５）に示したように、Ｋ_ｒを設定すると、δ_ｗ、δ_ｐは式（１６）、式（１７）のように示すことができ、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーが等しいとすると、式（１８）が得られる。式（１８）を変形することで、接線方向の特性方程式である式（１９）が得られる。ただし、接線方向の接触モデルには線形バネモデルを使用した。このように、接触モデルによって弾性エネルギーの計算式が異なるが、必要に応じて特性方程式を変更する等により、適切に弾性エネルギーを算出できる。

回転の運動方程式については第２の実施形態で詳述する。
（１－３）粗視化粒子の伝熱に関するパラメータについて
既述のＫ_ｒを用いて、粗視化粒子の伝熱に関するパラメータである熱伝導率を求めることもできる。

粒子の伝熱は熱伝導率を用いて、以下の式（２０）、式（２１）のようにあらわされる。

なお、式（２０）、式（２１）において、Ｑは熱流量を、ｈは熱流量係数を、ΔＴは壁面あるいは粒子同士の温度差を、ｋ_ｗは壁面１２あるいは外部粒子と粒子１１Ａとの熱伝導率を、ｋ_ｐは粒子群１１内部の熱伝導率を、ａは粒子群１１の接触半径をそれぞれ意味する。

粒子と壁面や、他の粒子との接触半径は、粒径によって変化するため、粒径が大きくなると粒子に与えられる熱流量が変化し、粒子の温度変化に影響する。

接触半径は、各粒子のオーバーラップ量に依存する。そして、既述のように粗視化前後でのオーバーラップ量は既述の垂直方向の特性方程式の解に関係づけられている。このことから、垂直方向の特性方程式の解Ｋ_ｒを用いて伝熱方程式を書き換えると以下の式（２２）、式（２３）のようになる。

式（２３）におけるａ´は粗視化粒子の接触半径を、ｒは粒子群１１を構成する粒子の半径をそれぞれ意味する。

そして、粗視化粒子の伝熱方程式は以下の式（２４）、式（２５）で表すことができる。

式（２４）、式（２５）におけるＱ_ｃは粗視化粒子の熱流量を、ｈ_ｃは粗視化粒子の熱流量係数を、ΔＴ_ｃは粗視化粒子２１と壁面１２あるいは外部粒子との温度差を、ｋ_ｗ´は壁面１２と粗視化粒子２１との熱伝導率を、ｋ_ｐ´は外部粒子と粗視化粒子との熱伝導率をそれぞれ意味する。

熱伝導率に注目すると、粗視化前後での熱伝導率を、Ｋ_ｒを用いて以下の式（２６）、式（２７）のように変換すると、粗視化前後で等価な熱伝導方程式を得ることができる。ただし、ΔＴ_ｃ＝αΔＴとした。

つまり粗視化する際に、熱伝導率をＫ_ｒ ^１／３α^１／２倍することで、粗視化前の粒子群に与えられる熱流量と、粗視化後の粒子に与えられる熱流量を等しくすることができ、その結果、温度の時間変化を粗視化前後で一致させることができる。

なお、ここでは粗視化前後での熱伝導率を例に説明したが、熱伝導率以外のパラメータについても同様にして、既述の特性方程式の解Ｋ_ｒを用いて粗視化後のパラメータを算出できる。例えば、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出することもできる。なお、計算に適用するモデルに応じて、これらの係数は調整可能であり、上述のように特性方程式を用いて算出、換算できる。
（２）シミュレーション装置
本実施形態のシミュレーション装置は、複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であり、以下の部材を有することができる。
粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析部。

第２パラメータ算出部は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出する。

図３に示したハードウェア構成図に示すように、本実施形態のシミュレーション装置３０は、例えば、情報処理装置（コンピュータ）で構成され、物理的には、演算処理部であるＣＰＵ（ＣｅｎｔｒａｌＰｒｏｃｅｓｓｉｎｇＵｎｉｔ：プロセッサ）３１と、主記憶装置であるＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）３２やＲＯＭ（ＲｅａｄＯｎｌｙＭｅｍｏｒｙ）３３と、補助記憶装置３４と、入出力インタフェース３５と、出力装置である表示装置３６等を含むコンピュータシステムとして構成することができる。これらは、バス３７で相互に接続されている。なお、補助記憶装置３４や表示装置３６は、外部に設けられていてもよい。

ＣＰＵ３１は、シミュレーション装置３０の全体の動作を制御し、各種の情報処理を行う。ＣＰＵ３１は、ＲＯＭ３３または補助記憶装置３４に格納された、例えば後述するシミュレーション方法や、プログラム（シュミレーションプログラム）を実行して、粗視化粒子についてのパラメータである第２パラメータの算出や、粗視化粒子の挙動の解析を行うことができる。

ＲＡＭ３２は、ＣＰＵ３１のワークエリアとして用いられ、主要な制御パラメータや情報を記憶する不揮発ＲＡＭを含んでもよい。

ＲＯＭ３３は、プログラム（シュミレーションプログラム）等を記憶することができる。

補助記憶装置３４は、ＳＳＤ（ＳｏｌｉｄＳｔａｔｅＤｒｉｖｅ）や、ＨＤＤ（ＨａｒｄＤｉｓｋＤｒｉｖｅ）等の記憶装置であり、シミュレーション装置の動作に必要な各種のデータ、ファイル等を格納できる。

入出力インタフェース３５は、タッチパネル、キーボード、表示画面、操作ボタン等のユーザインタフェースと、外部のデータ収録サーバ等からの情報を取り込み、他の電子機器に解析情報を出力する通信インタフェースとの双方を含む。

表示装置３６は、モニタディスプレイ等である。表示装置３６では、解析画面が表示され、入出力インタフェース３５を介した入出力操作に応じて画面が更新される。

図３に示したシミュレーション装置３０の各機能は、例えばＲＡＭ３２やＲＯＭ３３等の主記憶装置または補助記憶装置３４からプログラム（シミュレーションプログラム）等を読み込ませ、ＣＰＵ３１により実行することにより、ＲＡＭ３２等におけるデータの読み出しおよび書き込みを行うと共に、入出力インタフェース３５および表示装置３６を動作させることで実現できる。

図４に、本実施形態のシミュレーション装置３０の機能ブロック図を示す。

図４に示すように、シミュレーション装置３０は、受付部４１、処理装置４２、出力部４３を有することができる。これらの各部は、シミュレーション装置３０が有するＣＰＵ、記憶装置、各種インタフェース等を備えたパーソナルコンピュータ等の情報処理装置において、ＣＰＵが予め記憶されている例えば後述するシミュレーション方法や、プログラムを実行することでソフトウェアおよびハードウェアが協働して実現される。

各部の構成について以下に説明する。
（Ａ）受付部
受付部４１は、処理装置４２で実行される処理に関係するユーザーからのコマンドやデータの入力を受け付ける。受付部４１としてはユーザーが操作を行い、コマンド等を入力するキーボードやマウス、ネットワークを介して入力を行う通信装置、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ等の各種記憶媒体から入力を行う読み取り装置などが挙げられる。
（Ｂ）処理装置
処理装置４２は、第１パラメータ取得部４２１、第２パラメータ算出部４２２、粗視化粒子挙動解析部４２３を有することができる。なお、処理装置は、必要に応じてさらに任意の部材を有することができ、例えば初期設定部等を有することもできる。
（Ｂ－１）第１パラメータ取得部
第１パラメータ取得部４２１では、例えば解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。第１パラメータは、粉体に関連するパラメータ以外に、解析に要する各種パラメータを含むこともできる。第１パラメータは解析（シミュレーション）の内容に応じて選択できるため、その具体的な種類は特に限定されない。第１パラメータとしては、離散要素法計算に必要となる各種パラメータが挙げられ、具体的には例えば粒子径、粒子数、ヤング率、計算のＴｉｍｅｓｔｅｐ、ポワソン比、壁面との摩擦係数、粒子間の摩擦係数、転がり摩擦係数、密度等から選択された１種類以上が挙げられる。

第１パラメータは、データベース等に収録されているデータであってもよく、予め実験を行い求めた実験値であってもよい。また、第１パラメータは、実験結果からシミュレーション等によりフィッティングして算出した計算値であってもよい。
（Ｂ－２）第２パラメータ算出部
既述のように、本実施形態のシミュレーション装置３０では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。ただし、粗視化粒子は、粗視化前の粒子群を構成する個々の粒子とは質量等の各種パラメータが異なる。このため、粗視化粒子について、計算を行う際に必要となるパラメータの算出や、設定を行う必要がある。

第２パラメータ算出部４２２では、「（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて」の中で説明したように、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解であるＫ_ｒを用いて、第２パラメータを算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、既述の垂直方向のＫ_ｒの特性方程式である式（１３）を導出し、式（１３）から垂直方向のＫ_ｒを算出できる。そして、既述の式（１３）に示した特性方程式の解である垂直方向のＫ_ｒを用いて、第２パラメータを算出できる。また、既述の式（１９）に示した接線方向のＫ_ｒの特性方程式を用いて、接線方向のＫ_ｒを算出でき、係る接線方向のＫ_ｒを用いて第２パラメータを算出することもできる。

既述のように、Ｋ_ｒは、粗視化粒子の挙動を支配するパラメータであり、Ｋ_ｒを用いることで、粗視化粒子の挙動に関連する各種パラメータを算出できる。

後述する粗視化粒子挙動解析部で用いる第２パラメータの種類は解析の内容に応じて選択できるため、特に限定されない。例えば第２パラメータは、粗視化粒子の熱伝導率を含むこともできる。この場合、第２パラメータ算出部は、既述の特性方程式の解Ｋ_ｒを用いて、熱伝導率を算出できる。
（Ｂ－３）粗視化粒子挙動解析部
粗視化粒子挙動解析部４２３では、第１パラメータ取得部４２１により取得した第１パラメータ、および第２パラメータ算出部４２２で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。

なお、ここでいう挙動とは、粗視化粒子の運動による位置の変化だけではなく、温度変化等の状態変化も含む。
（Ｂ－４）初期設定部
図示しない初期設定部は、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析部４２３で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得部４２１により取得する場合には、初期設定部は設けなくてもよい。
（Ｃ）出力部
出力部４３は、ディスプレイ等を有することができる。粗視化粒子挙動解析部４２３で得られたシミュレーション結果を出力部４３に出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部４３に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部４３に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。

以上に説明した本実施形態のシミュレーション装置によれば、複数の粒子を含む粉体の挙動をシミュレーションすることができ、その用途等は特に限定されない。例えば、キルン等の回転体内での粉体の挙動のシミュレーションに好適に用いることができる。すなわち、本実施形態のシミュレーション装置によれば、回転体内での、粉体の挙動を解析することもできる。

以上に説明した本実施形態のシミュレーション装置によれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。

そして、粗視化粒子のパラメータを既述のパラメータＫ_ｒを用いて計算するため、精度よく計算を行うことができる。
［シミュレーション方法］
次に、本実施形態のシミュレーション方法について説明する。本実施形態のシミュレーション方法は、例えば既述のシミュレーション装置を用いて実施できる。このため、既に説明した事項の一部は説明を省略する。

本実施形態のシミュレーション方法は、複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法に関する。本実施形態のシミュレーション方法は、図５に示したフロー図に従って実施することができ、以下の工程を有することができる。

粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程（Ｓ１）。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程（Ｓ２）。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析工程（Ｓ３）。
そして、第２パラメータ算出工程（Ｓ２）は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出できる。

各工程について以下に説明する。
（１）第１パラメータ取得工程（Ｓ１）
第１パラメータ取得工程（Ｓ１）では、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。既述のシミュレーション装置を用いる場合、例えば第１パラメータ取得部４２１において、第１パラメータ取得工程を実施できる。

第１パラメータは解析の内容に応じて選択できるため、その具体的な種類は特に限定されない。第１パラメータとしては、離散要素法計算に必要となる各種パラメータが挙げられる。第１パラメータの具体的な例はシミュレーション装置で既に説明したため、ここでは説明を省略する。

第１パラメータは、データベース等に収録されているデータであってもよく、予め実験を行い求めた実験値であってもよい。また、第１パラメータは、実験結果からシミュレーション等によりフィッティングして算出した計算値であってもよい。
（２）第２パラメータ算出工程（Ｓ２）
本実施形態のシミュレーション方法では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。

そこで、第２パラメータ算出工程（Ｓ２）では、「（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて」の中で説明したように、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解であるＫ_ｒを用いて、第２パラメータを算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、既述の垂直方向のＫ_ｒの特性方程式である式（１３）を導出し、式（１３）から垂直方向のＫ_ｒを算出できる。そして、既述の式（１３）に示した特性方程式の解である垂直方向のＫ_ｒを用いて、第２パラメータを算出できる。また、既述の式（１９）に示した接線方向のＫ_ｒの特性方程式を用いて、接線方向のＫ_ｒを算出でき、係る接線方向のＫ_ｒを用いて第２パラメータを算出することもできる。

既述のシミュレーション装置を用いる場合、例えば第２パラメータ算出部４２２において、第２パラメータ算出工程を実施できる。

後述する粗視化粒子挙動解析工程で用いる第２パラメータの種類は解析の内容に応じて選択できるため、特に限定されない。例えば第２パラメータは、粗視化粒子の熱伝導率を含むこともできる。この場合、第２パラメータ算出工程は、既述の特性方程式の解Ｋ_ｒを用いて、熱伝導率を算出できる。
（３）粗視化粒子挙動解析工程（Ｓ３）
粗視化粒子挙動解析工程（Ｓ３）では、第１パラメータ取得工程（Ｓ１）で得た第１パラメータ、および第２パラメータ算出工程（Ｓ２）で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。

なお、ここでいう挙動とは、粗視化粒子の運動による位置の変化だけではなく、温度変化等の状態変化も含む。
（４）初期設定工程
本実施形態のシミュレーション方法は、例えばさらに初期設定工程を有することもできる。初期設定工程では、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析工程で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得工程により取得する場合には、初期設定工程は実施しなくてもよい。
（５）出力工程
本実施形態のシミュレーション方法は、例えばさらに出力工程を有することができる。出力工程では、例えば粗視化粒子挙動解析工程（Ｓ３）で得られたシミュレーション結果を、出力部へ出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。

以上に説明した本実施形態のシミュレーション方法によれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。

そして、粗視化粒子のパラメータを既述のパラメータＫ_ｒを用いて計算するため、精度よく計算を行うことができる。
［プログラム］
次に、本実施形態のプログラムについて説明する。

本実施形態のプログラムは、複数の粒子を含む粉体の挙動を解析するためのプログラムに関し、コンピュータを以下の各部として機能させることができる。

粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部。

複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部。

第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析部。

第２パラメータ算出部では、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出させることができる。

また、第２パラメータは、粗視化粒子の熱伝導率を含むことができ、この場合、第２パラメータ算出部では、既述の特性方程式の解を用いて、熱伝導率を算出できる。

本実施形態のプログラムは、例えば既述のシミュレーション装置のＲＡＭやＲＯＭ等の主記憶装置または補助記憶装置の各種記憶媒体に記憶させておくことができる。そして、係るプログラムを読み込ませ、ＣＰＵにより実行することにより、ＲＡＭ等におけるデータの読み出しおよび書き込みを行うと共に、入出力インタフェースおよび表示装置を動作させて実行できる。このため、シミュレーション装置で既に説明した事項については説明を省略する。

上述した本実施形態のプログラムは、インターネットなどのネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることで提供してもよい。また、本実施形態のプログラムをインターネットなどのネットワークを介して提供、配布するように構成してもよい。

本実施形態のプログラムは、ＣＤ－ＲＯＭ等の光ディスクや、半導体メモリ等の記録媒体に格納した状態で流通等させてもよい。

以上に説明した本実施形態のプログラムによれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。

そして、粗視化粒子のパラメータを既述のパラメータＫ_ｒを用いて計算するため、精度よく計算を行うことができる。
２．第２の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
第２の実施形態では、接線方向の運動方程式に関して、オーバーラップ量を算出する際に、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定する点がここまで説明した第１の実施形態の場合と異なる。また、回転方向の運動方程式に関して、上述のようにして求めたオーバーラップ量を用いてトルクを算出できる。これにより計算量を抑制したまま、粗視化後の粒子群の挙動についてより精度よく解析可能になる。
（１－１）オーバーラップ量について
通常、接線方向のオーバーラップ量は接触開始から接触終了までの速度の接線方向の成分ｖ_ｔ（以下の式中ではｖ_ｔの上に矢印で示されたもの）の時間積分を用いて以下のように求められる。ここでｔは時間、ベクトルｒ（以下の式中ではｒの上に矢印で示されたもの）は、粒子群１１を構成する粒子の中心から接触点までのベクトル、ベクトルω（以下の式中ではωの上に矢印で示されたもの）は粒子群１１の回転ベクトルを示す。

なお、以下の式中のハットｔ（以下の式中ではｔの上にハットを付して示されたもの）は、接線方向オーバーラップの単位ベクトルを示す。

このため、粒子と壁面との間での各オーバーラップ量は以下の式（２９）で表される。なお、以下の式中の添え字のｔは接線方向成分を意味する。

また、粒子間のオーバーラップ量は以下の式（３０）で表される。

一方、粗視化粒子２１のオーバーラップ量は以下の式（３１）で表される。なお、ベクトルω_ｃw（以下の式中ではω_ｃｗの上に矢印で示されたもの）は粗視化粒子２１の回転ベクトルである。

ここで、垂直方向における運動方程式である式（１１）の場合と同様に、粗視化前の粒子群１１と、粗視化粒子２１とで重心位置が一致しているとすると、以下の式（３２）の関係を満たすことになる。

ここで、粒子の回転を考慮する場合、重心移動に対して回転方向の自由度が残る。仮に重心位置が一致している場合であっても、式（３１）と式（３２）が一致しない場合がある。このため、粒子群１１の回転運動と粗視化粒子２１の回転運動を換算する必要がある。そこで、回転に関しては角運動量と回転エネルギーが粗視化前後で一致すると仮定する。この場合、角運動量は以下のような等式が成り立つ。

さらに回転運動エネルギーは以下のような等式が成り立つ。

上記式（３３）、式（３４）中、右辺の第１項、第２項は、それぞれ粗視化前の粒子群の自転運動成分（Ｓｐｉｎ）と、公転運動成分（Ｏｒｂｉｔ）とを意味する。ベクトルω_ｃｗ、ベクトルω_ｓ、ベクトルω_ｏ（それぞれの式（３３）、式（３４）中ではω_ｃｗ、ω_ｓ、ω_ｏの上に矢印で示されたもの）はそれぞれ粗視化粒子の角速度、粗視化前の粒子群の自転運動の角速度、粗視化前の粒子群の公転運動の角速度を意味する。また、Ｉ_ｃｗ、Ｉ_ｓ、Ｉ_ｏは粗視化粒子の慣性モーメント、粗視化前の粒子群の自転運動の慣性モーメント、粗視化前の粒子群の公転運動の慣性モーメントを意味する。

上記式（３３）、式（３４）を連立して公転運動成分を消去すると、以下の式（３５）になる。

これより、粗視化粒子の接線方向のオーバーラップ量を以下のように定義できる。

この定義の下で、垂直方向の場合と同様に以下の式（３７）が成り立つ。

ここでは計算の簡便化のため式（３３）、式（３４）、式（３５）としたが、粗視化粒子の形状によってはいくつかの形式が考えられる。例えば立方体状とした粗視化粒子の形状因子を考慮して慣性モーメントを算出しようとすると以下の式（３８）～式（４１）となる。

このため、角運動量保存とベクトルω_ｓ、ω_ｏについて、ω_ｓ＝ω_ｏを仮定すると式（３８）より以下の式（４２）、式（４３）となる。

また、あるいは同様に回転の運動エネルギーが保存することとベクトルω_ｓ、ω_ｏについて、ω_ｓ＝ω_ｏを仮定するなどしてベクトルω_ｓ、ベクトルω_ｃｗの関係を求めてもよい。計算対象とする現象に合わせてこれらを選択することができる。

このようにして、粗視化前の粒子群を構成する粒子のオーバーラップ量も求めることができる。
（１－２）回転方向の運動方程式について
ここまで説明した接線方向のオーバーラップ量から接線方向の力を算出し、トルクを計算できる。各粒子間や、粒子と壁間での転がり摩擦は粒子にかかる垂直抗力と接触半径の積に比例するトルクが発生するとした。それぞれの転がり摩擦は各粒子で発生するので、粗視化粒子全体での転がり摩擦抵抗は以下の式（４４）のように表せる。

式（４４）中、ベクトルＴ_{ｔｏｔ＿ｆｒｉｃ}が粗視化粒子全体の転がり摩擦抵抗を、ベクトルＴ_{ｗ，ｆｒｉｃ}が粒子群１１の壁面１２あるいは外部粒子との間の転がり摩擦抵抗を、ベクトルＴ_{ｐ，ｆｒｉｃ}が粒子群１１の内部粒子間の転がり摩擦抵抗を表している。各ベクトルは、以下の式中では文字の上に矢印を付して表記している。なお、粒子群１１では粒子－粒子間に作用反作用により逆向きの接触力が働くため、式（４４）に示した様に、ベクトルＴ_{ｐ，ｆｒｉｃ}は２倍となる。

そして、ベクトルＴ_{ｗ，ｆｒｉｃ}と、ベクトルＴ_{ｐ，ｆｒｉｃ}とは、式（４５）、式（４６）により表される。なお、式中のハットωは、回転方向の単位ベクトルを意味する。

式（４５）、式（４６）中の接触半径ｒ_ｗ、ｒ_ｐはオーバーラップ量から幾何学的に算出することができ、既述の接線方向の運動方程式の場合の定義に基づいて算出したオーバーラップ量を用いることで精度を高めることができる。また、μ_ｗは粒子群１１の壁面１２あるいは外部粒子との転がり摩擦係数を、μ_ｐは粒子群１１の内部粒子の転がり摩擦係数をそれぞれ意味する。

そして、回転方向の運動方程式は以下のよう表記できる。

（２）シミュレーション装置
本実施形態のシミュレーション装置においても、第２パラメータ算出部は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。

以上の点以外は第１の実施形態のシミュレーション装置の場合と同様に構成できるため、ここでは説明を省略する。
［シミュレーション方法］
本実施形態のシミュレーション方法においても、第２パラメータ算出工程では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。

以上の点以外は第１の実施形態のシミュレーション方法の場合と同様に構成できるため、ここでは説明を省略する。
［プログラム］
本実施形態のプログラムにおいても、第２パラメータ算出部では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。

以上の点以外は第１の実施形態のプログラムの場合と同様に構成できるため、ここでは説明を省略する。

以下に具体的な実施例を挙げて説明するが、本発明はこれらの実施例に限定されるものではない。
［実験例１］
［比較例１－１］
直方体の容器内に充填した粉体層の温度変化について、表１に示したパラメータを用いて、離散要素法計算により底部から加熱を行った場合の解析を行った。解析により求めた容器内の粉体層の平均温度の変化を図６に示す。なお、反発係数を０．１、摩擦係数を０．７、転がり摩擦係数を０．００１とした。

［比較例１－２］
粗視化倍率αを２とし、８個の粒子から構成される粒子群を１個の粗視化粒子とした点以外は比較例１－１と同様の条件で、直方体の容器内に充填した粉体層の温度変化について解析を行った。なお、粗視化粒子としたため、表１に示すように、粒子の粒径を変更しているが、粒径以外は比較例１－１と同じパラメータを用いて、解析を行った。解析により求めた容器内の粉体層の平均温度の変化を図６に示す。
［実施例１－１］
粗視化倍率αを２とし、既述のシミュレーション装置を用いて、直方体の容器内に充填した粉体層の温度変化について解析を行った。

具体的には、第１パラメータ取得部４２１により、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得した（第１パラメータ取得工程：Ｓ１）。

この際には、比較例１－２と同じパラメータを取得した。

次いで、第２パラメータ算出部４２２により、粗視化粒子についてのパラメータである第２パラメータを算出した（第２パラメータ算出工程：Ｓ２）。

この際、既述の式（１３）、式（１９）に示した特性方程式の解であるパラメータＫ_ｒを用いた式（２６）、式（２７）により、粗視化粒子に関する熱伝導率を算出した。また、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出した。既述のように、計算に適用するモデルに応じて、これらの係数は調整可能であり、上述のように特性方程式を用いて算出、換算できる。算出したパラメータＫ_ｒを表２に示す。

以上の取得、算出した第１パラメータ、および第２パラメータを表１に示す。そして、表１に示したパラメータに基づいて、粗視化粒子の挙動、具体的には温度変化を粗視化粒子挙動解析部４２３において解析した（粗視化粒子挙動解析工程：Ｓ３）。結果を図６に示す。

図６に示した結果によると、実施例１－１の結果は、粗視化を行っていない、比較例１－１の結果とほぼ一致しており、粗視化粒子のパラメータを適切に選択、設定し、解析を行えていることを確認できた。

実施例１－１では粗視化を行っているため、粒子の数が比較例１－１と比較して１／８倍になっており、比較例１－１と比較して計算量を抑制できている。
［実験例２］
［比較例２－１］
回転体であるキルン内の粉体の挙動について、表３に示したパラメータを用いて、離散要素法計算により解析を行った。なお、反発係数を０．７５、摩擦係数を０．３、転がり摩擦係数を０．５とした。
解析により求めたキルン内の粉体の動きを図８に、解析により求めたキルン内の粉体の温度分布を図１１に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。

図８（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図８（Ｂ）、図８（Ｃ）、図８（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。

図１１（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１１（Ｂ）、図１１（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。

［比較例２－２］
粗視化倍率αを４とし、６４個の粒子から構成される粒子群を１個の粗視化粒子とした点以外は比較例２－１と同様の条件で、回転体であるキルン内の粉体の挙動について解析を行った。なお、表３に示すように、粒径以外は比較例２－１と同じパラメータを用いて、解析を行った。解析により求めたキルン内の粉体の動きを図９に、解析により求めたキルン内の粉体の温度分布を図１２に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。

図９（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図９（Ｂ）、図９（Ｃ）、図９（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。

図１２（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１２（Ｂ）、図１２（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。
［実施例２－１］
粗視化倍率αを４とし、既述のシミュレーション装置を用いて、回転体であるキルン内の粉体の挙動について解析を行った。

この際には、比較例２－２と同じパラメータを取得した。

この際、既述の式（１３）、式（１９）に示した特性方程式の解であるパラメータＫ_ｒを用いた式（２６）、式（２７）により、粗視化粒子に関する熱伝導率を算出した。また、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出した。算出したパラメータＫ_ｒを表４に示す。

以上のようにして取得、算出した、第１パラメータ、および第２パラメータを表３に示す。そして、表３に示したパラメータに基づいて、粗視化粒子の挙動、具体的にはキルン内の動きや、温度変化を粗視化粒子挙動解析部４２３において解析した（粗視化粒子挙動解析工程：Ｓ３）。解析により求めたキルン内の粉体の動きを図７に、解析により求めたキルン内の粉体の温度分布を図１０に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。

図７（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図７（Ｂ）、図７（Ｃ）、図７（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。

図１０（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１０（Ｂ）、図１０（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。

図７と図８の比較や、図１０と図１１の比較から明らかなように、実施例２－１の結果は、粗視化を行っていない、比較例２－１の結果とほぼ一致しており、粗視化粒子のパラメータを適切に選択、設定し、解析を行えていることを確認できた。

図１３においても実施例２－１と比較例２－１とは変曲点を有するなど同様の傾向を示すことが確認できた。

そして、実施例２－１では粗視化を行っているため、粒子の数が比較例２－１と比較して１／６４倍になっており、比較例２－１と比較して計算量を抑制できている。

１１粒子群
１１Ａ、１１Ｂ粒子
２１粗視化粒子
３０シミュレーション装置
４２１第１パラメータ取得部
４２２第２パラメータ算出部
４２３粗視化粒子挙動解析部
Ｓ１第１パラメータ取得工程
Ｓ２第２パラメータ算出工程
Ｓ３粗視化粒子挙動解析工程

Claims

複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置。
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置。
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、
前記第２パラメータ算出部は、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項１または請求項２に記載のシミュレーション装置。
回転体内での、前記粉体の挙動を解析する請求項１から請求項３のいずれか１項に記載のシミュレーション装置。
複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析工程と、を有し、
前記第２パラメータ算出工程は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション方法。
複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析工程と、を有し、
前記第２パラメータ算出工程は、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション方法。
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、
前記第２パラメータ算出工程では、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項５または請求項６に記載のシミュレーション方法。
複数の粒子を含む粉体の挙動を解析するためのプログラムであって、
コンピュータを、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、して機能させ、
前記第２パラメータ算出部では、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出させるプログラム。
複数の粒子を含む粉体の挙動を解析するためのプログラムであって、
コンピュータを、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、して機能させ、
前記第２パラメータ算出部では、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出させるプログラム。
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、
前記第２パラメータ算出部は、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項８または請求項９に記載のプログラム。