(57)【要約】【課題】液相中に気泡が分散した混相流での現象をより高精度にシミュレーション可能なプロセス解析プログラムおよびプロセス解析装置を提供する。【解決手段】本実施の形態に従うプロセス解析は、せん断誘起乱流ＳＩによって生じる乱流応力τＬＳＩの影響の度合いを気泡誘起乱流ＢＩによって生じる乱流応力τＬＢＩに比較して小さくした乱流モデルを用いる。すなわち、以下のような乱流モデルとして、τＬＲｅ＝τＬＢＩ＋φτＬＳＩ（０＜φ＜１）を用いる。ここで、係数φの範囲は、検証実験によれば、０．０５≦φ≦０．２が好ましく、より好ましくは、０．０６≦φ≦０．１５である。さらに好ましくは、φ＝０．１である。【選択図】図４
2010040013.tif 000002
【特許請求の範囲】
【請求項１】
液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、
前記プロセス解析プログラムは、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、
前記計算モデルは、各相についての質量保存式および運動量保存式を含み、
前記状態値を時系列で計算するステップは、
液相と気泡との間の速度差、液相の体積率、気泡径に基づいて、気泡誘起乱流によって生じる乱流応力を計算するステップと、
液相の流れに伴う乱流エネルギーとエネルギー散逸とに関する輸送方程式に基づいて、せん断誘起乱流によって生じる乱流応力を計算するステップと、
気泡誘起乱流による影響に比較してせん断誘起乱流による影響の度合いを小さくして、液相乱流応力を計算するステップと、
計算された前記液相乱流応力を含む運動量保存式に基づいて、速度および圧力を計算するステップとを含む、プロセス解析プログラム。
【請求項２】
前記液相乱流応力を計算するステップでは、前記液相乱流応力τＬＲｅを、前記気泡誘起乱流によって生じる乱流応力τＬＢＩ、前記せん断誘起乱流によって生じる乱流応力τＬＳＩ、係数φ（０＜φ＜１）を用いて、
τＬＲｅ＝τＬＢＩ＋φτＬＳＩ
の式から計算する、請求項１に記載のプロセス解析プログラム。
【請求項３】
前記係数φは、０．０５≦φ≦０．２の範囲の値に設定される、請求項２に記載のプロセス解析プログラム。
【請求項４】
前記係数φは、０．１に設定される、請求項２に記載のプロセス解析プログラム。
【請求項５】
液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、
空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、
前記計算モデルは、各相についての質量保存式および運動量保存式を含み、
前記状態値を時系列で計算する手段は、
液相と気泡との間の速度差、液相の体積率、気泡径に基づいて、気泡誘起乱流によって生じる乱流応力を計算する手段と、
液相の流れに伴う乱流エネルギーとエネルギー散逸とに関する輸送方程式に基づいて、せん断誘起乱流によって生じる乱流応力を計算する手段と、
気泡誘起乱流による影響に比較してせん断誘起乱流による影響の度合いを小さくして、液相乱流応力を計算する手段と、
計算された前記液相乱流応力を含む運動量保存式に基づいて、速度および圧力を計算する手段とを含む、プロセス解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムおよびプロセス解析装置に関する。
【背景技術】
【０００２】
化学プラントの酸化、水素化、発酵、廃水処理などのさまざまなプロセスにおいて、気泡塔と呼ばれる装置が用いられている。この気泡塔は、液相中に気泡が分散した「気泡流」を利用した設備であり、他の気液反応装置に比較してより大きな気液界面積が利用でき、また構造を単純化できるため装置の管理が容易になるといった利点を有している。
【０００３】
このような気泡塔の解析では、気泡流構造が複雑であるので、実験的手法を用いて大規模な解析を行なうことは困難である。そのため、数値計算（シミュレーション）を用いて、気泡塔の解析を行なうことが注目されている。このようなシミュレーションにより気泡塔を解析する際には、対象の気泡塔についての計算モデルを高度化することが重要となっている。
【０００４】
このような状況下において、非特許文献１〜３には、非圧縮性流れ解法の枠内で各相密度の温度・圧力依存性を解析する計算モデルが提案されている。
【非特許文献１】Akio TOMIYAMA, Naoki SHIMADA, "(N+2)-FIELD MODELING FOR BUBBLY FLOW SIMULATION", Computational Fluid Dynamics JOURNAL, vol.9 no.4, January 2001, pp.418-426
【非特許文献２】島田他、「気泡塔内気泡流数値予測に必要となる相関式の感度解析」、社団法人  化学工学会、化学工学論文集、第２９巻、第６号、ｐｐ．７７８−７８５、２００３年
【非特許文献３】島田他、「化学反応・ガス吸収・熱輸送を伴う気泡塔内気泡流の数値解析」、社団法人  化学工学会、化学工学論文集、第３１巻、第６号、ｐｐ．３７７−３８７、２００５年
【非特許文献４】M.LOPEZ DE BERTONANO, R.T. LAHEY JR and O.C. JONES, "PHASE DISTRIBUTION IN BUBBLY TWO-PHASE FLOW IN VERTICAL DUCTS", Int. J. Multiphase Flow, Vol.20, No.5, pp.805-818, 1994
【非特許文献５】A.SOKOLICHIN and G. EIGENBERGER, "GAS-LIQUID FLOW BUBBLE COLUMNS AND LOOP REACTORS: PART I. DETAILED MODELLING AND NUMERICAL SIMULATION", Chemical Engineering Science, Vol.49, NO.24B, pp.5735-5746, 1994
【発明の開示】
【発明が解決しようとする課題】
【０００５】
上述のような気泡塔内での気泡を含む液相の流れをモデル化する際に重要なパラメータの１つとして、液相に生じる乱流応力（レイノルズ応力）がある。非特許文献２および４に記載されているように、この液相乱流応力は、理論的には、液相と反応器内部の壁面との間の摩擦によって生じるせん断誘起乱流と、液相と気泡との速度差によって生じる気泡誘起乱流とに依存して定まると考えることができる。
【０００６】
しかしながら、非特許文献２および５には、液相乱流応力にせん断誘起乱流を考慮すると、液相乱流応力を過大に評価してしまうという問題点が指摘されている。
【０００７】
この発明は、このような課題を解決するためになされたものであって、その目的は、液相中に気泡が分散した混相流での現象をより高精度にシミュレーション可能なプロセス解析プログラムおよびプロセス解析装置を提供することである。
【課題を解決するための手段】
【０００８】
この発明のある局面によれば、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムを提供する。プロセス解析プログラムは、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させる。計算モデルは、各相についての質量保存式および運動量保存式を含む。状態値を時系列で計算するステップは、液相と気泡との間の速度差、液相の体積率、気泡径に基づいて、気泡誘起乱流によって生じる乱流応力を計算するステップと、液相の流れに伴う乱流エネルギーとエネルギー散逸とに関する輸送方程式に基づいて、せん断誘起乱流によって生じる乱流応力を計算するステップと、気泡誘起乱流による影響に比較してせん断誘起乱流による影響の度合いを小さくして、液相乱流応力を計算するステップと、計算された液相乱流応力を含む運動量保存式に基づいて、速度および圧力を計算するステップとを含む。
【０００９】
好ましくは、液相乱流応力を計算するステップでは、液相乱流応力τＬＲｅを、気泡誘起乱流によって生じる乱流応力τＬＢＩ、せん断誘起乱流によって生じる乱流応力τＬＳＩ、係数φ（０＜φ＜１）を用いて、τＬＲｅ＝τＬＢＩ＋φτＬＳＩの式から計算する。
【００１０】
さらに好ましくは、係数φは、０．０５≦φ≦０．２の範囲の値に設定される。
また、さらに好ましくは、係数φは、０．１に設定される。
【００１１】
この発明の別の局面に従えば、液相中に気泡が分散した混相流を解析するためのプロセス解析装置を提供する。プロセス解析装置は、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を含む。計算モデルは、各相についての質量保存式および運動量保存式を含む。状態値を時系列で計算する手段は、液相と気泡との間の速度差、液相の体積率、気泡径に基づいて、気泡誘起乱流によって生じる乱流応力を計算する手段と、液相の流れに伴う乱流エネルギーとエネルギー散逸とに関する輸送方程式に基づいて、せん断誘起乱流によって生じる乱流応力を計算する手段と、気泡誘起乱流による影響に比較してせん断誘起乱流による影響の度合いを小さくして、液相乱流応力を計算する手段と、計算された液相乱流応力を含む運動量保存式に基づいて、速度および圧力を計算する手段とを含む。
【発明の効果】
【００１２】
この発明によれば、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションすることができる。
【発明を実施するための最良の形態】
【００１３】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【００１４】
＜解析対象＞
この発明の実施の形態に従うプロセス解析は、気泡塔などの液相中に気泡が分散した混相流を生じる系を解析対象とする。
【００１５】
図１は、この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【００１６】
図１を参照して、解析対象の一例として、反応器１０内に原料１（液相）を供給するとともに、原料１中に原料２（気相）が吹き込まれる系ＳＹＳを考える。反応器１０の底部には、原料１を供給するための供給ライン１２が接続されている。この供給ライン１２の経路上には、原料１を供給／遮断するための開閉弁１４が設けられている。なお、図１に示す系ＳＹＳでは、バッチ処理またはセミバッチ処理が行なわれるので、原料１は、当該バッチ処理毎に供給される。
【００１７】
また、反応器１０には、気相の原料２を供給するための供給ライン１６が接続されており、この供給ライン１６は、反応器１０内で水平方向に延びる複数の吹き出しライン１８と接続される。この吹き出しライン１８の各々には、原料１を液相ＬＱ内に吹き込むための孔が垂直上方向に向けて多数設けられている。これにより、気相の原料２が供給されると、原料１の液相中に原料２の気泡が分散し、この気泡が上昇することにより流れが生じる。また、気泡中の反応成分が液相中に溶解することにより、酸化、水素化、発酵、廃水処理などの化学反応が生じる。
【００１８】
反応器１０の底部には、化学反応によって生成された製品を取り出すための取り出しライン２０が接続されており、反応器１０内で生じた生成品は、この取り出しライン２０を通じて次工程へ運ばれる。また、この反応器１０の頂点部には、反応器１０内の圧力および液面を調整するため、開閉弁２６を設けた通気ライン２４が接続されている。
【００１９】
＜本発明の概要＞
図２は、液相乱流応力を生じる現象を説明するための模式図である。
【００２０】
図２を参照して、図１に示す気泡塔などの反応器１０では、気泡ＢＬが液相ＬＱに対して速度差をもって紙面上方向に流れる現象が生じる。なお、液相ＬＱ自身も何らかの速度で流れる。このような状況下では、この気泡ＢＬと液相ＬＱとの速度差に伴って、気泡ＢＬの下流側に気泡誘起乱流ＢＩが発生する。また、液相ＬＱの流れ、すなわち液相ＬＱと反応器１０内面との速度差に伴って、せん断誘起乱流ＳＩが発生する。液相ＬＱの乱流渦粘性は、これらの気泡誘起乱流ＢＩおよびせん断誘起乱流ＳＩに依存して定まると考えることができる。
【００２１】
図１に示すような気泡塔内に生じる気泡流では、非特許文献２に記載されているように、せん断誘起乱流ＳＩは気泡誘起乱流ＢＩに比較して小さいと考えられる。そこで、従来提案されていた計算モデルでは、気泡誘起乱流ＢＩのみを考慮した乱流応力が用いられていた（以下、「乱流モデルＡ」とも称す）。
【００２２】
あるいは、非特許文献４に記載されているように、せん断誘起乱流ＳＩおよび気泡誘起乱流ＢＩに起因するそれぞれの乱流応力を線形重畳させた乱流応力を用いることが提案されている（以下、「乱流モデルＢ」とも称す）。
【００２３】
すなわち、図２に示すような状況下で、液相ＬＱに生じる乱流応力τＬＲｅは、気泡誘起乱流ＢＩによって生じる乱流応力τＬＢＩおよびせん断誘起乱流ＳＩによって生じる乱流応力τＬＳＩを用いて、乱流モデルＡおよびＢのそれぞれについて、以下のように示すことができる。
【００２４】
乱流モデルＡ：τＬＲｅ＝τＬＢＩ
乱流モデルＢ：τＬＲｅ＝τＬＢＩ＋τＬＳＩ
しかしながら、乱流モデルＡでは、実際に生じていると考えられるせん断誘起乱流ＳＩが全く考慮されない。また、乱流モデルＢでは、非特許文献５においても指摘されているように、乱流エネルギーが過大に評価されてしまうという問題がある。
【００２５】
そこで、本実施の形態に従うプロセス解析では、せん断誘起乱流ＳＩによって生じる乱流応力τＬＳＩの影響の度合いを気泡誘起乱流ＢＩによって生じる乱流応力τＬＢＩに比較して小さくした乱流モデルを用いる。すなわち、本実施の形態では、以下のような乱流モデルを用いる。
【００２６】
本実施の形態に従う乱流モデル：τＬＲｅ＝τＬＢＩ＋φτＬＳＩ（０＜φ＜１）
ここで、係数φの範囲は、後述する検証実験によれば、０．０５≦φ≦０．２が好ましく、より好ましくは、０．０６≦φ≦０．１５である。さらに好ましくは、φ＝０．１である。
【００２７】
以下では、この乱流モデルを用いて、気泡塔内の液相中に気泡が分散した混相流に対するプロセス解析を行なう。
【００２８】
＜ハードウェア構成＞
本実施の形態に従うプロセス解析装置は、代表的にコンピュータによって実現される。
【００２９】
図３は、この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータ１００の概略構成図である。
【００３０】
図３を参照して、コンピュータ１００は、ＦＤ（Flexible Disk）駆動装置１１１およびＣＤ−ＲＯＭ（Compact Disk-Read Only Memory)駆動装置１１３を搭載したコンピュータ本体１０１と、モニタ１０２と、キーボード１０３と、マウス１０４とを含む。
【００３１】
コンピュータ本体１０１は、ＦＤ駆動装置１１１およびＣＤ−ＲＯＭ駆動装置１１３に加えて、相互にバスで接続された、演算装置であるＣＰＵ（Central Processing Unit）１０５と、メモリ１０６と、記憶装置である固定ディスク１０７と、通信インターフェース１０９とを含む。
【００３２】
本実施の形態に従うプロセス解析処理は、ＣＰＵ１０５がメモリ１０６などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、ＦＤ１１２やＣＤ−ＲＯＭ１１４などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、ＦＤ駆動装置１１１やＣＤ−ＲＯＭ駆動装置１１３などにより記録媒体から読取られて、または通信インターフェース１０９にて受信されて、固定ディスク１０７に格納される。さらに、このようなプログラムは、固定ディスク１０７からメモリ１０６に読出されて、ＣＰＵ１０５により実行される。
【００３３】
ＣＰＵ１０５は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ１０６は、ＣＰＵ１０５のプログラム実行に応じて各種の情報を記憶する。
【００３４】
モニタ１０２は、ＣＰＵ１０５が出力する情報を表示するための表示部であって、一例としてＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）などから構成される。すなわち、モニタ１０２には、プロセス制御の状態などが表示される。
【００３５】
マウス１０４は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード１０３は、入力されるキーに応じたユーザから指令を受付ける。
【００３６】
通信インターフェース１０９は、コンピュータ１００と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。
【００３７】
なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。
【００３８】
＜計算モデル＞
以下、本実施の形態に従うプロセス解析に用いられる計算モデルについて説明する。
【００３９】
（１．基礎式）
本実施の形態に従うプロセス解析に用いる計算モデルでは、気泡流に含まれる相を１つの連続気相（添え字Ｇ）、１つの連続液相（添え字Ｌ）、Ｎ群の気泡（Ｂｍ，Ｂｍ＝１，２，・・・，Ｎ）の合計（Ｎ＋２）種のグループに分類する。
【００４０】
本実施の形態に従うプロセス解析では、各相（（Ｎ＋２）種の各グループ）について、質量保存式および運動量保存式を考えることで、状態値の時間的変化を計算する。なお、
さらにエネルギー保存式を含むようにしてもよい。
【００４１】
以下の計算モデルでは、計算セルサイズΔｘに比べて、気液界面の長さスケールが大きな場合には、その界面で隔てられた二相を連続気相と連続液相とみなす。一方、計算セルサイズΔｘに比べて、小さな気泡は分散相とし、さらにその大きさに応じてＮ群に分類する。各相の体積率αは、体積率の定義から（１）式を満たす。
【００４２】
【数１】
2010040013.tif 000003
【００４３】
気液両相は相変化のない等温非圧縮性ニュートン流体と仮定し、気泡内部の分子および乱流拡散を無視すると、セル体積平均に基づいて、（２）〜（６）式に示す質量保存式および運動量保存式が成立する。
【００４４】
【数２】
2010040013.tif 000004
【００４５】
本実施の形態に従うプロセス解析で用いる計算モデルについても、２つの連続相に対して一流体近似を採用している。そのため、２つの連続相に対する運動量保存式としては、連続相混合体についての（６）式のみを用いる。（６）式における連続相混合体の体積率αＣは、αＣ＝αＧ＋αＬとして与え、混合体の密度ρＣを（７）式として与える。
【００４６】
【数３】
2010040013.tif 000005
【００４７】
液相と気泡との間の界面運動量輸送ＭＬＢｍには、（８）式に示す構成式を用いる。
【００４８】
【数４】
2010040013.tif 000006
【００４９】
ここで、粘性および乱流拡散項Ｆμは、（９）〜（１２）式を用いる。なお、（１２）式の連続液相の乱流応力τＬＲｅについては、後述する。
【００５０】
【数５】
2010040013.tif 000007
【００５１】
なお、本実施の形態に従うプロセス解析で用いる計算モデルでは、連続液相の粘性μＬとしては、分子粘性および渦粘性を考慮する。これについては、後述する。
【００５２】
表面張力ＦＳには、（１３）式を用いる。
【００５３】
【数６】
2010040013.tif 000008
【００５４】
（２．解法）
上述した式を用いて、対象プロセスにおける状態値の時間的変化を計算するための方法について説明する。
【００５５】
上述した（２）〜（４）式を時間的について離散化（時間刻み幅Δｔ）すると、それぞれ（１４）〜（１６）式のように示すことができる。
【００５６】
【数７】
2010040013.tif 000009
【００５７】
すべての相についての（１４）〜（１６）式の総和をとると、（１７）式が得られる。
【００５８】
【数８】
2010040013.tif 000010
【００５９】
（５）式および（６）式において、界面抗力ＭＤｍは半陰解法により離散化される。また、圧力Ｐの項を陰関数として取り扱い、他の項については陽関数として離散化する。その結果、（１８）式および（１９）式を得ることができる。
【００６０】
【数９】
2010040013.tif 000011
【００６１】
ここで、Ａｍ，Ｂｍ，Ｃｍ，Ｄ，Ｇ（ｎ）Ｂｍ，Ｇ（ｎ）Ｃは、時刻ｎで評価される定数またはベクトルを示す。さらに、（１８）式および（１９）式を行列表示すると、（２０）式のようになる。
【００６２】
【数１０】
2010040013.tif 000012
【００６３】
（２０）式の両辺に行列Ｍの逆行列Ｍ−１を乗じると、（２１）式および（２２）式が得られる。
【００６４】
【数１１】
2010040013.tif 000013
【００６５】
時刻（ｎ＋１）における速度および圧力の場については、公知の種々の数値解析方法を用いて、上述した（１７）式、（２１）式、（２２）式の方程式を連係させて解くことで容易に得ることができる。一例として、ＳＭＡＣ法を用いる場合の手順について以下に示す。
【００６６】
（２１）式および（２２）式の圧力Ｐ（ｎ＋１）を既知の圧力Ｐ（ｎ）で置き換えることで、速度の推定量を（２３）式および（２４）式のように定義する。
【００６７】
【数１２】
2010040013.tif 000014
【００６８】
（１７）式および（２１）〜（２４）式から導かれる（２５）式に示すポアソン方程式を用いて圧力修正量δＰを計算する。
【００６９】
【数１３】
2010040013.tif 000015
【００７０】
速度および圧力の場については、（２６）式〜（２８）式に従って、順次更新される。
【００７１】
【数１４】
2010040013.tif 000016
【００７２】
そして、更新された速度を（１４）式〜（１６）式に代入することで、体積率α（ｎ＋１）Ｂｍ，α（ｎ＋１）Ｇ，α（ｎ＋１）Ｌを計算する。連続した二層は、明確な界面で分離される必要があるので、α（ｎ＋１）Ｇおよびα（ｎ＋１）Ｌの計算には、界面追跡法を用いる。（１５）式および（１６）式の解法には、例えばＣＩＰ法（Cubic interpolated propagation）を用いることができる。
【００７３】
以上のように、場の計算に必要となるすべての式が導出されたので、これらの式を用いて、図４に示すような処理手順（後述する）に従って、シミュレーションが実行される。
【００７４】
＜液相乱流応力＞
次に、（１２）式に示す連続液相の乱流応力μＬＲｅについて説明する。
【００７５】
（乱流モデルＡ）
上述したように、従来の乱流モデルＡでは、乱流応力として、気泡誘起乱流ＢＩに起因するせん断力のみを考慮していた。すなわち、乱流モデルＡでは、乱流応力τＬＲｅは、（２９）式のように表わすことができる。
【００７６】
【数１５】
2010040013.tif 000017
【００７７】
そして、気泡誘起乱流ＢＩに起因する乱流応力μＬＢＩは、文献１（Sato, Y. and K. Sekoguchi; "Liquid Velocity Distribution in Two-Phase Bubbly Flow," Int. J. Multiphase Flow, 2, 79-95 (1975)）のモデルを用いると、（３０）式のように評価することができる。
【００７８】
【数１６】
2010040013.tif 000018
【００７９】
すなわち、気泡誘起乱流の度合いを示す値として、液相と気泡との間の速度差（ＶＢｍ−ＶＬ）、気泡の体積率（αＢｍ）、気泡径（ｄＢｍ）に基づいて、気泡誘起乱流ＢＩによって生じる乱流応力τＬＢＩが計算される。
【００８０】
（乱流モデルＢ）
上述したように、従来の乱流モデルＢでは、乱流応力として、気泡誘起乱流ＢＩおよびせん断誘起乱流ＳＩに起因するせん断力を同等に考慮していた。すなわち、乱流モデルＢでは、乱流応力τＬＲｅは、（３１）式のように表わすことができる。
【００８１】
【数１７】
2010040013.tif 000019
【００８２】
乱流モデルＡと同様に、文献１に記載されている乱流渦粘性モデルを用いると、気泡誘起乱流ＢＩによって生じる乱流応力τＬＢＩは、（３２）式のように表わすことができる。また、非特許文献４に記載されている乱流渦粘性モデルを用いて、せん断誘起乱流ＳＩによって生じる乱流応力τＬＳＩは、（３３）式のように表わすことができる。
【００８３】
【数１８】
2010040013.tif 000020
【００８４】
さらに、せん断誘起乱流エネルギーｋおよび乱流エネルギー散逸εの評価には、（３４）式〜（３６）式を用いる。
【００８５】
【数１９】
2010040013.tif 000021
【００８６】
なお、（３５）式〜（３７）式の各種係数には、単相乱流計算における推奨値、具体的には、Ｃμ＝０．０９，σκ＝１．０，σε＝１．３，Ｃε１＝１．４４，Ｃε２＝１．９２を用いる。
【００８７】
すなわち、せん断誘起乱流の度合いを示す値として、液相の流れに伴う乱流エネルギーとエネルギー散逸とに関する輸送方程式に基づいて従って、せん断誘起乱流ＳＩによって生じる乱流応力τＬＳＩが計算される。
【００８８】
（本実施の形態に従う乱流モデル）
上述したように、本実施の形態に従う乱流モデルは、乱流応力τＬＲｅに対する乱流応力τＬＳＩによる影響度を乱流応力τＬＢＩによる影響度に比較して小さくする。これにより、乱流エネルギーが過大に評価されることで生じる誤差を低減する。すなわち、本実施の形態に従う乱流モデルでは、乱流応力τＬＲｅは、１より小さい正の係数φを用いて、（３７）式のように表わすことができる。
【００８９】
【数２０】
2010040013.tif 000022
【００９０】
ここで、係数φの範囲は、後述する検証実験によれば、０．０５≦φ≦０．２が好ましく、より好ましくは、０．０６≦φ≦０．１５である。さらに好ましくは、φ＝０．１である。
【００９１】
＜各モデルのシミュレーション方法＞
図４は、この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。なお、図４に示す処理は、解析対象を空間的に分割して得られるそれぞれのセルについて実行される。
【００９２】
図４を参照して、まず、ユーザは、連続液相の粘性を設定するための係数φを設定する（ステップＳ１００）。すなわち、プロセス解析装置は、乱流応力τｋＲｅに対する係数φを受付ける。そして、ユーザがシミュレーションの開始を指示すると、ステップＳ１０２以降の処理が実行される。
【００９３】
ステップＳ１０２では、各相の体積率α，速度Ｖ，圧力Ｐ，せん断誘起乱流エネルギーｋ，乱流エネルギー散逸εが初期化される。すなわち、シミュレーション対象とするプロセスの初期状態における値を各変数にセットする。その後、処理はステップＳ１０４へ進む。
【００９４】
ステップＳ１０４では、（３２）式〜（３３）式に従って、乱流応力τｋＲｅを計算する。その後、処理はステップＳ１０６へ進む。
【００９５】
ステップＳ１０６では、（２３）式および（２４）式に従って、速度の推定量Ｖｋを計算する。その後、処理はステップＳ１０８へ進む。
【００９６】
ステップＳ１０８では、（２５）式に従って、圧力修正量δＰを計算する。その後、処理はステップＳ１１０へ進む。
【００９７】
ステップＳ１１０では、（２６）式〜（２８）式に従って、速度Ｖｋ（ｎ＋１）および圧力Ｐ（ｎ＋１）を計算する。その後、処理はステップＳ１１２へ進む。
【００９８】
ステップＳ１１２では、（１４）式〜（１６）式に従って、体積率αｋ（ｎ＋１）を計算する。その後、処理はステップＳ１１４へ進む。
【００９９】
ステップＳ１１４では、（３４）式および（３５）式に従って、（１７）式および（１８）式に従って、せん断誘起乱流エネルギーｋ（ｎ＋１），乱流エネルギー散逸ε（ｎ＋１）を計算する。その後、処理はステップＳ１１６へ進む。
【０１００】
ステップＳ１１６では、すべての計算処理が終了したか否かが判断される。すべての計算処理が終了していなければ（ステップＳ１１６においてＮＯ）、時刻ｎに１が加算され（ステップＳ１１８）、次の時刻における各パラメータを算出するために、処理はステップＳ１０４に戻る。そして、ステップＳ１０４以降の処理が再度実行される。
【０１０１】
一方、すべての計算処理が終了していれば（ステップＳ１１８においてＹＥＳ）、処理は終了する。
【０１０２】
このように、解析対象を空間的に分割して得られるそれぞれのセルについて、状態値が時系列に計算される。
【０１０３】
＜検証実験＞
本願発明者らは、模擬装置を用いて、上述した本実施の形態に従うプロセス解析方法について検証を行った。以下、この検証結果について説明する。
【０１０４】
（模擬装置）
まず、本検証実験に用いた模擬装置について説明する。
【０１０５】
図５は、この発明の実施の形態に従うプロセス解析方法の検証実験に用いた模擬装置２００の概略構成を示す図である。
【０１０６】
図５を参照して、模擬装置２００は、筒状の流通器２１０にトレーサを含む流体を一方向に流すとともに、この流通器２１０内の流体経路上に気泡を発生させるスパージャー２２１および２３１を設け、これにより流通器２１０内に液相中に気泡が分散した混相流を生じさせる。そして、トレーサ濃度の時間的変化を測定することで、流通器２１０内の混合状態を評価する。なお、本検証実験では、高さＨ＝１０００ｍｍの流通器２１０を用いた。
【０１０７】
本検証実験では、この模擬装置２００に対して、本実施の形態に従うプロセス解析を実行し、この解析結果と実験結果とを比較することで、本実施の形態に従うプロセス解析方法について評価した。
【０１０８】
模擬装置２００は、流通器２１０内に液相の流れを実現する構成として、流体タンク２１２と、トレーサタンク２１４と、トレーサ供給ポンプ２１６と、供給ポンプ２１８と、開閉弁２４１，２４２，２４３，２４４，２４７と、調整弁２４６と、流量計２４５と、濃度計２５２とを含む。
【０１０９】
流通器２１０は、その内部に流体を流すことのできる筒状（好ましくは、円筒状）の容器であり、流体を受け入れるための供給口２１０ａと、流体を排出するための排出口２１０ｂとを有している。
【０１１０】
この流通器２１０を流れる流体は、供給ポンプ２１８により、流体タンク２１２から流通器２１０の供給口２１０ａに供給される。この流体タンク２１２から供給口２１０ａまでの供給ライン２６２には、流体タンク２１２側から開閉弁２４３、供給ポンプ２１８、流量計２４５、調整弁２４６、開閉弁２４７の順で設けられる。流体タンク２１２から流通器２１０に供給される流体の単位時間当たりの量、すなわち流通器２１０内の流体の速度は、調整弁２４６によって調整可能である。また、流通器２１０に供給される流体の量は、流量計２４５によって測定される。
【０１１１】
なお、流体タンク２１２には、開閉弁２４１を開放することで、流体補充ライン２６１を通じて流体の補充が可能であるとともに、開閉弁２４２を開放することで、この補充される流体を供給ライン２６２に直接的に混合することも可能である。
【０１１２】
また、流通器２１０内の混合状態を評価するためのトレーサは、トレーサタンク２１４から供給される。このトレーサタンク２１４内のトレーサは、トレーサ供給ポンプ２１６によって吐出され、供給ポンプ２１８の前段で流体タンク２１２からの流体と混合する。そして、この混合した流体は、流通器２１０の供給口２１０ａに供給される。
【０１１３】
このトレーサを含む流体は、全体としては一定の速度（液空塔速度）で流通器２１０内を流れ、排出口２１０ｂに到達する。流通器の排出口２１０ｂには排出ライン２５１が接続されており、この排出ライン２５１上には、トレーサの濃度を測定するための濃度計２５２が設けられている。濃度計２５２は、排出口２１０ｂにおけるトレーサ濃度を、トレーサタンク２１４からトレーサが供給され始めてからの経過時間と対応付けて測定する。
【０１１４】
なお、トレーサの濃度の測定方法は、種々の方法が考えられるが、本検証実験には、トレーサとして電解質物質（塩化ナトリウム水溶液）を用い、ステップ応答方法により当該電解質物質を供給した。なお、ステップ応答方法を用いたのは、トレーサの濃度が大きく変化するため、トレーサ濃度の時間的変化をより正確に測定できるためである。
【０１１５】
なお、トレーサとしては、流体の流れに追従し、その存在量が検出できれものであればよく、電解質物質を用いる他に、色素や放射性同位元素を含む溶液を用いることができる。また、色素を用いた場合には、色素を含む流体を画像として捉え、この画像を解析することで濃度の時間的変化を測定することができる。
【０１１６】
また、本実施の形態に従うプロセス解析方法は、気泡塔などの液相中に気泡が分散した混相流が流れる反応器を対象とするので、流通器２１０を流れる流体に気泡を吹き込むための構成として、スパージャー２２１および２３１と、調整弁２２３および２３３と、圧力計２２４および２３４と、流量計２２５および２３５とを含む。
【０１１７】
スパージャー２２１および２３１は、流通器２１０内に垂直方向に並べて配置されており、供給ライン２７１から供給されるガス（空気）を吐出する。より具体的には、スパージャー２２１，２３１の各々は、リング状の導入配管２２０を含んで構成され、この導入配管２２０の垂直上方側には多数の孔が形成されている。なお、本検証実験では、流通器２１０の内径ｒ２＝３００ｍｍとし、導入配管２２０の外径ｒ１＝２５０ｍｍとした。
【０１１８】
スパージャー２２１および２３１から吐出されるガスの量は、それぞれ流量計２２５，２３５で測定されるとともに、その吐出圧は、それぞれ圧力計２２４および２３４で測定される。また、スパージャー２２１および２３１から吐出されるガス流量は、それぞれ調整弁２２３および２３３で調整される。なお、スパージャー２２１から吐出されるガス流量を下ガス流量と称し、スパージャー２３１から吐出されるガス流量を上ガス流量と称する。
【０１１９】
（混合拡散係数の推定）
図５に示す流通器２１０におけるトレーサ濃度変化は、（３８）式のように表わすことができる。
【０１２０】
【数２１】
2010040013.tif 000023
【０１２１】
なお、軸方向距離ｚは、流通器２１０の供給口２１０ａから排出口２１０ｂ方向への距離を意味する。
【０１２２】
（３８）式によれば、濃度計２５２で測定されたトレーサ濃度と、既知の液空塔速度および軸方向距離とに基づいて、混合拡散係数Ｄを推定することができる。すなわち、（３８）式に基づいて算出されるトレーサ濃度の時間的変化が、濃度計２５２で測定された測定値と一致するように、混合拡散係数Ｄについてのフィッティングを行なうことで、混合拡散係数Ｄを算出できる。
【０１２３】
図６は、図５に示す模擬装置２００において測定されたトレーサ濃度の時間的変化と、測定されたトレーサ濃度の時間的変化から推定された混合拡散係数Ｄを示す図である。
【０１２４】
混合拡散係数Ｄがゼロの場合には、流通器２１０内の流体の混合が生じていないことを意味しており、液空塔速度と同じ速度ですべてのトレーサが排出口２１０ｂに到達する。なお、このような流れはプラグフローとも称される。
【０１２５】
一方、混合拡散係数Ｄが大きくなると、排出口２１０ｂにおけるトレーサの濃度は、平均滞留時間よりも早い時間からゆっくりと増加する。
【０１２６】
このように混合拡散係数Ｄは、トレーサ濃度の時間的変化を示すグラフのカーブと近似したカーブを描く滞留時間分布関数を規定するものであり、この混合拡散係数Ｄを求めることにより、流通器２１０を流れる流体の混合状態を評価することができる。
【０１２７】
（検証結果）
上述したように、図５に示す模擬装置２００に対して、本実施の形態に従うプロセス解析を実行することで得られた解析結果と、実際の模擬装置２００を用いて測定された実験結果の比較を図７に示す。
【０１２８】
本検証実験は、連続液相の粘性を設定するための係数φを評価するためのものであるので、係数φを複数の値に変化させて、シミュレーションを実行した。なお、係数φ＝０の場合が従来の乱流モデルＡに対応し、係数φ＝１の場合が従来の乱流モデルＢに対応する。
【０１２９】
また、上述のシミュレーションでは、実験結果と比較するために、液体に追従するトレーサ粒子を流入口２１０ａから供給し、排出口２１０ｂにおけるトレーサ粒子数の時間的変化を算出した。さらに、その時間的変化に基づいて混合拡散係数Ｄを計算した。
【０１３０】
図７は、上ガス流量／総ガス流量と混合拡散係数Ｄとの関係をプロットした図である。なお、上ガス流量は、図５に示す流通器２１０の最上部に設けられたスパージャー２３１の流量を意味し、総ガス流量は、図５に示すスパージャー２２１および２３１に供給される流量の総和を意味する。したがって、上ガス流量／総ガス流量は、その値が大きいほど、流通器２１０の上部に気泡が相対的に多くなっていることを意味し、上ガス流量が多いほど、ガス空塔速度が液空塔速度に比較して速くなる。
【０１３１】
なお、液空塔速度Ｕは１．０ｃｍ／ｓとし、ガス空塔速度Ｕｇは１．５ｃｍ／ｓとして、解析および実験を行なった。
【０１３２】
図７（ａ）は、係数φ＝０（乱流モデルＡ）および係数φ＝１（乱流モデルＢ）と設定した場合にそれぞれ算出される混合拡散係数Ｄを検証実験により得られた混合拡散係数Ｄと比較したものである。図７（ａ）に示すように、乱流モデルＡおよび乱流モデルＢのいずれについても、上ガス流量の比率が小さい場合には、実験結果にやや近い値を示すものの、上ガス流量の比率が大きくなるにつれて実験結果に対する乖離は大きくなる。
【０１３３】
図７（ｂ）は、係数φ＝０．１と設定した場合に算出される混合拡散係数Ｄを検証実験により得られた混合拡散係数Ｄと比較したものである。図７（ｂ）に示すように、係数φ＝０．１と設定することで、上ガス流量の比率が変化した場合であっても、実験結果とほぼ同じ結果を示していることがわかる。すなわち、実際の混相流の挙動をより高い精度でシミュレーションできているといえる。
【０１３４】
図７（ｃ）は、係数φ＝０．０５，０．１，０．２と設定した場合にそれぞれ算出される混合拡散係数Ｄを検証実験により得られた混合拡散係数Ｄと比較したものである。図７（ｃ）に示すように、係数φ＝０．０５，０．２と設定した場合には、係数φ＝０．１と設定した場合に比較して、上ガス流量の比率が大きくなった場合に実験結果との乖離がやや大きくなっているが、上ガス流量の比率が相対的に小さい場合には、実験結果とほぼ同じ結果を示していることがわかる。すなわち、少なくとも係数φが０．０５≦φ≦０．２の範囲では、液相中に気泡が分散した混相流をより高精度に解析できることを意味する。
【０１３５】
図８は、他のパラメータと混合拡散係数Ｄとの関係をプロットした図である。
図８（ａ）は、液空塔速度Ｕと混合拡散係数Ｄとの関係が示されている。なお、図８（ａ）に示される混合拡散係数Ｄは、係数φ＝０．１と設定した場合に算出された値であり、また上ガス流量および下ガス流量は同一であるとする。
【０１３６】
図８（ａ）に示すように、液空塔速度が変化した場合であっても、シミュレーション結果と実験結果とが一致していることがわかる。
【０１３７】
図８（ｂ）は、ガス空塔速度Ｕｇと混合拡散係数Ｄとの関係が示されている。なお、図８（ｂ）に示される混合拡散係数Ｄは、係数φ＝０．１と設定した場合に算出された値であり、また上ガス流量および下ガス流量は同一であるとする。
【０１３８】
図８（ｂ）に示すように、ガス空塔速度が小さい場合には、若干の乖離があるものの、全体としては、ガス空塔速度が変化した場合であっても、シミュレーション結果と実験結果とが一致していることがわかる。なお、ガス空塔速度が小さい場合には、検証実験によって得られた混合拡散係数Ｄ自体もばらつきが大きいので、このような点を考慮すると、十分に実用に耐え得る精度でシミュレーションができているといえる。
【０１３９】
＜本実施の形態による効果＞
本実施の形態によれば、気泡誘起乱流によって生じる乱流応力およびせん断誘起乱流によって生じる乱流応力を適切な比率で反映した乱流応力を用いて、液相の特性の１つである液相乱流応力を算出するので、液相中に気泡が分散した混相流での現象をより高精度にプロセス解析することができる。
【０１４０】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は、上記した説明ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【図面の簡単な説明】
【０１４１】
【図１】この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【図２】液相乱流応力を生じる現象を説明するための模式図である。
【図３】この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図４】この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。
【図５】この発明の実施の形態に従うプロセス解析方法の検証実験に用いた模擬装置２００の概略構成を示す図である。
【図６】図５に示す模擬装置において測定されたトレーサ濃度の時間的変化と、測定されたトレーサ濃度の時間的変化から推定された混合拡散係数Ｄを示す図である。
【図７】上ガス流量／総ガス流量と混合拡散係数Ｄとの関係をプロットした図である。
【図８】他のパラメータと混合拡散係数Ｄとの関係をプロットした図である。
【符号の説明】
【０１４２】
１０  反応器、１２，１６  供給ライン、１４，２２，２６  開閉弁、１８  吹き出しライン、２０  取り出しライン、２４  通気ライン、１００  コンピュータ、１０１  コンピュータ本体、１０２  モニタ、１０３  キーボード、１０４  マウス、１０６  メモリ、１０７  固定ディスク、１０９  通信インターフェース、１１１  ＦＤ駆動装置、１１３  ＣＤ−ＲＯＭ駆動装置、２００  模擬装置、２１０  流通器、２１０ａ  供給口、２１０ｂ  排出口、２１２  流体タンク、２１４  トレーサタンク、２１６  トレーサ供給ポンプ、２１８  供給ポンプ、２２０  導入配管、２２１，２３１  スパージャー、２２３，２４６  調整弁、２２４  圧力計、２２５，２３５，２４５  流量計、２４１，２４２，２４３，２４４，２４７  開閉弁、２５１  排出ライン、２５２  濃度計、２６１  流体補充ライン、２６２，２７１  供給ライン。
【図１】
2010040013.tif 000024
【図２】
2010040013.tif 000025
【図３】
2010040013.tif 000026
【図４】
2010040013.tif 000027
【図５】
2010040013.tif 000028
【図６】
2010040013.tif 000029
【図７】
2010040013.tif 000030
【図８】
2010040013.tif 000031