特許 7122209 シミュレーション装置、シミュレーション方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2022-08-10

(45)【発行日】2022-08-19

(54)【発明の名称】シミュレーション装置、シミュレーション方法、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20220812BHJP

   G16C 10/00 20190101ALI20220812BHJP

【ＦＩ】

G16Z99/00

G16C10/00

【請求項の数】 8

(21)【出願番号】P 2018186520

(22)【出願日】2018-10-01

(65)【公開番号】P2020057135

(43)【公開日】2020-04-09

【審査請求日】2021-06-16

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】石原 定典

【審査官】宮地 匡人

(56)【参考文献】

【文献】特開２０１０－１０８１８３（ＪＰ，Ａ）

【文献】特開２０１２－１１８５７９（ＪＰ，Ａ）

【文献】廣岡 信行，振動状態におけるトナー微細粒子群の流動解析，計算工学講演会論文集［ＣＤ－ＲＯＭ］，2017年05月，Vol.22，pp.1-5

【文献】山井 三亀夫，スケール則に基づく DEM シミュレーション，粉体工学会誌，2018年02月10日，Vol.55 No.2，pp.95-103

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｚ ９９／００

Ｇ１６Ｃ １０／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

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

【請求項2】

前記処理装置は、
粗視化後の粉粒体の粒子の粒径に応じて、色相、彩度、明度の少なくとも１つを異ならせて、シミュレーションによって求められた粒子の挙動を前記出力装置に表示する請求項１に記載のシミュレーション装置。

【請求項3】

前記粒径分布を規定するパラメータは、粉粒体の粒径の分布を表すロジン－ラムラ分布の長さの次元を持つパラメータ、及び指数パラメータである請求項１または２に記載のシミュレーション装置。

【請求項4】

粗視化係数の入力された値をＫとしたとき、
前記処理装置は、
長さの次元を持つパラメータとして、入力された値のＫ倍の値を用い、指数パラメータとして、入力された値と同一の値を用いたロジン－ラムラ分布に基づいて粗視化された複数の粒子の粒径を設定し、粗視化された粉粒体の挙動を、離散要素法を用いてシミュレーションする請求項３に記載のシミュレーション装置。

【請求項5】

シミュレーション対象の粉粒体の物性値、物理量、粒径の分布、及び粒子を粗視化する基準となる粗視化係数の値を決定し、
決定された物性値、物理量、粒径の分布、及び粗視化係数の値に基づいて、粗視化後の物性値、物理量、及び粒径分布を算出し、
粗視化後の物性値、物理量、及び粒径分布に基づいて、粗視化後の粉粒体の粒子の挙動を、離散要素法を用いてシミュレーションするシミュレーション方法。

【請求項6】

シミュレーション対象の粉粒体の粒径の分布を決定する際に、ロジン－ラムラ分布の長さの次元を持つパラメータの値と指数パラメータの値とを決定し、
粗視化後の粒径分布を算出する際に、決定された値に基づいて、粗視化後の粉粒体の粒径の分布を規定するロジン－ラムラ分布の長さの次元を持つパラメータの値と指数パラメータの値とを決定する請求項５に記載のシミュレーション方法。

【請求項7】

シミュレーション対象の粉粒体を粗視化するときの粗視化係数の値をＫとしたとき、
粗視化後の粉粒体の粒径の分布を規定するロジン－ラムラ分布の長さの次元を持つパラメータ及び指数パラメータとして、シミュレーション対象の粉粒体の粒径の分布を規定するロジン－ラムラ分布の長さの次元を持つパラメータの値のＫ倍の値、及び指数パラメータと同一の値を用いる請求項６に記載のシミュレーション方法。

【請求項8】

シミュレーション対象の粉粒体の物性値、物理量、粒径分布を規定するパラメータの値、及び粉粒体を粗視化する基準となる粗視化係数の値を取得する機能と、
決定された物性値、物理量、粒径の分布、及び粗視化係数の値に基づいて、粗視化後の物性値、物理量、及び粒径分布を算出し、
粗視化後の物性値、物理量、及び粒径分布に基づいて、粗視化後の粉粒体の粒子の挙動を、離散要素法を用いてシミュレーションする機能と、
シミュレーションによって求められた粒子の挙動と、取得された粗視化係数の値とを関連付けて出力する機能と
をコンピュータに実現させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

