(57)【要約】【課題】気流中に粒子が混在する固気二相流を対象として，解析領域を計算格子（グリッド）に分割することなく精度よく解析できる，グリッドフリー手法を用いた固気二相流シミュレーションプログラム及びそれを記憶した記憶媒体並びに固気二相流シミュレーション装置を提供する．【解決手段】本発明は，微小な固体粒子が気流中に分散し，気流と相互作用を及ぼし合いながら流れる固気二相流の渦度場に着目して複数の渦要素で離散化し，気流に及ぼす粒子の力を考慮した運動量保存関係から渦度方程式を導出し，この渦度方程式のラグランジュ計算により渦要素の挙動を求めるものであり，解析領域を計算格子（グリッド）に分割する必要がない，グリッドフリー手法による解析方法である．【選択図】  図９
2009193110.tif 000002
【特許請求の範囲】
【請求項１】
気流中に微小な固体粒子を含む固気二相流を解析するための，グリッドフリー手法を用いた固気二相流シミュレーションプログラムであって，コンピュータに，次のステップ（ａ）〜（ｉ）を実行させて，解析領域における粒子の位置と速度，および前記解析領域内に予め設定した評価点における気流の速度を解析する固気二相流シミュレーションプログラム．
（ａ）前記解析領域内に存在する各粒子の速度をラグランジュ的に計算するステップ．
（ｂ）前記各粒子の位置をラグランジュ的に計算するステップ．
（ｃ）前記各粒子の周りに４つの正方形のセルを配置するステップ．
（ｄ）前記各粒子に作用する流体抗力を，前記ステップ（ａ）で算出された速度に基づいて計算するステップ．
（ｅ）前記ステップ（ｄ）で算出された流体抗力に基づいて，前記各セルにおける循環変化量を，気流の運動量保存方程式から導出される渦度方程式と，レイノルズの輸送定理とに基づいて計算するステップ．
（ｆ）前記ステップ（ｅ）で算出された前記各セルにおける循環変化量に基づいて，当該セル内の各渦要素の循環を計算し，各渦要素の位置と循環を記憶手段（１）に保存するステップ．
（ｇ）各渦要素のコア半径をラグランジュ的に計算して記憶手段（２）に保存するステップ．
（ｈ）前記ステップ（ｆ）および前記ステップ（ｇ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，各渦要素の移流すなわち位置をラグランジュ的に計算し，各渦要素の位置，循環およびコア半径を記憶手段（３）に保存するステップ．
（ｉ）前記ステップ（ｈ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，前記評価点における気体の速度をビオサバールの式を用いて計算し，記憶手段（４）に保存するステップ．
【請求項２】
請求項１に記載のグリッドフリー手法を用いた固気二相流シミュレーションプログラムを記憶した記憶媒体．
【請求項３】
気流中に微小な固体粒子を含む固気二相流を解析するための固気二相流シミュレーション装置であって，次の手段（ａ）〜（ｉ）を備え，これらの手段を実行させて，解析領域における粒子の位置と速度，および前記解析領域内に予め設定した評価点における気流の速度を解析する固気二相流シミュレーション装置．
（ａ）前記解析領域内に存在する各粒子の速度をラグランジュ的に計算する手段．
（ｂ）前記各粒子の位置をラグランジュ的に計算する手段．
（ｃ）前記各粒子の周りに４つの正方形のセルを配置する手段．
（ｄ）前記各粒子に作用する流体抗力を，前記手段（ａ）で算出された速度に基づいて計算する手段．
（ｅ）前記手段（ｄ）で算出された流体抗力に基づいて，前記各セルにおける循環変化量を，気流の運動量保存方程式から導出される渦度方程式と，レイノルズの輸送定理とに基づいて計算する手段．
（ｆ）前記手段（ｅ）で算出された前記各セルにおける循環変化量に基づいて，当該セル内の各渦要素の循環を計算し，各渦要素の位置と循環を記憶手段（１）に保存する手段．
（ｇ）各渦要素のコア半径をラグランジュ的に計算して記憶手段（２）に保存する手段．
（ｈ）前記手段（ｆ）および前記手段（ｇ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，各渦要素の移流すなわち位置をラグランジュ的に計算し，各渦要素の位置，循環およびコア半径を記憶手段（３）に保存する手段．
（ｉ）前記手段（ｈ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，前記評価点における気体の速度をビオサバールの式を用いて計算し，記憶手段（４）に保存する手段．
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は，微小な固体粒子が気流中に分散し，気流と相互作用を及ぼし合いながら流れる固気二相流の挙動の解析方法に関する，グリッドフリー手法を用いた固気二相流シミュレーションプログラム及びそれを記憶した記憶媒体並びに固気二相流シミュレーション装置に関する．
【背景技術】
【０００２】
微小な固体粒子が気流中に分散し，気流と相互作用を及ぼし合いながら流れる固気二相流は，微粉炭燃焼機器，固体ロケットエンジン，サンドブラスト装置など，様々な工業装置において観察される．このため，この固気二相流の挙動について，合理的な数値解析方法の確立が望まれている．
【０００３】
固気二相流の数値解析には，従来，差分法や有限体積法などのオイラー型解法が多用されている．これらの解析では，その前処理として，解析領域を多数の格子（グリッド）に分割する必要がある．格子は解析領域の全域をメッシュ状に覆うものであり，計算格子とよばれている．各格子の形状は，一般に四角形である．各格子の頂点で流れの変数（速度や圧力など）が定義され，上記解法により計算される．図１（ａ）は，液体を満たした円筒容器内に円柱が設置された場合に対する，円柱周りの液体流れ場（解析領域）を示す．図１（ｂ）は，格子分割の一例を示す．解析精度は格子分割の状態に依存するため，高い精度を確保するには，流れが急峻に変化するものと予測される領域は多数の細かな格子に分割する必要がある．また，格子の間には重なりや隙間があってはならないという幾何学的条件も満足しなければならない．したがって，複雑な幾何形状をもつ流れ場を扱う工学上の実務的シミュレーションでは，格子分割に細心の注意と多大な作業が要求され，解析を複雑化する原因となっていた．
【０００４】
また，前記の従来の解析方法では，流れの変数として速度と圧力が用いられ，それらを支配する複数の偏微分方程式が解析され，上記計算格子の頂点において流れの変数が計算されていた．これらの偏微分方程式は，微分を差分で近似表現することにより，代数方程式に置き換えられるが，非線形項が数値不安定性を有するため，固気二相流のレイノルズ数が高い場合には，計算過程で発散が生じるおそれがある．とくに，工業装置で見受けられる固気二相流は，レイノルズ数が高いものが多く，精度の高い解析ができないという問題があった． 
【０００５】
ところで，従来，気体あるいは液体のみの単相流の解析については渦法が有効に利用されている．ここで，渦法とは，渦度をもつ渦要素を追跡して渦度場の時間変化を求めるラグランジュ型解法であり，計算格子を用いることなく，大規模渦の形成や変形などの渦構造の発展過程を良好に計算できる特徴がある．下記非特許文献１〜８には，この渦法を利用して，固気二相流の挙動を解析する技術が開示されている．
【非特許文献１】Chung, J. N. and Troutt, T. R., Simulation of particle dispersion in an axisymmetric jet, J. Fluid Mech., 186(1988), 199-222.
【非特許文献２】内山知実・ほか2名，渦法による固気二相自由乱流の数値解法（数値モデルと二次元混合層への適用），日本機械学会論文集（B編），66-651(2000), 2853-2860.
【非特許文献３】内山知実・成瀬正章，Vortex in Cell法による固気二相自由乱流の数値解析，日本機械学会論文集（B編），69巻686号（2003），2200-2207.
【非特許文献４】内山知実・深瀬昭仁，三次元渦法による固気二相同軸円形噴流の数値解析，日本機械学会論文集（B編），70巻696号（2004），1957-1964．
【非特許文献５】内山知実・北野佳伸，自由落下粒子群が形成する粒子噴流の数値解析，日本機械学会論文集（B編），69巻684号（2003），1737-1745.
【非特許文献６】Uchiyama, T. and Yagami, H., Numerical Analysis of Gas-Particle Two-Phase Wake Flow by Vortex Method, Powder Technol., 149(2005), 112-120.
【非特許文献７】Uchiyama, T. and Fukase, A., Three-Dimensional Vortex Simulation of Particle-Laden Air Jet, Chem. Eng. Sci., 61(2006), 1767-1778.
【非特許文献８】内山知実・成瀬正章，自由落下粒子群が形成する粒子噴流の三次元渦法解析，日本機械学会論文集（B編），71巻707号（2005），1738-1745．
【特許文献１】特開２００４−１２６９２５公報
【特許文献２】特開２００６−８５６０３公報
【特許文献３】特開２００７−３０８９１６公報
【特許文献４】特開２００７−２８６９５５公報
【特許文献５】特開２００７−２８６５１４公報
【特許文献６】特開２００７−１７３８０公報
【特許文献７】特開２００６−２５９９１１公報
【特許文献８】特開２００６−７７２８８公報
【特許文献９】特開２０００−２１４１３４公報
【特許文献１０】特開平７−８４９９２公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら，上記非特許文献１に記載の技術は，まず粒子がない場合の気流を渦法により計算し，つぎに前記気流の中におかれた粒子の運動を計算するものであり，粒子が気流に及ぼす影響は無視されており，固気二相流の挙動を精度良く解析できるものではなかった．
【０００７】
上記非特許文献２および３に記載の技術は，本願が対象とする，渦法による固気二相流のシミュレーションに関するものである．しかし，粒子が気流に及ぼす影響を計算する際に解析領域を計算格子（グリッド）に分割する手法を用いている．よって，本願が対象とする計算格子を用いない手法（グリッドフリー手法）は導入されていない．本願のグリッドフリー手法は，上記非特許文献２および３に記載の方程式およびアルゴリズムとは全くことなるため，これらの手法を適用することはできない．
【０００８】
上記非特許文献４に記載の技術は，上記非特許文献２に記載の技術を三次元解析用に拡張したものであり，上記非特許文献２の技術と同様，解析領域を計算格子（グリッド）に分割する手法を用いている．よって，本願が対象とする計算格子を用いない手法（グリッドフリー手法）は導入されていない．また，本願のグリッドフリー手法は，上記非特許文献４に記載の方程式およびアルゴリズムとは全くことなるため，これらの手法を適用することはできない．
【０００９】
上記非特許文献５〜８は，上記非特許文献２〜４のいずれかに記載のシミュレーション技術を用いて，様々な固気二相流を解析した場合の解析結果を記述したものである．よって，本願のグリッドフリー手法とは関連しない内容である．
【００１０】
特許文献１の技術は，渦法を用いて流体の流れを解析するものであるが，水または空気などの単相流を解析の対象とするものであり，本願が解析対象とする固気二相流は解析できない．
特許文献２の技術は，液中に分散した気泡が液体と相互作用を及ぼし合いながら流れる気泡流に対する渦法による解析方法に関するものであり，本願が解析の対象とする固気二相流の解析には適用できない．
【００１１】
特許文献３〜１０には，固体粒子の運動をその支配方程式から計算によって解析しようとする技術が開示されているが，いずれも，本願が対象とする計算格子を使用しないグリッドフリー手法を用いた渦法に準拠した解法ではなく，かつ気流中を運動する粒子に関するものでもない．
【００１２】
すなわち，特許文献３の技術は，海底に投入された土砂を対象として，海浜流による移動現象をシミュレーションするものであり，固体粒子を含む気流すなわち固気二相流に関するものではない．
特許文献４の技術は，容器内の粒子の挙動シミュレーションに関するものであるが，渦法に基づく解法ではない．
【００１３】
特許文献５の技術は，離散粒子法による粒子挙動シミュレーション装置、粒子挙動シミュレーション方法、及びコンピュータプログラムを提供するものであり，粒子間の衝突や接触の計算に関するものであり，気流中の粒子挙動を対象としていない．
特許文献６の技術は，多数の粒子が存在する場合の粒子間衝突のためのシミュレーションに関するものであり，流体の計算は対象としていない．
【００１４】
特許文献７の技術は，個別要素法を用いて粒子間の影響を考慮する解析方法ならびにコンピュータプログラムに関するものであり，固気二相流としての扱いはしていない．
特許文献８の技術は，自熔炉内の精鉱粒子同士の衝突および成長現象を考慮した燃焼現象を予測する手法を提供するものであり，気体の流動については省略している．
【００１５】
特許文献９の技術は，粉粒体を収容した容器内での撹拌部材を用いた混合撹拌操作に伴って発生する静電気量を基にしたシミュレーションであり，流動に関するシミュレーションではない．
特許文献１０の技術は，回転する混合容器内の粒子挙動の予測に関するものであり，粒子間の衝突による運動を解析しており，流体運動に起因する粒子挙動までは扱ってはいない．
【００１６】
本発明は上記のような事情に基づいて完成されたものであって，気流中に粒子が混在する固気二相流を対象として，解析領域を計算格子（グリッド）に分割することなく精度よく解析できる，グリッドフリー手法を用いた固気二相流シミュレーションプログラム及びそれを記憶した記憶媒体並びに固気二相流シミュレーション装置を提供することを目的とする．
【課題を解決するための手段】
【００１７】
上記の目的を達成するための手段として，請求項１の発明に関するグリッドフリー手法を用いた固気二相流シミュレーションプログラムは，気流中に微小な固体粒子を含む固気二相流を解析するための固気二相流シミュレーションプログラムであって，コンピュータに，次のステップ（ａ）〜（ｉ）を実行させて，前記固気二相流の挙動を解析する．
（ａ）前記解析領域内に存在する各粒子の速度をラグランジュ的に計算するステップ．
（ｂ）前記各粒子の位置をラグランジュ的に計算するステップ．
（ｃ）前記各粒子の周りに４つの正方形のセルを配置するステップ．
（ｄ）前記各粒子に作用する流体抗力を，前記ステップ（ａ）で算出された速度に基づいて計算するステップ．
（ｅ）前記ステップ（ｄ）で算出された流体抗力に基づいて，前記各セルにおける循環変化量を，気流の運動量保存方程式から導出される渦度方程式と，レイノルズの輸送定理とに基づいて計算するステップ．
（ｆ）前記ステップ（ｅ）で算出された前記各セルにおける循環変化量に基づいて，当該セル内の各渦要素の循環を計算し，各渦要素の位置と循環を記憶手段（１）に保存するステップ．
（ｇ）各渦要素のコア半径をラグランジュ的に計算して記憶手段（２）に保存するステップ．
（ｈ）前記ステップ（ｆ）および前記ステップ（ｇ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，各渦要素の移流すなわち位置をラグランジュ的に計算し，各渦要素の位置，循環およびコア半径を記憶手段（３）に保存するステップ．
（ｉ）前記ステップ（ｈ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，前記評価点における気体の速度をビオサバールの式を用いて計算し，記憶手段（４）に保存するステップ．
【００１８】
請求項２の発明は，上記請求項１に記載のグリッドフリー手法を用いた固気二相流シミュレーションプログラムを記憶した記憶媒体である．なお，記憶媒体としては，フロッピー（登録商標）ディスク，ＣＤ−ＲＯＭ等が好適である．
【００１９】
請求項３に係る固気二相流シミュレーション装置は，気流中に微小な固体粒子を含む固気二相流を解析するためのグリッドフリー手法を用いた固気二相流シミュレーション装置であって，次の手段（ａ）〜（ｉ）を備え，これらの手段を実行させて，前記固気二相流の挙動を解析する．
（ａ）前記解析領域内に存在する各粒子の速度をラグランジュ的に計算する手段．
（ｂ）前記各粒子の位置をラグランジュ的に計算する手段．
（ｃ）前記各粒子の周りに４つの正方形のセルを配置する手段．
（ｄ）前記各粒子に作用する流体抗力を，前記手段（ａ）で算出された速度に基づいて計算する手段．
（ｅ）前記手段（ｄ）で算出された流体抗力に基づいて，前記各セルにおける循環変化量を，気流の運動量保存方程式から導出される渦度方程式と，レイノルズの輸送定理とに基づいて計算する手段．
（ｆ）前記手段（ｅ）で算出された前記各セルにおける循環変化量に基づいて，当該セル内の各渦要素の循環を計算し，各渦要素の位置と循環を記憶手段（１）に保存する手段．
（ｇ）各渦要素のコア半径をラグランジュ的に計算して記憶手段（２）に保存する手段．
（ｈ）前記手段（ｆ）および前記手段（ｇ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，各渦要素の移流すなわち位置をラグランジュ的に計算し，各渦要素の位置，循環およびコア半径を記憶手段（３）に保存する手段．
（ｉ）前記手段（ｈ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，前記評価点における気体の速度をビオサバールの式を用いて計算し，記憶手段（４）に保存する手段．
【発明の効果】
【００２０】
本発明は，微小な固体粒子が気流中に分散し，気流と相互作用を及ぼし合いながら流れる固気二相流の渦度場に着目して複数の渦要素で離散化し，気流に及ぼす粒子の力を考慮した運動量保存関係から渦度方程式を導出し，この渦度方程式のラグランジュ計算により渦要素の挙動を求めるものである．
解析領域を計算格子（グリッド）に分割する必要がない解析方法である．格子の寸法や配置は，解析精度に影響を及ぼすため，格子分割には細心の注意と多大な作業が要求される．よって，格子分割作業を省略できるこのグリッドフリー手法は，工業装置内など複雑な幾何形状をもつ流れ場の解析に極めて有用である．
【００２１】
また，渦度輸送方程式のラグランジュ計算には非線形項が現れないため，レイノルズ数に依存することなく，高レイノルズ数の固気二相流についても安定で高精度の解析が可能となる．また，渦度の高い領域には渦要素が能動的に密集する特徴があるため，高い解析精度を自動的に確保できる．
【００２２】
さらに，渦度を求めるに際して，渦要素のコア半径を時間の経過とともに増大させるコアスプレッディング法を適用することが望ましい．このようにすれば，粘性拡散による渦度の減衰が考慮され，適切な渦度を求めることができる．
【００２３】
＜理論＞
次の本発明の理論について説明する．
１．記号
本発明において下記に登場する各符号は次のように定義されている．
A：粒子の周りに設置するセルの一つあたりの面積 
CD：粒子の抵抗係数 
d：粒子直径 
FD：単位体積の気体が粒子から受ける力 
fD：粒子に作用する流体抗力 
g：重力加速度 
L：正方形セルの一辺の長さ 
Nv：解析領域内の渦要素数   
p：圧力 
t：時間 
u：速度 
vp：気体に対する粒子の相対速度 =up-ug
x, y：空間座標 
Γ：循環 
ν：気体の動粘度 
ρ：密度 
σ：渦要素のコア半径 
ω：渦度 =∇×ug
添字
g：気体 
p：粒子 
【００２４】
２．仮定
本発明は以下の仮定を前提としている．
（ａ）気体は非圧縮である．
（ｂ）粒子の密度は気体に比べ十分に大きい．
（ｃ）粒子は球形であり，均一な直径をもつ．
（ｄ）流れ場全体に及ぼす粒子間衝突の影響は小さく，無視できる．
【００２５】
３．基礎式
本発明では，所定の設定データ（後述する入力データ）に基づいて，計算に必要な条件を設定し，流れの基礎方程式に基づいて計算を実行する．以下，計算に用いる流れの基礎方程式について説明する．
【００２６】
（１）気体の支配方程式
気体の質量および運動量の保存方程式は，上記仮定（ａ）〜（ｄ）を前提とした場合，それぞれ次式で表される．
【数１】
2009193110.tif 000003
【数２】
2009193110.tif 000004
【００２７】
（２）粒子の支配方程式
粒子に作用する支配的な力は，上記仮定（ａ）〜（ｄ）を前提とした場合，流体抗力と重力であり，次式で与えられる．
【数３】
2009193110.tif 000005
ここで，mは粒子の質量であり，流体抗力fDは次式で与えられる．
【数４】
2009193110.tif 000006
ただし，dは粒子直径であり，抗力係数CDは次式で定められる．
【数５】
2009193110.tif 000007
ここで，Rep=d|ug-up|/νである．
【００２８】
（３）渦要素による渦度場の離散化
二次元流れ場を解析の対象とすれば，気体の渦度方程式は，上記数式２の回転をとり上記数式１を代入して得られる次式で表される．
【数６】
2009193110.tif 000008
一方，任意の位置xにおける気体の速度は次のビオサバールの式により求められる．
【数７】
2009193110.tif 000009
ここで，ug0は一様流またはポテンシャル流れの速度を表す．
【００２９】
つぎに，渦度場を多数の渦要素により離散化する．単相流解析に対する渦要素モデルを適用し，コア構造をもつ渦要素を導入する．渦要素αの循環をΓα，コア半径をσα，位置ベクトルをxα，とすれば，渦要素αによる位置xの渦度は次式で与えられる．
【数８】
2009193110.tif 000010
ここで，コア関数fは次式で与えられる．
【数９】
2009193110.tif 000011
【００３０】
渦度場をNv個の渦要素で離散化すれば，気体の速度は上記数式７と上記数式８から得られる次式で表される．
【数１０】
2009193110.tif 000012
ここで，kは渦度ベクトルと同一方向の単位ベクトルであり，関数gは次式で与えられる．
【数１１】
2009193110.tif 000013
【００３１】
また，渦要素は気流速度で移流するので，次のラグランジュ解析から計算できる．
【数１２】
2009193110.tif 000014
【００３２】
（４）粘性拡散によるコア半径の変化
渦度は粘性拡散により減衰する．ここで，液体のみまたは気体のみの単相流解析では，コア半径を時間の経過とともに増大させるコアスプレッディング（Core spreading）法により考慮される．そこで，本発明においても，これを準用し，次式をラグランジュ解析してコア半径の時間変化を求める．
【数１３】
2009193110.tif 000015
なお，Core spreading法は，Leonard, A., Vortex methods for flow simulation, J. Comput. Phys., 37(1980), 289-335. に記載されている．
【００３３】
（５）粒子運動による循環の変化
任意の閉曲線まわりの循環Γの時間変化率は，次式で表される．
【数１４】
2009193110.tif 000016
上記数式１４にレイノルズの輸送定理を適用したのち，上記数式１と上記数式２を代入し，変形すれば次式を得る．
【数１５】
2009193110.tif 000017
ここで，drは線素ベクトルである．また，粘性拡散項は，上記数式１３で考慮されているので，代入に際し無視している．
【００３４】
（６）粒子運動による循環変化の計算方法
図２に示すように，粒子Ｐの周囲に４つの正方形のセル６を設ける．セル６は，粒子の位置が既知であれば容易に配置できる．よって，その配置に要する作業は，解析領域を多数の格子（グリッド）に分割する従来技術に比べて著しく簡単である．とくに，複雑な幾何形状をもつ解析領域を対象とする場合には，特段の効果が現れる．セル６は座標軸xおよびyに沿っており，一辺の長さをLとする．気体Ｇに対する粒子の相対速度をvp（=up-ug）とすれば，-vpの方向に流体抗力fDが作用する．fDがx軸となす角度をθとする．この場合，粒子Ｐが位置するセル頂点Ｏに作用するFDのx方向成分FDxおよびy方向成分FDyは，それぞれ次式で与えられる．
【数１６】
2009193110.tif 000018
ここで，Δzはx-y平面に垂直方向（図２で紙面奥行方向）に粒子が等間隔で分布するもと仮定したときの分布間隔である．
セル頂点７の間でFDが線形変化するものと仮定し，４つのセル１，２，３，４に上記数式１５を適用すれば，セル１，２，３，４の循環変化ΔΓ1，ΔΓ2，ΔΓ3，ΔΓ4は次式で表される．
【数１７】
2009193110.tif 000019
【００３５】
各粒子Ｐのまわりに設置した４つのセルに対し，セルβ（β=1,2,3,4）にnv個の渦要素が存在する場合には，渦要素一つ当たりの循環をΔΓβ/nvとする．渦要素が存在しない場合には循環ΔΓβをもつ１つの渦要素をセル中央から発生させる．
【００３６】
いま，試みに，粒子運動により４つのセル６から渦要素５が一つずつ新たに発生する場合を考えてみる．θ=９０°のときには，図３（ａ）に示すように，気流に対する粒子Ｐの運動方向に関して，渦要素５の配置が左右対称となる．さらに，数式１７からわかるように循環の絶対値も対称となる．４つの渦要素５が誘起する粒子Ｐまわりの気流を計算すると図３（ｂ）のようになり，粒子運動方向に関して気流も左右対称となる．粒子は球形であり進行方向に対して形状は軸対称であるから，その運動に起因する流れは，二次元解析の場合には軸に関して線対称すなわち左右対称であらねばならない．よって，図３（ｂ）では合理的な速度場が求められている．θ=０°，１８０°，２７０°の場合にも同様な気流速度が求められる．
【００３７】
一方，θの値が０°，９０°，１８０°，２７０°でない場合には，図４（ａ）に示すように渦要素５の配置と循環の絶対値が粒子進行方向に関して対称でない．たとえばθ=６５°の場合の気流を求めると，図４（ｂ）に示すように非対称な速度分布となり，適切な結果は得られない． 
【００３８】
そこで，図５に示すように，粒子位置（セル頂点Ｏ）からセル頂点Ｎに向かうベクトルがつねにfDと一致するように４つのセル６を配置する．このセル配置は，図３（ａ）の配置に相当するため，粒子まわりの左右対称な気流を正確に計算できる［図３ （ｂ）参照］．
【００３９】
（７）解析手順
時刻t=tでの粒子および気体の流れが既知であれば，時刻t=t+Δtの流れを，例えば以下の手順で計算できる．なお，図６に示すように，解析領域Ｈには渦要素５と粒子Ｐが存在し，気体の速度を計算する評価点８が予め配置されている．
（ａ）粒子Ｐの運動を数式３から求める．
（ｂ）粒子Ｐの周囲に４つの正方形セル６を設置し，各セル６における循環変化量ΔΓ1〜ΔΓ4を数式１７においてθ=９０°として求める．
（ｃ）各セル６における循環変化量ΔΓ1〜ΔΓ4を当該セル６内の渦要素の循環に付与する．セル内に渦要素がない場合には，循環変化量に等しい循環をもつ１つの渦要素をセル中央から発生させる．
（ｄ）渦要素５のコア半径σを数式１３から求める．
（ｅ）渦要素５の位置を数式１２から求める．
（ｆ）評価点８における気流の速度ugを数式１０から求める．
【００４０】
このような手順を予め定めた時間間隔毎に繰り返し実行することで，計算格子を用いることなく（グリッドフリー），固気二相流を高精度で解析することができる．
【発明を実施するための最良の形態】
【００４１】
以下，本発明の一実施形態を前述した理論および図１〜図１０によって説明する．
【００４２】
本実施形態に係るグリッドフリー手法を用いた固気二相流シミュレーションプログラムは，後述する情報記憶装置１６に記録され，パーソナルコンピュータ（以下「コンピュータ１１」という）に導入することにより，そのコンピュータ１１に後述する固気二相流の挙動の解析ステップを実行させるためのものである．
【００４３】
１．コンピュータとその周辺機器の構成
図７には，コンピュータ１１とその周辺機器の構成を示した．コンピュータ１１には，演算および周辺機器の制御等を行うための中央処理部１２，この中央処理部１２と周辺機器との接続を行うためのＩ／Ｏインターフェイス１３，および周辺機器である表示装置１４，印刷装置１５，情報記憶装置１６および入力装置１７が接続されている．
【００４４】
なお，表示装置１４には，カソードレイチューブ（ＣＲＴ）の他に液晶装置，液晶プロジェクタ等が含まれる．また，印刷装置１５は，コンピュータ１１に直接に連結されている必要はなく，ＬＡＮ等のネットワークを介して連結されていてもよい．また，入力装置１７は，例えばキーボードやマウスなどのコンソールであって，ここでの操作によって解析領域Ｈの境界を設定する座標位置データや流体特性を設定する入力データなどを入力することができる．
【００４５】
中央処理部１２の詳細については図示しないが，中央演算処理装置（ＣＰＵ）のほかに，ＲＯＭ，各種ＲＡＭ，チップセット等から構成されている．情報記憶装置６は，ハードディスクやフロッピー（登録商標）ディスク，およびＭＯやＣＤ−ＲＯＭ等であり，コンピュータ１のプログラムや各種データを読み込みまたは記録することができる．
【００４６】
２．固気二相流の解析
さて，本実施形態の作用を上記中央処理部１２にて実行される固気二相流シミュレーションプログラムのフローチャートを参照しつつ説明する．
【００４７】
（１）計算に必要な入力データの入力および境界条件の設定
まず，上記入力装置１７の操作によって，計算に必要な入力データを入力する．ここで，入力データとしては，気体Ｇの密度ρgおよび動粘度νg，固体粒子Ｐの密度ρpおよび直径d，解析領域Ｈの入口における単位時間当たりの粒子Ｐの流入個数，解析領域Ｈの幾何形状および寸法，解析時間，解析時間タイミングΔtである．
【００４８】
（１）条件設定
本実施形態では，壁面で囲まれていない自由空間における，粒子Ｐを含む気流（固気二相流）に対する数値解析である．より，具体的には，図８に示すような，粒子Ｐが静止空気中に落下し，落下粒子が周囲に気体Ｇの流れ（気流）を誘起し，粒子Ｐと気体Ｇが相互作用を及ぼし合いながら流れる固気二相流を解析対象としている．
【００４９】
（２）解析処理
コンピュータ１１は，上記条件設定後，図９に示すフローチャートの処理を上記解析時間タイミングΔtで繰り返し実行し，固気二相流の挙動を数値解析する．
【００５０】
ステップＳ１で解析時間タイミングになったかどうかを判断し，解析時間タイミングになったとき（ステップＳ１で「Ｙ」），ステップＳ２で粒子Ｐの速度upを，数式３から計算する．すなわち，本発明の請求項に記載の「ステップ（ａ）」および「手段（ａ）」でいう「解析領域内に存在する各粒子の速度をラグランジュ的に計算」に相当する処理を実行する．より具体的には，数式３にオイラーの前進差分法を適用し，ラグランジュ解析により解析領域Ｈ内の各粒子Ｐを追跡するのである．この際，数式３のfDは数式４で与えられるが，数式４の右辺の気体Ｇと粒子Ｐの速度には，前回の解析時間タイミングで算出された速度を用いるのである．なお，気体Ｇの速度は，解析領域Ｈ内に予め配置した評価点８の速度の補間から定められるのである．補間には，評価したい位置の近傍点８を選択し，内挿補間法を利用すればよい．
【００５１】
ステップＳ３では，上記ステップＳ２で算出された粒子Ｐの速度upをもとにして，粒子Ｐの位置をラグランジュ法により計算する．すなわち，本発明の請求項に記載の「ステップ（ｂ）」および「手段（ｂ）」でいう「各粒子の位置をラグランジュ的に計算」に相当する処理を実行する．より具体的には，速度が既知であれば，その時間積分により位置が求められるので，上記粒子速度の場合と同様，オイラーの前進差分法により計算するのである．
【００５２】
ステップＳ４では，各粒子の周りに４つの正方形のセル６を配置する．すなわち，本発明の請求項に記載の「ステップ（ｃ）」および「手段（ｃ）」でいう「各粒子の周りに４つの正方形のセルを配置」に相当する処理を実行する．より具体的には，気体Ｇに対する粒子Ｐの相対速度vp（=up-ug）を計算し，図５に示したように，粒子位置（セル頂点Ｏ）からセル頂点Ｎに向かうベクトルがvpと一致するようにセル６を配置する．
【００５３】
ステップＳ５では，各粒子に作用する流体抗力fDを，前回の解析時間タイミングで算出された気体Ｇの速度ugと，ステップＳ１で計算した粒子Ｐの速度upを用いて，数式４から算出する．すなわち，本発明の請求項に記載の「ステップ（ｄ）」および「手段（ｄ）」でいう「各粒子に作用する流体抗力を，ステップ（ａ）あるいは手段（ａ）で算出された速度に基づいて計算」に相当する処理を実行する．なお，気体Ｇの速度ugは，解析領域Ｈ内に予め配置した評価点の速度の補間から定められるのである．
【００５４】
ステップＳ６では，各粒子のまわりの４つのセルにおける循環変化量ΔΓ1〜ΔΓ4を，上記ステップＳ５で算出された各粒子に作用する流体抗力に基づいて，数式１７から計算する．すなわち，本発明の請求項に記載の「ステップ（ｅ）」および「手段（ｅ）」でいう「ステップ（ｄ）あるいは手段（ｄ）で算出された流体抗力に基づいて，前記各セルにおける循環変化量を，気流の運動量保存方程式から導出される渦度方程式と，レイノルズの輸送定理とに基づいて計算」に相当する処理を実行する．なお，ステップＳ４におけるセルの配置からして，数式１７ではθ=90°として計算する．
【００５５】
ステップＳ７では，上記ステップＳ６で算出された各セルの循環変化量データを基にして，各セルに存在する渦要素の循環を計算し，各渦要素の位置と循環を，例えば上記情報記録装置１６の所定の保存領域あるいはＲＡＭの所定領域に保存する．すなわち，本発明の請求項に記載の「ステップ（ｆ）」および「手段（ｆ）」でいう「ステップ（ｅ）あるいは手段（ｅ）によって保存された前記各セルにおける循環変化量データを読み出して，当該セル内の各渦要素の循環を計算し，各渦要素の位置と循環を記憶手段（１）に保存」に相当する処理を実行する．具体的には，循環変化量がΔΓのセルの内部に渦要素がnv個存在する場合，各渦要素の循環の変化量をΔΓ/nvとする．渦要素が存在しない場合には，循環ΔΓをもつ渦要素を１つ新しく発生させる． 
【００５６】
ステップＳ８では，時間変化する渦要素のコア半径σを数式１３から計算し，例えば上記情報記録装置１６の所定の保存領域あるいはＲＡＭの所定領域に保存する．すなわち，本発明の請求項に記載の「ステップ（ｇ）」および「手段（ｇ）」でいう「各渦要素のコア半径をラグランジュ的に計算して記憶手段（２）に保存」に相当する処理を実行する．このラグランジュ計算には，オイラーの前進差分法が適用できる．
【００５７】
ステップＳ９では，渦要素５の移流を数式１２から計算し，例えば上記情報記録装置１６の所定の保存領域あるいはＲＡＭの所定領域に保存する．すなわち，本発明の請求項に記載の「ステップ（ｈ）」および「手段（ｈ）」でいう「ステップ（ｆ）およびステップ（ｇ）あるいは手段（ｆ）および手段（ｇ）によって保存された渦要素の位置，循環およびコア半径データを読み出して，各渦要素の移流すなわち位置をラグランジュ的に計算し，各渦要素の位置，循環およびコア半径を記憶手段（３）に保存」に相当する処理を実行する．より具体的には，上記粒子Ｐの運動計算の場合と同様，オイラーの前進差分法により数式１２を解析して，解析領域Ｈにおける渦要素１６の位置を追跡するのである．
【００５８】
ステップＳ１０では，ステップＳ７〜Ｓ９で計算した渦要素の位置，循環およびコア半径を数式１０に代入して，解析領域に予め設定した評価点における気体Ｇの速度を計算し，例えば上記情報記録装置１６の所定の保存領域内またはＲＡＭの所定領域に保存する．すなわち，本発明の請求項に記載の「ステップ（ｉ）」および「手段（ｉ）」でいう「ステップ（ｈ）あるいは手段（ｈ）によって保存された各渦要素の位置，循環およびコア半径データを読み出して，前記評価点における気体の速度をビオサバールの式を用いて計算し，記憶手段（４）に保存」に相当する処理を実行する．ステップ１２で解析継続の場合（ステップＳ１２で「Ｙ」），ステップＳ１に戻る．
【００５９】
以上のステップを繰り返し実行することにより，粒子Ｐと気体Ｇの流れ，すなわち固気二相流の時間変化をシミュレートすることができる．なお，以上の演算が終了すると，コンピュータ１１は，その解析データを例えば表示装置１４や印刷装置１５に引き渡す．解析データは，粒子Ｐおよび気体Ｇの速度upおよびug，粒子Ｐおよび渦要素５の位置座標であり，一般的に利用できる可視化ソフトウェアを用いることにより，ベクトル表示や等高線表示できる．
【００６０】
なお，図１０は計算結果の一例であり，粒子Ｐが空気中を自由落下する場合に発生する粒子，渦要素および気流の挙動であり，上記解析データをもとにして，粒子Ｐ，渦要素６，気流の速度ベクトルを表示してある．
【００６１】
３．本実施形態の効果
固気二相流の解析に従来用いられている方法は，例外なく，解析の前処理として解析領域を計算格子（グリッド）に分割する必要がある．格子は解析領域の全域をメッシュ状に覆うものであり，計算格子とよばれている．各格子の形状は，一般に四角形である．各格子の頂点で流れの変数（速度や圧力など）が定義され，上記解法により計算される．解析精度は格子分割の状態に依存するため，高い精度を確保するには，流れが急峻に変化するものと予測される領域は多数の細かな格子に分割する必要がある．また，格子の間には重なりや隙間があってはならないという幾何学的条件も満足しなければならない．したがって，工業装置の設計などの際に用いられる実務的シミュレーションでは，複雑な幾何形状をもつ流れ場を解析の対象とする場合が多く，格子分割に多大な作業が要求されてきた．本実施形態による固気二相流の解析方法では，解析領域を計算格子に分割する必要がない．このグリッドフリー技術は，解析の前処理を著しく簡単化するものであり，実務的シミュレーションの合理化をもたらす．
【００６２】
また，本実施形態による解析方法では，支配方程式の微分項を差分近似する必要がない．このため，従来の解析方法では数値不安定に起因して計算が破綻するような，高レイノルズ数の流れの解析にも適用可能である．
また，渦度の高い領域に渦要素が能動的に移流するため，高い解析精度を確保できる．
【００６３】
さらに，流体の速度や圧力を未知変数とする従来の解析方法では乱流モデルの利用が不可欠であったが，本実施形態の解析方法では乱流モデルを用いることなく，固気二相流を解析できる．
＜他の実施形態＞
本発明は上記記述および図面によって説明した実施形態に限定されるものではなく，例えば次のような実施形態も本発明の技術的範囲に含まれ，さらに，下記以外にも要旨を逸脱しない範囲内で種々変更して実施することができる．
（１）上記実施形態では，解析領域Ｈについて二次元解析を行った例を説明したが，本発明は三次元解析についても勿論適用することができる．
【００６４】
（２）上記実施形態では，粒子Ｐは球形としたが，これに限らず，高機能化を図る観点で，例えば楕円球などであってもよい．この場合，例えば入力装置１７にて粒子Ｐについての複数種の形状を選択可能とし，それらの複数種の粒子について計算してもよい．
【図面の簡単な説明】
【００６５】
【図１】流体解析における解析領域と計算格子を説明するための図
【図２】本発明の粒子のまわりに設置するセルを説明するための図
【図３】本発明のセルを適切に設置したときの図と得られる気体速度を説明するための図
【図４】本発明のセルを不適切に設置したときの図と得られる気体速度を説明するための図
【図５】本発明の粒子のまわりに設置するセルの配置を説明するための図
【図６】本発明の渦要素，固体粒子，気体速度評価点を説明するための図
【図７】本発明の一実施形態に係る固気二相流シミュレーションプログラムが導入されるコンピュータおよび周辺機器を示す図
【図８】本発明の計算が対象とする固気二相流を説明するための図
【図９】コンピュータによる処理内容を示すフローチャート
【図１０】本実施形態による解析結果を示す図
【符号の説明】
【００６６】
５・・・渦要素
６・・・セル
７・・・セル頂点
８・・・気体の速度の評価点
１１・・・コンピュータ
１２・・・中央処理部
１６・・・情報記録装置（記憶手段）
Ｇ・・・気体
Ｈ・・・解析領域
Ｐ・・・固体粒子
【図１】
2009193110.tif 000020
【図２】
2009193110.tif 000021
【図３】
2009193110.tif 000022
【図４】
2009193110.tif 000023
【図５】
2009193110.tif 000024
【図６】
2009193110.tif 000025
【図７】
2009193110.tif 000026
【図８】
2009193110.tif 000027
【図９】
2009193110.tif 000028
【図１０】
2009193110.tif 000029