ホーム > 特許ランキング > 住友金属鉱山株式会社
知財求人 - 知財ポータルサイト「IP Force」
▶ 住友金属鉱山株式会社の特許一覧 ▶ 公立大学法人大阪の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-07-07
(45)【発行日】2022-07-15
(54)【発明の名称】シミュレーション装置、シミュレーション方法、プログラム
(51)【国際特許分類】
G01N 15/00 20060101AFI20220708BHJP
G16Z 99/00 20190101ALI20220708BHJP
【ＦＩ】
G01N15/00 Z
G16Z99/00
【請求項の数】
10
(21)【出願番号】P 2020162719
(22)【出願日】2020-09-28
(65)【公開番号】P2021190060
(43)【公開日】2021-12-13
【審査請求日】2022-01-19
(31)【優先権主張番号】P 2020095620
(32)【優先日】2020-06-01
(33)【優先権主張国・地域又は機関】JP
【早期審査対象出願】
(73)【特許権者】
【識別番号】000183303
【氏名又は名称】住友金属鉱山株式会社
(73)【特許権者】
【識別番号】519135633
【氏名又は名称】公立大学法人大阪
(74)【代理人】
【識別番号】100107766
【弁理士】
【氏名又は名称】伊東 忠重
(74)【代理人】
【識別番号】100070150
【弁理士】
【氏名又は名称】伊東 忠彦
(72)【発明者】
【氏名】猿渡 元彬
(72)【発明者】
【氏名】仲村 英也
【審査官】北条 弥作子
(56)【参考文献】
【文献】特開２０２０－０５７１３５（ＪＰ，Ａ）
【文献】特開２０２０－０６４４５０（ＪＰ，Ａ）
【文献】特開２０００－３２２４０７（ＪＰ，Ａ）
【文献】特開２０１３－２１０９１８（ＪＰ，Ａ）
【文献】中国特許出願公開第１０４６３６５３８（ＣＮ，Ａ）
【文献】山井 三亀夫，スケール則に基づくDEMシミュレーション，J. Soc. Powder Technol，日本，2018年，55，pp.95-103，doi: 10.4164/sptj.55.95
【文献】A.V. Patil，Comparison of CFD-DEM heat transfer simulations with infrared/visual measurements，Chemical Engineering Journal，2015年，277，pp.388-401
【文献】瀬戸内 秀規，球要素間の回転剛性を導入した個別要素モデル，土木学会論文集Ａ２（応用力学），2012年，Vol.68, No.1，pp.18-29
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｎ １５／００－１５／１４
Ｇ１６Ｚ ９９／００
ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置。
【請求項2】
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置。
【請求項3】
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、前記第２パラメータ算出部は、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項１または請求項２に記載のシミュレーション装置。
【請求項4】
回転体内での、前記粉体の挙動を解析する請求項１から請求項３のいずれか１項に記載のシミュレーション装置。
【請求項5】
複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法であって、前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析工程と、を有し、
前記第２パラメータ算出工程は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション方法。
【請求項6】
複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析工程と、を有し、
前記第２パラメータ算出工程は、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション方法。
【請求項7】
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、前記第２パラメータ算出工程では、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項５または請求項６に記載のシミュレーション方法。
【請求項8】
複数の粒子を含む粉体の挙動を解析するためのプログラムであって、コンピュータを、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、して機能させ、
前記第２パラメータ算出部では、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出させるプログラム。
【請求項9】
複数の粒子を含む粉体の挙動を解析するためのプログラムであって、
コンピュータを、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、して機能させ、
前記第２パラメータ算出部では、前記粒子群と前記粗視化粒子とで、重心、角運動量、および回転エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出させるプログラム。
【請求項10】
前記第２パラメータは、前記粗視化粒子の熱伝導率を含んでおり、前記第２パラメータ算出部は、前記特性方程式の解を用いて、前記熱伝導率を算出する請求項８または請求項９に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、シミュレーション装置、シミュレーション方法、プログラム
に関する。
に関する。
【背景技術】
【0002】
特許文献１には、シミュレーション条件の入力を行う入力装置と、
シミュレーション結果を出力する出力装置と、
前記入力装置から入力されたシミュレーション条件に基づいて、大きさが異なる複数の
粒子を含む粉粒体の挙動を解析する処理装置とを有し、
前記処理装置は、
前記入力装置から入力されたシミュレーション対象の粉粒体の粒径分布を規定するパラメータの値、及び粒子を粗視化する基準となる粗視化係数の値に基づいて、粗視化された粉粒体の挙動をシミュレーションにより求め、
シミュレーションによって求められた粒子の挙動と、入力された粗視化係数の値とを関連付けて前記出力装置に出力するシミュレーション装置等に開示されているように、粉粒体についてのシミュレーション装置が開示されている。
シミュレーション結果を出力する出力装置と、
前記入力装置から入力されたシミュレーション条件に基づいて、大きさが異なる複数の
粒子を含む粉粒体の挙動を解析する処理装置とを有し、
前記処理装置は、
前記入力装置から入力されたシミュレーション対象の粉粒体の粒径分布を規定するパラメータの値、及び粒子を粗視化する基準となる粗視化係数の値に基づいて、粗視化された粉粒体の挙動をシミュレーションにより求め、
シミュレーションによって求められた粒子の挙動と、入力された粗視化係数の値とを関連付けて前記出力装置に出力するシミュレーション装置等に開示されているように、粉粒体についてのシミュレーション装置が開示されている。
【先行技術文献】
【特許文献】
【0003】
【文献】特開２０２０－５７１３５号公報
【発明の概要】
【発明が解決しようとする課題】
【0004】
工場等におけるプロセスの改善や、製造工程を検討する際の試験工数の削減等を目的として、複数の粒子を含む粉体（粉粒体）の挙動の解析を離散要素法（ＤＥＭ：Ｄｉｓｃｒｅｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ）計算等により行うことが従来からなされてきた。
【0005】
離散要素法計算は、個々の粒子について運動方程式を解くことで、粉体全体の運動を記述するシミュレーション技術である。
【0006】
しかしながら、離散要素法計算では、取り扱う粒子数が多くなるほど計算負荷が高くなる。このため、例えば工場で用いるプラント等のような大きな規模での粉体の挙動について解析を行う場合、計算量が膨大になり、現実的には計算を行うことが困難となる。
【0007】
そこで、複数個の粒子から構成される粒子群を１個の粗視化粒子とする粗視化の手法を用いたシミュレーション装置について従来から検討されていた（特許文献１を参照）。粗視化の手法を用いたシミュレーション装置では、正確な解析結果を得るため、計算で用いる粗視化粒子に関するパラメータを適切に選択することが求められる。このため、粗視化粒子に関するパラメータを新たな手法により選択、設定し、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置が求められていた。
【0008】
上記従来技術が有する問題に鑑み、本発明の一側面では、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【0009】
上記課題を解決するため本発明の一態様によれば、
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置を提供する。
複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であって、
前記粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部と、
複数個の前記粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、前記粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部と、
前記第１パラメータ、および前記第２パラメータに基づいて、前記粗視化粒子の挙動を解析する粗視化粒子挙動解析部と、を有し、
前記第２パラメータ算出部は、前記粒子群と前記粗視化粒子とで、重心および弾性エネルギーが一致しているとして、前記粒子群の弾性エネルギーと、前記粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、前記第２パラメータを算出するシミュレーション装置を提供する。
【発明の効果】
【0010】
本発明の一態様によれば、複数の粒子を含む粉体の挙動を解析できる新たなシミュレーション装置を提供することができる。
【図面の簡単な説明】
【0011】
【図1】複数個の粒子から構成される粒子群と、粒子群と壁面との衝突時の説明図である。
【図2】複数個の粒子から構成される粒子群を粗視化した粗視化粒子と、粗視化粒子と壁面との衝突時の説明図である。
【図3】本発明の実施形態に係るシミュレーション装置のハードウェア構成図である。
【図4】本発明の実施形態に係るシミュレーション装置の機能を示すブロック図である。
【図5】本発明の実施形態に係るシミュレーション方法を説明するフローチャートである。
【図6】実験例１における粉体の平均温度の変化を示すグラフである。
【図7】実施例２－１におけるキルン内での粉体の混合の状態を示す図である。
【図8】比較例２－１におけるキルン内での粉体の混合の状態を示す図である。
【図9】比較例２－２におけるキルン内での粉体の混合の状態を示す図である。
【図10】実施例２－１におけるキルン内の粉体の温度分布を示す図である。
【図11】比較例２－１におけるキルン内の粉体の温度分布を示す図である。
【図12】比較例２－２におけるキルン内の粉体の温度分布を示す図である。
【図13】実験例２における粉体の平均温度の変化を示すグラフである。
【発明を実施するための形態】
【0012】
本開示の一実施形態（以下「本実施形態」と記す）に係るシミュレーション装置、シミュレーション方法、プログラムの具体例を、以下に図面を参照しながら説明する。なお、本発明はこれらの例示に限定されるものではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内での全ての変更が含まれることが意図される。
１．第１の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
（１－１）粒子の粗視化について
本実施形態のシミュレーション装置の詳細について説明する前に、本実施形態のシミュレーション装置で用いることができる、複数個の粒子から構成される粒子群の粗視化、および粗視化した粒子である粗視化粒子に関連するパラメータの算出方法について以下に説明する。
１．第１の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
（１－１）粒子の粗視化について
本実施形態のシミュレーション装置の詳細について説明する前に、本実施形態のシミュレーション装置で用いることができる、複数個の粒子から構成される粒子群の粗視化、および粗視化した粒子である粗視化粒子に関連するパラメータの算出方法について以下に説明する。
【0013】
既述のように、離散要素法計算では、取り扱う粒子数が多くなるほど計算負荷が高くなる。このため、例えば工場で用いるプラント等のような大きな規模での粉体の挙動について解析を行う場合、計算量が膨大になり、現実的には計算を行うことが困難となる。
【0014】
そこで、多量の粒子を含む粉体について挙動を解析する場合、計算量を抑制するため、例えば図１（Ａ）に示した複数個の粒子から構成される粒子群１１を、図２（Ａ）に示すような１個の大きな粒子である粗視化粒子２１として取り扱う粗視化の技術が必要となる。
【0015】
ただし、粗視化前の個別の粒子と、粗視化粒子とでは比表面積等が異なるため、計算に要する一部のパラメータは変化することになる。このため、粗視化粒子のパラメータを適切に決定する必要がある。
（１－２）粗視化粒子の粒子挙動の計算に用いるパラメータについて
そこで、本発明の発明者は、粗視化粒子に関するパラメータを決定する方法について検討を行った。計算にあたって、図１（Ａ）に示した粗視化する前の複数個の粒子から構成される粒子群１１が壁面１２に衝突する場合と、図２（Ａ）に示した粗視化粒子２１が壁面１２に衝突する場合とをモデルに用いた。以下の説明では壁面と粒子とが衝突する場合を例に粗視化粒子に関するパラメータを決定する方法を記載するが、粒子同士が衝突する場合でも同様の議論となるため説明を省略する。
（１－２）粗視化粒子の粒子挙動の計算に用いるパラメータについて
そこで、本発明の発明者は、粗視化粒子に関するパラメータを決定する方法について検討を行った。計算にあたって、図１（Ａ）に示した粗視化する前の複数個の粒子から構成される粒子群１１が壁面１２に衝突する場合と、図２（Ａ）に示した粗視化粒子２１が壁面１２に衝突する場合とをモデルに用いた。以下の説明では壁面と粒子とが衝突する場合を例に粗視化粒子に関するパラメータを決定する方法を記載するが、粒子同士が衝突する場合でも同様の議論となるため説明を省略する。
【0016】
図１（Ａ）に示すように、複数個の粒子から構成される粒子群１１が、立方体状に、縦方向、横方向、高さ方向に２個ずつ、合計で２３個並んでいるとする。後述するように、係る８個の粒子をまとめて１個の粗視化粒子とする場合、一辺の方向に並べた粒子の数、すなわち２を粗視化倍率とする。
【0017】
図１（Ａ）に示した複数個の粒子１１Ａ、１１Ｂを含む粒子群１１が壁面１２に衝突する際に、粒子群１１のうち壁面１２側に位置する粒子１１Ａが壁面あるいは外部粒子から受ける力を図１（Ｂ）に示すようにＦｗとする。また、図１（Ｂ）に示すように粒子１１Ａの壁面１２あるいは外部粒子とのオーバーラップ量をδｗ、粒子１１Ｂの隣接する粒子１１Ａとのオーバーラップ量をδｐとする。なお、図１（Ｂ）は、粒子群１１が壁面１２に衝突する際の様子を側面側から見た図になる。
【0018】
この場合、粒子群１１に加わる力の大きさは、以下の式（１）で表すことができる。
【0019】
なお、式（１）におけるαは粗視化倍率であり、既述のように粒子群１１を１個の粗視化粒子とした場合に一辺の方向に並べた粒子の数を意味する。図１（Ａ）に示した粒子群１１を図２（Ａ）に示した１個の粗視化粒子２１とする場合、α＝２となる。
【0020】
また、ｍは各粒子１１Ａ、１１Ｂの質量を、ａＧは粒子群１１の重心の加速度を、ηｗは壁面１２あるいは外部粒子と粒子１１Ａとの反発係数から算出される粘性係数をそれぞれ意味している。作用反作用の法則により粒子間の接触力が打ち消されるため、壁面１２とは直接接しない粒子１１Ｂが、隣接する粒子１１Ａから受ける力Ｆｐは式（１）には表れないことになる。
【0021】
【数1】
【0022】
【数2】
【0023】
式（２）中のＦｃｗは、図２（Ｂ）に示すように粗視化粒子２１が壁面１２あるいは外部粒子から受ける力を、δｃｗは粗視化粒子２１の壁面１２あるいは外部粒子とのオーバーラップ量を、ηｃｗは粗視化粒子２１と壁面１２あるいは外部粒子との反発係数から算出される粘性係数をそれぞれ意味している。なお、図２（Ｂ）は、粗視化粒子２１が壁面１２に衝突する際の様子を側面側から見た図になる。
【0024】
【数3】
【0025】
従って、粒子群１１について計算を行った既述の式（１）と、粒子群を粗視化した粗視化粒子２１について計算を行った既述の式（２）とから、対応するパラメータが一致することを示す以下の式（３）、式（４）が導出される。
【0026】
【数4】
【0027】
【数5】
【0028】
【数6】
【0029】
【数7】
【0030】
【数8】
【0031】
【数9】
【0032】
【数10】
【0033】
【数11】
【0034】
【数12】
【0035】
そして、Ｋｒは、衝突中の粗視化前の粒子群１１の弾性エネルギーと、粗視化粒子２１の弾性エネルギーとの関係を用いた特性方程式により算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、特性方程式を作成してＫｒを算出できる。
【0036】
壁面１２と衝突中の、粒子群１１の弾性エネルギーと、粗視化粒子の弾性エネルギーは、既述の粒子群１１を構成する粒子１１Ａ、１１Ｂに加わる力や、粗視化粒子に加わる力を表した式（５）～式（７）をオーバーラップ量の距離で積分することで算出できる。
【0037】
従って、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーとを用いて、以下の式（１２）が得られる
【0038】
【数13】
【0039】
【数14】
【0040】
ここまで、壁面１２に対する垂直方向の運動方程式を用いて説明したが、接線方向の運動方程式や、回転の運動方程式についても同様のことが言える。
【0041】
具体的には、接線方向の運動方程式は、式（１４）で表すことができる。
【0042】
この場合も、式（１５）に示したように、Ｋｒを設定すると、δｗ、δｐは式（１６）、式（１７）のように示すことができ、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーが等しいとすると、式（１８）が得られる。式（１８）を変形することで、接線方向の特性方程式である式（１９）が得られる。ただし、接線方向の接触モデルには線形バネモデルを使用した。このように、接触モデルによって弾性エネルギーの計算式が異なるが、必要に応じて特性方程式を変更する等により、適切に弾性エネルギーを算出できる。
【0043】
【数15】
【0044】
【数16】
【0045】
【数17】
【0046】
【数18】
【0047】
【数19】
【0048】
【数20】
（１－３）粗視化粒子の伝熱に関するパラメータについて
既述のＫｒを用いて、粗視化粒子の伝熱に関するパラメータである熱伝導率を求めることもできる。
【0049】
粒子の伝熱は熱伝導率を用いて、以下の式（２０）、式（２１）のようにあらわされる。
【0050】
なお、式（２０）、式（２１）において、Ｑは熱流量を、ｈは熱流量係数を、ΔＴは壁面あるいは粒子同士の温度差を、ｋｗは壁面１２あるいは外部粒子と粒子１１Ａとの熱伝導率を、ｋｐは粒子群１１内部の熱伝導率を、ａは粒子群１１の接触半径をそれぞれ意味する。
【0051】
【数21】
【0052】
【数22】
【0053】
接触半径は、各粒子のオーバーラップ量に依存する。そして、既述のように粗視化前後でのオーバーラップ量は既述の垂直方向の特性方程式の解に関係づけられている。このことから、垂直方向の特性方程式の解Ｋｒを用いて伝熱方程式を書き換えると以下の式（２２）、式（２３）のようになる。
【0054】
式（２３）におけるａ´は粗視化粒子の接触半径を、ｒは粒子群１１を構成する粒子の半径をそれぞれ意味する。
【0055】
【数23】
【0056】
【数24】
【0057】
式（２４）、式（２５）におけるＱｃは粗視化粒子の熱流量を、ｈｃは粗視化粒子の熱流量係数を、ΔＴｃは粗視化粒子２１と壁面１２あるいは外部粒子との温度差を、ｋｗ´は壁面１２と粗視化粒子２１との熱伝導率を、ｋｐ´は外部粒子と粗視化粒子との熱伝導率をそれぞれ意味する。
【0058】
【数25】
【0059】
【数26】
【0060】
【数27】
【0061】
【数28】
【0062】
なお、ここでは粗視化前後での熱伝導率を例に説明したが、熱伝導率以外のパラメータについても同様にして、既述の特性方程式の解Ｋｒを用いて粗視化後のパラメータを算出できる。例えば、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出することもできる。なお、計算に適用するモデルに応じて、これらの係数は調整可能であり、上述のように特性方程式を用いて算出、換算できる。
（２）シミュレーション装置
本実施形態のシミュレーション装置は、複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であり、以下の部材を有することができる。
粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析部。
（２）シミュレーション装置
本実施形態のシミュレーション装置は、複数の粒子を含む粉体の挙動を解析するためのシミュレーション装置であり、以下の部材を有することができる。
粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析部。
【0063】
第２パラメータ算出部は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出する。
【0064】
図３に示したハードウェア構成図に示すように、本実施形態のシミュレーション装置３０は、例えば、情報処理装置（コンピュータ）で構成され、物理的には、演算処理部であるＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ：プロセッサ）３１と、主記憶装置であるＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）３２やＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）３３と、補助記憶装置３４と、入出力インタフェース３５と、出力装置である表示装置３６等を含むコンピュータシステムとして構成することができる。これらは、バス３７で相互に接続されている。なお、補助記憶装置３４や表示装置３６は、外部に設けられていてもよい。
【0065】
ＣＰＵ３１は、シミュレーション装置３０の全体の動作を制御し、各種の情報処理を行う。ＣＰＵ３１は、ＲＯＭ３３または補助記憶装置３４に格納された、例えば後述するシミュレーション方法や、プログラム（シュミレーションプログラム）を実行して、粗視化粒子についてのパラメータである第２パラメータの算出や、粗視化粒子の挙動の解析を行うことができる。
【0066】
ＲＡＭ３２は、ＣＰＵ３１のワークエリアとして用いられ、主要な制御パラメータや情報を記憶する不揮発ＲＡＭを含んでもよい。
【0067】
ＲＯＭ３３は、プログラム（シュミレーションプログラム）等を記憶することができる。
【0068】
補助記憶装置３４は、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）や、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）等の記憶装置であり、シミュレーション装置の動作に必要な各種のデータ、ファイル等を格納できる。
【0069】
入出力インタフェース３５は、タッチパネル、キーボード、表示画面、操作ボタン等のユーザインタフェースと、外部のデータ収録サーバ等からの情報を取り込み、他の電子機器に解析情報を出力する通信インタフェースとの双方を含む。
【0070】
表示装置３６は、モニタディスプレイ等である。表示装置３６では、解析画面が表示され、入出力インタフェース３５を介した入出力操作に応じて画面が更新される。
【0071】
図３に示したシミュレーション装置３０の各機能は、例えばＲＡＭ３２やＲＯＭ３３等の主記憶装置または補助記憶装置３４からプログラム（シミュレーションプログラム）等を読み込ませ、ＣＰＵ３１により実行することにより、ＲＡＭ３２等におけるデータの読み出しおよび書き込みを行うと共に、入出力インタフェース３５および表示装置３６を動作させることで実現できる。
【0072】
図４に、本実施形態のシミュレーション装置３０の機能ブロック図を示す。
【0073】
図４に示すように、シミュレーション装置３０は、受付部４１、処理装置４２、出力部４３を有することができる。これらの各部は、シミュレーション装置３０が有するＣＰＵ、記憶装置、各種インタフェース等を備えたパーソナルコンピュータ等の情報処理装置において、ＣＰＵが予め記憶されている例えば後述するシミュレーション方法や、プログラムを実行することでソフトウェアおよびハードウェアが協働して実現される。
【0074】
各部の構成について以下に説明する。
（Ａ）受付部
受付部４１は、処理装置４２で実行される処理に関係するユーザーからのコマンドやデータの入力を受け付ける。受付部４１としてはユーザーが操作を行い、コマンド等を入力するキーボードやマウス、ネットワークを介して入力を行う通信装置、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ等の各種記憶媒体から入力を行う読み取り装置などが挙げられる。
（Ｂ）処理装置
処理装置４２は、第１パラメータ取得部４２１、第２パラメータ算出部４２２、粗視化粒子挙動解析部４２３を有することができる。なお、処理装置は、必要に応じてさらに任意の部材を有することができ、例えば初期設定部等を有することもできる。
（Ｂ－１）第１パラメータ取得部
第１パラメータ取得部４２１では、例えば解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。第１パラメータは、粉体に関連するパラメータ以外に、解析に要する各種パラメータを含むこともできる。第１パラメータは解析（シミュレーション）の内容に応じて選択できるため、その具体的な種類は特に限定されない。第１パラメータとしては、離散要素法計算に必要となる各種パラメータが挙げられ、具体的には例えば粒子径、粒子数、ヤング率、計算のＴｉｍｅ ｓｔｅｐ、ポワソン比、壁面との摩擦係数、粒子間の摩擦係数、転がり摩擦係数、密度等から選択された１種類以上が挙げられる。
（Ａ）受付部
受付部４１は、処理装置４２で実行される処理に関係するユーザーからのコマンドやデータの入力を受け付ける。受付部４１としてはユーザーが操作を行い、コマンド等を入力するキーボードやマウス、ネットワークを介して入力を行う通信装置、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ等の各種記憶媒体から入力を行う読み取り装置などが挙げられる。
（Ｂ）処理装置
処理装置４２は、第１パラメータ取得部４２１、第２パラメータ算出部４２２、粗視化粒子挙動解析部４２３を有することができる。なお、処理装置は、必要に応じてさらに任意の部材を有することができ、例えば初期設定部等を有することもできる。
（Ｂ－１）第１パラメータ取得部
第１パラメータ取得部４２１では、例えば解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。第１パラメータは、粉体に関連するパラメータ以外に、解析に要する各種パラメータを含むこともできる。第１パラメータは解析（シミュレーション）の内容に応じて選択できるため、その具体的な種類は特に限定されない。第１パラメータとしては、離散要素法計算に必要となる各種パラメータが挙げられ、具体的には例えば粒子径、粒子数、ヤング率、計算のＴｉｍｅ ｓｔｅｐ、ポワソン比、壁面との摩擦係数、粒子間の摩擦係数、転がり摩擦係数、密度等から選択された１種類以上が挙げられる。
【0075】
第１パラメータは、データベース等に収録されているデータであってもよく、予め実験を行い求めた実験値であってもよい。また、第１パラメータは、実験結果からシミュレーション等によりフィッティングして算出した計算値であってもよい。
（Ｂ－２）第２パラメータ算出部
既述のように、本実施形態のシミュレーション装置３０では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。ただし、粗視化粒子は、粗視化前の粒子群を構成する個々の粒子とは質量等の各種パラメータが異なる。このため、粗視化粒子について、計算を行う際に必要となるパラメータの算出や、設定を行う必要がある。
（Ｂ－２）第２パラメータ算出部
既述のように、本実施形態のシミュレーション装置３０では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。ただし、粗視化粒子は、粗視化前の粒子群を構成する個々の粒子とは質量等の各種パラメータが異なる。このため、粗視化粒子について、計算を行う際に必要となるパラメータの算出や、設定を行う必要がある。
【0076】
第２パラメータ算出部４２２では、「（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて」の中で説明したように、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解であるＫｒを用いて、第２パラメータを算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、既述の垂直方向のＫｒの特性方程式である式（１３）を導出し、式（１３）から垂直方向のＫｒを算出できる。そして、既述の式（１３）に示した特性方程式の解である垂直方向のＫｒを用いて、第２パラメータを算出できる。また、既述の式（１９）に示した接線方向のＫｒの特性方程式を用いて、接線方向のＫｒを算出でき、係る接線方向のＫｒを用いて第２パラメータを算出することもできる。
【0077】
既述のように、Ｋｒは、粗視化粒子の挙動を支配するパラメータであり、Ｋｒを用いることで、粗視化粒子の挙動に関連する各種パラメータを算出できる。
【0078】
後述する粗視化粒子挙動解析部で用いる第２パラメータの種類は解析の内容に応じて選択できるため、特に限定されない。例えば第２パラメータは、粗視化粒子の熱伝導率を含むこともできる。この場合、第２パラメータ算出部は、既述の特性方程式の解Ｋｒを用いて、熱伝導率を算出できる。
（Ｂ－３）粗視化粒子挙動解析部
粗視化粒子挙動解析部４２３では、第１パラメータ取得部４２１により取得した第１パラメータ、および第２パラメータ算出部４２２で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。
（Ｂ－３）粗視化粒子挙動解析部
粗視化粒子挙動解析部４２３では、第１パラメータ取得部４２１により取得した第１パラメータ、および第２パラメータ算出部４２２で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。
【0079】
なお、ここでいう挙動とは、粗視化粒子の運動による位置の変化だけではなく、温度変化等の状態変化も含む。
（Ｂ－４）初期設定部
図示しない初期設定部は、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析部４２３で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得部４２１により取得する場合には、初期設定部は設けなくてもよい。
（Ｃ）出力部
出力部４３は、ディスプレイ等を有することができる。粗視化粒子挙動解析部４２３で得られたシミュレーション結果を出力部４３に出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部４３に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部４３に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。
（Ｂ－４）初期設定部
図示しない初期設定部は、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析部４２３で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得部４２１により取得する場合には、初期設定部は設けなくてもよい。
（Ｃ）出力部
出力部４３は、ディスプレイ等を有することができる。粗視化粒子挙動解析部４２３で得られたシミュレーション結果を出力部４３に出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部４３に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部４３に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。
【0080】
以上に説明した本実施形態のシミュレーション装置によれば、複数の粒子を含む粉体の挙動をシミュレーションすることができ、その用途等は特に限定されない。例えば、キルン等の回転体内での粉体の挙動のシミュレーションに好適に用いることができる。すなわち、本実施形態のシミュレーション装置によれば、回転体内での、粉体の挙動を解析することもできる。
【0081】
以上に説明した本実施形態のシミュレーション装置によれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。
【0082】
そして、粗視化粒子のパラメータを既述のパラメータＫｒを用いて計算するため、精度よく計算を行うことができる。
［シミュレーション方法］
次に、本実施形態のシミュレーション方法について説明する。本実施形態のシミュレーション方法は、例えば既述のシミュレーション装置を用いて実施できる。このため、既に説明した事項の一部は説明を省略する。
［シミュレーション方法］
次に、本実施形態のシミュレーション方法について説明する。本実施形態のシミュレーション方法は、例えば既述のシミュレーション装置を用いて実施できる。このため、既に説明した事項の一部は説明を省略する。
【0083】
本実施形態のシミュレーション方法は、複数の粒子を含む粉体の挙動を解析するためのシミュレーション方法に関する。本実施形態のシミュレーション方法は、図５に示したフロー図に従って実施することができ、以下の工程を有することができる。
【0084】
粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得工程（Ｓ１）。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程（Ｓ２）。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析工程（Ｓ３）。
そして、第２パラメータ算出工程（Ｓ２）は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出できる。
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出工程（Ｓ２）。
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析工程（Ｓ３）。
そして、第２パラメータ算出工程（Ｓ２）は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出できる。
【0085】
各工程について以下に説明する。
（１）第１パラメータ取得工程（Ｓ１）
第１パラメータ取得工程（Ｓ１）では、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。既述のシミュレーション装置を用いる場合、例えば第１パラメータ取得部４２１において、第１パラメータ取得工程を実施できる。
（１）第１パラメータ取得工程（Ｓ１）
第１パラメータ取得工程（Ｓ１）では、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得できる。既述のシミュレーション装置を用いる場合、例えば第１パラメータ取得部４２１において、第１パラメータ取得工程を実施できる。
【0086】
第１パラメータは解析の内容に応じて選択できるため、その具体的な種類は特に限定されない。第１パラメータとしては、離散要素法計算に必要となる各種パラメータが挙げられる。第１パラメータの具体的な例はシミュレーション装置で既に説明したため、ここでは説明を省略する。
【0087】
第１パラメータは、データベース等に収録されているデータであってもよく、予め実験を行い求めた実験値であってもよい。また、第１パラメータは、実験結果からシミュレーション等によりフィッティングして算出した計算値であってもよい。
（２）第２パラメータ算出工程（Ｓ２）
本実施形態のシミュレーション方法では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。
（２）第２パラメータ算出工程（Ｓ２）
本実施形態のシミュレーション方法では、計算量を抑制するため、粉体が有する複数個の粒子から構成される粒子群を、１個の粗視化粒子に粗視化し、粒子の数を少なくして計算を行うことができる。
【0088】
そこで、第２パラメータ算出工程（Ｓ２）では、「（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて」の中で説明したように、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解であるＫｒを用いて、第２パラメータを算出できる。具体的には例えば、粗視化前の粒子群１１全体の弾性エネルギーと、粗視化粒子全体の弾性エネルギーが等しくなると仮定して、既述の垂直方向のＫｒの特性方程式である式（１３）を導出し、式（１３）から垂直方向のＫｒを算出できる。そして、既述の式（１３）に示した特性方程式の解である垂直方向のＫｒを用いて、第２パラメータを算出できる。また、既述の式（１９）に示した接線方向のＫｒの特性方程式を用いて、接線方向のＫｒを算出でき、係る接線方向のＫｒを用いて第２パラメータを算出することもできる。
【0089】
既述のように、Ｋｒは、粗視化粒子の挙動を支配するパラメータであり、Ｋｒを用いることで、粗視化粒子の挙動に関連する各種パラメータを算出できる。
【0090】
既述のシミュレーション装置を用いる場合、例えば第２パラメータ算出部４２２において、第２パラメータ算出工程を実施できる。
【0091】
後述する粗視化粒子挙動解析工程で用いる第２パラメータの種類は解析の内容に応じて選択できるため、特に限定されない。例えば第２パラメータは、粗視化粒子の熱伝導率を含むこともできる。この場合、第２パラメータ算出工程は、既述の特性方程式の解Ｋｒを用いて、熱伝導率を算出できる。
（３）粗視化粒子挙動解析工程（Ｓ３）
粗視化粒子挙動解析工程（Ｓ３）では、第１パラメータ取得工程（Ｓ１）で得た第１パラメータ、および第２パラメータ算出工程（Ｓ２）で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。
（３）粗視化粒子挙動解析工程（Ｓ３）
粗視化粒子挙動解析工程（Ｓ３）では、第１パラメータ取得工程（Ｓ１）で得た第１パラメータ、および第２パラメータ算出工程（Ｓ２）で算出した第２パラメータを用いて、粗視化粒子の挙動を解析できる。具体的には離散要素法を用いて計算し、粗視化粒子の挙動を解析できる。粗視化粒子の挙動を解析することで、粉体の挙動を解析することができる。
【0092】
なお、ここでいう挙動とは、粗視化粒子の運動による位置の変化だけではなく、温度変化等の状態変化も含む。
（４）初期設定工程
本実施形態のシミュレーション方法は、例えばさらに初期設定工程を有することもできる。初期設定工程では、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析工程で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得工程により取得する場合には、初期設定工程は実施しなくてもよい。
（５）出力工程
本実施形態のシミュレーション方法は、例えばさらに出力工程を有することができる。出力工程では、例えば粗視化粒子挙動解析工程（Ｓ３）で得られたシミュレーション結果を、出力部へ出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。
（４）初期設定工程
本実施形態のシミュレーション方法は、例えばさらに初期設定工程を有することもできる。初期設定工程では、解析対象となる粉体を構成する粒子の位置を初期化するとともに、解析の条件、例えば必要に応じて粉体を配置する領域の温度等を設定できる。なお、例えば粗視化粒子挙動解析工程で粗視化粒子の挙動を解析する際に用いるプログラム等に予め初期条件が設定されている場合や、第１パラメータ取得工程により取得する場合には、初期設定工程は実施しなくてもよい。
（５）出力工程
本実施形態のシミュレーション方法は、例えばさらに出力工程を有することができる。出力工程では、例えば粗視化粒子挙動解析工程（Ｓ３）で得られたシミュレーション結果を、出力部へ出力できる。出力するシミュレーション結果の内容は特に限定されないが、例えば出力部に粗視化された粒子の位置を時系列で画像として出力し、表示することができる。また、例えば出力部に、粉体の温度分布の時系列変化を画像として出力し、表示することもできる。
【0093】
以上に説明した本実施形態のシミュレーション方法によれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。
【0094】
そして、粗視化粒子のパラメータを既述のパラメータＫｒを用いて計算するため、精度よく計算を行うことができる。
［プログラム］
次に、本実施形態のプログラムについて説明する。
［プログラム］
次に、本実施形態のプログラムについて説明する。
【0095】
本実施形態のプログラムは、複数の粒子を含む粉体の挙動を解析するためのプログラムに関し、コンピュータを以下の各部として機能させることができる。
【0096】
粉体に関連するパラメータを含む第１パラメータを取得する第１パラメータ取得部。
【0097】
複数個の粒子から構成される粒子群を粗視化し、１個の粗視化粒子とした場合の、粗視化粒子についてのパラメータである第２パラメータを算出する第２パラメータ算出部。
【0098】
第１パラメータ、および第２パラメータに基づいて、粗視化粒子の挙動を解析する粗視化粒子挙動解析部。
【0099】
第２パラメータ算出部では、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出させることができる。
【0100】
また、第２パラメータは、粗視化粒子の熱伝導率を含むことができ、この場合、第２パラメータ算出部では、既述の特性方程式の解を用いて、熱伝導率を算出できる。
【0101】
本実施形態のプログラムは、例えば既述のシミュレーション装置のＲＡＭやＲＯＭ等の主記憶装置または補助記憶装置の各種記憶媒体に記憶させておくことができる。そして、係るプログラムを読み込ませ、ＣＰＵにより実行することにより、ＲＡＭ等におけるデータの読み出しおよび書き込みを行うと共に、入出力インタフェースおよび表示装置を動作させて実行できる。このため、シミュレーション装置で既に説明した事項については説明を省略する。
【0102】
上述した本実施形態のプログラムは、インターネットなどのネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることで提供してもよい。また、本実施形態のプログラムをインターネットなどのネットワークを介して提供、配布するように構成してもよい。
【0103】
本実施形態のプログラムは、ＣＤ－ＲＯＭ等の光ディスクや、半導体メモリ等の記録媒体に格納した状態で流通等させてもよい。
【0104】
以上に説明した本実施形態のプログラムによれば、複数個の粒子からなる粒子群を１個の粗視化粒子とすることで、計算量を抑制できる。このため、工場で用いるプラント等のような大きな規模での粉体の挙動についても計算量を抑制し、効率的に計算を行うことができる。
【0105】
そして、粗視化粒子のパラメータを既述のパラメータＫｒを用いて計算するため、精度よく計算を行うことができる。
２．第２の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
第２の実施形態では、接線方向の運動方程式に関して、オーバーラップ量を算出する際に、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定する点がここまで説明した第１の実施形態の場合と異なる。また、回転方向の運動方程式に関して、上述のようにして求めたオーバーラップ量を用いてトルクを算出できる。これにより計算量を抑制したまま、粗視化後の粒子群の挙動についてより精度よく解析可能になる。
（１－１）オーバーラップ量について
通常、接線方向のオーバーラップ量は接触開始から接触終了までの速度の接線方向の成分ｖｔ（以下の式中ではｖｔの上に矢印で示されたもの）の時間積分を用いて以下のように求められる。ここでｔは時間、ベクトルｒ（以下の式中ではｒの上に矢印で示されたもの）は、粒子群１１を構成する粒子の中心から接触点までのベクトル、ベクトルω（以下の式中ではωの上に矢印で示されたもの）は粒子群１１の回転ベクトルを示す。
２．第２の実施形態
［シミュレーション装置］
（１）粒子の粗視化、および粗視化粒子の粒子挙動の計算に用いるパラメータについて
第２の実施形態では、接線方向の運動方程式に関して、オーバーラップ量を算出する際に、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定する点がここまで説明した第１の実施形態の場合と異なる。また、回転方向の運動方程式に関して、上述のようにして求めたオーバーラップ量を用いてトルクを算出できる。これにより計算量を抑制したまま、粗視化後の粒子群の挙動についてより精度よく解析可能になる。
（１－１）オーバーラップ量について
通常、接線方向のオーバーラップ量は接触開始から接触終了までの速度の接線方向の成分ｖｔ（以下の式中ではｖｔの上に矢印で示されたもの）の時間積分を用いて以下のように求められる。ここでｔは時間、ベクトルｒ（以下の式中ではｒの上に矢印で示されたもの）は、粒子群１１を構成する粒子の中心から接触点までのベクトル、ベクトルω（以下の式中ではωの上に矢印で示されたもの）は粒子群１１の回転ベクトルを示す。
【0106】
なお、以下の式中のハットｔ（以下の式中ではｔの上にハットを付して示されたもの）は、接線方向オーバーラップの単位ベクトルを示す。
【0107】
【数29】
【0108】
【数30】
【0109】
【数31】
【0110】
【数32】
【0111】
【数33】
【0112】
【数34】
【0113】
【数35】
【0114】
上記式（３３）、式（３４）を連立して公転運動成分を消去すると、以下の式（３５）になる。
【0115】
【数36】
【0116】
【数37】
【0117】
【数38】
【0118】
【数39】
【0119】
【数40】
【0120】
【数41】
【0121】
【数42】
【0122】
【数43】
【0123】
【数44】
【0124】
このようにして、粗視化前の粒子群を構成する粒子のオーバーラップ量も求めることができる。
（１－２）回転方向の運動方程式について
ここまで説明した接線方向のオーバーラップ量から接線方向の力を算出し、トルクを計算できる。各粒子間や、粒子と壁間での転がり摩擦は粒子にかかる垂直抗力と接触半径の積に比例するトルクが発生するとした。それぞれの転がり摩擦は各粒子で発生するので、粗視化粒子全体での転がり摩擦抵抗は以下の式（４４）のように表せる。
（１－２）回転方向の運動方程式について
ここまで説明した接線方向のオーバーラップ量から接線方向の力を算出し、トルクを計算できる。各粒子間や、粒子と壁間での転がり摩擦は粒子にかかる垂直抗力と接触半径の積に比例するトルクが発生するとした。それぞれの転がり摩擦は各粒子で発生するので、粗視化粒子全体での転がり摩擦抵抗は以下の式（４４）のように表せる。
【0125】
式（４４）中、ベクトルＴｔｏｔ＿ｆｒｉｃが粗視化粒子全体の転がり摩擦抵抗を、ベクトルＴｗ，ｆｒｉｃが粒子群１１の壁面１２あるいは外部粒子との間の転がり摩擦抵抗を、ベクトルＴｐ，ｆｒｉｃが粒子群１１の内部粒子間の転がり摩擦抵抗を表している。各ベクトルは、以下の式中では文字の上に矢印を付して表記している。なお、粒子群１１では粒子－粒子間に作用反作用により逆向きの接触力が働くため、式（４４）に示した様に、ベクトルＴｐ，ｆｒｉｃは２倍となる。
【0126】
そして、ベクトルＴｗ，ｆｒｉｃと、ベクトルＴｐ，ｆｒｉｃとは、式（４５）、式（４６）により表される。なお、式中のハットωは、回転方向の単位ベクトルを意味する。
【0127】
【数45】
【0128】
【数46】
【0129】
【数47】
【0130】
そして、回転方向の運動方程式は以下のよう表記できる。
【0131】
【数48】
本実施形態のシミュレーション装置においても、第２パラメータ算出部は、粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いた特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。
【0132】
以上の点以外は第１の実施形態のシミュレーション装置の場合と同様に構成できるため、ここでは説明を省略する。
［シミュレーション方法］
本実施形態のシミュレーション方法においても、第２パラメータ算出工程では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。
［シミュレーション方法］
本実施形態のシミュレーション方法においても、第２パラメータ算出工程では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。
【0133】
以上の点以外は第１の実施形態のシミュレーション方法の場合と同様に構成できるため、ここでは説明を省略する。
［プログラム］
本実施形態のプログラムにおいても、第２パラメータ算出部では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。
［プログラム］
本実施形態のプログラムにおいても、第２パラメータ算出部では、粗視化前の粒子群の弾性エネルギーと、粗視化粒子の弾性エネルギーとの関係を用いて導出した特性方程式の解を用いて、第２パラメータを算出できる。ただし、接線方向の運動方程式を用いる場合において、オーバーラップ量に関して、回転については角運動量と回転エネルギーが粗視化前後で一致すると仮定し算出できる。また、回転方向の運動方程式を用いる場合において、トルクの算出を行う際に、上記接線方向の運動方程式を用いる場合に算出したオーバーラップ量を用いることができる。
【0134】
以上の点以外は第１の実施形態のプログラムの場合と同様に構成できるため、ここでは説明を省略する。
【実施例】
【0135】
以下に具体的な実施例を挙げて説明するが、本発明はこれらの実施例に限定されるものではない。
［実験例１］
［比較例１－１］
直方体の容器内に充填した粉体層の温度変化について、表１に示したパラメータを用いて、離散要素法計算により底部から加熱を行った場合の解析を行った。解析により求めた容器内の粉体層の平均温度の変化を図６に示す。なお、反発係数を０．１、摩擦係数を０．７、転がり摩擦係数を０．００１とした。
［実験例１］
［比較例１－１］
直方体の容器内に充填した粉体層の温度変化について、表１に示したパラメータを用いて、離散要素法計算により底部から加熱を行った場合の解析を行った。解析により求めた容器内の粉体層の平均温度の変化を図６に示す。なお、反発係数を０．１、摩擦係数を０．７、転がり摩擦係数を０．００１とした。
【0136】
【表1】
粗視化倍率αを２とし、８個の粒子から構成される粒子群を１個の粗視化粒子とした点以外は比較例１－１と同様の条件で、直方体の容器内に充填した粉体層の温度変化について解析を行った。なお、粗視化粒子としたため、表１に示すように、粒子の粒径を変更しているが、粒径以外は比較例１－１と同じパラメータを用いて、解析を行った。解析により求めた容器内の粉体層の平均温度の変化を図６に示す。
［実施例１－１］
粗視化倍率αを２とし、既述のシミュレーション装置を用いて、直方体の容器内に充填した粉体層の温度変化について解析を行った。
【0137】
具体的には、第１パラメータ取得部４２１により、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得した（第１パラメータ取得工程：Ｓ１）。
【0138】
この際には、比較例１－２と同じパラメータを取得した。
【0139】
次いで、第２パラメータ算出部４２２により、粗視化粒子についてのパラメータである第２パラメータを算出した（第２パラメータ算出工程：Ｓ２）。
【0140】
この際、既述の式（１３）、式（１９）に示した特性方程式の解であるパラメータＫｒを用いた式（２６）、式（２７）により、粗視化粒子に関する熱伝導率を算出した。また、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出した。既述のように、計算に適用するモデルに応じて、これらの係数は調整可能であり、上述のように特性方程式を用いて算出、換算できる。算出したパラメータＫｒを表２に示す。
【0141】
【表2】
【0142】
図６に示した結果によると、実施例１－１の結果は、粗視化を行っていない、比較例１－１の結果とほぼ一致しており、粗視化粒子のパラメータを適切に選択、設定し、解析を行えていることを確認できた。
【0143】
実施例１－１では粗視化を行っているため、粒子の数が比較例１－１と比較して１／８倍になっており、比較例１－１と比較して計算量を抑制できている。
［実験例２］
［比較例２－１］
回転体であるキルン内の粉体の挙動について、表３に示したパラメータを用いて、離散要素法計算により解析を行った。なお、反発係数を０．７５、摩擦係数を０．３、転がり摩擦係数を０．５とした。
解析により求めたキルン内の粉体の動きを図８に、解析により求めたキルン内の粉体の温度分布を図１１に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。
［実験例２］
［比較例２－１］
回転体であるキルン内の粉体の挙動について、表３に示したパラメータを用いて、離散要素法計算により解析を行った。なお、反発係数を０．７５、摩擦係数を０．３、転がり摩擦係数を０．５とした。
解析により求めたキルン内の粉体の動きを図８に、解析により求めたキルン内の粉体の温度分布を図１１に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。
【0144】
図８（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図８（Ｂ）、図８（Ｃ）、図８（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。
【0145】
図１１（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１１（Ｂ）、図１１（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。
【0146】
【表3】
粗視化倍率αを４とし、６４個の粒子から構成される粒子群を１個の粗視化粒子とした点以外は比較例２－１と同様の条件で、回転体であるキルン内の粉体の挙動について解析を行った。なお、表３に示すように、粒径以外は比較例２－１と同じパラメータを用いて、解析を行った。解析により求めたキルン内の粉体の動きを図９に、解析により求めたキルン内の粉体の温度分布を図１２に、解析により求めたキルン内の粉体の平均温度を図１３にそれぞれ示す。
【0147】
図９（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図９（Ｂ）、図９（Ｃ）、図９（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。
【0148】
図１２（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１２（Ｂ）、図１２（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。
［実施例２－１］
粗視化倍率αを４とし、既述のシミュレーション装置を用いて、回転体であるキルン内の粉体の挙動について解析を行った。
［実施例２－１］
粗視化倍率αを４とし、既述のシミュレーション装置を用いて、回転体であるキルン内の粉体の挙動について解析を行った。
【0149】
具体的には、第１パラメータ取得部４２１により、解析の対象となる粉体に関連するパラメータを含む第１パラメータを取得した（第１パラメータ取得工程：Ｓ１）。
【0150】
この際には、比較例２－２と同じパラメータを取得した。
【0151】
次いで、第２パラメータ算出部４２２により、粗視化粒子についてのパラメータである第２パラメータを算出した（第２パラメータ算出工程：Ｓ２）。
【0152】
この際、既述の式（１３）、式（１９）に示した特性方程式の解であるパラメータＫｒを用いた式（２６）、式（２７）により、粗視化粒子に関する熱伝導率を算出した。また、特性方程式を用いて反発係数、摩擦係数、転がり摩擦係数を算出した。算出したパラメータＫｒを表４に示す。
【0153】
【表4】
【0154】
図７（Ａ）はキルン７０の回転を開始する前の状態を示しており、グループ分けした第１粉体群７１と、第２粉体群７２とが半分ずつ入っている。図７（Ｂ）、図７（Ｃ）、図７（Ｄ）は、それぞれキルンの回転開始後２秒経過時、４秒経過時、６秒経過時の状態を示しており、第１粉体群７１と、第２粉体群７２とが混合されている状態を示している。
【0155】
図１０（Ａ）はキルンの回転を開始する前の状態を示しており、温度が第１温度域１０１内にあり、均一であることを確認できる。図１０（Ｂ）、図１０（Ｃ）は、それぞれキルンの回転開始後３秒経過時、６秒経過時の状態を示しており、キルン７０の外壁側から加熱を行っているため、粉体の中心部側から順に第１温度域１０１、第２温度域１０２、第３温度域１０３が分布していることを確認できる。なお、第１温度域１０１、第２温度域１０２、第３温度域１０３の順に温度が高くなる。
【0156】
図７と図８の比較や、図１０と図１１の比較から明らかなように、実施例２－１の結果は、粗視化を行っていない、比較例２－１の結果とほぼ一致しており、粗視化粒子のパラメータを適切に選択、設定し、解析を行えていることを確認できた。
【0157】
図１３においても実施例２－１と比較例２－１とは変曲点を有するなど同様の傾向を示すことが確認できた。
【0158】
そして、実施例２－１では粗視化を行っているため、粒子の数が比較例２－１と比較して１／６４倍になっており、比較例２－１と比較して計算量を抑制できている。
【符号の説明】
【0159】
１１ 粒子群
１１Ａ、１１Ｂ 粒子
２１ 粗視化粒子
３０ シミュレーション装置
４２１ 第１パラメータ取得部
４２２ 第２パラメータ算出部
４２３ 粗視化粒子挙動解析部
Ｓ１ 第１パラメータ取得工程
Ｓ２ 第２パラメータ算出工程
Ｓ３ 粗視化粒子挙動解析工程
１１Ａ、１１Ｂ 粒子
２１ 粗視化粒子
３０ シミュレーション装置
４２１ 第１パラメータ取得部
４２２ 第２パラメータ算出部
４２３ 粗視化粒子挙動解析部
Ｓ１ 第１パラメータ取得工程
Ｓ２ 第２パラメータ算出工程
Ｓ３ 粗視化粒子挙動解析工程
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】