複数の粒子からなる粉粒体の挙動を解析する離散要素法（ＤＥＭ）と流体の流れ場を解析する数値流体力学（ＣＦＤ）とを連成させることで、固体粒子を流体中に浮遊させた状態の流動層の挙動を解析する手法が公知である（非特許文献１）。非特許文献１には、粒子数が増えたときの計算時間の増大を抑制するシミュレーション方法が提案されている。具体的には、粒子を拡大して粒子数を減らす処理（粗視化）を行い、粗視化の前後で支配方程式が同じになるように物性値や物理量を変換し、粗視化後の流動層についてシミュレーションを実行する。

【先行技術文献】

【非特許文献】

【0003】

【文献】鷲野 公彰、許 志宏、川口 寿裕、辻 裕、「流動層のＤＥＭ計算における相似則モデル」、粉体工学会誌 第４４巻第３号、２００７年、第１９８頁～第２０５頁

【発明の概要】

【発明が解決しようとする課題】

【0004】

従来のシミュレーション方法では、複数の粒子の粒径が均一である粉粒体が解析対象となる。例えばガラスビーズ等の工業製品のように粒径がほぼ均一とみなせる粉粒体を解析対象とすることが可能である。珪砂等のように粒径が均一ではない粉粒体を粗視化する手法は知られていない。粒径が均一ではない粉粒体を、粒径が均一であると仮定して粗視化を行い、粗視化後の粉粒体について解析を行っても、妥当な結果を得ることはできない。

【0005】

本発明の目的は、粒径が均一ではない粉粒体の挙動を解析することが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。

【課題を解決するための手段】

【0006】

本発明の一観点によると、
シミュレーション条件の入力を行う入力装置と、
シミュレーション結果を出力する出力装置と、
前記入力装置から入力されたシミュレーション条件に基づいて、大きさが異なる複数の粒子を含む粉粒体の挙動を解析する処理装置と
を有し、
前記処理装置は、
前記入力装置から入力されたシミュレーション対象の粉粒体の物性値、物理量、及び粒径分布を規定するパラメータの値、及び粒子を粗視化する基準となる粗視化係数の値に基づいて、粗視化後の物性値、物理量、及び粒径分布を算出し、
粗視化後の物性値、物理量、及び粒径分布に基づいて、離散要素法を用いて、粗視化された粉粒体の挙動をシミュレーションにより求め、
シミュレーションによって求められた粒子の挙動と、入力された粗視化係数の値とを関連付けて前記出力装置に出力するシミュレーション装置が提供される。

【0007】

本発明の他の観点によると、
シミュレーション対象の粉粒体の物性値、物理量、粒径の分布、及び粒子を粗視化する基準となる粗視化係数の値を決定し、
決定された物性値、物理量、粒径の分布、及び粗視化係数の値に基づいて、粗視化後の物性値、物理量、及び粒径分布を算出し、
粗視化後の物性値、物理量、及び粒径分布に基づいて、粗視化後の粉粒体の粒子の挙動を、離散要素法を用いてシミュレーションするシミュレーション方法が提供される。

【0008】

本発明のさらに他の観点によると、
シミュレーション対象の粉粒体の物性値、物理量、粒径分布を規定するパラメータの値、及び粉粒体を粗視化する基準となる粗視化係数の値を取得する機能と、
決定された物性値、物理量、粒径の分布、及び粗視化係数の値に基づいて、粗視化後の物性値、物理量、及び粒径分布を算出し、
粗視化後の物性値、物理量、及び粒径分布に基づいて、粗視化後の粉粒体の粒子の挙動を、離散要素法を用いてシミュレーションする機能と、
シミュレーションによって求められた粒子の挙動と、取得された粗視化係数の値とを関連付けて出力する機能と
をコンピュータに実現させるためのプログラムが提供される。

【発明の効果】

【0009】

粒径が均一ではない粉粒体の挙動を解析することが可能になる。さらに、シミュレーション結果の確認時に、シミュレーションを行う際に用いた粗視化係数を容易に把握することができる。

【図面の簡単な説明】

【0010】

【図1】図１Ａは、参考例によるシミュレーション対象の流動層の一例を示す模式図であり、図１Ｂは、シミュレーション対象の粗視化後の流動層の一例を示す模式図である。

