(57)【特許請求の範囲】
【請求項１】
２次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析手段と、
前記２次元領域を分割する複数の矩形領域を設定する矩形領域設定手段と、
前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出手段と、
前記複数の矩形領域に対応する複数のセルを設定するセル設定手段と、
前記セルの体積をそれぞれ算出するセル体積算出手段と、
前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出手段と、
前記セル体積算出手段により算出された、前記セルの体積と、前記粒子体積算出手段により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出手段と、を有し、
前記セルは、前記矩形領域の縦の辺である第１の辺と、前記矩形領域の横の辺である第２の辺と、当該第１の辺及び当該第２の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出手段により算出された平均径の長さを有する第３の辺とにより定められる平行六面体であり、前記第３の辺の中央における当該第３の辺に沿う方向の位置が、当該第３の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする粒子存在場解析装置。
【請求項２】
前記平均粒子径算出手段は、前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を算出することを特徴とする請求項１に記載の粒子存在場解析装置。
【請求項３】
前記空隙率算出手段により算出された粒子の空隙率のそれぞれに、一定値を加算又は減算する空隙率補正手段を更に有することを特徴とする請求項１又は２に記載の粒子存在場解析装置。
【請求項４】
前記粒子が存在する３次元領域における流体解析を、前記空隙率を用いて行う流体解析手段を有することを特徴とする請求項１〜３の何れか１項に記載の粒子存在場解析装置。
【請求項５】
前記粒子は、高炉に装入される粒子であることを特徴とする請求項１〜４の何れか１項に記載の粒子存在場解析装置。
【請求項６】
２次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析工程と、
前記２次元領域を分割する複数の矩形領域を設定する矩形領域設定工程と、
前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出工程と、
前記複数の矩形領域に対応する複数のセルを設定するセル設定工程と、
前記セルの体積をそれぞれ算出するセル体積算出工程と、
前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出工程と、
前記セル体積算出工程により算出された、前記セルの体積と、前記粒子体積算出工程により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出工程と、を有し、
前記セルは、前記矩形領域の縦の辺である第１の辺と、前記矩形領域の横の辺である第２の辺と、当該第１の辺及び当該第２の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出工程により算出された平均径の長さを有する第３の辺とにより定められる平行六面体であり、前記第３の辺の中央における当該第３の辺に沿う方向の位置が、当該第３の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする粒子存在場解析方法。
【請求項７】
前記平均粒子径算出工程は、前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を算出することを特徴とする請求項６に記載の粒子存在場解析方法。
【請求項８】
前記空隙率算出工程により算出された粒子の空隙率のそれぞれに、一定値を加算又は減算する空隙率補正工程を更に有することを特徴とする請求項６又は７に記載の粒子存在場解析方法。
【請求項９】
前記粒子が存在する３次元領域における流体解析を、前記空隙率を用いて行う流体解析工程を有することを特徴とする請求項６〜８の何れか１項に記載の粒子存在場解析方法。
【請求項１０】
前記粒子は、高炉に装入される粒子であることを特徴とする請求項６〜９の何れか１項に記載の粒子存在場解析方法。
【請求項１１】
請求項１〜５の何れか１項に記載の粒子存在場解析装置の各手段としてコンピュータを機能させることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子存在場解析装置、粒子存在場解析方法、及びコンピュータプログラムに関し、特に、粒子の空隙率を算出するために用いて好適なものである。
【背景技術】
【０００２】
例えば、鉄鋼業における製銑分野では、安定して高品質の溶銑を製造するために、高炉内における焼結鉱やコークス等の粒子の挙動を解析すると共に、高炉内におけるガスの流れを解析することを、コンピュータを使って解析することが行われている。
このような技術として、粒子の離散的な運動を記述する離散要素法（Distinct Element Method、以下の説明では必要に応じて「ＤＥＭ」と称する）と、気流の運動（ガスの流れ）を記述するNavier-Stokesの式とを連成して実高炉内の現象をシミュレーションする技術が提案されている（非特許文献１を参照）。この技術では、ＤＥＭにより解析した３次元領域における粒子の配置から粒子の空隙率を計算し、この空隙率をNavier-Stokesの式に与えることにより、粒子が装入されている高炉内におけるガスの流れを解析している。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】湯晋一、梅景俊彦、松崎眞六、門脇正具、国友和也、内藤誠章、廣川雄一、上原均共著，「ＣＯ2排出ミニマムを目指した実高炉内の多相（固気粉）流れの大規模シミュレーション」，平成20年度先端研究施設共用イノベーション創出事業（産業戦略利用）「地球シミュレータ産業戦略利用プログラム」利用成果報告書、独立行政法人海洋研究開発機構，2009年8月31日，ｐ．１５−２５
【非特許文献２】粉体工学会編、「粉体工学便覧」、日刊工業新聞社、初版、昭和61年2月28日、ｐ．１５
【発明の概要】
【発明が解決しようとする課題】
【０００４】
非特許文献１に記載の技術では、ＤＥＭによる解析を行って粒子の３次元領域での挙動を解析するようにしている。したがって、非特許文献１の技術では、実高炉の１／４のサイズの領域を計算領域とすると共に、実際には数億個存在する粒子の個数を、粒子の粒径を大きくして１０００万個の規模にし、且つ、地球シミュレータを用いて計算しても、ＤＥＭによる解析に、数か月のオーダの計算時間を要する（非特許文献１では、最大計算粒子数を合計で９００万個としている）。
【０００５】
このように、非特許文献１に記載の直接３次元をＤＥＭで扱う技術では、実高炉の内部における粒子の空隙率を算出するのに莫大な時間がかかるという問題点があった。そこで、２次元領域におけるＤＥＭによる解析を行って粒子の配置を求め、この粒子の配置から、２次元計算で空隙率を算出することも考えられる。
しかしながら、粒子の空隙率は３次元的な指標であるため、２次元領域における粒子の空隙率の値は、実際の空隙率の値と大きく異なることになる。
以上のように、従来の技術では、特に、計算対象となる粒子の数が多い場合には、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することが困難であるという問題点があった。
【０００６】
本発明は、以上のような問題点に鑑みてなされたものであり、計算対象となる粒子の数が多くなっても、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することを目的とする。
【課題を解決するための手段】
【０００７】
本発明の粒子存在場解析装置は、２次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析手段と、前記２次元領域を分割する複数の矩形領域を設定する矩形領域設定手段と、前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出手段と、前記複数の矩形領域に対応する複数のセルを設定するセル設定手段と、前記セルの体積をそれぞれ算出するセル体積算出手段と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出手段と、前記セル体積算出手段により算出された、前記セルの体積と、前記粒子体積算出手段により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出手段と、を有し、前記セルは、前記矩形領域の縦の辺である第１の辺と、前記矩形領域の横の辺である第２の辺と、当該第１の辺及び当該第２の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出手段により算出された平均径の長さを有する第３の辺とにより定められる平行六面体であり、前記第３の辺の中央における当該第３の辺に沿う方向の位置が、当該第３の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする。
【０００８】
本発明の粒子存在場解析方法は、２次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析工程と、前記２次元領域を分割する複数の矩形領域を設定する矩形領域設定工程と、前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出工程と、前記複数の矩形領域に対応する複数のセルを設定するセル設定工程と、前記セルの体積をそれぞれ算出するセル体積算出工程と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出工程と、前記セル体積算出工程により算出された、前記セルの体積と、前記粒子体積算出工程により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出工程と、を有し、前記セルは、前記矩形領域の縦の辺である第１の辺と、前記矩形領域の横の辺である第２の辺と、当該第１の辺及び当該第２の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出工程により算出された平均径の長さを有する第３の辺とにより定められる平行六面体であり、前記第３の辺の中央における当該第３の辺に沿う方向の位置が、当該第３の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする。
【０００９】
本発明のコンピュータプログラムは、前記粒子存在場解析装置の各手段としてコンピュータを機能させることを特徴とする。
【発明の効果】
【００１０】
本発明によれば、２次元領域における粒子の挙動を離散要素法により解析する。そして、当該２次元領域を分割する複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の平均径を個別に算出する。そして、当該複数の矩形領域に対応する複数のセルを設定する。このセルは、当該矩形領域の縦・横の辺である第１・第２の辺と、当該第１・第２の辺に垂直な辺であって、当該矩形領域に対して算出された平均径の長さを有する第３の辺とにより定められる平行六面体であり、第３の辺の中央における第３の辺に沿う方向の位置が、第３の辺に沿う方向の位置が夫々揃っている粒子の中心位置と揃う位置に設定される。そして、これらのセル毎に、粒子の空隙率を算出する。したがって、離散要素法により、３次元領域における粒子の挙動を解析しなくても、当該３次元領域における粒子の挙動を計算することにより得られる粒子の空隙率に近い値の空隙率を算出することができる。よって、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することができる。
【図面の簡単な説明】
【００１１】
【図１】粒子存在場解析装置の機能的な構成の一例を示す図である。
【図２】相互に接触している２つの粒子の様子の一例と、粒子の挙動の一例を概念的に示す図である。
【図３】各粒子が配置されている２次元領域に対して設定された矩形領域の一例を示す図である。
【図４】ｙ軸方向及びｚ軸方向から見たセルの一例を示す図である。
【図５】１つのセルの中に全ての領域が含まれていない粒子の体積の計算式のパラメータを示す図である。
【図６】ＤＥＭによる解析を行って計算対象の領域に配置された各粒子の一例を概念的に示す図である。
【図７】各区間Ａ、Ｂ、Ｃにおける各粒子の粒度構成を示す図である。
【図８】矩形領域に含まれる粒子の平均径の計算方法を異ならせて算出した空隙率を示す図である。
【図９】粒子存在場解析装置の処理の一例を説明するフローチャートである。
【図１０】図９のステップＳ９０６におけるセル内粒子体積算出処理を説明するフローチャートである。
【発明を実施するための形態】
【００１２】
以下、図面を参照しながら、本発明の一実施形態を説明する。
図１は、粒子存在場解析装置１００の機能的な構成の一例を示す図である。粒子存在場解析装置１００は、例えば、高炉のベルレス式炉頂装入装置用の旋回シュートから高炉内に装入される粒子（焼結鉱やコークス等）の挙動の解析を行い、その結果を利用して、高炉内におけるガスの流れの解析（流体解析）を行うものである。ただし、粒子存在場解析装置１００が解析を行う対象は、このようなものに限定されるものではない。
粒子存在場解析装置１００は、例えば、ＣＰＵ、ＲＯＭ、ＲＡＭ、ＨＤＤ、各種インターフェース、及びディスプレイ等を備えた情報処理装置（パーソナルコンピュータ）を用いることにより実現できる。
【００１３】
（粒子挙動解析部１０１）
粒子挙動解析部１０１は、球状の粒子の２次元領域における挙動をＤＥＭにより解析する。粒子挙動解析部１０１は、例えば、ユーザによるユーザインターフェースの操作に基づいて設定された、粒子の情報に基づいて、各粒子を所望の位置に配置する。そして、粒子挙動解析部１０１は、粒子間に発生している力のバランスを崩す条件を与えて、各粒子を動かし、各粒子が他の粒子や構造物と接触しているか否かを判定する。
粒子が他の粒子と接触している場合、粒子挙動解析部１０１は、法線方向及びせん断方向の反発力を求めて並進の運動方程式を解くことと、重心周りの力のモーメントを求めて回転の運動方程式を解くこと、を各粒子のそれぞれについて行い、その結果に基づいて各粒子を移動させる。一方、粒子が構造物と接触している場合、粒子挙動解析部１０１は、法線方向及びせん断方向の反発力を求めて並進の運動方程式を解くことと、重心周りの力のモーメントを求めて回転の運動方程式を解くこと、を粒子についてのみ行い、構造物については行わない。
粒子挙動解析部１０１は、このような接触判定処理と各粒子の移動処理とを所定の時間隔毎に行う。そして、その結果から各粒子の所定の期間における挙動を得る。
【００１４】
具体的に、粒子挙動解析部１０１は、例えば、以下の（１）式により、粒子ｉ、ｊの接触領域における法線方向の反発力Ｆn,ijを求め、以下の（２）式により、粒子ｉ、ｊの接触領域におけるせん断方向の反発力Ｆt,ijを求め、以下の（３）式の並進の運動方程式を解き、以下の（４）式により回転の運動方程式を解く。
【００１５】
【数１】
0005464160.tif 000002
【００１６】
図２は、相互に接触している２つの粒子２０１、２０２の様子の一例（図２（ａ））と、粒子２０３の挙動の一例（図２（ｂ））を概念的に示す図である。図２を参照しながら、（１）式から（４）式のパラメータを説明する。
（１）式において、Ｋnは、法線方向のばね定数であり、例えば図２（ａ）の距離Ｈやδnに応じて設定することができる。ただし、Ｋnは、予め設定された一定値にすることもできる。Δｕn,ijは、２つの粒子ｉ、ｊ間の重心の法線方向における相対的な並進変位である。ηnは、法線方向の粘性定数であり、予め設定されるものである。ｎijは、法線方向の単位ベクトルである。Δｔは、所定の時間隔であり、予め設定されるものである。
【００１７】
（２）式において、min[Ａ，Ｂ]は、ＡとＢのうち小さい方を採用することを意味する。μは、摩擦係数であり、予め設定された値にすることができる。ただし、μは、状態等に応じて異ならせることもできる。ｔijはせん断方向の単位ベクトルである。Ｋtは、せん断方向のばね定数であり、例えば図２（ａ）の距離Ｈやδn、δtに応じて設定することができる。ただし、Ｋtは、予め設定された一定値にすることもできる。Δｕt,ijは、２つの粒子ｉ、ｊ間の重心のせん断方向における相対的な並進変位である。Δφijは、粒子の回転に起因する接触領域の相対的な回転変位であり、例えば、図２（ａ）の粒子２０１の回転変位φiと粒子２０２の回転変位φjとの差により表すことができる。ηtは、せん断方向の粘性定数であり、予め設定されるものである。Δｔは、所定の時間隔であり、予め設定されるものである。
【００１８】
（３）式において、ｖは、粒子の速度である（図２を参照）。ｖの上に付されている・（ドット）は時間微分を表す。Ｆは、せん断方向及び法線方向の反発力の総和である。ｍは、粒子の質量であり、予め設定されるものである。ｇは、重力加速度である。尚、粒子が他の粒子や構造物と接触していない場合、（３）式の右辺の第１項は０（ゼロ）になる。
（４）式において、ωは、粒子の角速度である（図２を参照）。ωの上に付されている・は時間微分を表す。Ｍは、力のモーメントの総和であり、せん断方向及び法線方向の反発力Ｆn,ij、Ｆt,ijから求められる。Ｉは、慣性モーメントであり、予め設定されるものである。
【００１９】
本実施形態では、２次元領域における粒子の挙動をＤＥＭにより解析するに際し、当該粒子を球状とする。３次元の成分を有する指標（粒子の物性値等）を考慮して計算するためである。すなわち、２次元領域における粒子の挙動をＤＥＭにより解析するに際し、粒子に働く作用力の絶対値については３次元の情報として計算し、粒子に働く作用力の方向については２次元の情報として計算する。
【００２０】
また、高炉内に装入される粒子の２次元領域の挙動を解析する場合、高炉の半径方向と高さ方向とにより定まる２次元領域の挙動を解析することになる。高炉の奥行き方向への粒子の挙動は、解析の結果への影響が小さいためである。したがって、高炉内に装入される粒子の２次元領域の挙動を解析する場合、粒子に働く作用力の方向は、高炉を、その軸に沿って切った断面の面内方向となる。
以上のようにして、計算対象となる２次元領域における各粒子の配置が、所定の時間毎に得られる。
【００２１】
尚、粒子の挙動を解析する手法（粒子挙動解析部１０１が実行する解析手法）は、ＤＥＭにおける粒子の挙動についての種々の公知の解析手法を適用することができ、前述した（１）式〜（４）式を用いた手法に限定されるものではない。
粒子挙動解析部１０１は、例えば、ＣＰＵが、ＲＡＭ等から、粒子の情報を読み出して前述した演算処理を実行して、各時間における各粒子の２次元領域における座標を得てＲＡＭ等に記憶することにより実現される。
【００２２】
（矩形領域設定部１０２）
矩形領域設定部１０２は、粒子挙動解析部１０１により粒子の挙動が解析された２次元領域を分割する複数の矩形領域（所謂メッシュ）を設定する。前述したように、粒子挙動解析部１０１によって、２次元領域における所定時間毎の各粒子の配置が粒子の挙動として得られる。矩形領域設定部１０２は、この各粒子が配置されている２次元領域に対して、複数の矩形領域を設定する。高炉内に装入される粒子の２次元領域の挙動を解析する場合、矩形領域設定部１０２は、例えば、各粒子の動きが止まった時点において、各粒子が配置されている２次元領域の情報を取得し、当該２次元領域に対して、複数の矩形領域を設定する。ただし、必ずしも各粒子の動きが止まった時点における２次元領域に対して複数の矩形領域を設定する必要はなく、任意の時点における２次元領域に対して複数の矩形領域を設定することができる。
【００２３】
図３は、各粒子が配置されている２次元領域に対して設定される矩形領域の一例を示す図である。尚、各粒子が配置されている２次元領域は、例えば、後述する図６（ａ）に示すような領域である。
図３に示すように、本実施形態では、各粒子が配置されている２次元領域に、ｘ軸方向の長さがｄｘであり、ｚ軸方向の長さがｄｚである矩形領域が、ｘ軸方向及びｚ軸方向のそれぞれに連続的に設定される。図３に示すように、１つの矩形領域の中に全ての領域が含まれている粒子もあれば、複数の矩形領域に跨って存在する粒子もある。
尚、この２次元領域が高炉内の領域である場合、ｘ軸方向は、高炉の半径方向であり、ｚ軸方向は、高炉の高さ方向となる。
【００２４】
本実施形態では、全ての矩形領域の大きさは同じであるとする。また、１つの矩形領域に１０個以上の粒子が入るように、矩形領域の大きさを決めるようにするのが好ましい。さらに、粒子挙動解析部１０１により解析が行われる２次元領域の座標が分かれば矩形領域を設定することができるので、粒子挙動解析部１０１による解析が行われる前に矩形領域が設定されるようにしてもよい。
矩形領域設定部１０２は、例えば、ＣＰＵが、各矩形領域の２次元領域（ｘ−ｚ平面）における座標（ｘ軸座標、ｚ軸座標）を得てＲＡＭ等に記憶することにより実現される。
【００２５】
（平均粒子径算出部１０３）
平均粒子径算出部１０３は、矩形領域設定部１０２により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の平均径を算出する。本実施形態では、最も好ましい平均径として、長さ平均径を採用する。この理由については後述する。
ここで、１つの矩形領域の中に全ての領域が含まれていない粒子については、当該粒子の中心が含まれている矩形領域に属しているものとして長さ平均径を算出する。
平均粒子径算出部１０３は、例えば、ＣＰＵが、矩形領域に含まれる球状の粒子の長さ平均径を、全ての矩形領域について得てＲＡＭ等に記憶することにより実現される。
【００２６】
（セル設定部１０４）
セル設定部１０４は、「矩形領域の縦の辺であるｚ軸方向の辺（図３の長さｄｚの辺）」と、「矩形領域の横の辺であるｘ軸方向の辺（図３の長さｄｘの辺）」と、「当該ｚ軸方向の辺及び当該ｘ軸方向の辺に垂直なｙ軸方向の辺であって、当該矩形領域に対して平均粒子径算出部１０３により算出された長さ平均径の長さを有するｙ軸方向の辺」と、により定められる平行六面体であって、当該ｙ軸方向の辺の中央におけるｙ軸方向の位置が、球場の粒子のｙ軸方向の中心位置とそれぞれ揃えられた複数の平行六面体を設定する。以下の説明では、この「平行六面体」を必要に応じて「セル」と称する。
【００２７】
図４は、ｙ軸方向及びｚ軸方向から見たセルの一例を示す図である。尚、図４では、表記の都合上、各セル４０１ａ〜４０１ｃに含まれている粒子の数を１〜２個としているが、実際には、各セル４０１ａ〜４０１ｃには、これよりも多くの数の粒子が含まれる。
セル設定部１０４は、矩形領域設定部１０２により設定された矩形領域の「ｚ軸方向の辺（図３の長さｄｚの辺）と、ｘ軸方向の辺（図３の長さｄｘの辺）」の情報を取得すると共に、当該矩形領域に対して平均粒子径算出部１０３により算出された長さ平均径の情報を取得する。そして、セル設定部１０４は、「当該矩形領域のｚ軸方向の辺」と、「当該矩形領域のｘ軸方向の辺」と、「当該矩形領域の頂点を中央点とし、ｙ軸方向に当該矩形領域に対して設定された長さ平均径の長さ分だけ左右（ｙ軸正負の両方向に均等）に伸ばしたｙ軸方向の辺」とにより定まる平行六面体（直方体又は立方体）を、矩形領域設定部１０２により設定された全ての矩形領域について定める。
【００２８】
図４に示す例では、「当該矩形領域の頂点を中央点とし、ｙ軸方向に当該矩形領域に対して設定された長さ平均径の長さ分だけ左右に伸ばしたｙ軸方向の辺」は、長さがｄｙ1、ｄｙ2、ｄｙ3の辺である。
このようにして得られた平行六面体が、図４に示すセル４０１ａ〜４０１ｃになり、各セル４０１ａ〜４０１ｃのｙ軸方向の中央の位置は、ｘ軸及びｚ軸に平行な面４０２に位置する。
【００２９】
各矩形領域に含まれる粒子の大きさ及び数は矩形領域によって異なるので、各矩形領域に含まれる粒子の長さ平均径も矩形領域によって異なる。よって、図４に示すように、各セル４０１ａ〜４０１ｃのｙ軸方向の長さｄｙ1〜ｄｙ3もセル４０１ａ〜４０１ｃによって異なることになる（ただし、同じになる場合もある）。
尚、図４では、セル４０１の数を３つしか示していないが、セル設定部１０４により設定されるセル４０１の数は、矩形領域設定部１０２により設定された矩形領域の数と同じ数になる。
セル設定部１０４は、例えば、ＣＰＵが、セル４０１の３次元座標（ｘ軸座標、ｙ軸座標、ｚ軸座標）を得てＲＡＭ等に記憶することにより実現される。
【００３０】
（セル体積算出部１０５）
セル体積算出部１０５は、セル設定部１０４によって設定されたセル４０１の体積をそれぞれ算出する。前述したように、セルの形状は平行六面体（直方体又は立方体）であるので、各辺の長さ（ｄｘ、ｄｙ、ｄｚ）を掛け合わせることにより、各セル４０１の体積を得ることができる。すなわち、各セル４０１の体積Ｖcは、以下の（５）式により得られる。
【００３１】
【数２】
0005464160.tif 000003
【００３２】
セル体積算出部１０５は、例えば、ＣＰＵが、各セル４０１の体積を得てＲＡＭ等に記憶することにより実現される。
【００３３】
（粒子体積算出部１０６）
粒子体積算出部１０６は、セル４０１の中に含まれる粒子の体積を算出する。
粒子体積算出部１０６は、図４に示すように、各セル４０１ａ〜４０１ｃのｙ軸方向の中央の位置（座標）と、各粒子４１１〜４１５の中心４１１ａ〜４１５ａのｙ軸方向の位置とが、それぞれ、ｘ軸及びｚ軸に平行な面４０２に位置している。
【００３４】
そして、粒子体積算出部１０６は、各セル４０１ａ〜４０１ｃの中に含まれている粒子４１１〜４１５の体積を算出する。
【００３５】
各粒子４１１〜４１５の形状は球である。そこで、粒子体積算出部１０６は、セル４０１の中に全ての領域が含まれている粒子については、当該セル４０１の中に含まれている粒子の体積Ｖpを、以下の（６）式で算出する。
【００３６】
【数３】
0005464160.tif 000004
【００３７】
（６）式において、Ｄpは粒子（球）の直径である。図４に示す例では、粒子４１３〜４１５が、このような粒子に該当する。
また、粒子体積算出部１０６は、１つのセル４０１の中に全ての領域が含まれていない粒子については、当該粒子の領域のうち、当該セル４０１の中に含まれている領域の体積を、当該セル４０１の中に含まれている粒子の体積Ｖpとして算出する。
【００３８】
図５は、１つのセル４０１の中に全ての領域が含まれていない粒子の体積の計算式のパラメータを示す図である。
図５において、Ｄpは、粒子４１６の直径である。ａは、粒子４１６の中心４１６ａと、粒子４１６が跨いでいる「セル４０１ｄのｘ−ｚ平面上の辺」との間の最短距離である。これらＤp、ａを用いると、粒子体積算出部１０６は、セル４０１の中に一部の領域が含まれている粒子４１６については、当該粒子４１６の領域のうち、セル４０１ｄの中に含まれている領域の体積を、当該セル４０１ｄの中に含まれている粒子の体積Ｖpとして、以下の（７）式により算出する。
【００３９】
【数４】
0005464160.tif 000005
【００４０】
また、粒子体積算出部１０６は、セル４０１の中に一部の領域が含まれている粒子４１６については、当該粒子４１６の領域のうち、セル４０１ｄに含まれていない領域の体積を、当該領域が属しているセル４０１ｅの中に含まれている粒子の体積Ｖpとして、以下の（８）式により算出する。
【００４１】
【数５】
0005464160.tif 000006
【００４２】
図４に示す例では、粒子４１１、４１２が、このような粒子に該当する。
尚、１つのセル４０１の中に全ての領域が含まれていない粒子の領域のうち、当該セル４０１の中に含まれている領域の体積を算出する方法（計算式）は、前述した方法（（７）式、（８）式）に限定されず、粒子とセル４０１との位置関係に応じて適宜設定することができる。
そして、粒子体積算出部１０６は、１つのセル４０１の中に少なくとも一部の領域が含まれている粒子の全てについて（７）式又は（８）式により計算した「セル４０１の中に含まれている粒子の体積Ｖp」の総和を計算する。
粒子体積算出部１０６は、例えば、ＣＰＵが、各セル４０１の中に含まれている粒子の体積Ｖpを得て、ＲＡＭ等に記憶することにより実現される。
【００４３】
（空隙率算出部１０７）
空隙率算出部１０７は、セル体積算出部１０５で算出された「各セル４０１の体積Ｖc」と、粒子体積算出部１０６により算出された「各セル４０１の中に含まれている粒子の体積Ｖpの総和」とに基づいて、粒子の空隙率を、セル設定部１０４により設定された複数のセル４０１のそれぞれについて算出する。
【００４４】
空隙率算出部１０７は、以下の（９）式により、各矩形領域における粒子の仮空隙率ε´［−］を個別に算出する。
【００４５】
【数６】
0005464160.tif 000007
【００４６】
（９）式において、ｎは、粒子体積算出部１０６により、同一のセル４０１の中に含まれているとして計算された粒子の体積Ｖpの数である。例えば、図４に示す例では、セル４０１ａの中に含まれているとして計算された粒子の体積Ｖpの数は「２」であり、セル４０１ｂの中に含まれているとして計算された粒子の体積Ｖpの数は「３」であり、セル４０１ｃの中に含まれているとして計算された粒子の体積Ｖpの数は「１」である。
【００４７】
そして、空隙率算出部１０７は、算出した各矩形領域における粒子の仮空隙率ε´のそれぞれに定数（＝０．０２）を加算した値を、各矩形領域における粒子の空隙率εとして算出する。定数（＝０．０２）を加算して仮空隙率ε´を補正する理由については、後述する。
空隙率算出部１０７は、例えば、ＣＰＵが、各セル４０１における粒子の空隙率εを得て、ＲＡＭやＨＤＤ等に記憶することにより実現される。
【００４８】
（流体解析部１０８）
流体解析部１０８は、空隙率算出部１０７により算出された「各セル４０１における粒子の空隙率ε」を用いて、粒子が配置された３次元領域におけるガスの流れの解析（流体解析）を行い、その結果をディスプレイに表示する。流体解析の手法は、例えば、非特許文献１に記載されているので、ここでは、その詳細な説明を省略する。
流体解析部１０８は、例えば、ＣＰＵが、粒子が配置された３次元領域におけるガスの流れの解析結果を得て、そのデータをディスプレイに出力することにより実現される。
【００４９】
（長さ平均径を採用する理由と空隙率を調整する理由）
本実施形態では、以下に示す種々の粒度構成を有する条件で、２次元領域（２次元の計算場）における各粒子の挙動をＤＥＭにより計算し、その計算結果から、各種の平均径を用いて前述したようにして空隙率の値を求めると共に、３次元領域（３次元の計算場）における各粒子の挙動をＤＥＭにより計算し、それぞれの計算結果から空隙率の値を求めた。そして、これらの空隙率の値を比較した結果、一定値を用いていた従来法に比較して、平均粒子径でセル体積を設定する方法の方が、種々の粒度構成を有する条件でも３次元値に近い値を示し、かつ試みたいずれの平均径の場合も高炉内ガス流れ解析での使用に際して十分な精度を有していることが分かった。さらに、各種の平均径のうち、長さ平均径を採用することが好ましいということと、長さ平均径を用いて求めた空隙率の値を調整すると、３次元領域（３次元の計算場）で求めた空隙率の値に非常に近い空隙率が得られることを見出した。以下に、これらについて詳述する。
【００５０】
図６は、ＤＥＭによる解析を行って計算対象の領域に配置された各粒子の一例を概念的に示す図である。図６（ａ）は、２次元領域（２次元の計算場）における各粒子の挙動を解析して得られた各粒子の配置を示す図である。図６（ｂ）は、図６（ａ）（本実施形態）と比較するための図であり、３次元領域（３次元の計算場）における各粒子の挙動を解析して得られた各粒子の配置を示す図である。
【００５１】
図６（ａ）に示すように、ＤＥＭによって、２次元領域における各粒子の挙動を解析すると、２次元領域（ｘ−ｚ平面の領域）と直交する方向（ｙ軸方向）に積み重ならないように、つまり球場の粒子のｙ軸方向の中心位置が全て揃った球状の粒子が、ｘ−ｚ平面に並べられる。これに対し、図６（ｂ）に示すように、ＤＥＭによって、３次元領域における各粒子の挙動を解析すると、球状の粒子は、３次元領域（ｘ−ｙ−ｚ空間）に並べられるので、ｙ軸方向にも積み重なる。
【００５２】
図７は、図６に示す各区間Ａ、Ｂ、Ｃにおける各粒子の粒度構成を示す図である。具体的に、図７（ａ）は、図６（ａ）に示した各区間Ａ、Ｂ、Ｃにおける各粒子の面積比と、図６（ｂ）に示した区間Ａ、Ｂ、Ｃにおける各粒子の体積比とを示す。一方、図７（ｂ）は、図６（ａ）及び図６（ｂ）に示した各区間Ａ、Ｂ、Ｃにおける粒子個数比を示す。
図６は、各区間Ａ、Ｂ、Ｃのそれぞれにおいて、それぞれが接触しないように粒子を初期配置し、上方から（図３に示す矢印の方向に）各粒子を重力沈降させ、動きが止まったときの各粒子の配置を示す。尚、図６では、表記の都合上（空隙があることを示すために）、粒子と粒子との間に隙間があるが、実際には、図６に示す状態よりも各粒子は密に配置されている。図７に示すように、各区間Ａ、Ｂ、Ｃに落下させる粒子は、直径（粒子径）が２０．０［ｍｍ］、３７．５［ｍｍ］、５５．０［ｍｍ］の３種類の粒子である。
【００５３】
ここで、図６（ａ）に示す区間Ａ、Ｂ、Ｃは、２次元領域であるのに対し、図６（ｂ）に示す区間Ａ、Ｂ、Ｃは、３次元領域である。そこで、区間Ａ、Ｂ、Ｃのそれぞれにおいて、これら３種類の粒子の割合が、図６（ａ）に示す区間Ａ、Ｂ、Ｃと、図６（ｂ）に示す区間Ａ、Ｂ、Ｃとで同じになるようにしている。
【００５４】
具体的に、図７（ａ）に示すように、各粒子の図６（ａ）に示す各区間Ａ、Ｂ、Ｃにおける面積比と、各粒子の図６（ｂ）に示す各区間Ａ、Ｂ、Ｃにおける体積比とが同じになるようにしている。
例えば、図７（ａ）に示すように、図６（ａ）に示す区間Ａにおける「粒子径が２０．０［ｍｍ］の粒子の面積の総和」と、図６（ａ）に示す区間Ｂにおける「粒子径が２０．０［ｍｍ］の粒子の面積の総和」と、図６（ａ）に示す区間Ｃにおける「粒子径が２０．０［ｍｍ］の粒子の面積の総和」との比が「１：１：５」になるようにしている。
そこで、図７（ｂ）に示すように、図６（ｂ）に示す区間Ａにおける「粒子径が２０．０［ｍｍ］の粒子の体積の総和」と、図６（ｂ）に示す区間Ｂにおける「粒子径が２０．０［ｍｍ］の粒子の体積の総和」と、図６（ｂ）に示す区間Ｃにおける「粒子径が２０．０［ｍｍ］の粒子の体積の総和」との比が「１：１：５」となるようにしている。
【００５５】
また、各粒子の図６（ｂ）に示す各区間Ａ、Ｂ、Ｃにおける個数比（粒子個数比）と、各粒子の図６（ｂ）に示す各区間Ａ、Ｂ、Ｃにおける個数比（粒子個数比）は、図７に示す通りである。
例えば、図７（ａ）に示すように、図６（ａ）に示す区間Ａにおける「粒子径が２０．０［ｍｍ］の粒子の数」と、図６（ａ）に示す区間Ｂにおける「粒子径が２０．０［ｍｍ］の粒子の数」と、図６（ａ）に示す区間Ｃにおける「粒子径が２０．０［ｍｍ］の粒子の数」との比は、「８：８：４０＝１：１：５」となる。一方、図６（ｂ）に示す区間Ａにおける「粒子径が２０．０［ｍｍ］の粒子の数」と、図６（ｂ）に示す区間Ｂにおける「粒子径が２０．０［ｍｍ］の粒子の数」と、図６（ｂ）に示す区間Ｃにおける「粒子径が２０．０［ｍｍ］の粒子の数」との比も、「２１：２１：１０５＝１：１：５」となる。
【００５６】
図８は、矩形領域に含まれる粒子の平均径の計算方法を異ならせて算出した空隙率を示す図である。
図８に示す各点は、図３に示した区間Ａ、Ｂ、Ｃに、図４に示した割合で粒子を重力沈降させて、粒子の動きが止まった時点での各矩形領域の空隙率の値を示す。また、各区間Ａ、Ｂ、Ｃに対して、それぞれ６つの矩形領域を設定した（このことは、図８の各区間Ａ、Ｂ、Ｃに６つの点があることと対応する）。
【００５７】
また、図８において、「３Ｄ」は、ＤＥＭにより、３次元領域における粒子の挙動を解析し、その結果から算出した空隙率の値を示す（図３（ｂ）を参照）。この値が最も信頼できる値となる。
また、「平均径定数固定」とは、セルのｙ軸方向の長さを矩形領域によらずに３７．５［ｍｍ］に固定して算出した空隙率の値を示す（セルのｙ軸方向の長さを固定する以外は、本実施形態と同じ処理を行って得られた空隙率の値である（ただし、空隙率は（９）式で計算されるものであり、定数を加算していない））。
【００５８】
また、「個数平均径」、「面積平均径」、「体積平均径」、「平均表面積径」、「平均体積径」は、「長さ平均径」の代わりに、これらの平均径を使って算出した空隙率の値を示す（平均径の計算方法以外は、本実施形態と同じ処理を行って得られた空隙率の値である（ただし、空隙率は（９）式で計算されるものであり、定数を加算していない））。
【００５９】
また、「長さ平均径」は、定数を加算する前の空隙率（（９）式の仮空隙率ε´）である。
また、「長さ平均径＋０．０２」は、本実施形態と同じ処理を行って得られた空隙率である。
尚、「個数平均径」、「長さ平均径」、「面積平均径」、「体積平均径」、「平均表面積径」、「平均体積径」の計算方法については、非特許文献２に記載されている。また、「個数平均径」、「長さ平均径」、「面積平均径」、「体積平均径」については、それぞれ「重みつき個数平均径」、「重みつき長さ平均径」、「重みつき面積平均径」、「重みつき体積平均径」と称されることもある。
【００６０】
図８に示すように、程度の差はあるが、「個数平均径」、「長さ平均径」、「面積平均径」、「体積平均径」、「平均表面積径」、「平均体積径」のグラフの方が、「平均径定数固定」のグラフよりも、「３Ｄ」のグラフと、概形が似ている。区間Ａ、Ｂ、Ｃでは、３種類の粒子径の構成割合を極端に異ならせている（図４を参照）。よって、区間Ａ、Ｂ、Ｃの全てで、「３Ｄ」のグラフに対し、形と値とが近くなっていれば、粒子径がどのような条件であっても、精度よく空隙率を算出することができると言ってよい。したがって、粒子の平均径を矩形領域毎に算出し、各セルのｙ軸方向の長さを、算出した平均径に基づいて決定すると、精度よく空隙率を算出することができることが分かる。
【００６１】
特に、「長さ平均径」のグラフは、「３Ｄ」のグラフと、概形だけでなく、値も近くなっており、特に好ましいことも判る。そこで、本実施形態では、これらの平均径のうち、「長さ平均径」を採用することとした。
また、「長さ平均径」を用いて算出した空隙率に一定値を加算した空隙率の値と（図８の「長さ平均径」）、ＤＥＭにより３次元領域における粒子の挙動を解析した結果から算出した空隙率の値（図８の「３Ｄ」）との差を最小にするためには、当該一定値をいくつにすればよいのかを、最小二乗法により計算した。その結果、当該一定値は「０．０２」となった。よって、この例では、「長さ平均径」を用いて算出した空隙率の値に、当該一定値として「０．０２」を加算すれば、「３Ｄ」のグラフの値との差が非常に小さい空隙率を得ることができる。そこで、本実施形態では、「長さ平均径」を用いて算出した空隙率（仮空隙率ε´）に一定値（＝０．０２）を加算した値を各矩形領域の空隙率εとして採用することとした。
【００６２】
ただし、当該一定値を加算又は減算する前の空隙率の値よりも、当該一定値を加算又は減算した後の空隙率の値の方が、ＤＥＭにより３次元領域における粒子の挙動を解析した結果から算出した空隙率の値との差が小さくなるようにしていれば、必ずしも、「３Ｄ」のグラフの値との差が最小となるように当該一定値を計算する必要はない。また、「３Ｄ」のグラフの値との差が最小となるように当該一定値を計算するために用いる方法は、最小二乗法に限定されるものではない。
【００６３】
尚、「長さ平均径」以外の平均径を用いて算出した空隙率でも、当該空隙率に対し一定値を加算又は減算することにより、当該一定値を加算又は減算した空隙率の値と、「３Ｄ」のグラフの値との差を小さくすることができるので、このようにしてもよい。特に、「３Ｄ」のグラフの値との差が最小となるように、当該空隙率の値に対して加算又は減算する一定値を決定するのが好ましい。「３Ｄ」のグラフの値との差が最小となる一定値を得るための計算は、前述したように、例えば、最小二乗法を用いることにより実現することができる。
【００６４】
（動作フローチャート）
次に、図９のフローチャートを参照しながら、粒子存在場解析装置１００の処理の一例を説明する。
まず、ステップＳ９０１において、粒子挙動解析部１０１は、球状の粒子の２次元領域における挙動をＤＥＭにより解析する。これにより、所定の時間毎に、２次元領域における粒子の配置が得られる（図３（ａ）を参照）。
次に、ステップＳ９０２において、矩形領域設定部１０２は、ステップＳ９０１で粒子の挙動が解析された２次元領域を分割する複数の矩形領域を設定する（図３を参照）。
【００６５】
次に、ステップＳ９０３において、平均粒子径算出部１０３は、ステップＳ９０２で設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を個別に算出する。
次に、ステップＳ９０４において、セル設定部１０４は、ステップＳ９０２で設定された「矩形領域」の各辺をｘ軸、ｚ軸方向の辺とし、これらｘ軸、ｚ軸方向の辺と、ステップＳ９０３で算出された「当該矩形領域に含まれる球状の粒子の長さ平均径」の長さを有するｙ軸方向の辺とにより定まる平行六面体であって、当該ｙ軸方向の辺の中央の位置が夫々揃えられた複数の平行六面体を、セル４０１として設定する（図４を参照）。
【００６６】
次に、ステップＳ９０５において、セル体積算出部１０５は、ステップＳ９０４で設定された各セルの体積を算出する（（５）式を参照）。
次に、ステップＳ９０６において、粒子体積算出部１０６は、セル内粒子体積算出処理を実行して、セル４０１の中に含まれている粒子の体積Ｖpの総和を算出する（（６）式〜（８）式を参照）。尚、セル内粒子体積算出処理の詳細については、図１０を参照しながら後述する。
【００６７】
次に、ステップＳ９０７において、空隙率算出部１０７は、ステップＳ９０５で算出された「各セル４０１の体積Ｖc」と、ステップＳ９０６で算出された「各セル４０１の中に含まれている粒子の体積Ｖp」とに基づいて、ステップＳ９０４設定された複数のセル４０１のそれぞれにおける「粒子の空隙率ε」を算出する（粒子の空隙率εは（９）式の仮空隙率ε´に定数を加算した値となる）。
次に、ステップＳ９０８において、流体解析部１０８は、ステップＳ９０７で算出された「複数のセル４０１のそれぞれにおける『粒子の空隙率ε』」を用いて、粒子が配置された３次元領域におけるガスの流れの解析（流体解析）を行い、その結果をディスプレイに表示する。そして、図９のフローチャートによる処理を終了する。
【００６８】
次に、図１０のフローチャートを参照しながら、図９のステップＳ９０６におけるセル内粒子体積算出処理の詳細な説明を行う。尚、図１０のフローチャートによる処理は、ステップＳ９０４により設定された全てのセルについて行われる。
【００６９】
まず、ステップＳ１００１において、粒子体積算出部１０６は、計算対象となるセル４０１の中に少なくとも一部の領域が含まれている粒子のうち、未選択の粒子を１つ選択する。
次に、ステップＳ１００２において、粒子体積算出部１０６は、ステップＳ１００１で選択した粒子の全ての領域が、計算対象となるセル４０１の中に含まれているか否かを判定する。
【００７０】
この判定の結果、ステップＳ１００１で選択した粒子の全ての領域が、計算対象となるセル４０１の中に含まれている場合には、ステップＳ１００３に進む。
ステップＳ１００３に進むと、粒子体積算出部１０６は、ステップＳ１００１で選択した粒子の全体の体積を、計算対象となるセル４０１の中に含まれている粒子の体積Ｖpとして、（６）式により算出する。
【００７１】
一方、ステップＳ１００１で選択した粒子の全ての領域が、計算対象となるセル４０１の中に含まれていない場合には、ステップＳ１００４に進む。
ステップＳ１００４に進むと、粒子体積算出部１０６は、ステップＳ１００１で選択した粒子の領域のうち、当該セル４０１に含まれている領域の体積を、計算対象となるセル４０１の中に含まれている粒子の体積Ｖpとして、例えば（７）式又は（８）式により計算する。
【００７２】
以上のようにして、ステップＳ１００１で選択した粒子について、計算対象となるセル４０１の中に含まれている粒子の体積Ｖpが算出されると、ステップＳ１００５に進む。ステップＳ１００５に進むと、粒子体積算出部１０６は、計算対象となるセル４０１の中に少なくとも一部の領域が含まれている粒子を全て選択したか否かを判定する。この判定の結果、計算対象となるセル４０１の中に少なくとも一部の領域が含まれている粒子を全て選択していない場合には、ステップＳ１００１の処理に戻る。そして、計算対象となるセル４０１の中に少なくとも一部の領域が含まれている粒子の全てについて、計算対象となるセル４０１の中に含まれている粒子の体積Ｖpを算出するまで、ステップＳ１００１〜Ｓ１００５の処理を繰り返し行う。その後、粒子体積算出部１０６は、計算対象となるセル４０１の中に含まれている粒子の体積Ｖpの総和を計算する。
【００７３】
（まとめ）
以上のように本実施形態では、球状の粒子の２次元領域における挙動をＤＥＭにより解析した後、当該２次元領域を分割する複数の矩形領域を設定する。複数の矩形領域のそれぞれについて、当該矩形領域に含まれる粒子の長さ平均径を個別に算出する。そして、そして、当該複数の矩形領域に対応する複数のセル４０１を設定する。このセル４０１は、矩形領域のｘ軸、ｚ軸方向の辺と、当該矩形領域における長さ平均径の長さを有するｙ軸方向の辺と、により定まる平行六面体であり、当該ｙ軸方向の辺の中央の位置が、当該ｙ軸方向の位置が夫々揃っている球状の粒子の中心位置と揃う位置に設定される。そして、各セル４０１に含まれる粒子の割合を求めて、各矩形領域における粒子の空隙率εを算出する。したがって、ＤＥＭによって、３次元領域における粒子の挙動を計算しなくても、当該３次元領域における粒子の挙動を計算することにより得られる粒子の空隙率εに近い値の空隙率εを算出することができる。例えば、高炉の炉頂から羽口の間の荷下がりを、３次元計算場でＤＥＭにより計算すると、計算時間が約２ヶ月であったのに対し、２次元計算場でＤＥＭにより計算すると、計算時間が約５日であった。また、これらにおける空隙率εの差は、概ね図８に示した結果と同じになった。よって、計算対象となる粒子の数が多くなっても、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することができる。
【００７４】
＜変形例＞
本実施形態では、流体解析を行うようにしたが、粒子の空隙率εをセル４０１毎に算出するようにしていれば、必ずしも流体解析を行う必要はない。また、粒子の空隙率εをディスプレイに表示してユーザに提示するようにしてもよい。
【００７５】
尚、以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、プログラムをコンピュータに供給するための手段、例えばかかるプログラムを記録したＣＤ−ＲＯＭ等のコンピュータ読み取り可能な記録媒体、又はかかるプログラムを伝送する伝送媒体も本発明の実施の形態として適用することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体などのプログラムプロダクトも本発明の実施の形態として適用することができる。前記のプログラム、コンピュータ読み取り可能な記録媒体、伝送媒体及びプログラムプロダクトは、本発明の範疇に含まれる。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
【符号の説明】
【００７６】
１００    粒子存在場解析装置
１０１    粒子挙動解析部
１０２    矩形領域設定部
１０３    平均粒子径算出部
１０４    セル設定部
１０５    セル体積算出部
１０６    粒子体積算出部
１０７    空隙率算出部
１０８    流体解析部
４０１    セル
４１１〜４１６    粒子
【図１】
0005464160.tif 000008
【図２】
0005464160.tif 000009
【図３】
0005464160.tif 000010
【図４】
0005464160.tif 000011
【図５】
0005464160.tif 000012
【図６】
0005464160.tif 000013
【図７】
0005464160.tif 000014
【図８】
0005464160.tif 000015
【図９】
0005464160.tif 000016
【図１０】
0005464160.tif 000017