数値流体計算方法及び数値流体計算装置
【課題】 Solid pressure項に関わる問題点を解消して分散相と連続相の混相流を精度良く解析することができる数値流体計算方法の提供。
【解決手段】 本発明は、連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法において、分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする。好ましい実施例では、繰り込み分として、分散相の真の体積分率と、見かけの体積分率または嵩と、の差を用いる。また、分散相の繰り込み分として、圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻す操作が付加されてもよい。
【解決手段】 本発明は、連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法において、分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする。好ましい実施例では、繰り込み分として、分散相の真の体積分率と、見かけの体積分率または嵩と、の差を用いる。また、分散相の繰り込み分として、圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻す操作が付加されてもよい。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算方法及び数値流体計算装置に関する。
【背景技術】
【０００２】
粒子を入れ物に詰めた場合、粒子1個の体積と個数の積である粒子の真の体積と、粒子間の隙間を含め粒子が存在する空間の体積（見かけの体積または嵩）とは異なり、球形粒子では六方最密充填構造の時に充填率が７４％であり、この最大充填率以上にならないことが一般的に知られている。
【０００３】
一方、流体計算の粒子相の運動量輸送方程式には付加的な項として、Solid Pressure項があり、非特許文献１および２などに記載されているようにGidaspowらは、粒子相の空隙率αSの関数として表現している。すなわち、
【０００４】
【数１】
ここで、G0は代表弾性率、cは圧縮係数、αsmは粒子相の最大空隙率である。この関連技術によれば、運動量の連続の式で粒子相の空隙率すなわち、粒子の充填率が設定した値に近づくにつれて大きな圧力勾配が生じ、粒子の充填率が設定した値に近づくような粒子相の流れを阻害することで、粒子が設定した充填率以上にならない状態を間接的にモデル化している。この圧力勾配の項は、あたかも粒子同士が反発力を持つようであり、Solid pressure項と呼ばれる。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Gidaspow, D., Multiphase Flow and Fluidization , Academic Press, (1994).
【非特許文献２】Bouillard, J. X., Lyczcowski, R. W. and Gidaspow, D., AIChE, J. , 35, 908(1996).
【非特許文献３】Wen, C.Y. and Yu, Y.H. (1996).CEP. Symp. Series, 62, 100-111.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上記記載した方法では、計算した充填率が、設定した最大充填率の値に近くなったときに、Solid pressure項の値を指数関数的に急激に大きくする必要がある。このため、数値的な安定性に問題があった。具体的には短い時間ステップで計算する必要があるとともに、最大充填率に達した計算セルを含む粒子の運動を計算する場合に、粒子相の運動量輸送方程式に付加したSolid pressure項の値の計算で、数値誤差が発生するなどという不具合があった。この影響に起因して、設定した最大充填率を超える計算結果となったり、計算値が発散して計算できないという欠点があった。
【０００７】
そこで、本発明は、Solid pressure項に関わる問題点を解消して分散相と連続相の混相流を精度良く解析することができる数値流体計算方法及び数値流体計算装置の提供を目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するため、本発明の一局面によれば、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算方法において、
分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする数値流体計算方法が提供される。
【０００９】
本発明のその他の一局面によれば、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算装置であって、
分散相の体積分率に連続相の体積分率の一部が繰り込まれるように、分散相と連続相をモデル化するモデル化手段と、
前記モデル化された分散相と連続相の混相流を解析する解析手段とを備えることを特徴とする、数値流体計算装置が提供される。
【発明の効果】
【００１０】
本発明によれば、Solid pressure項に関わる問題点を解消して分散相と連続相の混相流を精度良く解析することができる数値流体計算方法及び数値流体計算装置が得られる。
【図面の簡単な説明】
【００１１】
【図１】本発明による数値計算方法の基本概念の説明図である。
【図２】本発明の実施態様１による数値計算装置１の要部構成の一例を示す図である。
【図３】本発明の実施態様１による数値計算方法の要部の一例の流れを示すフローチャートである。
【図４】本発明の実施態様２による数値計算装置２の要部構成の一例を示す図である。
【図５】図５（A）は、沈殿物底からの高さと圧力の関係を示すグラフであり、図５（B）は、沈殿物底からの高さと設定充填率の関係を示すグラフである。
【図６】本発明の実施態様２による数値計算方法の要部の一例の流れを示すフローチャートである。
【図７】実施例１による体積分率最大値の時間変化を示す図である。
【図８】実施例１による計算モデル断面での密度を示す図である。
【図９】実施例１による個々の計算セル中での真の分散相の体積分率の分布を示す図である。
【図１０】実施例２による体積分率最大値の時間変化を示す図である。
【図１１】実施例２による計算モデル断面での密度を示す図である。
【図１２】実施例２による個々の計算セル中での真の分散相の体積分率の分布を示す図である。
【図１３】Solid Pressureを用いて分散相の充填率が上限値を持つようにモデル化する計算手順（比較例）を示すフローチャートである。
【図１４】比較例１による計算結果を示す図である。
【図１５】比較例２による計算結果を示す図である。
【発明を実施するための形態】
【００１２】
以下、図面を参照して、本発明を実施するための最良の形態の説明を行う。
【００１３】
本発明による数値計算方法は、連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法に関する。
【００１４】
ここで取り扱う流体の運動方程式は、例えば次のような式で表される。
【００１５】
【数２】
これを、別の書き方で表現すると、次の通りである。
【００１６】
【数３】
この式で、Kは外力（重力加速度ｇ等）であり、ρは密度であり、ｐは圧力分布である。
【００１７】
Eulerian-Euleria法では、固体粒子も粒ではなく連続体として考えてあたかも２つの流体の流れを計算するかのように、多相流の流れ場を計算する手法である。尚、ANSYS社の汎用流体計算ソフト「CFX」等においても、固体粒子を連続体として扱っており、かかる手法が一般的である。
【００１８】
従って、Eulerian-Euleria法では、連続相と分散相のそれぞれの相に対して、数２の運動方程式が適用される。この際、各相の運動方程式おける外力の項Ｋには、互いの相から受ける抵抗（相間の運動量交換量）が含まれることになる。
【００１９】
ところで、上述の如く、粒子を入れ物に詰めた場合、粒子1個の体積と個数の積である粒子の真の体積と、粒子間の隙間を含め粒子が存在する空間の体積（見かけの体積または嵩）とは異なり、例えば球形粒子では六方最密充填構造の時に充填率が７４％であり、この最大充填率以上にならないことが一般的に知られている（図１（Ａ）及び図１（Ｂ）参照）。
【００２０】
このような粒子の充填特徴を考慮するために、本発明による数値計算方法は、従来使用されてきたSolid pressure項を使用せずに、それに代えて、分散相の体積分率に連続相の体積分率の一部を繰り込むことを特徴とする。
【００２１】
図１は、このような本発明による数値計算方法の基本概念の説明図である。
【００２２】
図１（Ａ）は、球形の粒子が液体底に沈殿している様子を模式的に示す図である。図１（Ｂ）は、図１（Ａ）の粒子の実際の体積分率を模式的に示す図である。図１（Ｃ）は、分散相の体積分率に連続相の体積分率の一部を繰り込む方法を概念的に示す図である。
【００２３】
概念的には固体は球形の粒であるため、図１（Ａ）に示すように、固体を完全に充填することはできず、固体間の隙間に液体が存在する。従って、図１（Ｂ）に示すように、分散相の体積分率は、嵩よりも小さい。
【００２４】
そこで、本発明では、図１（Ｃ）に示すように、分散相の体積分率に連続相の体積分率の一部を繰り込むことで、図１（Ａ）に示すような固体充填時の状態をモデル化する。即ち、図１（Ｃ）に示すように、固体成分の内訳と嵩との差を考慮して、球形の固体の隙間の液体を分散相に繰り込む。概念的には、分散相を球形の粒子ではなく、ブロック状の粒子と考え、球形であった場合に隙間に存在する連続相の液体（図１（Ｃ）の符号Ａで示す部分）を、分散相の一部として繰り込むというモデルである。尚、このモデルでは、ブロック状の粒子間には、Solid pressure項を付加することで発生するような反発力は生じない。
【００２５】
このように本発明によれば、分散相の体積分率に連続相の体積分率の一部を繰り込むことで、Solid pressure項を使用せずに、固体成分の内訳と嵩との差を補償することができる。これにより、Solid pressure項の問題点を解消して分散相と連続相の混相流を精度良く解析することができる。
【００２６】
尚、図１（Ｃ）を参照したモデルは、固体の大きさは、隙間の液体分だけ大きくした繰り込み固体の大きさで計算するため、個体の大きさに関連するパラメータ（例えば抵抗係数等）は、好ましくは、繰り込み前の元の大きさを基準にして（換算して）計算する。この換算方法については後述する。
【００２７】
次に、以上の本発明の基本概念を実現するための好ましい実施態様について説明する。
［実施態様１］
図２は、本発明の実施態様１による数値計算装置１の要部構成の一例を示す図である。
【００２８】
本実施態様１の数値計算装置１は、図２に示すように、入力パラメータ生成部１０と、繰り込みモデル化部１２と、数値解析部１４と、パラメータ補正部１６とを含む。尚、これらの各部１０，１２，１４，１６は、ソフトウェアとハードウェアの組み合わせにより実現される。この際、これらの各部１０，１２，１４，１６は、共通のソフトウェアで実現されてもよいし、複数のソフトウェアにより協動して実現されてもよい。また、同様に、これらの各部１０，１２，１４，１６は、共通のハードウェアで実現されてもよいし、複数のハードウェアにより協動して実現されてもよい。
【００２９】
入力パラメータ生成部１０は、外部からの入力データ（典型的には、ユーザからの入力データ）に基づいて、数値解析部１４における解析に必要な入力パラメータ（解析条件）を生成（設定）する。入力パラメータは、例えば連続相の特性（水といった媒体の特定を含む）、分散相の粒子の密度、最大充填率等を含む。
【００３０】
繰り込みモデル化部１２は、図１（Ｃ）を参照して説明したモデル化を行う。即ち、繰り込みモデル化部１２は、入力パラメータ生成部１０と協動して、分散相の体積分率に連続相の体積分率の一部が繰り込まれるような入力パラメータを生成する。
【００３１】
具体的には、繰り込みモデル化部１２は、連続相の一部を分散相に繰り込むために、分散相を粒子モデルとし、粒子径として真の粒子径Dp,realと見かけの粒子径Dpとの関係を(2)式で定義する。
【００３２】
【数４】
ここで、εmaxは、最大充填率である。尚、繰り込みモデル化部１２の処理の詳細については後述する。
【００３３】
数値解析部１４は、入力パラメータ生成部１０で生成された入力パラメータと、後述のパラメータ補正部１６から供給されるパラメータとを用いて、連続体モデル（Eulerian-Euleria法）による分散相と連続相の混相流を解析する。数値解析部１４における流体解析方法は、連続体モデルを用いるものであればその詳細は任意である。例えば、数値解析部１４は、流体ソルバーと称される汎用流体計算ソフトとＣＰＵ等のハードウェアにより実現されてもよい。かかる汎用流体計算ソフトとしては、ANSYS社の汎用流体計算ソフトCFXが知られている。数値解析部１４は、遠隔位置にあるハードウェア資源（例えばスーパーコンピューター）と協動して実現されてもよい。
【００３４】
ここで、数値解析部１４において考慮される相間の運動量交換量について説明する。
【００３５】
相間の運動量交換量は、単位体積当たりの抵抗Ｄd,cが、１つの分散相粒子が運動するときに周りの連続相から受ける抵抗Ｄsと単位体積当りの粒子数ｎpの積であるから、式（３）のようになる。
【００３６】
【数５】
ここで、CDは抵抗係数、ρCは連続相の密度、Ud及びUcはそれぞれ分散相と連続相の速度、rpは分散相の体積分率である。ここで、分散相の真の体積分率をrp,realとし、分散相に式（２）の繰り込みを行うと以下の式（３'）に示すように、抵抗係数CDを、真の抵抗係数CD,realで置き換えたのみであり、同じ形の式に書くことができる。
【００３７】
【数６】
すなわち、真の抵抗係数CD,realを用いることにより、数値解析部１４内は標準的な式（３）に従った計算を行うが、式（３'）のように、実質的には真の粒子径Dp,realと分散相の真の体積分率rp,realを用いた計算結果を得ることが出来る。
なお、通常の流体ソルバーを用いる数値解析部１４では、相間の運動量交換に関する抵抗係数として標準的な値以外の定数などを設定でき、後述のパラメータ補正部１６により導出される抵抗係数（真の抵抗係数CD,real）を設定することできる。
【００３８】
パラメータ補正部１６は、数値解析部１４の計算ステップ毎に、次の計算ステップで使用する入力パラメータを補正する。この補正は、例えば、個体の大きさに関連するパラメータ（例えば抵抗係数や密度等）に対して行われる。以下、これらについて順に説明する。
【００３９】
先ず、抵抗係数の補正について、具体的に説明する。抵抗係数として、WenYuモデル（非特許文献３）を用いると、WenYuの標準の抵抗係数CDは式（３"）となる。
【００４０】
【数７】
ここで、rcは連続相の真体積分率、Repは分散相のレイノルズ係数である。
連続相の真体積分率rc,realは、分散相の真の体積分率rp,realと式（４）の関係であるから、これを用いて分散相のレイノルズ係数を式（５）のように表す。
【００４１】
【数８】
ここで、Vrelは連続相と分散相の相対速度、μcは連続相の粘性係数である。これらから、Wen Yuの抵抗係数を分散相と連続相の真の値を用いて、式（６）で与えると、連続相を繰り込む前の真の粒径に対する抵抗係数で計算することができる。
【００４２】
【数９】
なお、式（６）は分散相の密度が真の分散相の密度である場合の式であるため、補正を必要とする。この補正は、「見かけの分散相が取り込んだ連続相」と「見かけの分散相の周囲に存在する連続相」の「密度の違い」対する重力の効果であるから、式の意味に従えば、連続相の密度が定数である限り、補正の必要はない。しかしながら、式（６）は実験式であり、実際の抵抗係数に関する物理現象を表現しているものではないため、誤差を含む。このため、実際には実験で得られる沈降速度と本モデルによる計算結果が合うように、式（６）での連続相の体積分率の指数部パラメータ「−１．６５」を調整すると良い。
【００４３】
次に、分散相の密度の補正について説明する。分散相は従来の方法では真の分散相のみ単成分であるため、分散相の密度ρpは分散相の物性値そのものである。しかし、連続相の一部を取り込んだ状態の分散相は２成分からなり、見かけの分散相に対する平均密度ρpの空間分布は、連続相の密度ρWと真の分散相の密度ρp,realと最大充填率εmaxから式（７）で求める値を用いて補正する。
【００４４】
【数１０】
尚、分散相の密度の補正については、式（７）から明らかなように、本実施態様１のように最大充填率εmaxが一定であるとき（尚、後述の実施態様２では、最大充填率εmaxが変化される）、補正値も一定である。従って、分散相の密度の補正は、計算ステップ毎に実行される必要はなく、分散相の密度の補正は、入力パラメータ生成部１０にて実現されてもよい。実際に、アンシス社製の流体ソルバーCFXなどを用いる数値解析部１４では、分散相が２成分からなることを入力パラメータ生成部１０を介して設定するのみで、計算セル内での分散相の平均の密度ρpについては、数値解析部１４内で個々の計算セルでの各成分の割合（例えば質量分率）から都度計算される。このため、ユーザにとっては、平均密度ρpの空間分布は、一般的には、各々の計算時間ステップでの計算結果として認識される。
【００４５】
尚、この分散相が多成分からなることを設定できる機能は、通常、固溶体や擬似的な化合物などを多成分からなる混合物（Mixture）として取り扱うことで、たとえば、2種類の分散相の間での化学反応などに伴う物質移動を取り扱うために用意されている。本実施態様を実施するためには、「多成分からなる混合物を扱うことができ、連続相に指定した物質を、分散相の成分の１つとして設定すること」が可能であれば付加的なプログラムを用意する必要はなく、実施例で用いたアンシス社製の流体ソルバーCFX以外でも、利用することができる。
【００４６】
以上説明したように、連続相を繰り込まないときの真の分散相の粒径に対する抵抗係数で計算することにより、数値解析部１４では分散相の見かけの体積分率で計算しているにもかかわらず、粒子としての振る舞いは連続相を繰り込まないとき真の分散相として振舞うという状態を実現することができる。
【００４７】
ところで、アンシス社製の流体ソルバーCFXなどを用いる数値解析部１４が直接計算するのは、見かけの分散相の体積分率であるが、既に述べたようにセル内での平均密度の空間分布などを都度計算するため、分散相中の真の分散相の質量分率MFpと、分散相中の水の質量分率MFWも解いている。このため、計算結果として利用できるこれらの値から以下の式によって、次の式により、流体全体に占める真の分散相の体積分率を求めることができる。
【００４８】
【数１１】
ここで、ρp,realは真の分散相の密度、ρWは水の密度である。
【００４９】
図３は、本発明の実施態様１による数値計算方法の要部の一例の流れを示すフローチャートである。即ち、図３は、連続相を分散相に一部繰り込むことにより分散相の充填率が上限値を持つようにモデル化する計算手順を示す。
【００５０】
ステップ３００では、入力パラメータ生成部１０において、粒子の最大充填率εmaxが設定される。本実施態様１では、最大充填率εmaxは、一定であり、例えばユーザからの入力データに基づいて設定される。
【００５１】
ステップ３０２では、繰り込みモデル化部１２において、上記ステップ３００で設定された最大充填率εmaxに基づいて、体積分率の繰り込みが設定される。分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分をαとすると、εmax＝１−αである。従って、（連続相の体積分率）＝（真の連続相の体積分率）−αとして設定される。また、分散相は、真の分散相と、繰り込み分αの連続相との２成分からなることが設定される。
【００５２】
ステップ３０４では、繰り込みモデル化部１２において、上記ステップ３００で設定された最大充填率εmaxに基づいて、上記の(2)式に従い、（分散相の粒子径）＝（真の分散相の粒子径）×εmax（−１／３）と設定される。
【００５３】
ステップ３０６では、入力パラメータ生成部１０において設定されたその他の初期条件に従って、数値解析部１４において計算領域全体の速度と圧力が仮定される。
【００５４】
ステップ３０８では、数値解析部１４（流体ソルバー）により場所（計算セル）毎に２相間の相対速度が算出される。
【００５５】
ステップ３１０では、パラメータ補正部１６において、上記ステップ３０８で場所毎に算出された相対速度に基づいて、場所毎に真の分散相の抵抗係数CD,realが算出される（上記式（６）参照）。このようにして算出された抵抗係数CD,realは、数値解析部１４（流体ソルバー）に提供される。
【００５６】
ステップ３１２では、数値解析部１４において、上記ステップ３１０で場所毎に算出された真の分散相の抵抗係数CD,realに基づいて、連続相と分散相との間の相間の運動量交換量が算出される（式（３）'参照）。
【００５７】
ステップ３１４では、数値解析部１４において、上記ステップ３１２の算出結果に基づいて、計算領域全体の速度と圧力が変更される。
【００５８】
ステップ３１６では、数値解析部１４において、時間刻みに従って計算が実行される。
【００５９】
ステップ３１８では、数値解析部１４において、計算終了条件が判断され、計算終了条件が満たされた場合は、終了し、計算終了条件が満たされない場合は、ステップ３０８に戻り、計算終了条件が満たされるまで、ステップ３０８から３１６の処理が時間ステップ毎に繰り返される。
［実施態様２］
図４は、本発明の実施態様２による数値計算装置２の要部構成の一例を示す図である。
【００６０】
本実施態様２の数値計算装置２は、図３に示すように、入力パラメータ生成部１０と、繰り込みモデル化部１２と、数値解析部１４と、パラメータ補正部１６と、圧密反映部１８とを含む。
【００６１】
実施態様２では、主に圧密反映部１８が追加された点が上述の実施態様１と異なる。以下では、実施態様２に特有の構成についてのみ説明し、他の構成については上述の実施態様１と同様であってよい。
【００６２】
圧密反映部１８は、「分散相の繰り込み分として、圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻す」というモデル化を行う。上述の実施態様１では、式（２）の最大充填率（分散相中の連続相の体積割合）εmaxを定数としていたが、この値は、例えば堆積した沈殿物では、その上部より下部のほうが大きな値となる。これは「圧密」と呼ばれ、堆積物自身の重みにより堆積物の隙間の水がしみ出すことにより、堆積物の「嵩」が低下する現象である。
【００６３】
実施態様２においては、このような「圧密」の現象を、例えばεmaxを静水圧Pの関数として与えることとしてもよい。この場合、「圧密」の現象を、容易に計算することが出来る。
【００６４】
具体的には、式（８）により、各場所でVRatio（見かけの分散相に占める真の分散相の体積分率）の値が求まるので、この値が式（９）のVRatio,specとなるように、分散相中の連続相と周りの連続相との間で物質移動を行うように設定する。
【００６５】
【数１２】
式（９）で、Pは圧力であるが、分散相の自重のみによる圧力の関数とするため、連続相の密度を基準とした静水圧とする。静水圧の例である図５（A）において、高さ約0.03[m]以上で、値が０となっているが、これは、分散相の自重のみによる圧力を示しているためである。分散相の自重のみによる圧力とは、あたかも、連続相が空気のように密度が０に近い状態での流体全体（分散相＋連続相）の静水圧であり、通常の静水圧から連続相のみの静水圧を差し引いた値である。Cprsはモデルパラメータであり、例えば図５（A）に示す圧力の時に図５（B）に示す充填率とするためには、0.035となる。HTotalは、計算モデルの高さであり（この例では0.16[m]）、VRatio,iniは圧密前の最大充填率、VFp,iniは圧密を考慮しない場合の沈殿物の嵩位置（見かけ）の体積分率であり（この例では0.15）、ρp,realは真の分散相の密度であり（この例では1608[kg/m3]）、ρWは連続相の密度であり（この例では997[kg/m3]）、gは重力加速度である。
【００６６】
図６は、本発明の実施態様２による数値計算方法の要部の一例の流れを示すフローチャートである。即ち、図６は、本発明の連続相を分散相の充填率が上限値を持つようにモデル化する場合に、圧密の影響による分散相の嵩の減少をさらに考慮する場合の手順を示す。
【００６７】
ステップ６００では、入力パラメータ生成部１０において、粒子の最大充填率εmaxの初期値が設定される。本実施態様１では、最大充填率εmaxは、変化するが、その初期値は、例えばユーザからの入力データに基づいて設定される。
【００６８】
ステップ６０２では、繰り込みモデル化部１２において、上記ステップ６００で設定された最大充填率εmaxに基づいて、体積分率の繰り込みが設定される。分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分をαとすると、εmax＝１−αである。従って、（連続相の体積分率）＝（真の連続相の体積分率）−αとして設定される。また、分散相は、真の分散相と、繰り込み分αの連続相との２成分からなることが設定される。
【００６９】
ステップ６０４では、繰り込みモデル化部１２において、上記ステップ６００で設定された最大充填率εmaxに基づいて、上記の(2)式に従い、（分散相の粒子径）＝（真の分散相の粒子径）×εmax（−１／３）と設定される。
【００７０】
ステップ６０６では、入力パラメータ生成部１０において設定されたその他の初期条件に従って、数値解析部１４において計算領域全体の速度と圧力が仮定される。
【００７１】
ステップ６０８では、数値解析部１４（流体ソルバー）により場所（計算セル）毎に２相間の相対速度が算出される。
【００７２】
ステップ６１０では、パラメータ補正部１６において、上記ステップ６０８で場所毎に算出された相対速度に基づいて、場所毎に真の分散相の抵抗係数CD,realが算出される（上記式（６）参照）。このようにして算出された抵抗係数CD,realは、数値解析部１４（流体ソルバー）に提供される。
【００７３】
ステップ６１２では、圧密反映部１８において、充填率εmaxが各場所での圧力の関数として算出され、数値解析部１４（流体ソルバー）に提供される。
【００７４】
ステップ６１４では、数値解析部１４において、上記ステップ６１２で算出された充填率εmaxに基づいて繰り込み分を物質移動により調整される（上記式（９）参照）。
【００７５】
ステップ６１６では、数値解析部１４において、上記ステップ６１０で場所毎に算出された真の分散相の抵抗係数CD,realに基づいて、連続相と分散相との間の相間の運動量交換量が算出される（式（３）'参照）。
【００７６】
ステップ６１８では、数値解析部１４において、上記ステップ６１４、６１６の算出結果に基づいて、計算領域全体の速度と圧力が変更される。
【００７７】
ステップ６２０では、数値解析部１４において、時間刻みに従って計算が実行される。
【００７８】
ステップ６２２では、数値解析部１４において、計算終了条件が判断され、計算終了条件が満たされた場合は、終了し、計算終了条件が満たされない場合は、ステップ６０８に戻り、計算終了条件が満たされるまで、ステップ６０８から６２０の処理が時間ステップ毎に繰り返される。
【００７９】
次に、本願発明者が行った実施例とそれに対する比較例について説明する。
【実施例１】
【００８０】
実施例１は、上述の実施態様１に対応した実施例である。
【００８１】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、最大充填率を0.76に設定して本発明の繰り込みを用い(Solid pressure項を使用しない)アンシス社製の流体ソルバーCFXを用いて計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、初期状態で上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ50[ms]とし、200秒まで計算した。
【００８２】
体積分率最大値の時間変化の一例を図６に示すが、沈殿物が計算モデルの下部に堆積し始めた10秒以降計算終了まで、計算が発散することはなかった。この図で、真の分散相の体積分率とは、上述の式（８）によって、ソルバーの通常の計算結果から求めたものである。
【００８３】
図８に、200秒後の計算結果の例として、計算モデル断面での密度を示す。この図で計算モデル下部の分散相が堆積した部分の密度が1461[kg/m3]であり、最大充填率0.76の設定が反映されていることが分かる。また、図９に200秒後の計算結果の例として、個々の計算セル中での真の分散相（密度1608[kg/m3]の粒子）の体積分率の分布を示す。この値は、図６と同様に式（８）によってソルバーの計算結果の変数から求めたものである。
【００８４】
以上、示したように、本発明を用いることにより、従来のSolid Pressure項を用いずに、分散相が沈降して計算モデルの下部に堆積する状況を安定して計算できることが分かる。
【実施例２】
【００８５】
実施例２は、上述の実施態様２に対応した実施例である。
【００８６】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、圧密の影響を考慮する前の最大充填率を0.76に設定して本発明の繰り込みを用い(Solid pressure項を使用しない)、圧密の影響により分散相中に繰り込まれた水が計算モデル中の沈殿物の最下部で約0.79となる分だけ粒子分の圧力に応じて連続相に戻す設定とし、アンシス社製の流体ソルバーCFXを用いて計算した。具体的には、最大充填率を真の分散相の量に相当する圧力（基準の密度を連続相として計算した静水圧）の関数として設定した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、初期状態で上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ100[ms]とし、200秒まで計算した。なお、この計算では、計算の安定性を確認するため圧力に対する分散相の体積分率変化を大きな値とした。
【００８７】
体積分率最大値の時間変化の一例を図１０に示すが、沈殿物が計算モデルの下部に堆積し始めた10秒以降計算終了まで、計算が発散することはなかった。
【００８８】
図１１に、200秒後の計算結果の例として、計算モデル断面での密度を示す。この図で計算モデル下部の分散相が堆積した部分の密度が1481[kg/m3]であり、分散相中に繰り込まれた水が連続相に一部戻されていることが分かる。また、図１２に200秒後の計算結果の例として、見かけの分散相中に対する真の分散相（密度1608[kg/m3]の粒子）の体積分率の分布を示す。この値は、式（８）のVratioの値をソルバーの計算結果の変数から求めて表示したものである。この図から、最大充填率は圧密の影響を考慮することにより、0.76から約0.79に上昇し、圧力の高い下部ほど高い充填率になっていることがわかる。
【００８９】
以上、示したように、本発明を用いることにより、分散相が沈降して計算モデルの下部に堆積し、圧密の影響により分散相中の水が連続相に染み出す状況を安定して計算できることが分かる。
【比較例１】
【００９０】
比較例１は、Solid Pressureを用いた実施例である。図１３は、この比較例１に対応する計算手順、即ち、Solid Pressureを用いて分散相の充填率が上限値を持つようにモデル化する計算手順を示す。
【００９１】
比較例１では、水を連続相とし、密度が2300[kg/m3]の粒子を分散相とし、アンシス社製の流体ソルバーCFXを用いてSolid Pressureで最大充填率を0.76に設定して計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、実施例１と同様の設定で、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ100[ms]として計算した。
【００９２】
図１４に示すように、計算ステップの途中までは計算できたが、１００秒を越えたあたりから最大充填率が設定した0.76を越え、途中で発散して計算が停止したために140秒(1400ステップ)以降の計算結果は得られなかった。
【比較例２】
【００９３】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、アンシス社製の流体ソルバーCFXを用いてSolid Pressureで最大充填率を0.76に設定して計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、実施例１と同一の初期状態（上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０）に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ50[ms]として計算した。
【００９４】
図１５に示すように、計算ステップの途中までは計算できたが、最大充填率は、１３０秒を越えたあたりから設定した0.76を越え、１６０秒あたりから不安定に推移するとともに、発散して計算が停止したために170秒(3400ステップ)以降の計算結果が得られなかった。
【００９５】
以上のように、Solid Pressureを用いた比較例１，２では充填率の上限値を運動量保存式中の圧力勾配という間接的な物理量でモデル化しているのに比べ、本発明による実施例１，２では、体積分率の繰り込みを用いているため、より直接的に、分散相の充填率上限値を設定している。かかる点が、本発明による実施例１，２において安定して計算できる理由の一つである。
【００９６】
以上のように本発明の実施例によれば、分散相に連続相の一部を繰り込むために、分散相の体積分率は真の粒子隙間の体積分率との合計で表すことができ、分散相の体積分率をみかけの体積分率または嵩の値として計算することができる。すなわち、流体計算ソルバーで計算する見かけの体積分率が１．０という状態は、真の体積分率が設定した充填率の最大値となった状態を表している。さらに、体積分率が１．０という状態は、数値的にも安定である。したがって、本発明の実施例によれば、数値的に不安定なSolid pressure項を用いて粒子相が一定の充填率以上にならない状態をモデル化する比較例に比べ、格段に安定した計算が可能となる。
【００９７】
以上、本発明の好ましい実施例について詳説したが、本発明は、上述した実施例に制限されることはなく、本発明の範囲を逸脱することなく、上述した実施例に種々の変形及び置換を加えることができる。
【符号の説明】
【００９８】
１，２ 数値計算装置
１０ 入力パラメータ生成部
１２ 繰り込みモデル化部
１４ 数値解析部
１６ パラメータ補正部
１８ 圧密反映部
【技術分野】
【０００１】
本発明は、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算方法及び数値流体計算装置に関する。
【背景技術】
【０００２】
粒子を入れ物に詰めた場合、粒子1個の体積と個数の積である粒子の真の体積と、粒子間の隙間を含め粒子が存在する空間の体積（見かけの体積または嵩）とは異なり、球形粒子では六方最密充填構造の時に充填率が７４％であり、この最大充填率以上にならないことが一般的に知られている。
【０００３】
一方、流体計算の粒子相の運動量輸送方程式には付加的な項として、Solid Pressure項があり、非特許文献１および２などに記載されているようにGidaspowらは、粒子相の空隙率αSの関数として表現している。すなわち、
【０００４】
【数１】
ここで、G0は代表弾性率、cは圧縮係数、αsmは粒子相の最大空隙率である。この関連技術によれば、運動量の連続の式で粒子相の空隙率すなわち、粒子の充填率が設定した値に近づくにつれて大きな圧力勾配が生じ、粒子の充填率が設定した値に近づくような粒子相の流れを阻害することで、粒子が設定した充填率以上にならない状態を間接的にモデル化している。この圧力勾配の項は、あたかも粒子同士が反発力を持つようであり、Solid pressure項と呼ばれる。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Gidaspow, D., Multiphase Flow and Fluidization , Academic Press, (1994).
【非特許文献２】Bouillard, J. X., Lyczcowski, R. W. and Gidaspow, D., AIChE, J. , 35, 908(1996).
【非特許文献３】Wen, C.Y. and Yu, Y.H. (1996).CEP. Symp. Series, 62, 100-111.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上記記載した方法では、計算した充填率が、設定した最大充填率の値に近くなったときに、Solid pressure項の値を指数関数的に急激に大きくする必要がある。このため、数値的な安定性に問題があった。具体的には短い時間ステップで計算する必要があるとともに、最大充填率に達した計算セルを含む粒子の運動を計算する場合に、粒子相の運動量輸送方程式に付加したSolid pressure項の値の計算で、数値誤差が発生するなどという不具合があった。この影響に起因して、設定した最大充填率を超える計算結果となったり、計算値が発散して計算できないという欠点があった。
【０００７】
そこで、本発明は、Solid pressure項に関わる問題点を解消して分散相と連続相の混相流を精度良く解析することができる数値流体計算方法及び数値流体計算装置の提供を目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するため、本発明の一局面によれば、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算方法において、
分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする数値流体計算方法が提供される。
【０００９】
本発明のその他の一局面によれば、連続体モデルを用いて分散相と連続相の混相流を解析する数値流体計算装置であって、
分散相の体積分率に連続相の体積分率の一部が繰り込まれるように、分散相と連続相をモデル化するモデル化手段と、
前記モデル化された分散相と連続相の混相流を解析する解析手段とを備えることを特徴とする、数値流体計算装置が提供される。
【発明の効果】
【００１０】
本発明によれば、Solid pressure項に関わる問題点を解消して分散相と連続相の混相流を精度良く解析することができる数値流体計算方法及び数値流体計算装置が得られる。
【図面の簡単な説明】
【００１１】
【図１】本発明による数値計算方法の基本概念の説明図である。
【図２】本発明の実施態様１による数値計算装置１の要部構成の一例を示す図である。
【図３】本発明の実施態様１による数値計算方法の要部の一例の流れを示すフローチャートである。
【図４】本発明の実施態様２による数値計算装置２の要部構成の一例を示す図である。
【図５】図５（A）は、沈殿物底からの高さと圧力の関係を示すグラフであり、図５（B）は、沈殿物底からの高さと設定充填率の関係を示すグラフである。
【図６】本発明の実施態様２による数値計算方法の要部の一例の流れを示すフローチャートである。
【図７】実施例１による体積分率最大値の時間変化を示す図である。
【図８】実施例１による計算モデル断面での密度を示す図である。
【図９】実施例１による個々の計算セル中での真の分散相の体積分率の分布を示す図である。
【図１０】実施例２による体積分率最大値の時間変化を示す図である。
【図１１】実施例２による計算モデル断面での密度を示す図である。
【図１２】実施例２による個々の計算セル中での真の分散相の体積分率の分布を示す図である。
【図１３】Solid Pressureを用いて分散相の充填率が上限値を持つようにモデル化する計算手順（比較例）を示すフローチャートである。
【図１４】比較例１による計算結果を示す図である。
【図１５】比較例２による計算結果を示す図である。
【発明を実施するための形態】
【００１２】
以下、図面を参照して、本発明を実施するための最良の形態の説明を行う。
【００１３】
本発明による数値計算方法は、連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法に関する。
【００１４】
ここで取り扱う流体の運動方程式は、例えば次のような式で表される。
【００１５】
【数２】
これを、別の書き方で表現すると、次の通りである。
【００１６】
【数３】
この式で、Kは外力（重力加速度ｇ等）であり、ρは密度であり、ｐは圧力分布である。
【００１７】
Eulerian-Euleria法では、固体粒子も粒ではなく連続体として考えてあたかも２つの流体の流れを計算するかのように、多相流の流れ場を計算する手法である。尚、ANSYS社の汎用流体計算ソフト「CFX」等においても、固体粒子を連続体として扱っており、かかる手法が一般的である。
【００１８】
従って、Eulerian-Euleria法では、連続相と分散相のそれぞれの相に対して、数２の運動方程式が適用される。この際、各相の運動方程式おける外力の項Ｋには、互いの相から受ける抵抗（相間の運動量交換量）が含まれることになる。
【００１９】
ところで、上述の如く、粒子を入れ物に詰めた場合、粒子1個の体積と個数の積である粒子の真の体積と、粒子間の隙間を含め粒子が存在する空間の体積（見かけの体積または嵩）とは異なり、例えば球形粒子では六方最密充填構造の時に充填率が７４％であり、この最大充填率以上にならないことが一般的に知られている（図１（Ａ）及び図１（Ｂ）参照）。
【００２０】
このような粒子の充填特徴を考慮するために、本発明による数値計算方法は、従来使用されてきたSolid pressure項を使用せずに、それに代えて、分散相の体積分率に連続相の体積分率の一部を繰り込むことを特徴とする。
【００２１】
図１は、このような本発明による数値計算方法の基本概念の説明図である。
【００２２】
図１（Ａ）は、球形の粒子が液体底に沈殿している様子を模式的に示す図である。図１（Ｂ）は、図１（Ａ）の粒子の実際の体積分率を模式的に示す図である。図１（Ｃ）は、分散相の体積分率に連続相の体積分率の一部を繰り込む方法を概念的に示す図である。
【００２３】
概念的には固体は球形の粒であるため、図１（Ａ）に示すように、固体を完全に充填することはできず、固体間の隙間に液体が存在する。従って、図１（Ｂ）に示すように、分散相の体積分率は、嵩よりも小さい。
【００２４】
そこで、本発明では、図１（Ｃ）に示すように、分散相の体積分率に連続相の体積分率の一部を繰り込むことで、図１（Ａ）に示すような固体充填時の状態をモデル化する。即ち、図１（Ｃ）に示すように、固体成分の内訳と嵩との差を考慮して、球形の固体の隙間の液体を分散相に繰り込む。概念的には、分散相を球形の粒子ではなく、ブロック状の粒子と考え、球形であった場合に隙間に存在する連続相の液体（図１（Ｃ）の符号Ａで示す部分）を、分散相の一部として繰り込むというモデルである。尚、このモデルでは、ブロック状の粒子間には、Solid pressure項を付加することで発生するような反発力は生じない。
【００２５】
このように本発明によれば、分散相の体積分率に連続相の体積分率の一部を繰り込むことで、Solid pressure項を使用せずに、固体成分の内訳と嵩との差を補償することができる。これにより、Solid pressure項の問題点を解消して分散相と連続相の混相流を精度良く解析することができる。
【００２６】
尚、図１（Ｃ）を参照したモデルは、固体の大きさは、隙間の液体分だけ大きくした繰り込み固体の大きさで計算するため、個体の大きさに関連するパラメータ（例えば抵抗係数等）は、好ましくは、繰り込み前の元の大きさを基準にして（換算して）計算する。この換算方法については後述する。
【００２７】
次に、以上の本発明の基本概念を実現するための好ましい実施態様について説明する。
［実施態様１］
図２は、本発明の実施態様１による数値計算装置１の要部構成の一例を示す図である。
【００２８】
本実施態様１の数値計算装置１は、図２に示すように、入力パラメータ生成部１０と、繰り込みモデル化部１２と、数値解析部１４と、パラメータ補正部１６とを含む。尚、これらの各部１０，１２，１４，１６は、ソフトウェアとハードウェアの組み合わせにより実現される。この際、これらの各部１０，１２，１４，１６は、共通のソフトウェアで実現されてもよいし、複数のソフトウェアにより協動して実現されてもよい。また、同様に、これらの各部１０，１２，１４，１６は、共通のハードウェアで実現されてもよいし、複数のハードウェアにより協動して実現されてもよい。
【００２９】
入力パラメータ生成部１０は、外部からの入力データ（典型的には、ユーザからの入力データ）に基づいて、数値解析部１４における解析に必要な入力パラメータ（解析条件）を生成（設定）する。入力パラメータは、例えば連続相の特性（水といった媒体の特定を含む）、分散相の粒子の密度、最大充填率等を含む。
【００３０】
繰り込みモデル化部１２は、図１（Ｃ）を参照して説明したモデル化を行う。即ち、繰り込みモデル化部１２は、入力パラメータ生成部１０と協動して、分散相の体積分率に連続相の体積分率の一部が繰り込まれるような入力パラメータを生成する。
【００３１】
具体的には、繰り込みモデル化部１２は、連続相の一部を分散相に繰り込むために、分散相を粒子モデルとし、粒子径として真の粒子径Dp,realと見かけの粒子径Dpとの関係を(2)式で定義する。
【００３２】
【数４】
ここで、εmaxは、最大充填率である。尚、繰り込みモデル化部１２の処理の詳細については後述する。
【００３３】
数値解析部１４は、入力パラメータ生成部１０で生成された入力パラメータと、後述のパラメータ補正部１６から供給されるパラメータとを用いて、連続体モデル（Eulerian-Euleria法）による分散相と連続相の混相流を解析する。数値解析部１４における流体解析方法は、連続体モデルを用いるものであればその詳細は任意である。例えば、数値解析部１４は、流体ソルバーと称される汎用流体計算ソフトとＣＰＵ等のハードウェアにより実現されてもよい。かかる汎用流体計算ソフトとしては、ANSYS社の汎用流体計算ソフトCFXが知られている。数値解析部１４は、遠隔位置にあるハードウェア資源（例えばスーパーコンピューター）と協動して実現されてもよい。
【００３４】
ここで、数値解析部１４において考慮される相間の運動量交換量について説明する。
【００３５】
相間の運動量交換量は、単位体積当たりの抵抗Ｄd,cが、１つの分散相粒子が運動するときに周りの連続相から受ける抵抗Ｄsと単位体積当りの粒子数ｎpの積であるから、式（３）のようになる。
【００３６】
【数５】
ここで、CDは抵抗係数、ρCは連続相の密度、Ud及びUcはそれぞれ分散相と連続相の速度、rpは分散相の体積分率である。ここで、分散相の真の体積分率をrp,realとし、分散相に式（２）の繰り込みを行うと以下の式（３'）に示すように、抵抗係数CDを、真の抵抗係数CD,realで置き換えたのみであり、同じ形の式に書くことができる。
【００３７】
【数６】
すなわち、真の抵抗係数CD,realを用いることにより、数値解析部１４内は標準的な式（３）に従った計算を行うが、式（３'）のように、実質的には真の粒子径Dp,realと分散相の真の体積分率rp,realを用いた計算結果を得ることが出来る。
なお、通常の流体ソルバーを用いる数値解析部１４では、相間の運動量交換に関する抵抗係数として標準的な値以外の定数などを設定でき、後述のパラメータ補正部１６により導出される抵抗係数（真の抵抗係数CD,real）を設定することできる。
【００３８】
パラメータ補正部１６は、数値解析部１４の計算ステップ毎に、次の計算ステップで使用する入力パラメータを補正する。この補正は、例えば、個体の大きさに関連するパラメータ（例えば抵抗係数や密度等）に対して行われる。以下、これらについて順に説明する。
【００３９】
先ず、抵抗係数の補正について、具体的に説明する。抵抗係数として、WenYuモデル（非特許文献３）を用いると、WenYuの標準の抵抗係数CDは式（３"）となる。
【００４０】
【数７】
ここで、rcは連続相の真体積分率、Repは分散相のレイノルズ係数である。
連続相の真体積分率rc,realは、分散相の真の体積分率rp,realと式（４）の関係であるから、これを用いて分散相のレイノルズ係数を式（５）のように表す。
【００４１】
【数８】
ここで、Vrelは連続相と分散相の相対速度、μcは連続相の粘性係数である。これらから、Wen Yuの抵抗係数を分散相と連続相の真の値を用いて、式（６）で与えると、連続相を繰り込む前の真の粒径に対する抵抗係数で計算することができる。
【００４２】
【数９】
なお、式（６）は分散相の密度が真の分散相の密度である場合の式であるため、補正を必要とする。この補正は、「見かけの分散相が取り込んだ連続相」と「見かけの分散相の周囲に存在する連続相」の「密度の違い」対する重力の効果であるから、式の意味に従えば、連続相の密度が定数である限り、補正の必要はない。しかしながら、式（６）は実験式であり、実際の抵抗係数に関する物理現象を表現しているものではないため、誤差を含む。このため、実際には実験で得られる沈降速度と本モデルによる計算結果が合うように、式（６）での連続相の体積分率の指数部パラメータ「−１．６５」を調整すると良い。
【００４３】
次に、分散相の密度の補正について説明する。分散相は従来の方法では真の分散相のみ単成分であるため、分散相の密度ρpは分散相の物性値そのものである。しかし、連続相の一部を取り込んだ状態の分散相は２成分からなり、見かけの分散相に対する平均密度ρpの空間分布は、連続相の密度ρWと真の分散相の密度ρp,realと最大充填率εmaxから式（７）で求める値を用いて補正する。
【００４４】
【数１０】
尚、分散相の密度の補正については、式（７）から明らかなように、本実施態様１のように最大充填率εmaxが一定であるとき（尚、後述の実施態様２では、最大充填率εmaxが変化される）、補正値も一定である。従って、分散相の密度の補正は、計算ステップ毎に実行される必要はなく、分散相の密度の補正は、入力パラメータ生成部１０にて実現されてもよい。実際に、アンシス社製の流体ソルバーCFXなどを用いる数値解析部１４では、分散相が２成分からなることを入力パラメータ生成部１０を介して設定するのみで、計算セル内での分散相の平均の密度ρpについては、数値解析部１４内で個々の計算セルでの各成分の割合（例えば質量分率）から都度計算される。このため、ユーザにとっては、平均密度ρpの空間分布は、一般的には、各々の計算時間ステップでの計算結果として認識される。
【００４５】
尚、この分散相が多成分からなることを設定できる機能は、通常、固溶体や擬似的な化合物などを多成分からなる混合物（Mixture）として取り扱うことで、たとえば、2種類の分散相の間での化学反応などに伴う物質移動を取り扱うために用意されている。本実施態様を実施するためには、「多成分からなる混合物を扱うことができ、連続相に指定した物質を、分散相の成分の１つとして設定すること」が可能であれば付加的なプログラムを用意する必要はなく、実施例で用いたアンシス社製の流体ソルバーCFX以外でも、利用することができる。
【００４６】
以上説明したように、連続相を繰り込まないときの真の分散相の粒径に対する抵抗係数で計算することにより、数値解析部１４では分散相の見かけの体積分率で計算しているにもかかわらず、粒子としての振る舞いは連続相を繰り込まないとき真の分散相として振舞うという状態を実現することができる。
【００４７】
ところで、アンシス社製の流体ソルバーCFXなどを用いる数値解析部１４が直接計算するのは、見かけの分散相の体積分率であるが、既に述べたようにセル内での平均密度の空間分布などを都度計算するため、分散相中の真の分散相の質量分率MFpと、分散相中の水の質量分率MFWも解いている。このため、計算結果として利用できるこれらの値から以下の式によって、次の式により、流体全体に占める真の分散相の体積分率を求めることができる。
【００４８】
【数１１】
ここで、ρp,realは真の分散相の密度、ρWは水の密度である。
【００４９】
図３は、本発明の実施態様１による数値計算方法の要部の一例の流れを示すフローチャートである。即ち、図３は、連続相を分散相に一部繰り込むことにより分散相の充填率が上限値を持つようにモデル化する計算手順を示す。
【００５０】
ステップ３００では、入力パラメータ生成部１０において、粒子の最大充填率εmaxが設定される。本実施態様１では、最大充填率εmaxは、一定であり、例えばユーザからの入力データに基づいて設定される。
【００５１】
ステップ３０２では、繰り込みモデル化部１２において、上記ステップ３００で設定された最大充填率εmaxに基づいて、体積分率の繰り込みが設定される。分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分をαとすると、εmax＝１−αである。従って、（連続相の体積分率）＝（真の連続相の体積分率）−αとして設定される。また、分散相は、真の分散相と、繰り込み分αの連続相との２成分からなることが設定される。
【００５２】
ステップ３０４では、繰り込みモデル化部１２において、上記ステップ３００で設定された最大充填率εmaxに基づいて、上記の(2)式に従い、（分散相の粒子径）＝（真の分散相の粒子径）×εmax（−１／３）と設定される。
【００５３】
ステップ３０６では、入力パラメータ生成部１０において設定されたその他の初期条件に従って、数値解析部１４において計算領域全体の速度と圧力が仮定される。
【００５４】
ステップ３０８では、数値解析部１４（流体ソルバー）により場所（計算セル）毎に２相間の相対速度が算出される。
【００５５】
ステップ３１０では、パラメータ補正部１６において、上記ステップ３０８で場所毎に算出された相対速度に基づいて、場所毎に真の分散相の抵抗係数CD,realが算出される（上記式（６）参照）。このようにして算出された抵抗係数CD,realは、数値解析部１４（流体ソルバー）に提供される。
【００５６】
ステップ３１２では、数値解析部１４において、上記ステップ３１０で場所毎に算出された真の分散相の抵抗係数CD,realに基づいて、連続相と分散相との間の相間の運動量交換量が算出される（式（３）'参照）。
【００５７】
ステップ３１４では、数値解析部１４において、上記ステップ３１２の算出結果に基づいて、計算領域全体の速度と圧力が変更される。
【００５８】
ステップ３１６では、数値解析部１４において、時間刻みに従って計算が実行される。
【００５９】
ステップ３１８では、数値解析部１４において、計算終了条件が判断され、計算終了条件が満たされた場合は、終了し、計算終了条件が満たされない場合は、ステップ３０８に戻り、計算終了条件が満たされるまで、ステップ３０８から３１６の処理が時間ステップ毎に繰り返される。
［実施態様２］
図４は、本発明の実施態様２による数値計算装置２の要部構成の一例を示す図である。
【００６０】
本実施態様２の数値計算装置２は、図３に示すように、入力パラメータ生成部１０と、繰り込みモデル化部１２と、数値解析部１４と、パラメータ補正部１６と、圧密反映部１８とを含む。
【００６１】
実施態様２では、主に圧密反映部１８が追加された点が上述の実施態様１と異なる。以下では、実施態様２に特有の構成についてのみ説明し、他の構成については上述の実施態様１と同様であってよい。
【００６２】
圧密反映部１８は、「分散相の繰り込み分として、圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻す」というモデル化を行う。上述の実施態様１では、式（２）の最大充填率（分散相中の連続相の体積割合）εmaxを定数としていたが、この値は、例えば堆積した沈殿物では、その上部より下部のほうが大きな値となる。これは「圧密」と呼ばれ、堆積物自身の重みにより堆積物の隙間の水がしみ出すことにより、堆積物の「嵩」が低下する現象である。
【００６３】
実施態様２においては、このような「圧密」の現象を、例えばεmaxを静水圧Pの関数として与えることとしてもよい。この場合、「圧密」の現象を、容易に計算することが出来る。
【００６４】
具体的には、式（８）により、各場所でVRatio（見かけの分散相に占める真の分散相の体積分率）の値が求まるので、この値が式（９）のVRatio,specとなるように、分散相中の連続相と周りの連続相との間で物質移動を行うように設定する。
【００６５】
【数１２】
式（９）で、Pは圧力であるが、分散相の自重のみによる圧力の関数とするため、連続相の密度を基準とした静水圧とする。静水圧の例である図５（A）において、高さ約0.03[m]以上で、値が０となっているが、これは、分散相の自重のみによる圧力を示しているためである。分散相の自重のみによる圧力とは、あたかも、連続相が空気のように密度が０に近い状態での流体全体（分散相＋連続相）の静水圧であり、通常の静水圧から連続相のみの静水圧を差し引いた値である。Cprsはモデルパラメータであり、例えば図５（A）に示す圧力の時に図５（B）に示す充填率とするためには、0.035となる。HTotalは、計算モデルの高さであり（この例では0.16[m]）、VRatio,iniは圧密前の最大充填率、VFp,iniは圧密を考慮しない場合の沈殿物の嵩位置（見かけ）の体積分率であり（この例では0.15）、ρp,realは真の分散相の密度であり（この例では1608[kg/m3]）、ρWは連続相の密度であり（この例では997[kg/m3]）、gは重力加速度である。
【００６６】
図６は、本発明の実施態様２による数値計算方法の要部の一例の流れを示すフローチャートである。即ち、図６は、本発明の連続相を分散相の充填率が上限値を持つようにモデル化する場合に、圧密の影響による分散相の嵩の減少をさらに考慮する場合の手順を示す。
【００６７】
ステップ６００では、入力パラメータ生成部１０において、粒子の最大充填率εmaxの初期値が設定される。本実施態様１では、最大充填率εmaxは、変化するが、その初期値は、例えばユーザからの入力データに基づいて設定される。
【００６８】
ステップ６０２では、繰り込みモデル化部１２において、上記ステップ６００で設定された最大充填率εmaxに基づいて、体積分率の繰り込みが設定される。分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分をαとすると、εmax＝１−αである。従って、（連続相の体積分率）＝（真の連続相の体積分率）−αとして設定される。また、分散相は、真の分散相と、繰り込み分αの連続相との２成分からなることが設定される。
【００６９】
ステップ６０４では、繰り込みモデル化部１２において、上記ステップ６００で設定された最大充填率εmaxに基づいて、上記の(2)式に従い、（分散相の粒子径）＝（真の分散相の粒子径）×εmax（−１／３）と設定される。
【００７０】
ステップ６０６では、入力パラメータ生成部１０において設定されたその他の初期条件に従って、数値解析部１４において計算領域全体の速度と圧力が仮定される。
【００７１】
ステップ６０８では、数値解析部１４（流体ソルバー）により場所（計算セル）毎に２相間の相対速度が算出される。
【００７２】
ステップ６１０では、パラメータ補正部１６において、上記ステップ６０８で場所毎に算出された相対速度に基づいて、場所毎に真の分散相の抵抗係数CD,realが算出される（上記式（６）参照）。このようにして算出された抵抗係数CD,realは、数値解析部１４（流体ソルバー）に提供される。
【００７３】
ステップ６１２では、圧密反映部１８において、充填率εmaxが各場所での圧力の関数として算出され、数値解析部１４（流体ソルバー）に提供される。
【００７４】
ステップ６１４では、数値解析部１４において、上記ステップ６１２で算出された充填率εmaxに基づいて繰り込み分を物質移動により調整される（上記式（９）参照）。
【００７５】
ステップ６１６では、数値解析部１４において、上記ステップ６１０で場所毎に算出された真の分散相の抵抗係数CD,realに基づいて、連続相と分散相との間の相間の運動量交換量が算出される（式（３）'参照）。
【００７６】
ステップ６１８では、数値解析部１４において、上記ステップ６１４、６１６の算出結果に基づいて、計算領域全体の速度と圧力が変更される。
【００７７】
ステップ６２０では、数値解析部１４において、時間刻みに従って計算が実行される。
【００７８】
ステップ６２２では、数値解析部１４において、計算終了条件が判断され、計算終了条件が満たされた場合は、終了し、計算終了条件が満たされない場合は、ステップ６０８に戻り、計算終了条件が満たされるまで、ステップ６０８から６２０の処理が時間ステップ毎に繰り返される。
【００７９】
次に、本願発明者が行った実施例とそれに対する比較例について説明する。
【実施例１】
【００８０】
実施例１は、上述の実施態様１に対応した実施例である。
【００８１】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、最大充填率を0.76に設定して本発明の繰り込みを用い(Solid pressure項を使用しない)アンシス社製の流体ソルバーCFXを用いて計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、初期状態で上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ50[ms]とし、200秒まで計算した。
【００８２】
体積分率最大値の時間変化の一例を図６に示すが、沈殿物が計算モデルの下部に堆積し始めた10秒以降計算終了まで、計算が発散することはなかった。この図で、真の分散相の体積分率とは、上述の式（８）によって、ソルバーの通常の計算結果から求めたものである。
【００８３】
図８に、200秒後の計算結果の例として、計算モデル断面での密度を示す。この図で計算モデル下部の分散相が堆積した部分の密度が1461[kg/m3]であり、最大充填率0.76の設定が反映されていることが分かる。また、図９に200秒後の計算結果の例として、個々の計算セル中での真の分散相（密度1608[kg/m3]の粒子）の体積分率の分布を示す。この値は、図６と同様に式（８）によってソルバーの計算結果の変数から求めたものである。
【００８４】
以上、示したように、本発明を用いることにより、従来のSolid Pressure項を用いずに、分散相が沈降して計算モデルの下部に堆積する状況を安定して計算できることが分かる。
【実施例２】
【００８５】
実施例２は、上述の実施態様２に対応した実施例である。
【００８６】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、圧密の影響を考慮する前の最大充填率を0.76に設定して本発明の繰り込みを用い(Solid pressure項を使用しない)、圧密の影響により分散相中に繰り込まれた水が計算モデル中の沈殿物の最下部で約0.79となる分だけ粒子分の圧力に応じて連続相に戻す設定とし、アンシス社製の流体ソルバーCFXを用いて計算した。具体的には、最大充填率を真の分散相の量に相当する圧力（基準の密度を連続相として計算した静水圧）の関数として設定した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、初期状態で上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ100[ms]とし、200秒まで計算した。なお、この計算では、計算の安定性を確認するため圧力に対する分散相の体積分率変化を大きな値とした。
【００８７】
体積分率最大値の時間変化の一例を図１０に示すが、沈殿物が計算モデルの下部に堆積し始めた10秒以降計算終了まで、計算が発散することはなかった。
【００８８】
図１１に、200秒後の計算結果の例として、計算モデル断面での密度を示す。この図で計算モデル下部の分散相が堆積した部分の密度が1481[kg/m3]であり、分散相中に繰り込まれた水が連続相に一部戻されていることが分かる。また、図１２に200秒後の計算結果の例として、見かけの分散相中に対する真の分散相（密度1608[kg/m3]の粒子）の体積分率の分布を示す。この値は、式（８）のVratioの値をソルバーの計算結果の変数から求めて表示したものである。この図から、最大充填率は圧密の影響を考慮することにより、0.76から約0.79に上昇し、圧力の高い下部ほど高い充填率になっていることがわかる。
【００８９】
以上、示したように、本発明を用いることにより、分散相が沈降して計算モデルの下部に堆積し、圧密の影響により分散相中の水が連続相に染み出す状況を安定して計算できることが分かる。
【比較例１】
【００９０】
比較例１は、Solid Pressureを用いた実施例である。図１３は、この比較例１に対応する計算手順、即ち、Solid Pressureを用いて分散相の充填率が上限値を持つようにモデル化する計算手順を示す。
【００９１】
比較例１では、水を連続相とし、密度が2300[kg/m3]の粒子を分散相とし、アンシス社製の流体ソルバーCFXを用いてSolid Pressureで最大充填率を0.76に設定して計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、実施例１と同様の設定で、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ100[ms]として計算した。
【００９２】
図１４に示すように、計算ステップの途中までは計算できたが、１００秒を越えたあたりから最大充填率が設定した0.76を越え、途中で発散して計算が停止したために140秒(1400ステップ)以降の計算結果は得られなかった。
【比較例２】
【００９３】
水を連続相とし、密度が1608[kg/m3]の粒子を分散相とし、アンシス社製の流体ソルバーCFXを用いてSolid Pressureで最大充填率を0.76に設定して計算した。計算モデルは、直径約３[cm]、体積100[ml]の円柱とし、実施例１と同一の初期状態（上半分の分散相体積分率を一様な値で0.3、下半分をほぼ０）に設定し、分散相が時間と共に計算モデルの下部に体積する様子の時間進展を計算した。時間ステップ50[ms]として計算した。
【００９４】
図１５に示すように、計算ステップの途中までは計算できたが、最大充填率は、１３０秒を越えたあたりから設定した0.76を越え、１６０秒あたりから不安定に推移するとともに、発散して計算が停止したために170秒(3400ステップ)以降の計算結果が得られなかった。
【００９５】
以上のように、Solid Pressureを用いた比較例１，２では充填率の上限値を運動量保存式中の圧力勾配という間接的な物理量でモデル化しているのに比べ、本発明による実施例１，２では、体積分率の繰り込みを用いているため、より直接的に、分散相の充填率上限値を設定している。かかる点が、本発明による実施例１，２において安定して計算できる理由の一つである。
【００９６】
以上のように本発明の実施例によれば、分散相に連続相の一部を繰り込むために、分散相の体積分率は真の粒子隙間の体積分率との合計で表すことができ、分散相の体積分率をみかけの体積分率または嵩の値として計算することができる。すなわち、流体計算ソルバーで計算する見かけの体積分率が１．０という状態は、真の体積分率が設定した充填率の最大値となった状態を表している。さらに、体積分率が１．０という状態は、数値的にも安定である。したがって、本発明の実施例によれば、数値的に不安定なSolid pressure項を用いて粒子相が一定の充填率以上にならない状態をモデル化する比較例に比べ、格段に安定した計算が可能となる。
【００９７】
以上、本発明の好ましい実施例について詳説したが、本発明は、上述した実施例に制限されることはなく、本発明の範囲を逸脱することなく、上述した実施例に種々の変形及び置換を加えることができる。
【符号の説明】
【００９８】
１，２ 数値計算装置
１０ 入力パラメータ生成部
１２ 繰り込みモデル化部
１４ 数値解析部
１６ パラメータ補正部
１８ 圧密反映部
【特許請求の範囲】
【請求項１】
連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法において、
分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする数値流体計算方法。
【請求項２】
分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分として、分散相の真の体積分率と、見かけの体積分率または嵩と、の差を用いることを特徴とする、請求項１に記載の数値流体計算方法。
【請求項３】
圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻すことを特徴とする、請求項１又は２に記載の数値流体計算方法。
【請求項４】
Solid pressure項を使用しないことを特徴とする、請求項１〜３のうちのいずれか１項に記載の数値流体計算方法。
【請求項５】
分散相の粒子の大きさに関係するパラメータは、前記連続相の体積分率の一部を繰り込む前の同粒子の元の大きさに基づいて設定する、請求項１〜４のうちのいずれか１項に記載の数値流体計算方法。
【請求項６】
前記パラメータは、抵抗係数を含む、請求項５に記載の数値流体計算方法。
【請求項７】
連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算装置であって、
分散相の体積分率に連続相の体積分率の一部が繰り込まれるように、分散相と連続相をモデル化するモデル化手段と、
前記モデル化された分散相と連続相の混相流を解析する解析手段とを備えることを特徴とする、数値流体計算装置。
【請求項８】
前記分散相の体積分率に繰り込まれる連続相の体積分率の一部は、分散相の真の体積分率と見かけの体積分率または嵩との差に基づいて設定される、請求項７に記載の数値流体計算装置。
【請求項９】
圧密の影響による分散相の嵩の減少をモデル化する圧密反映手段を更に備える、請求項７又は８に記載の数値流体計算装置。
【請求項１０】
分散相の粒子の大きさに関係するパラメータであって、前記解析手段の解析に用いるパラメータを、前記連続相の体積分率の一部を繰り込む前の元の大きさに基づいて設定するパラメータ補正手段を更に備える、請求項７〜９のうちのいずれか１項に記載の数値流体計算装置。
【請求項１１】
前記パラメータは、抵抗係数を含む、請求項１０に記載の数値流体計算装置。
【請求項１】
連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算方法において、
分散相の体積分率に連続相の体積分率の一部を繰り込むこと特徴とする数値流体計算方法。
【請求項２】
分散相の体積分率に連続相の体積分率の一部を繰り込む際の繰り込み分として、分散相の真の体積分率と、見かけの体積分率または嵩と、の差を用いることを特徴とする、請求項１に記載の数値流体計算方法。
【請求項３】
圧密の影響による分散相の嵩の減少を表現するために分散相に繰り込んだ連続相の体積分率の一部を連続相に戻すことを特徴とする、請求項１又は２に記載の数値流体計算方法。
【請求項４】
Solid pressure項を使用しないことを特徴とする、請求項１〜３のうちのいずれか１項に記載の数値流体計算方法。
【請求項５】
分散相の粒子の大きさに関係するパラメータは、前記連続相の体積分率の一部を繰り込む前の同粒子の元の大きさに基づいて設定する、請求項１〜４のうちのいずれか１項に記載の数値流体計算方法。
【請求項６】
前記パラメータは、抵抗係数を含む、請求項５に記載の数値流体計算方法。
【請求項７】
連続体モデル（Eulerian-Euleria法）を用いて分散相と連続相の混相流を解析する数値流体計算装置であって、
分散相の体積分率に連続相の体積分率の一部が繰り込まれるように、分散相と連続相をモデル化するモデル化手段と、
前記モデル化された分散相と連続相の混相流を解析する解析手段とを備えることを特徴とする、数値流体計算装置。
【請求項８】
前記分散相の体積分率に繰り込まれる連続相の体積分率の一部は、分散相の真の体積分率と見かけの体積分率または嵩との差に基づいて設定される、請求項７に記載の数値流体計算装置。
【請求項９】
圧密の影響による分散相の嵩の減少をモデル化する圧密反映手段を更に備える、請求項７又は８に記載の数値流体計算装置。
【請求項１０】
分散相の粒子の大きさに関係するパラメータであって、前記解析手段の解析に用いるパラメータを、前記連続相の体積分率の一部を繰り込む前の元の大きさに基づいて設定するパラメータ補正手段を更に備える、請求項７〜９のうちのいずれか１項に記載の数値流体計算装置。
【請求項１１】
前記パラメータは、抵抗係数を含む、請求項１０に記載の数値流体計算装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１３】
【図１４】
【図１５】
【図８】
【図９】
【図１１】
【図１２】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１３】
【図１４】
【図１５】
【図８】
【図９】
【図１１】
【図１２】
【公開番号】特開２０１１−１２３８２６（Ｐ２０１１−１２３８２６Ａ）
【公開日】平成２３年６月２３日（２０１１．６．２３）
【国際特許分類】
【出願番号】特願２００９−２８３１２３（Ｐ２００９−２８３１２３）
【出願日】平成２１年１２月１４日（２００９．１２．１４）
【出願人】（０００１８３３０３）住友金属鉱山株式会社 (2,015)
【Ｆターム（参考）】
【公開日】平成２３年６月２３日（２０１１．６．２３）
【国際特許分類】
【出願日】平成２１年１２月１４日（２００９．１２．１４）
【出願人】（０００１８３３０３）住友金属鉱山株式会社 (2,015)
【Ｆターム（参考）】
[ Back to top ]