【図2】図２は、粒子及びガスの物性値、粒子及びガスに関して定義される種々の物理量について、本明細書で用いる記号及び粗視化の係数の一覧を示す図表である。

【図3】図３は、本参考例によるシミュレーション装置のブロック図である。

【図4】図４は、本参考例によるシミュレーション方法のフローチャートである。

【図5】図５は、参考例による手法を用いて実際に行ったシミュレーションのシミュレーション領域を示す斜視図である。

【図6】図６は、粗視化した流動層の粗視化粒子の位置及び温度のシミュレーション結果を、時系列で示す図である。

【図7】図７Ａ及び図７Ｂは、参考例による手法のシミュレーション結果から求めた粒子の平均温度の時間変化を示すグラフである。

【図8】図８は、粗視化前及び粗視化後の粉粒体の粒径の質量分率の一例を示すグラフである。

【図9】図９は、実施例によるシミュレーション装置を用いたシミュレーション方法のフローチャートである。

【図10】図１０Ａは、初期充填領域に充填されている状態の粉粒体の側面図であり、図１０Ｂ及び図１０Ｃは、それぞれ粗視化前の粉粒体、及び粗視化係数Ｋ＝３０で粗視化した粉粒体の、十分時間が経過した後の分布を示す図である。

【図11】図１１は、実際に粒径が均一ではない珪砂を自然落下させて形成された珪砂の山を撮影した図である。

【発明を実施するための形態】

【0011】

［参考例］
本願の発明の実施例について説明する前に、図１Ａ～図７Ｂを参照して、実施例に関連する参考例について説明する。

【0012】

図１Ａ～図７Ｂを参照して、参考例によるシミュレーション方法及び装置について説明する。
図１Ａは、シミュレーション対象の流動層の一例を示す模式図である。シミュレーション対象の領域１０内に複数の粒子１１を配置し、下方から上方に向かって領域１０内にガス１２を導入することにより形成される流動層の挙動をシミュレーションする。粒子１１の直径をＤ_ｐ１で表す。本参考例では、粒子１１の各々を拡大して、その個数を減らす（以下、粗視化という。）ことにより、計算負荷を軽減させる。

【0013】

図１Ｂは、シミュレーション対象の粗視化後の流動層の一例を示す模式図である。粒子１１を拡大して仮想的な粒子２１を得る。仮想的な粒子２１は、シミュレーション対象の領域２０内に配置される。粗視化後の領域２０の寸法は、粗視化前の領域１０の領域の寸法と同一である。仮想的な粒子２１の直径をＤ_ｐ２で表す。拡大率（粗視化係数）Ｋを、粗視化前の粒子１１の直径に対する粗視化後の仮想的な粒子２１の直径の比と定義する。拡大率Ｋは以下の式で定義される。

【数1】

【0014】

粗視化後の粒子２１が配置された領域２０内に下方から上方に向かってガス２２を導入することにより形成される粗視化後の流動層について、数値流体力学（ＣＦＤ）と離散要素法（ＤＥＭ）とを連成させた解析を行う。粗視化後の仮想的な流動層と、粗視化前の実際の流動層とが相似則を満たすように、粗視化に際し、粒子１１及びガス１２の物性値及び種々の物理量を変換する。

【0015】

次に、図２を参照して、粒子１１及びガス１２の物性値及び種々の物理量の変換則について説明する。

【0016】

図２は、粒子及びガスの物性値、粒子及びガスに関して定義される種々の物理量について、本明細書で用いる記号及び粗視化の係数の一覧を示す図表である。粗視化前の実際の物性値及び物理量に粗視化の係数を乗ずることにより、粗視化後の流動層に関する物性値及び物理量が得られる。本明細書において、例えば式（１）に示したように、粗視化前の物性値及び物理量を表す記号には、下付きの添え字「１」を付し、粗視化後の物性値及び物理量を表す記号には、下付きの添え字「２」を付す。

【0017】

流動層の流れに関する無次元量として、粒子レイノルズ数Ｒｅ_ｐ、アルキメデス数Ａｒ_ｐ、及びフルード数Ｆｒが挙げられる。これらの無次元量は以下の式で定義される。

【数2】

