"2次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析手段と、前記2次元領域を分割する複数の矩形領域を設定する矩形領域設定手段と、前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出手段と、前記複数の矩形領域に対応する複数のセルを設定するセル設定手段と、前記セルの体積をそれぞれ算出するセル体積算出手段と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出手段と、前記セル体積算出手段により算出された、前記セルの体積と、前記粒子体積算出手段により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出手段と、を有し、前記セルは、前記矩形領域の縦の辺である第1の辺と、前記矩形領域の横の辺である第2の辺と、当該第1の辺及び当該第2の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出手段により算出された平均径の長さを有する第3の辺とにより定められる平行六面体であり、前記第3の辺の中央における当該第3の辺に沿う方向の位置が、当該第3の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする粒子存在場解析装置。前記平均粒子径算出手段は、前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を算出することを特徴とする請求項1に記載の粒子存在場解析装置。前記空隙率算出手段により算出された粒子の空隙率のそれぞれに、一定値を加算又は減算する空隙率補正手段を更に有することを特徴とする請求項1又は2に記載の粒子存在場解析装置。前記粒子が存在する3次元領域における流体解析を、前記空隙率を用いて行う流体解析手段を有することを特徴とする請求項1~3の何れか1項に記載の粒子存在場解析装置。前記粒子は、高炉に装入される粒子であることを特徴とする請求項1~4の何れか1項に記載の粒子存在場解析装置。2次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析工程と、前記2次元領域を分割する複数の矩形領域を設定する矩形領域設定工程と、前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出工程と、前記複数の矩形領域に対応する複数のセルを設定するセル設定工程と、前記セルの体積をそれぞれ算出するセル体積算出工程と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出工程と、前記セル体積算出工程により算出された、前記セルの体積と、前記粒子体積算出工程により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出工程と、を有し、前記セルは、前記矩形領域の縦の辺である第1の辺と、前記矩形領域の横の辺である第2の辺と、当該第1の辺及び当該第2の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出工程により算出された平均径の長さを有する第3の辺とにより定められる平行六面体であり、前記第3の辺の中央における当該第3の辺に沿う方向の位置が、当該第3の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする粒子存在場解析方法。前記平均粒子径算出工程は、前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を算出することを特徴とする請求項6に記載の粒子存在場解析方法。前記空隙率算出工程により算出された粒子の空隙率のそれぞれに、一定値を加算又は減算する空隙率補正工程を更に有することを特徴とする請求項6又は7に記載の粒子存在場解析方法。前記粒子が存在する3次元領域における流体解析を、前記空隙率を用いて行う流体解析工程を有することを特徴とする請求項6~8の何れか1項に記載の粒子存在場解析方法。前記粒子は、高炉に装入される粒子であることを特徴とする請求項6~9の何れか1項に記載の粒子存在場解析方法。請求項1~5の何れか1項に記載の粒子存在場解析装置の各手段としてコンピュータを機能させることを特徴とするコンピュータプログラム。非特許文献1に記載の技術では、DEMによる解析を行って粒子の3次元領域での挙動を解析するようにしている。したがって、非特許文献1の技術では、実高炉の1/4のサイズの領域を計算領域とすると共に、実際には数億個存在する粒子の個数を、粒子の粒径を大きくして1000万個の規模にし、且つ、地球シミュレータを用いて計算しても、DEMによる解析に、数か月のオーダの計算時間を要する非特許文献1では、最大計算粒子数を合計で900万個としている。このように、非特許文献1に記載の直接3次元をDEMで扱う技術では、実高炉の内部における粒子の空隙率を算出するのに莫大な時間がかかるという問題点があった。そこで、2次元領域におけるDEMによる解析を行って粒子の配置を求め、この粒子の配置から、2次元計算で空隙率を算出することも考えられる。しかしながら、粒子の空隙率は3次元的な指標であるため、2次元領域における粒子の空隙率の値は、実際の空隙率の値と大きく異なることになる。以上のように、従来の技術では、特に、計算対象となる粒子の数が多い場合には、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することが困難であるという問題点があった。本発明は、以上のような問題点に鑑みてなされたものであり、計算対象となる粒子の数が多くなっても、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することを目的とする。本発明の粒子存在場解析装置は、2次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析手段と、前記2次元領域を分割する複数の矩形領域を設定する矩形領域設定手段と、前記矩形領域設定手段により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出手段と、前記複数の矩形領域に対応する複数のセルを設定するセル設定手段と、前記セルの体積をそれぞれ算出するセル体積算出手段と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出手段と、前記セル体積算出手段により算出された、前記セルの体積と、前記粒子体積算出手段により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出手段と、を有し、前記セルは、前記矩形領域の縦の辺である第1の辺と、前記矩形領域の横の辺である第2の辺と、当該第1の辺及び当該第2の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出手段により算出された平均径の長さを有する第3の辺とにより定められる平行六面体であり、前記第3の辺の中央における当該第3の辺に沿う方向の位置が、当該第3の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする。本発明の粒子存在場解析方法は、2次元領域における粒子の挙動を、当該粒子の形状を球状として、離散要素法により解析する粒子挙動解析工程と、前記2次元領域を分割する複数の矩形領域を設定する矩形領域設定工程と、前記矩形領域設定工程により設定された複数の矩形領域のそれぞれについて、当該矩形領域の中に全ての領域が入っている球状の粒子と、当該矩形領域の中に一部の領域が入っている球状の粒子であって、中心の位置が当該矩形領域の中に含まれている球状の粒子と、の平均径を算出する平均粒子径算出工程と、前記複数の矩形領域に対応する複数のセルを設定するセル設定工程と、前記セルの体積をそれぞれ算出するセル体積算出工程と、前記セルの中に全ての領域が含まれている粒子の体積と、前記セルの中に一部の領域が含まれている粒子の領域のうち当該セルの中に含まれている領域の体積との総和を算出する粒子体積算出工程と、前記セル体積算出工程により算出された、前記セルの体積と、前記粒子体積算出工程により算出された、前記体積の総和とに基づいて、粒子の空隙率を、前記セルのそれぞれについて算出する空隙率算出工程と、を有し、前記セルは、前記矩形領域の縦の辺である第1の辺と、前記矩形領域の横の辺である第2の辺と、当該第1の辺及び当該第2の辺に垂直な辺であって、当該矩形領域に対して前記平均粒子径算出工程により算出された平均径の長さを有する第3の辺とにより定められる平行六面体であり、前記第3の辺の中央における当該第3の辺に沿う方向の位置が、当該第3の辺に沿う方向の位置が夫々揃っている前記粒子の中心位置と揃う位置に設定されることを特徴とする。本発明のコンピュータプログラムは、前記粒子存在場解析装置の各手段としてコンピュータを機能させることを特徴とする。本発明によれば、2次元領域における粒子の挙動を離散要素法により解析する。そして、当該2次元領域を分割する複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の平均径を個別に算出する。そして、当該複数の矩形領域に対応する複数のセルを設定する。このセルは、当該矩形領域の縦横の辺である第1第2の辺と、当該第1第2の辺に垂直な辺であって、当該矩形領域に対して算出された平均径の長さを有する第3の辺とにより定められる平行六面体であり、第3の辺の中央における第3の辺に沿う方向の位置が、第3の辺に沿う方向の位置が夫々揃っている粒子の中心位置と揃う位置に設定される。そして、これらのセル毎に、粒子の空隙率を算出する。したがって、離散要素法により、3次元領域における粒子の挙動を解析しなくても、当該3次元領域における粒子の挙動を計算することにより得られる粒子の空隙率に近い値の空隙率を算出することができる。よって、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することができる。以下、図面を参照しながら、本発明の一実施形態を説明する。図1は、粒子存在場解析装置100の機能的な構成の一例を示す図である。粒子存在場解析装置100は、例えば、高炉のベルレス式炉頂装入装置用の旋回シュートから高炉内に装入される粒子焼結鉱やコークス等の挙動の解析を行い、その結果を利用して、高炉内におけるガスの流れの解析流体解析を行うものである。ただし、粒子存在場解析装置100が解析を行う対象は、このようなものに限定されるものではない。粒子存在場解析装置100は、例えば、CPU、ROM、RAM、HDD、各種インターフェース、及びディスプレイ等を備えた情報処理装置パーソナルコンピュータを用いることにより実現できる。粒子挙動解析部101粒子挙動解析部101は、球状の粒子の2次元領域における挙動をDEMにより解析する。粒子挙動解析部101は、例えば、ユーザによるユーザインターフェースの操作に基づいて設定された、粒子の情報に基づいて、各粒子を所望の位置に配置する。そして、粒子挙動解析部101は、粒子間に発生している力のバランスを崩す条件を与えて、各粒子を動かし、各粒子が他の粒子や構造物と接触しているか否かを判定する。粒子が他の粒子と接触している場合、粒子挙動解析部101は、法線方向及びせん断方向の反発力を求めて並進の運動方程式を解くことと、重心周りの力のモーメントを求めて回転の運動方程式を解くこと、を各粒子のそれぞれについて行い、その結果に基づいて各粒子を移動させる。一方、粒子が構造物と接触している場合、粒子挙動解析部101は、法線方向及びせん断方向の反発力を求めて並進の運動方程式を解くことと、重心周りの力のモーメントを求めて回転の運動方程式を解くこと、を粒子についてのみ行い、構造物については行わない。粒子挙動解析部101は、このような接触判定処理と各粒子の移動処理とを所定の時間隔毎に行う。そして、その結果から各粒子の所定の期間における挙動を得る。具体的に、粒子挙動解析部101は、例えば、以下の1式により、粒子i、jの接触領域における法線方向の反発力Fnijを求め、以下の2式により、粒子i、jの接触領域におけるせん断方向の反発力Ftijを求め、以下の3式の並進の運動方程式を解き、以下の4式により回転の運動方程式を解く。図2は、相互に接触している2つの粒子201、202の様子の一例図2aと、粒子203の挙動の一例図2bを概念的に示す図である。図2を参照しながら、1式から4式のパラメータを説明する。1式において、Knは、法線方向のばね定数であり、例えば図2aの距離Hやδnに応じて設定することができる。ただし、Knは、予め設定された一定値にすることもできる。Δunijは、2つの粒子i、j間の重心の法線方向における相対的な並進変位である。ηnは、法線方向の粘性定数であり、予め設定されるものである。nijは、法線方向の単位ベクトルである。Δtは、所定の時間隔であり、予め設定されるものである。2式において、minA,Bは、AとBのうち小さい方を採用することを意味する。μは、摩擦係数であり、予め設定された値にすることができる。ただし、μは、状態等に応じて異ならせることもできる。tijはせん断方向の単位ベクトルである。Ktは、せん断方向のばね定数であり、例えば図2aの距離Hやδn、δtに応じて設定することができる。ただし、Ktは、予め設定された一定値にすることもできる。Δutijは、2つの粒子i、j間の重心のせん断方向における相対的な並進変位である。Δφijは、粒子の回転に起因する接触領域の相対的な回転変位であり、例えば、図2aの粒子201の回転変位φiと粒子202の回転変位φjとの差により表すことができる。ηtは、せん断方向の粘性定数であり、予め設定されるものである。Δtは、所定の時間隔であり、予め設定されるものである。3式において、vは、粒子の速度である図2を参照。vの上に付されているドットは時間微分を表す。Fは、せん断方向及び法線方向の反発力の総和である。mは、粒子の質量であり、予め設定されるものである。gは、重力加速度である。尚、粒子が他の粒子や構造物と接触していない場合、3式の右辺の第1項は0ゼロになる。4式において、ωは、粒子の角速度である図2を参照。ωの上に付されているは時間微分を表す。Mは、力のモーメントの総和であり、せん断方向及び法線方向の反発力Fnij、Ftijから求められる。Iは、慣性モーメントであり、予め設定されるものである。本実施形態では、2次元領域における粒子の挙動をDEMにより解析するに際し、当該粒子を球状とする。3次元の成分を有する指標粒子の物性値等を考慮して計算するためである。すなわち、2次元領域における粒子の挙動をDEMにより解析するに際し、粒子に働く作用力の絶対値については3次元の情報として計算し、粒子に働く作用力の方向については2次元の情報として計算する。また、高炉内に装入される粒子の2次元領域の挙動を解析する場合、高炉の半径方向と高さ方向とにより定まる2次元領域の挙動を解析することになる。高炉の奥行き方向への粒子の挙動は、解析の結果への影響が小さいためである。したがって、高炉内に装入される粒子の2次元領域の挙動を解析する場合、粒子に働く作用力の方向は、高炉を、その軸に沿って切った断面の面内方向となる。以上のようにして、計算対象となる2次元領域における各粒子の配置が、所定の時間毎に得られる。尚、粒子の挙動を解析する手法粒子挙動解析部101が実行する解析手法は、DEMにおける粒子の挙動についての種々の公知の解析手法を適用することができ、前述した1式~4式を用いた手法に限定されるものではない。粒子挙動解析部101は、例えば、CPUが、RAM等から、粒子の情報を読み出して前述した演算処理を実行して、各時間における各粒子の2次元領域における座標を得てRAM等に記憶することにより実現される。矩形領域設定部102矩形領域設定部102は、粒子挙動解析部101により粒子の挙動が解析された2次元領域を分割する複数の矩形領域所謂メッシュを設定する。前述したように、粒子挙動解析部101によって、2次元領域における所定時間毎の各粒子の配置が粒子の挙動として得られる。矩形領域設定部102は、この各粒子が配置されている2次元領域に対して、複数の矩形領域を設定する。高炉内に装入される粒子の2次元領域の挙動を解析する場合、矩形領域設定部102は、例えば、各粒子の動きが止まった時点において、各粒子が配置されている2次元領域の情報を取得し、当該2次元領域に対して、複数の矩形領域を設定する。ただし、必ずしも各粒子の動きが止まった時点における2次元領域に対して複数の矩形領域を設定する必要はなく、任意の時点における2次元領域に対して複数の矩形領域を設定することができる。図3は、各粒子が配置されている2次元領域に対して設定される矩形領域の一例を示す図である。尚、各粒子が配置されている2次元領域は、例えば、後述する図6aに示すような領域である。図3に示すように、本実施形態では、各粒子が配置されている2次元領域に、x軸方向の長さがdxであり、z軸方向の長さがdzである矩形領域が、x軸方向及びz軸方向のそれぞれに連続的に設定される。図3に示すように、1つの矩形領域の中に全ての領域が含まれている粒子もあれば、複数の矩形領域に跨って存在する粒子もある。尚、この2次元領域が高炉内の領域である場合、x軸方向は、高炉の半径方向であり、z軸方向は、高炉の高さ方向となる。本実施形態では、全ての矩形領域の大きさは同じであるとする。また、1つの矩形領域に10個以上の粒子が入るように、矩形領域の大きさを決めるようにするのが好ましい。さらに、粒子挙動解析部101により解析が行われる2次元領域の座標が分かれば矩形領域を設定することができるので、粒子挙動解析部101による解析が行われる前に矩形領域が設定されるようにしてもよい。矩形領域設定部102は、例えば、CPUが、各矩形領域の2次元領域x-z平面における座標x軸座標、z軸座標を得てRAM等に記憶することにより実現される。平均粒子径算出部103平均粒子径算出部103は、矩形領域設定部102により設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の平均径を算出する。本実施形態では、最も好ましい平均径として、長さ平均径を採用する。この理由については後述する。ここで、1つの矩形領域の中に全ての領域が含まれていない粒子については、当該粒子の中心が含まれている矩形領域に属しているものとして長さ平均径を算出する。平均粒子径算出部103は、例えば、CPUが、矩形領域に含まれる球状の粒子の長さ平均径を、全ての矩形領域について得てRAM等に記憶することにより実現される。セル設定部104セル設定部104は、矩形領域の縦の辺であるz軸方向の辺図3の長さdzの辺と、矩形領域の横の辺であるx軸方向の辺図3の長さdxの辺と、当該z軸方向の辺及び当該x軸方向の辺に垂直なy軸方向の辺であって、当該矩形領域に対して平均粒子径算出部103により算出された長さ平均径の長さを有するy軸方向の辺と、により定められる平行六面体であって、当該y軸方向の辺の中央におけるy軸方向の位置が、球場の粒子のy軸方向の中心位置とそれぞれ揃えられた複数の平行六面体を設定する。以下の説明では、この平行六面体を必要に応じてセルと称する。図4は、y軸方向及びz軸方向から見たセルの一例を示す図である。尚、図4では、表記の都合上、各セル401a~401cに含まれている粒子の数を1~2個としているが、実際には、各セル401a~401cには、これよりも多くの数の粒子が含まれる。セル設定部104は、矩形領域設定部102により設定された矩形領域のz軸方向の辺図3の長さdzの辺と、x軸方向の辺図3の長さdxの辺の情報を取得すると共に、当該矩形領域に対して平均粒子径算出部103により算出された長さ平均径の情報を取得する。そして、セル設定部104は、当該矩形領域のz軸方向の辺と、当該矩形領域のx軸方向の辺と、当該矩形領域の頂点を中央点とし、y軸方向に当該矩形領域に対して設定された長さ平均径の長さ分だけ左右y軸正負の両方向に均等に伸ばしたy軸方向の辺とにより定まる平行六面体直方体又は立方体を、矩形領域設定部102により設定された全ての矩形領域について定める。図4に示す例では、当該矩形領域の頂点を中央点とし、y軸方向に当該矩形領域に対して設定された長さ平均径の長さ分だけ左右に伸ばしたy軸方向の辺は、長さがdy1、dy2、dy3の辺である。このようにして得られた平行六面体が、図4に示すセル401a~401cになり、各セル401a~401cのy軸方向の中央の位置は、x軸及びz軸に平行な面402に位置する。各矩形領域に含まれる粒子の大きさ及び数は矩形領域によって異なるので、各矩形領域に含まれる粒子の長さ平均径も矩形領域によって異なる。よって、図4に示すように、各セル401a~401cのy軸方向の長さdy1~dy3もセル401a~401cによって異なることになるただし、同じになる場合もある。尚、図4では、セル401の数を3つしか示していないが、セル設定部104により設定されるセル401の数は、矩形領域設定部102により設定された矩形領域の数と同じ数になる。セル設定部104は、例えば、CPUが、セル401の3次元座標x軸座標、y軸座標、z軸座標を得てRAM等に記憶することにより実現される。セル体積算出部105セル体積算出部105は、セル設定部104によって設定されたセル401の体積をそれぞれ算出する。前述したように、セルの形状は平行六面体直方体又は立方体であるので、各辺の長さdx、dy、dzを掛け合わせることにより、各セル401の体積を得ることができる。すなわち、各セル401の体積Vcは、以下の5式により得られる。セル体積算出部105は、例えば、CPUが、各セル401の体積を得てRAM等に記憶することにより実現される。粒子体積算出部106粒子体積算出部106は、セル401の中に含まれる粒子の体積を算出する。粒子体積算出部106は、図4に示すように、各セル401a~401cのy軸方向の中央の位置座標と、各粒子411~415の中心411a~415aのy軸方向の位置とが、それぞれ、x軸及びz軸に平行な面402に位置している。そして、粒子体積算出部106は、各セル401a~401cの中に含まれている粒子411~415の体積を算出する。各粒子411~415の形状は球である。そこで、粒子体積算出部106は、セル401の中に全ての領域が含まれている粒子については、当該セル401の中に含まれている粒子の体積Vpを、以下の6式で算出する。6式において、Dpは粒子球の直径である。図4に示す例では、粒子413~415が、このような粒子に該当する。また、粒子体積算出部106は、1つのセル401の中に全ての領域が含まれていない粒子については、当該粒子の領域のうち、当該セル401の中に含まれている領域の体積を、当該セル401の中に含まれている粒子の体積Vpとして算出する。図5は、1つのセル401の中に全ての領域が含まれていない粒子の体積の計算式のパラメータを示す図である。図5において、Dpは、粒子416の直径である。aは、粒子416の中心416aと、粒子416が跨いでいるセル401dのx-z平面上の辺との間の最短距離である。これらDp、aを用いると、粒子体積算出部106は、セル401の中に一部の領域が含まれている粒子416については、当該粒子416の領域のうち、セル401dの中に含まれている領域の体積を、当該セル401dの中に含まれている粒子の体積Vpとして、以下の7式により算出する。また、粒子体積算出部106は、セル401の中に一部の領域が含まれている粒子416については、当該粒子416の領域のうち、セル401dに含まれていない領域の体積を、当該領域が属しているセル401eの中に含まれている粒子の体積Vpとして、以下の8式により算出する。図4に示す例では、粒子411、412が、このような粒子に該当する。尚、1つのセル401の中に全ての領域が含まれていない粒子の領域のうち、当該セル401の中に含まれている領域の体積を算出する方法計算式は、前述した方法7式、8式に限定されず、粒子とセル401との位置関係に応じて適宜設定することができる。そして、粒子体積算出部106は、1つのセル401の中に少なくとも一部の領域が含まれている粒子の全てについて7式又は8式により計算したセル401の中に含まれている粒子の体積Vpの総和を計算する。粒子体積算出部106は、例えば、CPUが、各セル401の中に含まれている粒子の体積Vpを得て、RAM等に記憶することにより実現される。空隙率算出部107空隙率算出部107は、セル体積算出部105で算出された各セル401の体積Vcと、粒子体積算出部106により算出された各セル401の中に含まれている粒子の体積Vpの総和とに基づいて、粒子の空隙率を、セル設定部104により設定された複数のセル401のそれぞれについて算出する。空隙率算出部107は、以下の9式により、各矩形領域における粒子の仮空隙率ε ́[-]を個別に算出する。9式において、nは、粒子体積算出部106により、同一のセル401の中に含まれているとして計算された粒子の体積Vpの数である。例えば、図4に示す例では、セル401aの中に含まれているとして計算された粒子の体積Vpの数は2であり、セル401bの中に含まれているとして計算された粒子の体積Vpの数は3であり、セル401cの中に含まれているとして計算された粒子の体積Vpの数は1である。そして、空隙率算出部107は、算出した各矩形領域における粒子の仮空隙率ε ́のそれぞれに定数=0.02を加算した値を、各矩形領域における粒子の空隙率εとして算出する。定数=0.02を加算して仮空隙率ε ́を補正する理由については、後述する。空隙率算出部107は、例えば、CPUが、各セル401における粒子の空隙率εを得て、RAMやHDD等に記憶することにより実現される。流体解析部108流体解析部108は、空隙率算出部107により算出された各セル401における粒子の空隙率εを用いて、粒子が配置された3次元領域におけるガスの流れの解析流体解析を行い、その結果をディスプレイに表示する。流体解析の手法は、例えば、非特許文献1に記載されているので、ここでは、その詳細な説明を省略する。流体解析部108は、例えば、CPUが、粒子が配置された3次元領域におけるガスの流れの解析結果を得て、そのデータをディスプレイに出力することにより実現される。長さ平均径を採用する理由と空隙率を調整する理由本実施形態では、以下に示す種々の粒度構成を有する条件で、2次元領域2次元の計算場における各粒子の挙動をDEMにより計算し、その計算結果から、各種の平均径を用いて前述したようにして空隙率の値を求めると共に、3次元領域3次元の計算場における各粒子の挙動をDEMにより計算し、それぞれの計算結果から空隙率の値を求めた。そして、これらの空隙率の値を比較した結果、一定値を用いていた従来法に比較して、平均粒子径でセル体積を設定する方法の方が、種々の粒度構成を有する条件でも3次元値に近い値を示し、かつ試みたいずれの平均径の場合も高炉内ガス流れ解析での使用に際して十分な精度を有していることが分かった。さらに、各種の平均径のうち、長さ平均径を採用することが好ましいということと、長さ平均径を用いて求めた空隙率の値を調整すると、3次元領域3次元の計算場で求めた空隙率の値に非常に近い空隙率が得られることを見出した。以下に、これらについて詳述する。図6は、DEMによる解析を行って計算対象の領域に配置された各粒子の一例を概念的に示す図である。図6aは、2次元領域2次元の計算場における各粒子の挙動を解析して得られた各粒子の配置を示す図である。図6bは、図6a本実施形態と比較するための図であり、3次元領域3次元の計算場における各粒子の挙動を解析して得られた各粒子の配置を示す図である。図6aに示すように、DEMによって、2次元領域における各粒子の挙動を解析すると、2次元領域x-z平面の領域と直交する方向y軸方向に積み重ならないように、つまり球場の粒子のy軸方向の中心位置が全て揃った球状の粒子が、x-z平面に並べられる。これに対し、図6bに示すように、DEMによって、3次元領域における各粒子の挙動を解析すると、球状の粒子は、3次元領域x-y-z空間に並べられるので、y軸方向にも積み重なる。図7は、図6に示す各区間A、B、Cにおける各粒子の粒度構成を示す図である。具体的に、図7aは、図6aに示した各区間A、B、Cにおける各粒子の面積比と、図6bに示した区間A、B、Cにおける各粒子の体積比とを示す。一方、図7bは、図6a及び図6bに示した各区間A、B、Cにおける粒子個数比を示す。図6は、各区間A、B、Cのそれぞれにおいて、それぞれが接触しないように粒子を初期配置し、上方から図3に示す矢印の方向に各粒子を重力沈降させ、動きが止まったときの各粒子の配置を示す。尚、図6では、表記の都合上空隙があることを示すために、粒子と粒子との間に隙間があるが、実際には、図6に示す状態よりも各粒子は密に配置されている。図7に示すように、各区間A、B、Cに落下させる粒子は、直径粒子径が20.0[mm]、37.5[mm]、55.0[mm]の3種類の粒子である。ここで、図6aに示す区間A、B、Cは、2次元領域であるのに対し、図6bに示す区間A、B、Cは、3次元領域である。そこで、区間A、B、Cのそれぞれにおいて、これら3種類の粒子の割合が、図6aに示す区間A、B、Cと、図6bに示す区間A、B、Cとで同じになるようにしている。具体的に、図7aに示すように、各粒子の図6aに示す各区間A、B、Cにおける面積比と、各粒子の図6bに示す各区間A、B、Cにおける体積比とが同じになるようにしている。例えば、図7aに示すように、図6aに示す区間Aにおける粒子径が20.0[mm]の粒子の面積の総和と、図6aに示す区間Bにおける粒子径が20.0[mm]の粒子の面積の総和と、図6aに示す区間Cにおける粒子径が20.0[mm]の粒子の面積の総和との比が1:1:5になるようにしている。そこで、図7bに示すように、図6bに示す区間Aにおける粒子径が20.0[mm]の粒子の体積の総和と、図6bに示す区間Bにおける粒子径が20.0[mm]の粒子の体積の総和と、図6bに示す区間Cにおける粒子径が20.0[mm]の粒子の体積の総和との比が1:1:5となるようにしている。また、各粒子の図6bに示す各区間A、B、Cにおける個数比粒子個数比と、各粒子の図6bに示す各区間A、B、Cにおける個数比粒子個数比は、図7に示す通りである。例えば、図7aに示すように、図6aに示す区間Aにおける粒子径が20.0[mm]の粒子の数と、図6aに示す区間Bにおける粒子径が20.0[mm]の粒子の数と、図6aに示す区間Cにおける粒子径が20.0[mm]の粒子の数との比は、8:8:40=1:1:5となる。一方、図6bに示す区間Aにおける粒子径が20.0[mm]の粒子の数と、図6bに示す区間Bにおける粒子径が20.0[mm]の粒子の数と、図6bに示す区間Cにおける粒子径が20.0[mm]の粒子の数との比も、21:21:105=1:1:5となる。図8は、矩形領域に含まれる粒子の平均径の計算方法を異ならせて算出した空隙率を示す図である。図8に示す各点は、図3に示した区間A、B、Cに、図4に示した割合で粒子を重力沈降させて、粒子の動きが止まった時点での各矩形領域の空隙率の値を示す。また、各区間A、B、Cに対して、それぞれ6つの矩形領域を設定したこのことは、図8の各区間A、B、Cに6つの点があることと対応する。また、図8において、3Dは、DEMにより、3次元領域における粒子の挙動を解析し、その結果から算出した空隙率の値を示す図3bを参照。この値が最も信頼できる値となる。また、平均径定数固定とは、セルのy軸方向の長さを矩形領域によらずに37.5[mm]に固定して算出した空隙率の値を示すセルのy軸方向の長さを固定する以外は、本実施形態と同じ処理を行って得られた空隙率の値であるただし、空隙率は9式で計算されるものであり、定数を加算していない。また、個数平均径、面積平均径、体積平均径、平均表面積径、平均体積径は、長さ平均径の代わりに、これらの平均径を使って算出した空隙率の値を示す平均径の計算方法以外は、本実施形態と同じ処理を行って得られた空隙率の値であるただし、空隙率は9式で計算されるものであり、定数を加算していない。また、長さ平均径は、定数を加算する前の空隙率9式の仮空隙率ε ́である。また、長さ平均径0.02は、本実施形態と同じ処理を行って得られた空隙率である。尚、個数平均径、長さ平均径、面積平均径、体積平均径、平均表面積径、平均体積径の計算方法については、非特許文献2に記載されている。また、個数平均径、長さ平均径、面積平均径、体積平均径については、それぞれ重みつき個数平均径、重みつき長さ平均径、重みつき面積平均径、重みつき体積平均径と称されることもある。図8に示すように、程度の差はあるが、個数平均径、長さ平均径、面積平均径、体積平均径、平均表面積径、平均体積径のグラフの方が、平均径定数固定のグラフよりも、3Dのグラフと、概形が似ている。区間A、B、Cでは、3種類の粒子径の構成割合を極端に異ならせている図4を参照。よって、区間A、B、Cの全てで、3Dのグラフに対し、形と値とが近くなっていれば、粒子径がどのような条件であっても、精度よく空隙率を算出することができると言ってよい。したがって、粒子の平均径を矩形領域毎に算出し、各セルのy軸方向の長さを、算出した平均径に基づいて決定すると、精度よく空隙率を算出することができることが分かる。特に、長さ平均径のグラフは、3Dのグラフと、概形だけでなく、値も近くなっており、特に好ましいことも判る。そこで、本実施形態では、これらの平均径のうち、長さ平均径を採用することとした。また、長さ平均径を用いて算出した空隙率に一定値を加算した空隙率の値と図8の長さ平均径、DEMにより3次元領域における粒子の挙動を解析した結果から算出した空隙率の値図8の3Dとの差を最小にするためには、当該一定値をいくつにすればよいのかを、最小二乗法により計算した。その結果、当該一定値は0.02となった。よって、この例では、長さ平均径を用いて算出した空隙率の値に、当該一定値として0.02を加算すれば、3Dのグラフの値との差が非常に小さい空隙率を得ることができる。そこで、本実施形態では、長さ平均径を用いて算出した空隙率仮空隙率ε ́に一定値=0.02を加算した値を各矩形領域の空隙率εとして採用することとした。ただし、当該一定値を加算又は減算する前の空隙率の値よりも、当該一定値を加算又は減算した後の空隙率の値の方が、DEMにより3次元領域における粒子の挙動を解析した結果から算出した空隙率の値との差が小さくなるようにしていれば、必ずしも、3Dのグラフの値との差が最小となるように当該一定値を計算する必要はない。また、3Dのグラフの値との差が最小となるように当該一定値を計算するために用いる方法は、最小二乗法に限定されるものではない。尚、長さ平均径以外の平均径を用いて算出した空隙率でも、当該空隙率に対し一定値を加算又は減算することにより、当該一定値を加算又は減算した空隙率の値と、3Dのグラフの値との差を小さくすることができるので、このようにしてもよい。特に、3Dのグラフの値との差が最小となるように、当該空隙率の値に対して加算又は減算する一定値を決定するのが好ましい。3Dのグラフの値との差が最小となる一定値を得るための計算は、前述したように、例えば、最小二乗法を用いることにより実現することができる。動作フローチャート次に、図9のフローチャートを参照しながら、粒子存在場解析装置100の処理の一例を説明する。まず、ステップS901において、粒子挙動解析部101は、球状の粒子の2次元領域における挙動をDEMにより解析する。これにより、所定の時間毎に、2次元領域における粒子の配置が得られる図3aを参照。次に、ステップS902において、矩形領域設定部102は、ステップS901で粒子の挙動が解析された2次元領域を分割する複数の矩形領域を設定する図3を参照。次に、ステップS903において、平均粒子径算出部103は、ステップS902で設定された複数の矩形領域のそれぞれについて、当該矩形領域に含まれる球状の粒子の長さ平均径を個別に算出する。次に、ステップS904において、セル設定部104は、ステップS902で設定された矩形領域の各辺をx軸、z軸方向の辺とし、これらx軸、z軸方向の辺と、ステップS903で算出された当該矩形領域に含まれる球状の粒子の長さ平均径の長さを有するy軸方向の辺とにより定まる平行六面体であって、当該y軸方向の辺の中央の位置が夫々揃えられた複数の平行六面体を、セル401として設定する図4を参照。次に、ステップS905において、セル体積算出部105は、ステップS904で設定された各セルの体積を算出する5式を参照。次に、ステップS906において、粒子体積算出部106は、セル内粒子体積算出処理を実行して、セル401の中に含まれている粒子の体積Vpの総和を算出する6式~8式を参照。尚、セル内粒子体積算出処理の詳細については、図10を参照しながら後述する。次に、ステップS907において、空隙率算出部107は、ステップS905で算出された各セル401の体積Vcと、ステップS906で算出された各セル401の中に含まれている粒子の体積Vpとに基づいて、ステップS904設定された複数のセル401のそれぞれにおける粒子の空隙率εを算出する粒子の空隙率εは9式の仮空隙率ε ́に定数を加算した値となる。次に、ステップS908において、流体解析部108は、ステップS907で算出された複数のセル401のそれぞれにおける粒子の空隙率εを用いて、粒子が配置された3次元領域におけるガスの流れの解析流体解析を行い、その結果をディスプレイに表示する。そして、図9のフローチャートによる処理を終了する。次に、図10のフローチャートを参照しながら、図9のステップS906におけるセル内粒子体積算出処理の詳細な説明を行う。尚、図10のフローチャートによる処理は、ステップS904により設定された全てのセルについて行われる。まず、ステップS1001において、粒子体積算出部106は、計算対象となるセル401の中に少なくとも一部の領域が含まれている粒子のうち、未選択の粒子を1つ選択する。次に、ステップS1002において、粒子体積算出部106は、ステップS1001で選択した粒子の全ての領域が、計算対象となるセル401の中に含まれているか否かを判定する。この判定の結果、ステップS1001で選択した粒子の全ての領域が、計算対象となるセル401の中に含まれている場合には、ステップS1003に進む。ステップS1003に進むと、粒子体積算出部106は、ステップS1001で選択した粒子の全体の体積を、計算対象となるセル401の中に含まれている粒子の体積Vpとして、6式により算出する。一方、ステップS1001で選択した粒子の全ての領域が、計算対象となるセル401の中に含まれていない場合には、ステップS1004に進む。ステップS1004に進むと、粒子体積算出部106は、ステップS1001で選択した粒子の領域のうち、当該セル401に含まれている領域の体積を、計算対象となるセル401の中に含まれている粒子の体積Vpとして、例えば7式又は8式により計算する。以上のようにして、ステップS1001で選択した粒子について、計算対象となるセル401の中に含まれている粒子の体積Vpが算出されると、ステップS1005に進む。ステップS1005に進むと、粒子体積算出部106は、計算対象となるセル401の中に少なくとも一部の領域が含まれている粒子を全て選択したか否かを判定する。この判定の結果、計算対象となるセル401の中に少なくとも一部の領域が含まれている粒子を全て選択していない場合には、ステップS1001の処理に戻る。そして、計算対象となるセル401の中に少なくとも一部の領域が含まれている粒子の全てについて、計算対象となるセル401の中に含まれている粒子の体積Vpを算出するまで、ステップS1001~S1005の処理を繰り返し行う。その後、粒子体積算出部106は、計算対象となるセル401の中に含まれている粒子の体積Vpの総和を計算する。まとめ以上のように本実施形態では、球状の粒子の2次元領域における挙動をDEMにより解析した後、当該2次元領域を分割する複数の矩形領域を設定する。複数の矩形領域のそれぞれについて、当該矩形領域に含まれる粒子の長さ平均径を個別に算出する。そして、そして、当該複数の矩形領域に対応する複数のセル401を設定する。このセル401は、矩形領域のx軸、z軸方向の辺と、当該矩形領域における長さ平均径の長さを有するy軸方向の辺と、により定まる平行六面体であり、当該y軸方向の辺の中央の位置が、当該y軸方向の位置が夫々揃っている球状の粒子の中心位置と揃う位置に設定される。そして、各セル401に含まれる粒子の割合を求めて、各矩形領域における粒子の空隙率εを算出する。したがって、DEMによって、3次元領域における粒子の挙動を計算しなくても、当該3次元領域における粒子の挙動を計算することにより得られる粒子の空隙率εに近い値の空隙率εを算出することができる。例えば、高炉の炉頂から羽口の間の荷下がりを、3次元計算場でDEMにより計算すると、計算時間が約2ヶ月であったのに対し、2次元計算場でDEMにより計算すると、計算時間が約5日であった。また、これらにおける空隙率εの差は、概ね図8に示した結果と同じになった。よって、計算対象となる粒子の数が多くなっても、粒子の空隙率を、その算出精度を大きく落とすことなく高速に算出することができる。<変形例>本実施形態では、流体解析を行うようにしたが、粒子の空隙率εをセル401毎に算出するようにしていれば、必ずしも流体解析を行う必要はない。また、粒子の空隙率εをディスプレイに表示してユーザに提示するようにしてもよい。尚、以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、プログラムをコンピュータに供給するための手段、例えばかかるプログラムを記録したCD-ROM等のコンピュータ読み取り可能な記録媒体、又はかかるプログラムを伝送する伝送媒体も本発明の実施の形態として適用することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体などのプログラムプロダクトも本発明の実施の形態として適用することができる。前記のプログラム、コンピュータ読み取り可能な記録媒体、伝送媒体及びプログラムプロダクトは、本発明の範疇に含まれる。また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。"