"液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、前記状態値を時系列で計算するステップは、前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算するステップと、前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算するステップと、前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算するステップと、前記化学種保存式に基づいてモル濃度を計算するステップとを含み、前記液相の体積率を計算するステップは、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、前記モル濃度を計算するステップは、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析プログラム。前記第1液相領域の体積率は、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる、請求項1に記載のプロセス解析プログラム。前記パラメータは、0.01~1.0の範囲の値に設定される、請求項2に記載のプロセス解析プログラム。液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、前記状態値を時系列で計算する手段は、前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算する手段と、前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算する手段と、前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算する手段と、前記化学種保存式に基づいてモル濃度を計算する手段とを含み、前記液相の体積率を計算する手段は、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、前記モル濃度を計算する手段は、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析装置。しかしながら、従来の計算モデルによるプロセス解析では、系のサイズによっては、シミュレーション結果が実際の反応系を精度よく再現できないという課題があった。具体的には、たとえば上記したような気泡塔が、直径0.4m以上、高さ4m以上のものとなった場合に、シミュレーション結果の実際の反応系からの乖離が顕著に現われていた。本発明は係る実情に鑑み考え出されたものであり、その目的は、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようにすることである。本発明に従ったプロセス解析プログラムは、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算するステップは、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算するステップと、質量保存式および状態保存式に基づいて、液相の体積率を計算するステップと、質量保存式および状態保存式に基づいて、気相の体積率を計算するステップと、化学種保存式に基づいてモル濃度を計算するステップとを含み、液相の体積率を計算するステップは、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算するステップは、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。好ましくは、第1液相領域の体積率は、気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる。好ましくは、パラメータは、0.01~1.0の範囲の値に設定される。本発明に従ったプロセス解析装置は、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算する手段は、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算する手段と、質量保存式および状態保存式に基づいて、液相の体積率を計算する手段と、質量保存式および状態保存式に基づいて、気相の体積率を計算する手段と、化学種保存式に基づいてモル濃度を計算する手段とを含み、液相の体積率を計算する手段は、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算する手段は、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。本発明によれば、液相中に気泡が分散した混相流を解析する際に、液相において、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算される。これにより、気泡によって液体が流れを変えられる作用を考慮でき、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようになる。本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。<解析対象>この発明の実施の形態に従うプロセス解析は、気泡塔などの液相中に気泡が分散した混相流を生じる系を解析対象とする。図1は、この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。図1を参照して、解析対象の一例として、反応器10内に原料1液相を供給するとともに、原料1中に原料2気相が吹き込まれる系SYSを考える。反応器10の底部には、原料1を供給するための供給ライン12が接続されている。この供給ライン12の経路上には、原料1を供給/遮断するための開閉弁14が設けられている。なお、図1に示す系SYSでは、バッチ処理またはセミバッチ処理が行なわれるので、原料1は、当該バッチ処理毎に供給される。また、反応器10には、気相の原料2を供給するための供給ライン16が接続されており、この供給ライン16は、反応器10内で水平方向に延びる吹き出しライン18と接続される。吹き出しライン18には、原料1を液相L内に吹き込むための孔が垂直上方向に向けて多数設けられている。これにより、気相の原料2が供給されると、原料1の液相中に原料2の気泡が分散し、この気泡が上昇することにより流れが生じる。また、気泡中の反応成分が液相中に溶解することにより、酸化、水素化、発酵、廃水処理などの化学反応が生じる。反応器10の底部には、化学反応によって生成された製品を取り出すための取り出しライン20が接続されており、反応器10内で生じた生成品は、この取り出しライン20を通じて次工程へ運ばれる。また、この反応器10の頂点部には、反応器10内の圧力および液面を調整するため、開閉弁26を設けた通気ライン24が接続されている。<本発明の概要>図1に示したような気泡塔内では、たとえば以下のa~eに示すような現象が起こっていると考えられる。a底部から液体が供給されるb塔内部に空気が供給されるc液体液相中を、空気が気泡として分散しながら流れるd気泡から液体液相中にガスが溶解するe溶解したガスと液体が反応する本明細書では、これだけの複雑な過程を含む反応系について、まず、ラボの規模で反応系を検討し、そして、ミニプラント、実機へと順にスケールアップさせることとした。つまり、ラボやミニプラントで得た知見を実機の合理化設計へと適切に展開するために、上記現象を様々なスケールで捉え、スケールごとに重要となる因子を議論することにより、スケールアップのあり方を調整することとした。そこで、本明細書では、次の表1に挙げられた、ミクロ、メゾ、マクロの各スケールにおいて検討を行なう。なお、表1では、生じると考えられる主な現象によって、スケールが定義されている。ミクロスケールは、反応機構レベルで現象を捉えるスケールである。マクロスケールは、反応塔内における液混合全体的に現象を捉えるスケールである。メゾは、ミクロとマクロの中間のスケールであり、ガスの溶解や、その際に各気泡が液を同伴して移動するなど、相間の界面での現象を捉えるスケールである。ミクロスケール同じ反応機構を持っていた場合であっても、気泡塔内の液体の流れによって反応成績が変化する。その原因について、以下説明する。反応容器における流体の流れとして、供給した液体が瞬時にして容器全体に混合する流れを完全混合と、供給した液体が流出部へと一次元的に移動する押し出し流れとがある。ラボで使用される攪拌槽反応器における流れは、完全混合とみなすことができる。一方、押し出し流れは、反応容器内の液体が全く混合していない状態と捉えることができる。たとえば、管内の流れなどは、押し出し流れの一例と言える。そして、通常の反応容器内の流れは、完全混合と押し出し流れの間に位置づけられる。ここで、次のEQ式に示すような、出発物質Aが中間体Rを経て、最終生成物Sへと変化するA→R→Sという逐次反応を対象として、Rの選択率を検討する。なお、k1,k2は、それぞれの反応の反応速度定数である。供給されたAは、反応容器の流出部に進むにつれてRへと変化する。Rの一部は、反応が進んでSへと変化する。反応容器において混合があると、供給したAが反応せずに流出部へと流れやすくなる一方で、Aが反応して生成したRが入口に逆流して、Sへと変化する割合が増える。反応液に限らず、流体の流れの状態を表す指標として混合段数混合槽数という値が知られている。この混合段数とは、筒状の胴体を有する容器以下、適宜流通器と称するの中を一方向に流れる流体が、流通器を流れる間に混合される度合いを示す値であり、上記流通器を、直列に連結された互いに同一体積の複数の完全混合槽にモデル化した場合の、当該完全混合槽の数である。混合段数については、文献A橋本健治、反応工学、改訂版、培風館、1993年3月20日、p.179-187、195―197および文献BOctaveLevenspielChemicalreactionengineeringFifthprintingUSAJohnWileyandSonsINCMay1967p261265にその説明が記載されている。そして、この混合段数が1最小の場合が、押し出し流れの状態であり、混合段数が無限大の場合が、完全混合の状態である。EQ式の反応における、A成分、O成分、R成分、および、S成分の各濃度についての、従来の技術に従ったシミュレーション結果について、説明する。図2は、k1とk2の値が同じである場合の、各成分の濃度の時間変化を示す図である。なお、図2Aは、反応容器内が均一の場合つまり、完全混合の場合を示し、図2Bは、反応容器内の不均一の場合つまり、押し出し流れの場合を示す。また、各図では、A成分の濃度が実線で、O成分の濃度が破線で、R成分の濃度が一点鎖線で、そして、S成分の濃度が点線で、記載されている。図2Bの不均一の場合には、反応することなく残留するA成分およびR成分がほとんど見られないのに対し、図2Aの均一の場合には、或る程度の量のA成分およびR成分が反応することなく残留している。図3は、k1がk2よりも大きい場合の、各成分の濃度の時間変化を示す図である。なお、図3Aは、反応容器内が均一の場合つまり、完全混合の場合を示し、図3Bは、反応容器内の不均一の場合つまり、押し出し流れの場合を示す。図3では、各成分の濃度変化が、図2と同様の線の種類で記載されている。また、図4は、k1がk2よりも小さい場合の、各成分の濃度の時間変化を示す図である。なお、図4Aは、反応容器内が均一の場合つまり、完全混合の場合を示し、図4Bは、反応容器内の不均一の場合つまり、押し出し流れの場合を示す。図4では、各成分の濃度変化が、図2と同様の線の種類で記載されている。図3と図4に示されたいずれの場合においても、反応容器内が均一である場合と不均一である場合とでは、各成分の濃度変化の挙動に差異を生じている。つまり、図2~図4から理解されるように、反応系における各成分の濃度の時間変化は、k1とk2の間の大小関係に関わらず、反応容器内が均一であるか不均一であるかによって相違することとなる。図5は、EQ式の中間生成物であるR成分の、生成量溶液中の濃度の時間変化を示す図である。なお、図5Aは、図2A,図2B,図3A,図3B,図4A,図4Bのそれぞれに示されたR成分についての濃度変化を示す。そして、図5Bは、図5Aに示されたそれぞれの状況下でのR成分の濃度変化の、反応初期における変化を拡大してつまり、図5Aの各グラフを横軸方向に拡大して、示すものである。図5Aおよび図5Bに示されるように、中間生成物の生成量は、反応条件によって、大きく異なるといえる。ここでいう反応条件とは、k1とk2の間の大小関係と、反応容器における流体の流れ完全混合/押し出し流れとの組合せである。また、図5Aおよび図5Bにおいて、R成分の生成量の挙動の差異は、反応初期だけでなく、反応中期においても見られている。なお、k1とk2の大小関係のいずれであってもk1=k2、k1>k2、または、k1<k2、R成分の濃度は、反応初期には、完全混合よりも押し出し流れの方が高くなっているが、反応中期以降は、完全混合の方が高くなる傾向にある。つまり、同じ反応率において完全混合よりも押し出し流れの方が、Sを生成する反応についての選択率が高いと言える。以上のことから、反応系における生成物の生成量を予測するためには、k1とk2の間の大小関係とともに、反応容器における流体の流れについても検討することは重要であるといえる。メゾスケール上記したように、反応における選択率は反応容器内の流体の流れに影響を受け、反応機構が同じであれば、完全混合よりも選択率の低い流れは存在しないものと考えられる。しかしながら、実機における反応成績には、完全混合よりも選択率が劣る場合が観測されている。ここで、上記したような考察は、反応容器内では液体が単相で流れているという想定に基づいたものであった。これに対して、実際の気泡塔内では、気泡と液体が混在した混相流の状態になっている。このことから、上記したような、完全混合よりも選択率の低い系が存在したという結果は、気泡が液体を同伴して流れたことによる特別な作用に起因するものではないかと考えられる。図6は、気泡が液体を同伴して流れる状態の概念図である。図6では、液相L内を、気泡Bが、液体ELを同伴して、上方へ流れている状態が示されている。従来の計算モデルでは、計算セルサイズΔxに比べて、小さな気泡を分散相とし、また、その大きさに応じてN群に分類することにより、各相の体積率αを、1式で表していた。1式では、気泡流に含まれる相を1つの連続気相添え字G、1つの連続液相添え字L、そして、N群の気泡Bm:Bm=1,2,,Nの合計N2種のグループに分類がなされていた。本実施の形態では、液相において、気泡の影響を受ける可能性を考慮して、体積率の検討において領域の分類を図る。具体的には、図7に示されるように、液相において、気泡の影響を受ける液体ELを第1液相領域とし、気泡の影響を受けにくいと考えられるバルクの液相Lを第2液相領域とする。本実施の形態では、1式においてαLで記述されていた液相の体積率は、第1液相領域の体積率αL1と第2液相領域の体積率αL2とを用いて、次の式のように記述される。そして、第1液相領域の体積率αL1は、それを同伴する気泡の体積率αBに依存すると考えられる。したがって、各液体ELの体積率αL1は、各気泡の体積率αBを用いて、次の式のように記述される。ここで、βは、気泡の影響範囲を表すパラメータである。そして、液相の体積率は、次の式で示すように、第1液相領域の体積率αL1の総量と第2液相領域の体積率αL2の総量の和として記述される。上記したように、本実施の形態では、計算モデルにおける体積率について、気泡の影響の有無を考慮し、液相の領域を第1液相領域と第2液相領域に分割して検討している。そして、第1液相領域の体積率αL1を、気泡の体積率αBとパラメータβによって、つまり、系のスケールに依存しない値によって、言及している。これにより、本実施の形態では、系のスケールに依存することなく、液相における気泡の影響を考慮できるため、より実際の反応における挙動に沿ったシミュレーションが可能となる。なお、本明細書では、気泡Bと区別するために、反応系内のバルクの気相を適宜気相Gと表記する。<ハードウェア構成>本実施の形態に従うプロセス解析装置は、代表的にコンピュータによって実現される。図8は、この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータ100の概略構成図である。図8を参照して、コンピュータ100は、FDFlexibleDisk駆動装置111およびCD-ROMCompactDiskReadOnlyMemory駆動装置113を搭載したコンピュータ本体101と、モニタ102と、キーボード103と、マウス104とを含む。コンピュータ本体101は、FD駆動装置111およびCD-ROM駆動装置113に加えて、相互にバスで接続された、演算装置であるCPUCentralProcessingUnit105と、メモリ106と、記憶装置である固定ディスク107と、通信インターフェース109とを含む。本実施の形態に従うプロセス解析処理は、CPU105がメモリ106などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、FD112やCD-ROM114などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、FD駆動装置111やCD-ROM駆動装置113などにより記録媒体から読取られて、または通信インターフェース109にて受信されて、固定ディスク107に格納される。さらに、このようなプログラムは、固定ディスク107からメモリ106に読出されて、CPU105により実行される。CPU105は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ106は、CPU105のプログラム実行に応じて各種の情報を記憶する。モニタ102は、CPU105が出力する情報を表示するための表示部であって、一例としてLCDLiquidCrystalDisplayやCRTCathodeRayTubeなどから構成される。すなわち、モニタ102には、プロセス制御の状態などが表示される。マウス104は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード103は、入力されるキーに応じたユーザから指令を受付ける。通信インターフェース109は、コンピュータ100と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。また、コンピュータ100を本実施の形態に従うように機能させるプログラムが記録される記憶媒体は、上記したFDやCD-ROMに限定されない。記憶媒体の他の例としては、DVD-ROMDigitalVersatileDiskReadOnlyMemory、USBUniversalSerialBusメモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、MOMagneticOpticalDisc、MDMiniDisc、ICIntegratedCircuitカードメモリカードを除く、光カード、マスクROM、EPROM、EEPROMElectronicallyErasableProgrammableReadOnlyMemoryなどの、不揮発的にプログラムを格納する媒体が挙げられる。<計算モデル>以下、本実施の形態に従うプロセス解析に用いられる計算モデルについて説明する。1.基礎式本実施の形態に従うプロセス解析に用いる計算モデルでは、気泡流に含まれる相を1つの連続気相添え字G、液相領域添え字L、N群の気泡Bm,Bm=1,2,,Nの合計N2種のグループに分類する。本実施の形態に従うプロセス解析では、各相N2種の各グループについて、質量保存式、運動量保存式、状態方程式、および、化学種保存式を考えることで、状態値の時間的変化を計算する。以下の計算モデルでは、計算セルサイズΔxに比べて、気液界面の長さスケールが大きな場合には、その界面で隔てられた二相を連続気相と連続液相上記した第2液相領域とみなす。一方、計算セルサイズΔxに比べて、小さな気泡は分散相とし、さらにその大きさに応じてN群に分類する。各相の体積率αは、体積率の定義から上記1式を満たす。気液両相は相変化のない等温非圧縮性ニュートン流体と仮定し、気泡内部の分子および乱流拡散を無視すると、セル体積平均に基づいて、2式~6式に示す質量保存式および運動量保存式が成立する。ここで、下付添字のGは連続気相を、Lは連続液相を、Bmはグループmm=1,2,...,Nの気泡を意味する。また、tは時間、Vは速度、ρは密度、ΓCBは合体分裂による単位体積単位時間当りの質量生成量、ΓAEは吸収蒸発による単位体積単位時間当りの質量生成量を表す。なお、質量生成量Γの下付添字は輸送に関与する二相を意味し、たとえばLBmは気泡グループBmから連続液相Lへの輸送を表す。また、DBm/DtおよびDC/Dtは、以下の実質微分を意味する。本実施の形態に従うプロセス解析で用いる計算モデルについても、2つの連続相に対して一流体近似を採用している。そのため、2つの連続相に対する運動量保存式としては、連続相混合体についての6式のみを用いる。6式における連続相混合体の体積率αCは、αC=αGαLとして与え、混合体の密度ρCを次の7式として与える。また、本実施の形態に従うプロセス解析で用いる計算モデルでは、化学種の量を評価するパラメータとして、モル濃度モル分率分密度分圧質量分率絶対湿度があげられる。ここでは、質量分率Yとモル濃度Cを用いて、状態方程式を構成する場合を考える。グループmの気泡のみが存在する場合、Yの定義から次の8式が成立する。ここで、下付添字のiは、化学種を表す。また、ρBmiP,TBmは、純物質時の状態式である8式から、平均密度に関して次の9式が得られる。連続液相と連続気相の平均密度も同様の手順で得られる。本実施の形態では、各相Bm,L,Gを一括して表す下付添字kを用い、以下の10式のように表記する。10式では、添字kを用いて表記する便宜上、TcをTL,TG一流体近似からTc=TL=TGに置き換えている。一方、モル濃度Ckiと平均密度との関係は、次の11式を採用している。ここで、Mは分子量をあらわす。変数YとCは、10式を用いる方法と11式を用いる方法は等価である。これらの二式を比較した場合、変分が単純となる11式を用いる方が簡便である。また、多くの場合、反応速度はCに直接依存するため、多くの等温系液相流れの計算でCが使用されている。しかしながら、圧縮性流体を扱う場合はCに温度圧力依存性を加味する必要がある。一方、10式を用いる場合は、ρkiP,Tkを物性値推算法データベースなどを用いて設定することで、密度の温度圧力依存性を加味できる。このため、圧縮性非圧縮性の区別なくYの輸送計算ができる利点があり、ρkiP,Tkが明確なガス燃焼流れなどの計算に良く使用されている。工業的なガス吸収操作では、ガスを溶質として液相側で反応させる場合が極めて多い。本実施の形態では、気相Bm,Gの化学種を質量分率Y、液相Lの化学種をモル濃度Cで評価し、化学種保存式として、次の12式~14式を用いる。なお、13式において、Kは気泡BMまたは液相に溶解している気体Cを表す。これにより、13式では、気体成分について、気泡BMの状態および液相内の溶質の状態の運動が考慮されている。ここで、wは反応に伴う単位体積当たりのモル生成率、Dは拡散係数、BLBmiは気泡グループと連続液相における単位体積単位時間当たりの化学種iの生成量、ΦBmiとΦGiはΓGBmに伴う化学種iの生成量を表す。12式~14式の右辺第一項は生成項、12式の第二項と14式の第三項はΓに伴う生成項、13式と14式の第二項は拡散項、13式の第三項は温度圧力依存項である。なお、本実施の形態では、図7を参照して上記したように、液相を第1液相領域と第2液相領域に分割して検討を行なうため、13式が改良される。改良にあたり、検討を容易にするため、次の仮定を設ける。仮定1:連続気相から液相への物質移動はない仮定2:各モル濃度は圧力の影響を受けない仮定3:気相からの物質移動は第1液相領域への移動のみが生じ、その後、第1液相領域と第2液相領域の間で交換が生じる仮定1によれば次の15式が導かれ、仮定2によれば次の16式が導かれる。また、モル濃度に関する13式を第1液相領域と第2液相領域に分割して検討すると、次の17式および18式が導かれる。ここで、Ljiは液相の領域jにいる成分iを意味している。また、BL12iは第1液相領域と第2液相領域の間での物質交換量であり、BL1Bmiは気相からの物質移動についての物質交換量である。気相からの物質移動が無い成分については、の場合、BL1Bmi=0となる。したがって、この場合、17式と18式は次の19式と20式のように書ける。17式~20式を解くには、本実施の形態において新たに採用したパラメータβとBL12を指定する必要がある。βは、気泡の形状、上昇速度、密集度、乱れの状況など、多くの流れの作用を受ける複雑なパラメータであるが、既往の研究例えば、RCliftJRGraceMEWeberBubblesDropsandParticlesDoverPublications2005からb001~10の範囲が考えられる。BL12についても、気泡周りに関する詳細な知見から与える必要がある。これらの情報は、以下の研究内容から取得できる。レーザー誘起蛍光法例えば、文献MizutaKOdoYYoshimitsuSKajiHMurakamiMandMatsumotoTTransientBehaviorduringaGasDissolutionProcessaroundaStationarySingleBubbleJChemEngJapan417pp553556に記載の方法で溶解する物質に及ぼす気泡の影響を調べられる。シミュレーション例えば、界面追跡法、境界適合格子法による計算で予測する。2.構成式界面運動量輸送MLBmFには、次の21式で表される構成式が用いられる。また、粘性および乱流拡散項Fμには、次の22式で表される構成式が用いられる。ここで、上付添字Tは転置であり、μcは連続相混合体の粘度である。また、表面張力FSには、次の23式が用いられる。また、化学種輸送量Bおよびこれに伴う質量輸送ΓAEには、次の24式~27式が用いられる。3.解法上述した式を用いて、対象プロセスにおける状態値の時間的変化を計算するための方法について説明する。3-1.化学種保存式の解法まず、化学種を求めるための解法を提示する。12式~14式および15式~20式を参照して説明した化学種保存式には、対流項、反応項、相間輸送項、拡散項、温度圧力依存項、体積輸送に伴う項が含まれている。各項の特性時間が大きく異なる場合、同じ時間積分法で各項を計算する方法は非効率である。そこで、各項を相間輸送項上付数字INT、反応項RT、温度圧力依存項PT、それ以降の項OTに分類し別々に取り扱うことで、計算の効率向上を図る。YまたはCを一括してΨとすると、その時間変化は、次の28式で評価できる。ここで、上付添字n1は時刻n1Δtを表す。28式は、後述する圧力速度計算時に必要となる。圧力速度計算後、Ψの時間更新を行なう。時間更新時には、計算不安定化を回避するため、陰解法を用いて相間輸送項による計算を行なった後、次の29式を用いて相間輸送項以外の項で時間更新する。ここで、上付添字は、相間輸送項のみを用いて更新する値を意味する。以下に、相間輸送項、反応項、温度圧力依存項、それ以外の項の評価法を整理する。3-1-1.相間輸送項気泡塔内気泡流では、気泡から液相へのガス吸収が最も支配的な相間輸送項である。そこで、簡便のため、以下では相間の化学種iの輸送を気泡から液相へのガス吸収に限定する。12式~14式等を参照して説明したような相間輸送項に後述のモル流束に関する相関式を代入すると、以下の30式および31式のように整理できる。ここで、ξ,Ψ,ζ,ηは代入操作によって生じるYBmi,CLiの各係数をまとめたものである。時間変化項を評価する際には、30式~31式の右辺におけるYBmi,CLiを時刻nで評価した、次の32式および33式を用いる。一方、時間更新時には本項による計算不安定化を回避するため、30式~31式の右辺を陰的に扱った、次の34式および35式を用いる。34式と35式を行列表示に置換え、さらに、YBmi,CLiについて解くと、次の36式が得られる。3-1-2.反応項ここでは、液相内での反応のみを取扱う。化学反応速度式は非線形連立方程式で、特性時間が運動量輸送式離散化的に用いる時間刻み幅Δtより短い場合が多い。計算不安定化と計算誤差を回避するため、反応項の離散化には半陰解法を用いる。さらに、ΔtをJ分割したΔtJ=Δt/Jを定義し、l=1,...,Jとして、次の37式による計算をCLiに関して繰返す。ここで、上付添字のnは時刻nΔtを、lは時刻nΔtlΔtJを表しており、l=0は時刻nΔtを、l=Jは時刻n1Δtを意味する。37式の計算過程においてCLiが平衡状態に達したとみなされる場合は、J回以下でも繰返し計算を終了できる。反応項による時間変化は、次の38式で評価する。3-1-3.温度圧力依存項連続液相のモル濃度に対して温度圧力に関する変分を取ると、次の39式が得られる。δCLi/δTCn,δCLj/δPは物性値推算法データベースなどから取得する。DCTC/DtおよびDCP/Dtの評価法に関しては、後述する。Yを使用する場合は、本項は存在しない。3-1-4.上記以外の項上記以外の項として、拡散項および質量輸送に伴う項がある。これらの項では、次の40式を用い、各変数を陽的に扱う。3-2.密度の化学種温度圧力依存性を考慮した解法上記の説明により、化学種の時間変化および時間更新計算法が構築できたので、これを質量運動量エネルギ保存式に組込み、全体の計算手順を整理する。密度の化学種温度圧力依存性を考慮して、2式~4式で示した質量保存式を変形する。15式で表した状態方程式から、密度の変分に関して、次の41式が成り立つ。そして、41式により、密度のラグランジュ微分は、次の42式のように書くことができる。ここでは、Tcと同様に、Dc/Dtの、DL/Dt,DG/Dtに置換えられている。そして、モル濃度を用いる場合は、11式から、次の43式のように記述することができる。また、2式で表される気泡の質量保存式と上記42式により、次の44式が得られる。そして、44式を質量保存式として、解法を構成する。次に、各項を評価する時刻を、45式に基づいて定める。次に、運動量およびエネルギ保存式の各項を評価する時刻を、46式および47式のように定める。そして、連続液相と連続気相の基礎式に関しても同様の操作を行なう。ただし、連続液相ではモル濃度を用いているため、45式に対応する式は、次の48式となる。Tomiyamaetalは、温度が速い過渡変化を伴う場合、温度および密度を評価する時刻が数値解に有意な差を生じさせることを指摘している。45式~48式においても密度変化の影響が各式に可能な限り反映されるように、各項を評価する時刻を選択している。圧力のラグランジュ微分項であるDkP/Dtは、次の49式で評価する。ShimadaandTomiyamaでは、大きな静水圧差下の気泡流計算を通して、49式の右辺第1項を無視する近似の妥当性を確認している。28式と47式,48式により、時刻nの値のみを用いてDkYki/Dt,DkP/Dtが求められる。このため、質量保存式気泡の場合44式の右辺は既知量となり、エネルギ保存式および化学種保存式を用いずに質量運動量保存式のみを用いて新しい時刻の速度Vk圧力Pを計算できる。1式を考慮して45式に相当する式をすべての相に関して辺々を加え合せると、以下の混合体に対する質量保存式である50式を得る。次に、21式を運動量保存式気泡の場合46式に代入し、さらに、VBmおよびVGについて整理すると、次の51式および52式が得られる。ここで、Am,Bm,Cm,Dは、時刻nで評価される項をまとめたものである。そして、39式と40式を行列表示すると、次の53式が得られる。そして、53式をVBm,VCについて解き、整理すると、次の54式が得られる。ここで、ρの上付添字~およびGの上付添字~は、41式の右辺に逆行列を乗じた項を置換えたものである。50式および54式からVk,Pを求める手順には、種々の非圧縮性流れの数値解法を利用できる。ここでは、SMAC法SimplifiedMarkerandCellmethodAmsdenandHallow1970に基づく求解手順を整理しておく。まず、圧力修正量δPを次の55式で定義する。次に、54式における圧力を既知の圧力Pで置換えた次の56式で、速度の推定値Vkの上付添字~を定義する。50式に、54式~56式を代入すると、以下のδPに関するポアソン方程式である57式が得られる。化学種と同様に温度を時間更新する際も、計算安定性確保のため界面エネルギ輸送項を陰的に評価する必要がある。温度計算の詳細については、文献ShimadaandTomiyama2005などで開示されている技術を採用することができる。<各モデルのシミュレーション方法>図9は、この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。なお、図9に示す処理は、解析対策を空間的に分割して得られるそれぞれのセルについて実行される。図9を参照して、まず、プロセス解析装置は、ステップS100で、各相の密度ρ,体積率α,速度V,圧力P,モル濃度C,質量分率Yを初期化する。すなわち、シミュレーション対象となるプロセスの初期状態における値を各変数にセットする。その後、処理は、ステップS102へ進められる。ステップS102では、10式および11式に従って、各化学種についてのδρ/δPnを算出する。その後、処理はステップS104に進められる。ステップS104では、28式、47式および48式に従って、DkYki/Dt,DkP/Dtを計算する。その後、処理はステップS106へ進められる。ステップS106では、56式に従って、Vkの上付添字~を計算する。その後、処理はステップS108へ進められる。ステップS108では、57式に従って、圧力修正量δPを計算する。その後、処理はステップS110へ進められる。ステップS110では、54式および55式に従って、速度Vkおよび圧力Pを計算する。その後、処理はステップS112へ進められる。ステップS112~ステップS116では、気泡の体積率αBmn1と、気相の体積率αGn1と、液相の体積率αLn1が、それぞれ計算される。なお、これらの体積率は、1式、45式、および、48式に従って求められる。そして、処理はステップS118へ進められる。ステップS118では、図6および図7を参照して説明した第1液相領域ELと第2液相領域のそれぞれの体積率αLjn1が、以下の58式および59式に従って計算される。そして、処理はステップS120へ進められる。ステップS120では、17式~20式に従って第1液相領域と第2液相領域それぞれのモル濃度CLjin1が計算される。そして、ステップS122へ処理が進められる。ステップS122では、各化学種の質量分率Ykn1が計算される。そして、ステップS124へ処理が進められる。ステップS124では、10式および11式に従って、密度ρkn1が計算される。その後、処理はステップS126へ進められる。ステップS126では、すべての計算処理が終了したか否かが判断される。すべての計算処理が終了していなければステップS126においてNO、時刻nに1が加算されステップS128、次の時刻における各パラメータを算出するために処理はステップS102に戻される。そして、ステップS102以降の処理が再度実行される。一方、すべての計算処理が終了していればステップS128においてYES、処理は終了する。このように、解析対象を空間的に分離して得られるそれぞれのセルについて、状態値が時系列に計算される。<検証実験>本願発明者らは、模擬装置を用いて、本実施の形態におけるパラメータβ,BL12を取得するためのシミュレーションについての検証を行なった。検証は、水のSO2ガス吸収操作を対象としている。図10は、気泡流入開始後1秒におけるpHの分布を示す図である。なお、図10Aは実験結果であり、図10Bは、本実施の形態に従うプロセス解析方法による計算結果である。そして、図10Aおよび図10Bでは、pHの値がハッチングの違いで表現されている。具体的には、図中、pH4.4からpH5.2までのpHの変化が、P1~P6の6段階で示されている。これらは、P1がpH4.4に対応し、P6がpH5.2に対応する。また、図11は、気泡流入開始後1秒における速度ベクトルの分布を示す図である。なお、図11Aは実験結果であり、図11Bは、本実施の形態に従うプロセス解析方法による計算結果である。図10および図11を参照して、気泡流入開始後1秒では、pHと速度ベクトルのいずれについても、反射の影響により、気泡周辺でのパラメータの計算結果が実験値に沿ったものとはなっていない。図12および図13に、気泡流入開始後2秒における実験結果および計算結果を示す。なお、図12は、pHの分布を示す図である。図12Aは実験結果であり、図12Bは、本実施の形態に従うプロセス解析方法による計算結果である。図12Aおよび図12Bにおいても、pHの変化は、図10Aおよび図10Bと同様に、6段階のハッチングで示されている。また、図13は、速度ベクトルの分布を示す図である。なお、図13Aは実験結果であり、図13Bは、本実施の形態に従うプロセス解析方法による計算結果である。図12および図13を参照して、気泡流入開始後2秒では、実験結果において、気泡周辺に速度の高い領域が見られ、この傾向は、計算結果でも見られている。このことから、計算結果は、実験結果に沿ったものと言える。<実施の形態の効果>たとえば、混相流におけるA→R→Sの逐次反応を対象としてRの選択率を考えた場合、反応容器に供給されたAは、当該反応容器の流出部に進むにつれてRへと変化する。そして、Rの一部は反応が進んでSへと変化する。このような反応系において、反応容器内で混合があると、供給したAが反応せずRへと変化することなくに流出部へと流れやすくなる。これと同時に、Aが反応することによって生じたRのうち、反応容器の入口に逆流して、Sへと変化する割合が高くなる。このように、反応容器反応塔塔内の液体の流れによって、反応成績は異なると考えられる。以上説明した本実施の形態では、液相中に気泡が分散した混相流の解析において、液相に関し、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算された。これにより、液相が気泡に同伴することによる、液相において生じる流れの作用が考慮するためのパラメータβ,BL12を取得できる。したがって、反応をメゾスケールで捉えて高精度にシミュレーションできるようになる。"