ここで、ｇは重力加速度である。太字のＶ及びＵは、ベクトルであることを意味する。ボイド率εは、充填された粒子の全質量をＭ、粒子が充填された領域の見かけの体積をＶ_Ａとして、以下の式で定義される。

【数3】

【0018】

粗視化の前後で、流動層の流れに関する無次元量である粒子レイノルズ数Ｒｅ_ｐ、アルキメデス数Ａｒ_ｐ、及びフルード数Ｆｒが変化しないという条件を設定する。さらに、ボイド率εが変化しないという条件、及びガス粘性係数μが変化しないという条件の下で、粗視化前後の物性値及び物理量の変換則を求めると、以下の変換則が得られる。

【数4】

【0019】

ガス密度ρ_ｆ２の変換則から、ガス圧力ｐについて、以下の変換則が得られる。

【数5】

【0020】

粗視化の前後で粒子が充填された領域の見かけの体積Ｖ_Ａが変化せず、粒子の個数が粗視化によって１／Ｋ^３に減少すると仮定すると、以下の変換則が得られる。

【数6】

【0021】

粒子質量流量ｍ_ｐドットは、流路面積をＡとして以下の式で定義される。

【数7】

この式から、以下の変換則が導出される。

【数8】

【0022】

さらに、熱輸送に関する無次元量についても、粗視化の前後で変化しないという条件を付す。熱輸送に関する無次元量として、プラントル数Ｐｒ、粒子ヌセルト数Ｎｕ_ｐ、ビオ数Ｂｉが挙げられる。プラントル数Ｐｒ、粒子ヌセルト数Ｎｕ_ｐ、ビオ数Ｂｉは以下の式で定義される。

【数9】

ここで、Ｌ_ｐは粒子の特徴長さであり、Ｌ_ｐ＝Ｄ_ｐ／６で定義することができる。

【0023】

物性値の温度依存性を簡単化するために、粗視化の前後で粒子温度Ｔ_ｐ及びガス温度Ｔが変化しないと仮定する。さらに、粒子熱伝達係数ｈも、粗視化の前後で変化しないと仮定する。この仮定の下で、以下の変換則が得られる。

【数10】

【0024】

上述の仮定のみでは、粒子比熱ｃの変換則が定まらない。本参考例では、粒子比熱ｃの変換則を決定するために、粗視化の前後で粒子全体の顕熱Ｑ_{ｐ，ａｌｌ}が変化しないという仮定を導入する。粒子全体の顕熱Ｑ_{ｐ，ａｌｌ}は、粒子の個数をＮ_ｐ、流動層に導入するガス温度Ｔに対する粒子の初期温度の差をΔＴ_ｐとして、以下の式で定義される。

【数11】

【0025】

粒子の個数Ｎ_ｐは、粗視化によって約１／Ｋ^３に減少するため、粒子全体の顕熱Ｑ_{ｐ，ａｌｌ}が粗視化の前後で不変と仮定すると、以下の変換則が得られる。

【数12】

【0026】

粒子表面の伝熱量Ｑドットは、以下の式で定義される。

【数13】

この定義から、伝熱量Ｑドットについて以下の変換則が得られる。

【数14】

粒子表面の熱流束ｑドットについては、以下の変換則が得られる。

【数15】

【0027】

図３は、本参考例によるシミュレーション装置のブロック図である。本参考例によるシミュレーション装置は、処理装置３０、入力装置３８、及び出力装置３９を含む。処理装置３０は、シミュレーション条件取得部３１、拡大率取得部３２、演算部３３、及び出力制御部３４を含む。

【0028】

図３に示す各ブロックは、ハードウェア的には、コンピュータの中央処理ユニット（ＣＰＵ）をはじめとする素子や機械装置で実現することができ、ソフトウェア的にはコンピュータプログラム等によって実現することができる。図３では、ハードウェア及びソフトウェアの連携によって実現される機能ブロックが示されている。従って、これらの機能ブロックは、ハードウェア及びソフトウェアの組み合わせによって、種々の態様で実現することが可能である。

【0029】

処理装置３０は入力装置３８及び出力装置３９と接続される。入力装置３８は、処理装置３０で実行される処理に関係するユーザからのコマンド及びデータの入力を受ける。入力装置３８として、例えばユーザが操作を行うことにより入力を行うキーボードやマウス、インターネット等のネットワークを介して入力を行う通信装置、ＣＤ、ＤＶＤ等の記録媒体から入力を行う読取装置等を用いることができる。

【0030】

シミュレーション条件取得部３１は、入力装置３８を介してシミュレーション条件を取得する。シミュレーション条件には、シミュレーションに必要な種々の情報が含まれる。例えば、シミュレーション対象の粒子及びガスの物性値、粒子及びガスに関する物理量の初期条件、境界条件等が含まれる。拡大率取得部３２は、入力装置３８を介して拡大率Ｋ（図２）を取得する。

【0031】

演算部３３は、シミュレーション条件及び拡大率Ｋに基づいて、粗視化前の物性値及び物理量に粗視化の係数（図２）を乗じることにより、粒子及びガスの粗視化後の物性値及び物理量の初期条件を算出する。粗視化後の物性値及び物理量の初期条件に基づき、ＣＦＤとＤＥＭとを連成させた流動層のシミュレーションを行う。

【0032】

出力制御部３４は、シミュレーション結果を出力装置３９に出力する。例えば、粒子の位置及び温度の変動、ガスの温度分布の変動を、出力装置３９の表示画面に図形で表示する。

【0033】

図４は、本参考例によるシミュレーション方法のフローチャートである。まず、シミュレーション条件取得部３１（図３）がシミュレーション条件を取得し（ステップＳ１）、拡大率取得部３２（図３）が拡大率Ｋ（図２）を取得する（ステップＳ２）。

【0034】

その後、演算部３３（図３）は、シミュレーション条件として入力された物性値及び物理量の初期値を、粗視化後の値に変換する（ステップＳ３）。さらに、変換後の物性値及び物理量に基づいてシミュレーションを実行する（ステップＳ４）。シミュレーションが終了すると、出力制御部３４（図３）がシミュレーション結果を出力する（ステップＳ５）。

【0035】

次に、図５～図７Ｂを参照して、本参考例によるシミュレーション方法を用いて実際にシミュレーションを行った結果について説明する。このシミュレーションの対象は、非特許文献１に記載されているものと同一である。

【0036】

図５は、シミュレーション領域４０を示す斜視図である。シミュレーション領域４０は、幅８ｃｍ、厚さ１．５ｃｍ、高さ２５ｃｍの直方体である。シミュレーション領域４０内に、直径１ｍｍの複数のガラス粒子を充填し、シミュレーション領域４０の底面からシミュレーション領域４０内にガスを導入する。粒子密度ρ_ｐを２５００ｋｇ／ｍ^３とした。粒子比熱ｃを８４０Ｊ／ｋｇ／Ｋとし、ガス定圧比熱ｃ_ｐ，ｆを１０１０Ｊ／ｋｇ／Ｋとし、ガス粘性係数μを２．０×１０^－５Ｐａ・ｓとした。シミュレーション領域４０内に充填する粒子の質量の合計を７５ｇとした。粒子の初期温度よりも低温のガスをシミュレーション領域４０内に導入した。ガスの流速を１．２０ｍ／ｓとした場合（流速が遅い場合）と、１．５４ｍ／ｓとした場合（流速が速い場合）とについてシミュレーションを行った。

【0037】

拡大率Ｋを２にして粗視化した流動層と、元の流動層との２つについてシミュレーションを行った。

【0038】

図６は、粗視化した流動層のシミュレーションにより求めた粗視化粒子の位置及び温度を、時系列で示す図である。図６の左から１番目、２番目、３番目、及び４番目の図は、それぞれ冷却開始時点、冷却開始からの経過時間がｔ、２ｔ、及び３ｔの流動層の状態を示す。各粒子の濃さは粒子の温度を表しており、温度が高いほど濃く表されている。ガスの流入によって粒子が流動し、時間の経過とともに粒子の温度が低下していることがわかる。

【0039】

図７Ａ及び図７Ｂは、シミュレーション結果から求めた粒子の平均温度の時間変化を示すグラフである。横軸は冷却開始からの経過時間を任意単位で表し、縦軸は粒子の平均温度を、初期温度を基準とした相対値で表す。図７Ａは、ガス流速が遅い場合を示し、図７Ｂは、ガス流速が速い場合を示している。グラフ中の破線は粗視化前の流動層のシミュレーション結果を示し、実線は粗視化後の流動層のシミュレーション結果を示している。参考のために、非特許文献１に示された実験結果による粒子の温度変化を丸記号で示している。

【0040】

図７Ａ及び図７Ｂに示したシミュレーション結果から、本参考例による方法で粗視化してシミュレーションを行っても、シミュレーション結果は実験結果とよく一致していることが確認できる。ガスの流速を速くすると、粒子の温度低下が速くなることも確認できる。このように、本参考例による粗視化の手法は、温度変化を伴う流動層の挙動のシミュレーションに適用することが可能である。

【0041】

［実施例］
次に、図８～図１１を参照して、実施例によるシミュレーション装置及びシミュレーション方法について説明する。

【0042】

上述の参考例では、流動層を構成する複数の粒子の粒径を均一と仮定した。以下に説明する実施例では、粒径の異なる複数の粒子からなる粉粒体をシミュレーション対象とする。また、上述の参考例では、粉粒体の挙動を解析する離散要素法（ＤＥＭ）と流体の流れ場を解析する数値流体力学（ＣＦＤ）とを連成させて流動層のシミュレーションを行ったが、以下に説明する実施例では、数値流体力学との連成を行うことなく、離散要素法を用いて粉粒体の挙動を解析する。なお、以下の実施例で説明する離散要素法を用いた解析手法を、流動層の解析に応用することも可能である。

【0043】

粉粒体の粒径の分布を表す手法として、ロジン－ラムラ分布が広く用いられている。ロジン－ラムラ分布Ｑ_１は、以下の式で表される。

【数16】

ここで、ｄ_ｐ１は粉粒体の粒子の粒径を表している。Ｑ（ｄ_ｐ１）は、粒径がｄ_ｐ１以下の粒子の質量の、全質量に対する比率（質量分率）であり、ｄ_ｅ１は長さの次元を持つパラメータ（以下、基準粒径パラメータという。）であり、ｎ_１は指数パラメータである。指数パラメータｎ_１は分布の広がり具合を表す。なお、粒径分布の範囲を決めるために、最小粒径ｄ_{ｐ１，ｍｉｎ}及び最大粒径ｄ_{ｐ１，ｍａｘ}を定義する。粒径ｄ_ｐ１を決めると、式（１６）の分布から、粒径ｄ_ｐ１以下の粒子の質量分率が求まる。

【0044】

粗視化後の粉粒体の粒子の粒径をｄ_ｐ２、粗視化後の粉粒体に適用されるロジン－ラムラ分布の基準粒径パラメータをｄ_ｅ２、指数パラメータをｎ_２で表す。粗視化後の粉粒体に適用されるロジン－ラムラ分布Ｑ_２は、以下の式で表される。

【数17】

【0045】

粗視化後のロジン－ラムラ分布の基準粒径パラメータｄ_ｅ２、最小粒径ｄ_{ｐ２，ｍｉｎ}及び最大粒径ｄ_{ｐ２，ｍａｘ}として、粗視化前の値に粗視化の係数を乗じた値を用い、指数パラメータｎ_２として、粗視化前の値を使用する。粗視化係数をＫで表すと、粗視化の変換則は、以下の式で表される。

【数18】

【0046】

式（１８）を式（１７）に代入すると、以下の式が得られる。

【数19】

【0047】

粗視化前の粉粒体に関するパラメータｄ_ｅ１、ｄ_{ｐ１，ｍｉｎ}、ｄ_{ｐ１,ｍａｘ}、ｎ_１は、シミュレーション対象の粉粒体の粒径の分布を調べることにより決定することができる。従って、粗視化後の粒子の粒径分布として、式（１９）のロジン－ラムラ分布が決まる。粗視化後の粉粒体の挙動をシミュレーションする際には、式（１９）の分布を用いればよい。

【0048】

図８は、粗視化前及び粗視化後の粉粒体の粒径の質量分率の一例を示すグラフである。横軸は、粒径を単位「μｍ」とした対数目盛で表し、縦軸はある粒径範囲に含まれる粒子の質量分率を単位「％」で表す。図８のグラフ中の破線及び実線は、それぞれ粗視化前の粉粒体及び粗視化後の粉粒体の粒径分布を示す。図８に示したグラフは、ロジン－ラムラ分布を粒径で微分したものに相当する。粗視化後の粒径分布を示すグラフは、粗視化前の粒径分布を示すグラフを、対数目盛上で平行移動したものに等しい。

【0049】

粉粒体の粒子の密度の変換則として、図２に示した変換則を用いる。すなわち、粗視化前の粉粒体の粒子の密度をρ_ｐ１で表し、粗視化後の粉粒体の粒子の密度をρ_ｐ２で表すと、変換則は以下の式で表される。

【数20】

【0050】

次に、図９を参照して実施例によるシミュレーション方法について説明する。
図９は、実施例によるシミュレーション装置を用いたシミュレーション方法のフローチャートである。シミュレーション装置は、図３に示したように、処理装置３０、入力装置３８、及び出力装置３９を含む。

【0051】

まず、シミュレーション対象の粉粒体の粒径の分布を決定する（ステップＳＡ１）。粒径分布の測定には、例えば画像法、ふるい分け法、沈降法、光回折法等を用いることができる。測定された粒径分布に基づいて、ロジン－ラムラ分布の基準粒径パラメータｄ_ｅ１及び指数パラメータｎ_１を決定することができる。さらに、最小粒径ｄ_{ｐ１，ｍｉｎ}及び最大粒径ｄ_{ｐ１，ｍａｘ}を決定することができる。

【0052】

粒径の分布を決定した後、粗視化係数Ｋを決定する（ステップＳＡ２）。粗視化係数Ｋは、シミュレーション対象の粉粒体の粒子数、シミュレーションに用いるコンピュータの処理速度等に基づいて決定するとよい。

【0053】

粗視化係数Ｋを決定したら、シミュレーション対象の粉粒体の粒径の分布（すなわち、基準粒径パラメータｄ_ｅ１と指数パラメータｎ_１）、粗視化係数Ｋ、その他のシミュレーション条件を、シミュレーション装置の入力装置３８（図３）を操作してシミュレーション装置に入力する（ステップＳＡ３）。その他のシミュレーション条件には、例えば、粒子の密度、最小粒径ｄ_{ｐ１，ｍｉｎ}、最大粒径ｄ_{ｐ１，ｍａｘ}、初期条件、境界条件、粒子に作用する外力、解析終了条件等が含まれる。

【0054】

これらのシミュレーション条件が入力されたら、処理装置３０（図３）は、粗視化後の粉粒体の粒径の分布を算出する。具体的には、式（１８）を適用して、粗視化後の粉粒体に適用されるロジン－ラムラ分布の基準粒径パラメータｄ_ｅ２、指数パラメータｎ_２、最小粒径ｄ_{ｐ２，ｍｉｎ}、最大粒径ｄ_{ｐ２，ｍａｘ}等を決定する。さらに、式（２０）を適用して、粗視化後の粒子の密度を算出する。

【0055】

処理装置３０は、粗視化後の粉粒体に適用されるロジン－ラムラ分布に基づいて、粗視化後の粉粒体の挙動を、離散要素法を用いてシミュレーションする（ステップＳＡ５）。例えば、粗視化後の粉粒体の複数の粒子の大きさを定義する。さらに、入力された初期条件に基づいて、粗視化後の粒子の初期状態を決定する。その後、この初期状態からの粒子の挙動を解析する。

【0056】

シミュレーションの演算が終了すると、処理装置３０は、シミュレーション結果と粗視化係数Ｋとを関連付けて出力装置３９（図３）に出力する。例えば図６に示すように、粗視化された粒子の位置を時系列で画像として表示するとともに、表示画面内に粗視化係数Ｋの値を数値で表示する。粒子の位置を表示する画像において、粗視化後の粉粒体の粒子の粒径に応じて、色相、彩度、明度の少なくとも１つを異ならせるとよい。

【0057】

次に、図１０Ａ～図１０Ｃを参照して、実際に実施例による方法でシミュレーションを行った結果について説明する。このシミュレーションにおいては、初期充填領域として円筒を採用し、水平な平坦面の上に粉粒体を自然落下させたときの粉粒体の挙動を求めた。粗視化前の粉粒体、粗視化係数Ｋ＝２０で粗視化した粉粒体、及び粗視化係数Ｋ＝３０で粗視化した粉粒体についてシミュレーションを行った。粗視化前の粉粒体の粒径分布を表すロジン－ラムラ分布の基準粒径パラメータｄ_ｅ１＝０．１１８１ｍｍ、指数パラメータｎ_１＝２．３２１、最小粒径ｄ_{ｐ１，ｍｉｎ}＝０．０５ｍｍ、最大粒径ｄ_{ｐ１，ｍａｘ}＝０．２ｍｍとした。

【0058】

実際に粗視化してシミュレーションを行う場合には、粗視化によって粉粒体の粒子数が減少するが、本シミュレーションにおいては、粗視化の前後で粒子の個数がほぼ等しくなるように、初期充填領域の直径及び高さを粗視化係数Ｋに応じて変更している。具体的には、粗視化前の粉粒体の初期充填領域は、半径３ｍｍ、高さ５ｍｍの円筒である。粗視化係数Ｋ＝２０の粉粒体の初期充填領域は、半径６０ｍｍ、高さ１００ｍｍの円筒である。粗視化係数Ｋ＝３０の粉粒体の初期充填領域は、半径９０ｍｍ、高さ１５０ｍｍの円筒である。

【0059】

図１０Ａは、初期充填領域に充填されている状態における粗視化前の粉粒体の側面図である。図１０Ａでは、粒子の大きさに応じて、各粒子の明度を異ならせて示している。

【0060】

図１０Ｂ及び図１０Ｃは、それぞれ粗視化前の粉粒体、及び粗視化係数Ｋ＝３０で粗視化した粉粒体の、十分時間が経過した後の分布を示す図である。初期状態の粉粒体が自然落下することにより、水平面上に粉粒体の山が形成されている。図１０Ｂ及び図１０Ｃにおいても、図１０Ａと同様に、粒子の大きさに応じて、各粒子の明度を異ならせて示している。なお、図１０Ｂと図１０Ｃとでは、縮尺が異なっている。図１０Ｂは図１０Ｃと比べて、相対的に大きな縮尺で表している。

【0061】

図１０Ｂ及び図１０Ｃに示した粉粒体の分布から安息角を求めた。粗視化前の粉粒体（図１０Ｂ）からなる山の安息角は２７．９°であり、粗視化係数Ｋ＝３０で粗視化した粉粒体（図１０Ｃ）からなる山の安息角は２９．１°であった。なお、粗視化係数Ｋ＝２０で粗視化した粉粒体からなる山の安息角は２８．７°であった。このように、粗視化前、及び粗視化後の粉粒体について、シミュレーションに求めた安息角はほぼ等しい値になった。このシミュレーション結果から、実施例による粗視化手法が、粒径が均一でない粉粒体の挙動のシミュレーションに適用可能であることが確認された。

【0062】

図１１は、実際に粒径が均一ではない珪砂を自然落下させて形成した珪砂の山を撮影した図である。この珪砂からなる粉粒体を粗視化してシミュレーションを行って求めた安息角は、実際の珪砂の山の安息角の測定結果とほぼ等しい値であった。この評価実験によっても、上記実施例による粗視化手法が有効であることが確認された。

【0063】

次に、上記実施例の変形例について説明する。上述の図１０Ａ～図１０Ｃに示した例では、粉粒体を自然落下させたときの粉粒体の挙動を実際にシミュレーションした。上記実施例による粗視化手法は、その他に、複数の固体粒子からなる粉粒体の挙動が重要となる種々の装置、例えばスクリューフィーダ、コークス炉、バグフィルタ、流動層装置における粉粒体のシミュレーションに適用することができる。

【0064】

また、上記実施例では、粗視化による粒子の物性値の変換則として、式（２０）に示した密度の変換則を用いたが、粒子の物性値の変換則として、その他の変換則を適用してもよい。

【0065】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

【符号の説明】

【0066】

１０ シミュレーション対象領域
１１ 粒子
１２ 流体
２０ 粗視化後のシミュレーション対象領域
２１ 粗視化後の粒子
２２ 粗視化後の流体
３０ 処理装置
３１ シミュレーション条件取得部
３２ 拡大率取得部
３３ 演算部
３４ 出力制御部
３８ 入力装置
３９ 出力装置
４０ シミュレーション領域

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】