知財求人 - 知財ポータルサイト「IP Force」
▶ 住友化学株式会社の特許一覧 ▶ 国立大学法人 鹿児島大学の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】5791105
(24)【登録日】2015年8月14日
(45)【発行日】2015年10月7日
(54)【発明の名称】プロセス解析プログラム
(51)【国際特許分類】
G06F 17/50 20060101AFI20150917BHJP
G06F 19/00 20110101ALI20150917BHJP
【ＦＩ】
G06F17/50 612G
G06F17/50 680Z
G06F19/00 110
【請求項の数】3
【全頁数】31
(21)【出願番号】特願2011-185059(P2011-185059)
(22)【出願日】2011年8月26日
(65)【公開番号】特開2013-45423(P2013-45423A)
(43)【公開日】2013年3月4日
【審査請求日】2014年4月28日
(73)【特許権者】
【識別番号】000002093
【氏名又は名称】住友化学株式会社
(73)【特許権者】
【識別番号】504258527
【氏名又は名称】国立大学法人 鹿児島大学
(74)【代理人】
【識別番号】110001195
【氏名又は名称】特許業務法人深見特許事務所
(72)【発明者】
【氏名】島田 直樹
(72)【発明者】
【氏名】斉木 理奈
(72)【発明者】
【氏名】ダール アビナーブ
(72)【発明者】
【氏名】水田 敬
【審査官】
合田 幸裕
(56)【参考文献】
【文献】
特開２０１１−１５０６８０（ＪＰ，Ａ）
【文献】
特開２０１１−０６０２０５（ＪＰ，Ａ）
【文献】
特開２００６−０８５６０３（ＪＰ，Ａ）
【文献】
特開２０１０−０４００１３（ＪＰ，Ａ）
【文献】
坂口 忠司, 赤對 秀明, 高橋 広光, 冨山 明男, 南川 久人，垂直管内固気液三相気泡流の体積率に関する研究，日本機械学会論文集 B編，日本，社団法人日本機械学会，１９９１年１０月，第57巻／第542号，第3325-3332頁
【文献】
前川 宗則, 島田 直樹, 宋 明良, 冨山 明男，多流体モデルに基づく非均質気泡流の数値予測，日本機械学会論文集 B編，日本，社団法人日本機械学会，２００６年１１月，第72巻／第723号，第2680-2686頁
(58)【調査した分野】（Int.Cl.，ＤＢ名）
Ｇ０６Ｆ １７／５０
Ｇ０６Ｆ １９／００
ＩＥＥＥ Ｘｐｌｏｒｅ
ＣｉＮｉｉ
ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算するステップは、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算するステップと、
前記化学種保存式に基づいてモル濃度を計算するステップとを含み、
前記液相の体積率を計算するステップは、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる前記液相内で前記気泡に同伴される第１液相領域の体積率と、前記液相内の前記第１液相領域以外の領域である第２液相領域の体積率とを計算し、
前記モル濃度を計算するステップは、前記第１液相領域と前記第２液相領域との間での物質交換量を示すパラメータを採用して改良した化学種保存式に基づいて前記第１液相領域のモル濃度と前記第２液相領域のモル濃度とを計算する、プロセス解析プログラム。
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算するステップは、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算するステップと、
前記化学種保存式に基づいてモル濃度を計算するステップとを含み、
前記液相の体積率を計算するステップは、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる前記液相内で前記気泡に同伴される第１液相領域の体積率と、前記液相内の前記第１液相領域以外の領域である第２液相領域の体積率とを計算し、
前記モル濃度を計算するステップは、前記第１液相領域と前記第２液相領域との間での物質交換量を示すパラメータを採用して改良した化学種保存式に基づいて前記第１液相領域のモル濃度と前記第２液相領域のモル濃度とを計算する、プロセス解析プログラム。
【請求項2】
前記パラメータは、０．０１〜１．０の範囲の値に設定される、請求項１に記載のプロセス解析プログラム。
【請求項3】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、
空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算する手段は、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算する手段と、
前記化学種保存式に基づいてモル濃度を計算する手段とを含み、
前記液相の体積率を計算する手段は、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる前記液相内で前記気泡に同伴される第１液相領域の体積率と、前記液相内の前記第１液相領域以外の領域である第２液相領域の体積率とを計算し、
前記モル濃度を計算する手段は、前記第１液相領域と前記第２液相領域との間での物質交換量を示すパラメータを採用して改良した化学種保存式に基づいて前記第１液相領域のモル濃度と前記第２液相領域のモル濃度とを計算する、プロセス解析装置。
空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算する手段は、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算する手段と、
前記化学種保存式に基づいてモル濃度を計算する手段とを含み、
前記液相の体積率を計算する手段は、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる前記液相内で前記気泡に同伴される第１液相領域の体積率と、前記液相内の前記第１液相領域以外の領域である第２液相領域の体積率とを計算し、
前記モル濃度を計算する手段は、前記第１液相領域と前記第２液相領域との間での物質交換量を示すパラメータを採用して改良した化学種保存式に基づいて前記第１液相領域のモル濃度と前記第２液相領域のモル濃度とを計算する、プロセス解析装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、気液混合相におけるプロセス解析プログラムに関する。
【背景技術】
【0002】
化学プラントの酸化、水素化、発酵、廃水処理などのさまざまなプロセスにおいて、気泡塔と呼ばれる装置が用いられている。この気泡塔は、液相中に気泡が分散した「気泡流」を利用した設備であり、他の気液反応装置に比較してより大きな気液界面積が利用でき、また構造を単純化できるため装置の管理が容易になるといった利点を有している。
【0003】
このような気泡塔の解析では、気泡流構造が複雑であるので、実験的手法を用いて大規模な解析を行なうことは困難である。そのため、数値計算（シミュレーション）を用いて、気泡塔の解析を行なうことが注目されている。このようなシミュレーションにより気泡塔を解析する際には、対象の気泡塔についての計算モデルを高度化することが重要となっている。
【0004】
このような状況下において、たとえば非特許文献１〜３や特許文献１には、二流体モデルを用いた種々の計算モデルが提案されている。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特開２０１０−４００１３号公報
【非特許文献】
【0006】
【非特許文献1】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
【非特許文献2】島田他、「気泡塔内気泡流数値予測に必要となる相関式の感度解析」、社団法人 化学工学会、化学工学論文集、第２９巻、第６号、ｐｐ．７７８−７８５、２００３年
【非特許文献3】島田他、「化学反応・ガス吸収・熱輸送を伴う気泡塔内気泡流の数値解法」、社団法人 化学工学会、化学工学論文集、第３１巻、第６号、ｐｐ．３７７−３８７、２００５年
【非特許文献4】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
【非特許文献5】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
【発明の開示】
【発明が解決しようとする課題】
【0007】
しかしながら、従来の計算モデルによるプロセス解析では、系のサイズによっては、シミュレーション結果が実際の反応系を精度よく再現できないという課題があった。
【0008】
具体的には、たとえば上記したような気泡塔が、直径０．４ｍ以上、高さ４ｍ以上のものとなった場合に、シミュレーション結果の実際の反応系からの乖離が顕著に現われていた。
【0009】
本発明は係る実情に鑑み考え出されたものであり、その目的は、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようにすることである。
【課題を解決するための手段】
【0010】
本発明に従ったプロセス解析プログラムは、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算するステップは、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算するステップと、質量保存式および状態保存式に基づいて、液相の体積率を計算するステップと、質量保存式および状態保存式に基づいて、気相の体積率を計算するステップと、化学種保存式に基づいてモル濃度を計算するステップとを含み、液相の体積率を計算するステップは、液相内で気泡の影響を受ける第１液相領域の体積率と、液相内の第１液相領域以外の領域である第２液相領域の体積率とを計算し、モル濃度を計算するステップは、第１液相領域のモル濃度と第２液相領域のモル濃度とを計算する。
【0011】
好ましくは、第１液相領域の体積率は、気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる。
【0012】
好ましくは、パラメータは、０．０１〜１．０の範囲の値に設定される。本発明に従ったプロセス解析装置は、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算する手段は、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算する手段と、質量保存式および状態保存式に基づいて、液相の体積率を計算する手段と、質量保存式および状態保存式に基づいて、気相の体積率を計算する手段と、化学種保存式に基づいてモル濃度を計算する手段とを含み、液相の体積率を計算する手段は、液相内で気泡の影響を受ける第１液相領域の体積率と、液相内の第１液相領域以外の領域である第２液相領域の体積率とを計算し、モル濃度を計算する手段は、第１液相領域のモル濃度と第２液相領域のモル濃度とを計算する。
【発明の効果】
【0013】
本発明によれば、液相中に気泡が分散した混相流を解析する際に、液相において、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算される。
【0014】
これにより、気泡によって液体が流れを変えられる作用を考慮でき、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようになる。
【図面の簡単な説明】
【0015】
【図1】本発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【図2】反応における各成分の濃度の時間変化の一例を示す図である。
【図3】反応における各成分の濃度の時間変化の他の例を示す図である。
【図4】反応における各成分の濃度の時間変化のさらに他の例を示す図である。
【図5】反応の中間生成物の、生成量（溶液中の濃度）の時間変化を示す図である。
【図6】気泡が液体を同伴して流れる状態の概念図である。
【図7】液相の第１液相領域と第２液相領域を説明するための図である。
【図8】本発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図9】本発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。
【図10】気泡流入開始後１秒におけるｐＨの分布を示す図である。
【図11】気泡流入開始後１秒における速度ベクトルの分布を示す図である。
【図12】ｐＨの分布を示す図である。
【図13】速度ベクトルの分布を示す図である。
【発明を実施するための形態】
【0016】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【0017】
＜解析対象＞この発明の実施の形態に従うプロセス解析は、気泡塔などの液相中に気泡が分散した混相流を生じる系を解析対象とする。
【0018】
図１は、この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【0019】
図１を参照して、解析対象の一例として、反応器１０内に原料１（液相）を供給するとともに、原料１中に原料２（気相）が吹き込まれる系ＳＹＳを考える。反応器１０の底部には、原料１を供給するための供給ライン１２が接続されている。この供給ライン１２の経路上には、原料１を供給／遮断するための開閉弁１４が設けられている。なお、図１に示す系ＳＹＳでは、バッチ処理またはセミバッチ処理が行なわれるので、原料１は、当該バッチ処理毎に供給される。
【0020】
また、反応器１０には、気相の原料２を供給するための供給ライン１６が接続されており、この供給ライン１６は、反応器１０内で水平方向に延びる吹き出しライン１８と接続される。吹き出しライン１８には、原料１を液相Ｌ内に吹き込むための孔が垂直上方向に向けて多数設けられている。これにより、気相の原料２が供給されると、原料１の液相中に原料２の気泡が分散し、この気泡が上昇することにより流れが生じる。また、気泡中の反応成分が液相中に溶解することにより、酸化、水素化、発酵、廃水処理などの化学反応が生じる。
【0021】
反応器１０の底部には、化学反応によって生成された製品を取り出すための取り出しライン２０が接続されており、反応器１０内で生じた生成品は、この取り出しライン２０を通じて次工程へ運ばれる。また、この反応器１０の頂点部には、反応器１０内の圧力および液面を調整するため、開閉弁２６を設けた通気ライン２４が接続されている。
【0023】
ａ） 底部から液体が供給されるｂ） 塔内部に空気が供給される
ｃ） 液体（液相）中を、空気が気泡として分散しながら流れる
ｄ） 気泡から液体（液相）中にガスが溶解する
ｅ） 溶解したガスと液体が反応する
本明細書では、これだけの複雑な過程を含む反応系について、まず、ラボの規模で反応系を検討し、そして、ミニプラント、実機へと順にスケールアップさせることとした。つまり、ラボやミニプラントで得た知見を実機の合理化設計へと適切に展開するために、上記現象を様々なスケールで捉え、スケールごとに重要となる因子を議論することにより、スケールアップのあり方を調整することとした。
【0024】
そこで、本明細書では、次の表１に挙げられた、ミクロ、メゾ、マクロの各スケールにおいて検討を行なう。なお、表１では、生じると考えられる主な現象によって、スケールが定義されている。ミクロスケールは、反応機構レベルで現象を捉えるスケールである。マクロスケールは、反応塔内における液混合全体的に現象を捉えるスケールである。メゾは、ミクロとマクロの中間のスケールであり、ガスの溶解や、その際に各気泡が液を同伴して移動するなど、相間の界面での現象を捉えるスケールである。
【0025】
【表1】
【0026】
（ミクロスケール）同じ反応機構を持っていた場合であっても、気泡塔内の液体の流れによって反応成績が変化する。その原因について、以下説明する。
【0027】
反応容器における流体の流れとして、供給した液体が瞬時にして容器全体に混合する流れを「完全混合」と、供給した液体が流出部へと一次元的に移動する「押し出し流れ」とがある。
【0028】
ラボで使用される攪拌槽反応器における流れは、完全混合とみなすことができる。一方、押し出し流れは、反応容器内の液体が全く混合していない状態と捉えることができる。たとえば、管内の流れなどは、押し出し流れの一例と言える。
【0029】
そして、通常の反応容器内の流れは、完全混合と押し出し流れの間に位置づけられる。ここで、次の（ＥＱ）式に示すような、出発物質Ａが中間体Ｒを経て、最終生成物Ｓへと変化する「Ａ→Ｒ→Ｓ」という逐次反応を対象として、Ｒの選択率を検討する。なお、ｋ１，ｋ２は、それぞれの反応の反応速度定数である。
【0030】
【数1】
【0031】
供給されたＡは、反応容器の流出部に進むにつれてＲへと変化する。Ｒの一部は、反応が進んでＳへと変化する。
【0032】
反応容器において混合があると、供給したＡが反応せずに流出部へと流れやすくなる一方で、Ａが反応して生成したＲが入口に逆流して、Ｓへと変化する割合が増える。
【0033】
反応液に限らず、流体の流れの状態を表す指標として混合段数（混合槽数）という値が知られている。この混合段数とは、筒状の胴体を有する容器（以下、適宜「流通器」と称する）の中を一方向に流れる流体が、流通器を流れる間に混合される度合いを示す値であり、上記流通器を、直列に連結された互いに同一体積の複数の完全混合槽にモデル化した場合の、当該完全混合槽の数である。混合段数については、文献Ａ（橋本健治、「反応工学」、改訂版、培風館、１９９３年３月２０日、ｐ．１７９−１８７、１９５―１９７）および文献Ｂ（Octave Levenspiel, “Chemical reaction engineering”, Fifth printing, USA, John Wiley and Sons, INC., May, 1967 p.261-265）にその説明が記載されている。そして、この混合段数が１（最小）の場合が、押し出し流れの状態であり、混合段数が無限大の場合が、完全混合の状態である。
【0034】
（ＥＱ）式の反応における、Ａ成分、Ｏ成分、Ｒ成分、および、Ｓ成分の各濃度についての、従来の技術に従ったシミュレーション結果について、説明する。
【0035】
図２は、ｋ１とｋ２の値が同じである場合の、各成分の濃度の時間変化を示す図である。なお、図２（Ａ）は、反応容器内が均一の場合（つまり、完全混合の場合）を示し、図２（Ｂ）は、反応容器内の不均一の場合（つまり、押し出し流れの場合）を示す。また、各図では、Ａ成分の濃度が実線で、Ｏ成分の濃度が破線で、Ｒ成分の濃度が一点鎖線で、そして、Ｓ成分の濃度が点線で、記載されている。
【0036】
図２（Ｂ）の不均一の場合には、反応することなく残留するＡ成分およびＲ成分がほとんど見られないのに対し、図２（Ａ）の均一の場合には、或る程度の量のＡ成分およびＲ成分が反応することなく残留している。
【0037】
図３は、ｋ１がｋ２よりも大きい場合の、各成分の濃度の時間変化を示す図である。なお、図３（Ａ）は、反応容器内が均一の場合（つまり、完全混合の場合）を示し、図３（Ｂ）は、反応容器内の不均一の場合（つまり、押し出し流れの場合）を示す。図３では、各成分の濃度変化が、図２と同様の線の種類で記載されている。
【0038】
また、図４は、ｋ１がｋ２よりも小さい場合の、各成分の濃度の時間変化を示す図である。なお、図４（Ａ）は、反応容器内が均一の場合（つまり、完全混合の場合）を示し、図４（Ｂ）は、反応容器内の不均一の場合（つまり、押し出し流れの場合）を示す。図４では、各成分の濃度変化が、図２と同様の線の種類で記載されている。
【0041】
図５は、（ＥＱ）式の中間生成物であるＲ成分の、生成量（溶液中の濃度）の時間変化を示す図である。なお、図５（Ａ）は、図２（Ａ），図２（Ｂ），図３（Ａ），図３（Ｂ），図４（Ａ），図４（Ｂ）のそれぞれに示されたＲ成分についての濃度変化を示す。そして、図５（Ｂ）は、図５（Ａ）に示されたそれぞれの状況下でのＲ成分の濃度変化の、反応初期における変化を拡大して（つまり、図５（Ａ）の各グラフを横軸方向に拡大して）、示すものである。
【0042】
図５（Ａ）および図５（Ｂ）に示されるように、中間生成物の生成量は、反応条件によって、大きく異なるといえる。ここでいう反応条件とは、ｋ１とｋ２の間の大小関係と、反応容器における流体の流れ（完全混合／押し出し流れ）との組合せである。また、図５（Ａ）および図５（Ｂ）において、Ｒ成分の生成量の挙動の差異は、反応初期だけでなく、反応中期においても見られている。
【0043】
なお、ｋ１とｋ２の大小関係のいずれであっても（ｋ１＝ｋ２、ｋ１＞ｋ２、または、ｋ１＜ｋ２）、Ｒ成分の濃度は、反応初期には、完全混合よりも押し出し流れの方が高くなっているが、反応中期以降は、完全混合の方が高くなる傾向にある。つまり、同じ反応率において完全混合よりも押し出し流れの方が、Ｓを生成する反応についての選択率が高いと言える。
【0044】
以上のことから、反応系における生成物の生成量を予測するためには、ｋ１とｋ２の間の大小関係とともに、反応容器における流体の流れについても検討することは重要であるといえる。
【0045】
（メゾスケール）上記したように、反応における選択率は反応容器内の流体の流れに影響を受け、反応機構が同じであれば、完全混合よりも選択率の低い流れは存在しないものと考えられる。しかしながら、実機における反応成績には、完全混合よりも選択率が劣る場合が観測されている。
【0046】
ここで、上記したような考察は、反応容器内では液体が単相で流れているという想定に基づいたものであった。これに対して、実際の気泡塔内では、気泡と液体が混在した混相流の状態になっている。このことから、上記したような、完全混合よりも選択率の低い系が存在したという結果は、気泡が液体を同伴して流れたことによる特別な作用に起因するものではないかと考えられる。
【0048】
従来の計算モデルでは、計算セルサイズΔｘに比べて、小さな気泡を分散相とし、また、その大きさに応じてＮ群に分類することにより、各相の体積率αを、（１）式で表していた。
【0049】
【数2】
【0050】
（１）式では、気泡流に含まれる相を１つの連続気相（添え字Ｇ）、１つの連続液相（添え字Ｌ）、そして、Ｎ群の気泡（Ｂｍ：Ｂｍ＝１，２，・・・，Ｎ）の合計（Ｎ＋２）種のグループに分類がなされていた。
【0051】
本実施の形態では、液相において、気泡の影響を受ける可能性を考慮して、体積率の検討において領域の分類を図る。具体的には、図７に示されるように、液相において、気泡の影響を受ける液体ＥＬを第１液相領域とし、気泡の影響を受けにくいと考えられるバルクの液相Ｌを第２液相領域とする。本実施の形態では、（１）式においてαＬで記述されていた液相の体積率は、第１液相領域の体積率αＬ１と第２液相領域の体積率αＬ２とを用いて、次の式のように記述される。
【0052】
【数3】
【0053】
そして、第１液相領域の体積率αＬ１は、それを同伴する気泡の体積率αＢに依存すると考えられる。したがって、各液体ＥＬの体積率αＬ１は、各気泡の体積率αＢを用いて、次の式のように記述される。
【0054】
【数4】
【0055】
ここで、βは、気泡の影響範囲を表すパラメータである。そして、液相の体積率は、次の式で示すように、第１液相領域の体積率αＬ１の総量と第２液相領域の体積率αＬ２の総量の和として記述される。
【0056】
【数5】
【0057】
上記したように、本実施の形態では、計算モデルにおける体積率について、気泡の影響の有無を考慮し、液相の領域を第１液相領域と第２液相領域に分割して検討している。そして、第１液相領域の体積率αＬ１を、気泡の体積率αＢとパラメータβによって、つまり、系のスケールに依存しない値によって、言及している。
【0058】
これにより、本実施の形態では、系のスケールに依存することなく、液相における気泡の影響を考慮できるため、より実際の反応における挙動に沿ったシミュレーションが可能となる。
【0059】
なお、本明細書では、気泡Ｂと区別するために、反応系内のバルクの気相を適宜気相Ｇと表記する。
【0060】
＜ハードウェア構成＞本実施の形態に従うプロセス解析装置は、代表的にコンピュータによって実現される。
【0061】
図８は、この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータ１００の概略構成図である。
【0062】
図８を参照して、コンピュータ１００は、ＦＤ（Flexible Disk）駆動装置１１１およびＣＤ−ＲＯＭ（Compact Disk-Read Only Memory)駆動装置１１３を搭載したコンピュータ本体１０１と、モニタ１０２と、キーボード１０３と、マウス１０４とを含む。
【0063】
コンピュータ本体１０１は、ＦＤ駆動装置１１１およびＣＤ−ＲＯＭ駆動装置１１３に加えて、相互にバスで接続された、演算装置であるＣＰＵ（Central Processing Unit）１０５と、メモリ１０６と、記憶装置である固定ディスク１０７と、通信インターフェース１０９とを含む。
【0064】
本実施の形態に従うプロセス解析処理は、ＣＰＵ１０５がメモリ１０６などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、ＦＤ１１２やＣＤ−ＲＯＭ１１４などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、ＦＤ駆動装置１１１やＣＤ−ＲＯＭ駆動装置１１３などにより記録媒体から読取られて、または通信インターフェース１０９にて受信されて、固定ディスク１０７に格納される。さらに、このようなプログラムは、固定ディスク１０７からメモリ１０６に読出されて、ＣＰＵ１０５により実行される。
【0065】
ＣＰＵ１０５は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ１０６は、ＣＰＵ１０５のプログラム実行に応じて各種の情報を記憶する。
【0066】
モニタ１０２は、ＣＰＵ１０５が出力する情報を表示するための表示部であって、一例としてＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）などから構成される。すなわち、モニタ１０２には、プロセス制御の状態などが表示される。
【0067】
マウス１０４は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード１０３は、入力されるキーに応じたユーザから指令を受付ける。
【0068】
通信インターフェース１０９は、コンピュータ１００と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。
【0069】
なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。
【0070】
また、コンピュータ１００を本実施の形態に従うように機能させるプログラムが記録される記憶媒体は、上記したＦＤやＣＤ−ＲＯＭに限定されない。記憶媒体の他の例としては、ＤＶＤ−ＲＯＭ（Digital Versatile Disk - Read Only Memory）、ＵＳＢ（Universal Serial Bus）メモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、ＭＯ（Magnetic Optical Disc）、ＭＤ（Mini Disc）、ＩＣ（Integrated Circuit）カード（メモリカードを除く）、光カード、マスクＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ（Electronically Erasable Programmable Read-Only Memory）などの、不揮発的にプログラムを格納する媒体が挙げられる。
【0071】
＜計算モデル＞以下、本実施の形態に従うプロセス解析に用いられる計算モデルについて説明する。
【0072】
（１．基礎式）本実施の形態に従うプロセス解析に用いる計算モデルでは、気泡流に含まれる相を１つの連続気相（添え字Ｇ）、液相領域（添え字Ｌ）、Ｎ群の気泡（Ｂｍ，Ｂｍ＝１，２，・・・，Ｎ）の合計（Ｎ＋２）種のグループに分類する。
【0073】
本実施の形態に従うプロセス解析では、各相（（Ｎ＋２）種の各グループ）について、質量保存式、運動量保存式、状態方程式、および、化学種保存式を考えることで、状態値の時間的変化を計算する。
【0074】
以下の計算モデルでは、計算セルサイズΔｘに比べて、気液界面の長さスケールが大きな場合には、その界面で隔てられた二相を連続気相と連続液相（上記した第２液相領域）とみなす。一方、計算セルサイズΔｘに比べて、小さな気泡は分散相とし、さらにその大きさに応じてＮ群に分類する。各相の体積率αは、体積率の定義から上記（１）式を満たす。
【0075】
気液両相は相変化のない等温非圧縮性ニュートン流体と仮定し、気泡内部の分子および乱流拡散を無視すると、セル体積平均に基づいて、（２）式〜（６）式に示す質量保存式および運動量保存式が成立する。
【0076】
【数6】
【0077】
ここで、下付添字のＧは連続気相を、Ｌは連続液相を、Ｂｍはグループｍ（ｍ＝１，２，…，Ｎ）の気泡を意味する。また、ｔは時間、Ｖは速度、ρは密度、ΓＣＢは合体・分裂による単位体積・単位時間当りの質量生成量、ΓＡＥは吸収・蒸発による単位体積・単位時間当りの質量生成量を表す。なお、質量生成量Γの下付添字は輸送に関与する二相を意味し、たとえばＬＢｍは気泡グループ Ｂｍから連続液相Ｌへの輸送を表す。また、ＤＢｍ／ＤｔおよびＤＣ／Ｄｔは、以下の実質微分を意味する。
【0078】
【数7】
【0079】
本実施の形態に従うプロセス解析で用いる計算モデルについても、２つの連続相に対して一流体近似を採用している。そのため、２つの連続相に対する運動量保存式としては、連続相混合体についての（６）式のみを用いる。（６）式における連続相混合体の体積率αＣは、αＣ＝αＧ＋αＬとして与え、混合体の密度ρＣを次の（７）式として与える。
【0080】
【数8】
【0081】
また、本実施の形態に従うプロセス解析で用いる計算モデルでは、化学種の量を評価するパラメータとして、モル濃度・モル分率・分密度・分圧・質量分率・絶対湿度があげられる。ここでは、質量分率Ｙとモル濃度Ｃを用いて、状態方程式を構成する場合を考える。グループｍの気泡のみが存在する場合、Ｙの定義から次の（８）式が成立する。
【0082】
【数9】
【0083】
ここで、下付添字のｉは、化学種を表す。また、ρＢｍｉ（Ｐ，ＴＢｍ）は、純物質時の状態式である（８）式から、平均密度に関して次の（９）式が得られる。
【0084】
【数10】
【0085】
連続液相と連続気相の平均密度も同様の手順で得られる。本実施の形態では、各相（Ｂｍ，Ｌ，Ｇ）を一括して表す下付添字ｋを用い、以下の（１０）式のように表記する。
【0086】
【数11】
【0087】
（１０）式では、添字ｋを用いて表記する便宜上、ＴｃをＴＬ，ＴＧ（一流体近似からＴｃ＝ＴＬ＝ＴＧ）に置き換えている。一方、モル濃度Ｃｋｉと平均密度との関係は、次の（１１）式を採用している。
【0088】
【数12】
【0089】
ここで、Ｍは分子量をあらわす。変数ＹとＣは、（１０）式を用いる方法と（１１）式を用いる方法は等価である。これらの二式を比較した場合、変分が単純となる（１１）式を用いる方が簡便である。また、多くの場合、反応速度はＣに直接依存するため、多くの等温系液相流れの計算でＣが使用されている。しかしながら、圧縮性流体を扱う場合はＣに温度・圧力依存性を加味する必要がある。一方、（１０）式を用いる場合は、ρｋｉ（Ｐ，Ｔｋ）を物性値推算法・データベースなどを用いて設定することで、密度の温度・圧力依存性を加味できる。このため、圧縮性・非圧縮性の区別なくＹの輸送計算ができる利点があり、ρｋｉ（Ｐ，Ｔｋ）が明確なガス燃焼流れなどの計算に良く使用されている。工業的なガス吸収操作では、ガスを溶質として液相側で反応させる場合が極めて多い。本実施の形態では、気相（Ｂｍ，Ｇ）の化学種を質量分率Ｙ、液相（Ｌ）の化学種をモル濃度Ｃで評価し、化学種保存式として、次の（１２）式〜（１４）式を用いる。
【0090】
【数13】
【0091】
なお、（１３）式において、Ｋは気泡（ＢＭ）または液相に溶解している気体（Ｃ）を表す。これにより、（１３）式では、気体成分について、気泡（ＢＭ）の状態および液相内の溶質の状態の運動が考慮されている。
【0092】
ここで、ｗは反応に伴う単位体積当たりのモル生成率、Ｄは拡散係数、ＢＬＢｍｉは気泡グループと連続液相における単位体積・単位時間当たりの化学種ｉの生成量、ΦＢｍｉとΦＧｉはΓＧＢｍに伴う化学種ｉの生成量を表す。（１２）式〜（１４）式の右辺第一項は生成項、（１２）式の第二項と（１４）式の第三項はΓに伴う生成項、（１３）式と（１４）式の第二項は拡散項、（１３）式の第三項は温度・圧力依存項である。
【0093】
なお、本実施の形態では、図７を参照して上記したように、液相を第１液相領域と第２液相領域に分割して検討を行なうため、（１３）式が改良される。改良にあたり、検討を容易にするため、次の仮定を設ける。
【0094】
仮定１：連続気相から液相への物質移動はない仮定２：各モル濃度は圧力の影響を受けない
仮定３：気相からの物質移動は第１液相領域への移動のみが生じ、その後、第１液相領域と第２液相領域の間で交換が生じる
仮定１によれば次の（１５）式が導かれ、仮定２によれば次の（１６）式が導かれる。
【0095】
【数14】
【0096】
また、モル濃度に関する（１３）式を第１液相領域と第２液相領域に分割して検討すると、次の（１７）式および（１８）式が導かれる。
【0097】
【数15】
【0098】
ここで、Ｌｊｉは液相の領域ｊにいる成分ｉを意味している。また、ＢＬ１２ｉは第１液相領域と第２液相領域の間での物質交換量であり、ＢＬ１Ｂｍｉは気相からの物質移動についての物質交換量である。気相からの物質移動が無い成分については、の場合、ＢＬ１Ｂｍｉ＝０となる。したがって、この場合、（１７）式と（１８）式は次の（１９）式と（２０）式のように書ける。
【0099】
【数16】
【0100】
（１７）式〜（２０）式を解くには、本実施の形態において新たに採用したパラメータβとＢＬ１２を指定する必要がある。
【0101】
βは、気泡の形状、上昇速度、密集度、乱れの状況など、多くの流れの作用を受ける複雑なパラメータであるが、既往の研究（例えば、R. Clift , J. R. Grace , M. E. Weber“Bubbles,Drops and Particles”, Dover Publications, 2005）からb = 0.01〜1.0の範囲が考えられる。
【0102】
ＢＬ１２についても、気泡周りに関する詳細な知見から与える必要がある。これらの情報は、以下の研究内容から取得できる。
【0103】
（レーザー誘起蛍光法）例えば、文献（Mizuta, K., Odo, Y. Yoshimitsu, S., Kaji, H., Murakami, M. and Matsumoto, T., “Transient Behavior during a Gas Dissolution Process around a Stationary Single Bubble”, J. Chem. Eng. Japan, 41, 7, pp. 553-556, (2008)）に記載の方法で溶解する物質に及ぼす気泡の影響を調べられる。
【0104】
（シミュレーション）例えば、界面追跡法、境界適合格子法による計算で予測する。
【0105】
（２．構成式）界面運動量輸送ＭＬＢｍＦには、次の（２１）式で表される構成式が用いられる。
【0106】
【数17】
【0107】
また、粘性および乱流拡散項Ｆμには、次の（２２）式で表される構成式が用いられる。
【0108】
【数18】
【0109】
ここで、上付添字Ｔは転置であり、μｃは連続相混合体の粘度である。また、表面張力ＦＳには、次の（２３）式が用いられる。
【0110】
【数19】
【0111】
また、化学種輸送量Ｂおよびこれに伴う質量輸送ΓＡＥには、次の（２４）式〜（２７）式が用いられる。
【0112】
【数20】
【0113】
（３．解法）上述した式を用いて、対象プロセスにおける状態値の時間的変化を計算するための方法について説明する。
【0114】
（３−１．化学種保存式の解法）まず、化学種を求めるための解法を提示する。
【0115】
（１２）式〜（１４）式および（１５）式〜（２０）式を参照して説明した化学種保存式には、対流項、反応項、相間輸送項、拡散項、温度・圧力依存項、体積輸送に伴う項が含まれている。各項の特性時間が大きく異なる場合、同じ時間積分法で各項を計算する方法は非効率である。そこで、各項を相間輸送項（上付数字ＩＮＴ）、反応項（ＲＴ）、温度・圧力依存項（ＰＴ）、それ以降の項（ＯＴ）に分類し別々に取り扱うことで、計算の効率向上を図る。
【0116】
ＹまたはＣを一括してΨとすると、その時間変化は、次の（２８）式で評価できる。
【0117】
【数21】
【0118】
ここで、上付添字（ｎ＋１）は時刻（ｎ＋１）Δｔを表す。（２８）式は、後述する圧力・速度計算時に必要となる。圧力・速度計算後、Ψの時間更新を行なう。時間更新時には、計算不安定化を回避するため、陰解法を用いて相間輸送項による計算を行なった後、次の（２９）式を用いて相間輸送項以外の項で時間更新する。
【0119】
【数22】
【0120】
ここで、上付添字（＊）は、相間輸送項のみを用いて更新する値を意味する。以下に、相間輸送項、反応項、温度・圧力依存項、それ以外の項の評価法を整理する。
【0121】
（３−１−１．相間輸送項）気泡塔内気泡流では、気泡から液相へのガス吸収が最も支配的な相間輸送項である。そこで、簡便のため、以下では相間の化学種ｉの輸送を気泡から液相へのガス吸収に限定する。（１２）式〜（１４）式等を参照して説明したような相間輸送項に後述のモル流束に関する相関式を代入すると、以下の（３０）式および（３１）式のように整理できる。
【0122】
【数23】
【0123】
ここで、ξ，Ψ，ζ，ηは代入操作によって生じるＹBmi，ＣLiの各係数をまとめたものである。時間変化項を評価する際には、（３０）式〜（３１）式の右辺におけるＹBmi，ＣLiを時刻ｎで評価した、次の（３２）式および（３３）式を用いる。
【0124】
【数24】
【0125】
一方、時間更新時には本項による計算不安定化を回避するため、（３０）式〜（３１）式の右辺を陰的に扱った、次の（３４）式および（３５）式を用いる。
【0126】
【数25】
【0127】
（３４）式と（３５）式を行列表示に置換え、さらに、ＹBmi (*)，ＣLi (*)について解くと、次の（３６）式が得られる。
【0128】
【数26】
【0129】
（３−１−２．反応項）ここでは、液相内での反応のみを取扱う。
【0130】
化学反応速度式は非線形連立方程式で、特性時間が運動量輸送式離散化的に用いる時間刻み幅Δｔより短い場合が多い。計算不安定化と計算誤差を回避するため、反応項の離散化には半陰解法を用いる。さらに、ΔｔをＪ分割したΔｔJ（＝Δｔ／Ｊ）を定義し、ｌ＝１，…，Ｊとして、次の（３７）式による計算をＣLi (l+1)に関して繰返す。
【0131】
【数27】
【0132】
ここで、上付添字の（ｎ）は時刻ｎΔｔを、（ｌ）は時刻ｎΔｔ＋ｌΔｔJを表しており、ｌ＝０は時刻ｎΔｔを、ｌ＝Ｊは時刻（ｎ＋１）Δｔを意味する。
【0133】
（３７）式の計算過程においてＣLiが平衡状態に達したとみなされる場合は、Ｊ回以下でも繰返し計算を終了できる。反応項による時間変化は、次の（３８）式で評価する。
【0134】
【数28】
【0135】
（３−１−３．温度・圧力依存項）連続液相のモル濃度に対して温度・圧力に関する変分を取ると、次の（３９）式が得られる。
【0136】
【数29】
【0137】
（δＣLi／δＴC）（n)，（δＣLj／δＰ）(n)は物性値推算法・データベースなどから取得する。（ＤCＴC／Ｄt）(n+1)および（ＤCＰ／Ｄt）(n+1)の評価法に関しては、後述する。Ｙを使用する場合は、本項は存在しない。
【0138】
（３−１−４．上記以外の項）上記以外の項として、拡散項および質量輸送に伴う項がある。これらの項では、次の（４０）式を用い、各変数を陽的に扱う。
【0139】
【数30】
【0140】
（３−２．密度の化学種・温度・圧力依存性を考慮した解法）上記の説明により、化学種の時間変化および時間更新計算法が構築できたので、これを質量・運動量・エネルギ保存式に組込み、全体の計算手順を整理する。
【0141】
密度の化学種・温度・圧力依存性を考慮して、（２）式〜（４）式で示した質量保存式を変形する。
【0142】
（１５）式で表した状態方程式から、密度の変分に関して、次の（４１）式が成り立つ。
【0143】
【数31】
【0144】
そして、（４１）式により、密度のラグランジュ微分は、次の（４２）式のように書くことができる。
【0145】
【数32】
【0146】
ここでは、Ｔcと同様に、Ｄc／Ｄｔの、ＤL／Ｄｔ，ＤG／Ｄｔに置換えられている。そして、モル濃度を用いる場合は、（１１）式から、次の（４３）式のように記述することができる。
【0147】
【数33】
【0148】
また、（２）式で表される気泡の質量保存式と上記（４２）式により、次の（４４）式が得られる。
【0149】
【数34】
【0150】
そして、（４４）式を質量保存式として、解法を構成する。次に、各項を評価する時刻を、（４５）式に基づいて定める。
【0151】
【数35】
【0152】
次に、運動量およびエネルギ保存式の各項を評価する時刻を、（４６）式および（４７）式のように定める。
【0153】
【数36】
【0154】
そして、連続液相と連続気相の基礎式に関しても同様の操作を行なう。ただし、連続液相ではモル濃度を用いているため、（４５）式に対応する式は、次の（４８）式となる。
【0155】
【数37】
【0156】
Tomiyama et al. (1991)は、温度が速い過渡変化を伴う場合、温度および密度を評価する時刻が数値解に有意な差を生じさせることを指摘している。（４５）式〜（４８）式においても密度変化の影響が各式に可能な限り反映されるように、各項を評価する時刻を選択している。
【0157】
圧力のラグランジュ微分項である（ＤkＰ／Ｄｔ）(n+1)は、次の（４９）式で評価する。
【0158】
【数38】
【0159】
Shimada and Tomiyama (2005)では、大きな静水圧差下の気泡流計算を通して、（４９）式の右辺第１項を無視する近似の妥当性を確認している。（２８）式と（４７）式，（４８）式により、時刻ｎの値のみを用いて（ＤkＹki／Ｄｔ）(n+1)，（ＤkＰ／Ｄｔ）(n+1)が求められる。このため、質量保存式（気泡の場合（４４）式）の右辺は既知量となり、エネルギ保存式および化学種保存式を用いずに質量・運動量保存式のみを用いて新しい時刻の速度Ｖk(n+1)・圧力Ｐ(n+1)を計算できる。（１）式を考慮して（４５）式に相当する式をすべての相に関して辺々を加え合せると、以下の混合体に対する質量保存式である（５０）式を得る。
【0160】
【数39】
【0161】
次に、（２１）式を運動量保存式（気泡の場合（４６）式）に代入し、さらに、ＶBm (n+1)およびＶG (n+1)について整理すると、次の（５１）式および（５２）式が得られる。
【0162】
【数40】
【0163】
ここで、Ａm，Ｂm，Ｃm，Ｄは、時刻ｎで評価される項をまとめたものである。そして、（３９）式と（４０）式を行列表示すると、次の（５３）式が得られる。
【0164】
【数41】
【0165】
そして、（５３）式をＶBm(n+1)，ＶC(n+1)について解き、整理すると、次の（５４）式が得られる。
【0166】
【数42】
【0167】
ここで、ρの上付添字（〜）およびＧの上付添字（〜）は、（４１）式の右辺に逆行列を乗じた項を置換えたものである。（５０）式および（５４）式からＶk(n+1)，Ｐ(n+1)を求める手順には、種々の非圧縮性流れの数値解法を利用できる。ここでは、ＳＭＡＣ法（Simplified Marker and Cell method）（Amsden and. Hallow, 1970）に基づく求解手順を整理しておく。
【0168】
まず、圧力修正量δＰを次の（５５）式で定義する。
【0169】
【数43】
【0170】
次に、（５４）式における圧力を既知の圧力Ｐ(n)で置換えた次の（５６）式で、速度の推定値Ｖkの上付添字（〜）を定義する。
【0171】
【数44】
【0172】
（５０）式に、（５４）式〜（５６）式を代入すると、以下のδＰに関するポアソン方程式である（５７）式が得られる。
【0173】
【数45】
【0174】
化学種と同様に温度を時間更新する際も、計算安定性確保のため界面エネルギ輸送項を陰的に評価する必要がある。温度計算の詳細については、文献（Shimada and Tomiyama, 2005）などで開示されている技術を採用することができる。
【0175】
＜各モデルのシミュレーション方法＞図９は、この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。なお、図９に示す処理は、解析対策を空間的に分割して得られるそれぞれのセルについて実行される。
【0176】
図９を参照して、まず、プロセス解析装置は、ステップＳ１００で、各相の密度ρ，体積率α，速度Ｖ，圧力Ｐ，モル濃度Ｃ，質量分率Ｙを初期化する。すなわち、シミュレーション対象となるプロセスの初期状態における値を各変数にセットする。その後、処理は、ステップＳ１０２へ進められる。
【0177】
ステップＳ１０２では、（１０）式および（１１）式に従って、各化学種についての（δρ／δＰ）nを算出する。その後、処理はステップＳ１０４に進められる。
【0178】
ステップＳ１０４では、（２８）式、（４７）式および（４８）式に従って、（ＤkＹki／Ｄｔ）(n+1)，（ＤkＰ／Ｄｔ）(n+1)を計算する。その後、処理はステップＳ１０６へ進められる。
【0179】
ステップＳ１０６では、（５６）式に従って、Ｖkの上付添字（〜）を計算する。その後、処理はステップＳ１０８へ進められる。
【0180】
ステップＳ１０８では、（５７）式に従って、圧力修正量δＰを計算する。その後、処理はステップＳ１１０へ進められる。
【0181】
ステップＳ１１０では、（５４）式および（５５）式に従って、速度Ｖk(n+1)および圧力Ｐ(n+1)を計算する。その後、処理はステップＳ１１２へ進められる。
【0182】
ステップＳ１１２〜ステップＳ１１６では、気泡の体積率αBmn+1と、気相の体積率αGn+1と、液相の体積率αLn+1が、それぞれ計算される。なお、これらの体積率は、（１）式、（４５）式、および、（４８）式に従って求められる。そして、処理はステップＳ１１８へ進められる。
【0183】
ステップＳ１１８では、図６および図７を参照して説明した第１液相領域ＥＬと第２液相領域のそれぞれの体積率αLjn+1が、以下の（５８）式および（５９）式に従って計算される。そして、処理はステップＳ１２０へ進められる。
【0184】
【数46】
【0185】
ステップＳ１２０では、（１７）式〜（２０）式に従って第１液相領域と第２液相領域それぞれのモル濃度ＣLjin+1が計算される。そして、ステップＳ１２２へ処理が進められる。
【0186】
ステップＳ１２２では、各化学種の質量分率Ｙkn+1が計算される。そして、ステップＳ１２４へ処理が進められる。
【0187】
ステップＳ１２４では、（１０）式および（１１）式に従って、密度ρkn+1が計算される。その後、処理はステップＳ１２６へ進められる。
【0188】
ステップＳ１２６では、すべての計算処理が終了したか否かが判断される。すべての計算処理が終了していなければ（ステップＳ１２６においてＮＯ）、時刻ｎに１が加算され（ステップＳ１２８）、次の時刻における各パラメータを算出するために処理はステップＳ１０２に戻される。そして、ステップＳ１０２以降の処理が再度実行される。一方、すべての計算処理が終了していれば（ステップＳ１２８においてＹＥＳ）、処理は終了する。
【0189】
このように、解析対象を空間的に分離して得られるそれぞれのセルについて、状態値が時系列に計算される。
【0190】
＜検証実験＞本願発明者らは、模擬装置を用いて、本実施の形態におけるパラメータβ，ＢＬ１２を取得するためのシミュレーションについての検証を行なった。
【0191】
検証は、水のＳＯ２ガス吸収操作を対象としている。図１０は、気泡流入開始後１秒におけるｐＨの分布を示す図である。なお、図１０（Ａ）は実験結果であり、図１０（Ｂ）は、本実施の形態に従うプロセス解析方法による計算結果である。そして、図１０（Ａ）および図１０（Ｂ）では、ｐＨの値がハッチングの違いで表現されている。具体的には、図中、ｐＨ４．４からｐＨ５．２までのｐＨの変化が、Ｐ１〜Ｐ６の６段階で示されている。これらは、Ｐ１がｐＨ４．４に対応し、Ｐ６がｐＨ５．２に対応する。また、図１１は、気泡流入開始後１秒における速度ベクトルの分布を示す図である。なお、図１１（Ａ）は実験結果であり、図１１（Ｂ）は、本実施の形態に従うプロセス解析方法による計算結果である。
【0193】
図１２および図１３に、気泡流入開始後２秒における実験結果および計算結果を示す。なお、図１２は、ｐＨの分布を示す図である。図１２（Ａ）は実験結果であり、図１２（Ｂ）は、本実施の形態に従うプロセス解析方法による計算結果である。図１２（Ａ）および図１２（Ｂ）においても、ｐＨの変化は、図１０（Ａ）および図１０（Ｂ）と同様に、６段階のハッチングで示されている。また、図１３は、速度ベクトルの分布を示す図である。なお、図１３（Ａ）は実験結果であり、図１３（Ｂ）は、本実施の形態に従うプロセス解析方法による計算結果である。
【0194】
図１２および図１３を参照して、気泡流入開始後２秒では、実験結果において、気泡周辺に速度の高い領域が見られ、この傾向は、計算結果でも見られている。このことから、計算結果は、実験結果に沿ったものと言える。
【0195】
＜実施の形態の効果＞たとえば、混相流における「Ａ→Ｒ→Ｓ」の逐次反応を対象としてＲの選択率を考えた場合、反応容器に供給されたＡは、当該反応容器の流出部に進むにつれてＲへと変化する。そして、Ｒの一部は反応が進んでＳへと変化する。
【0196】
このような反応系において、反応容器内で混合があると、供給したＡが反応せず（Ｒへと変化することなく）に流出部へと流れやすくなる。これと同時に、Ａが反応することによって生じたＲのうち、反応容器の入口に逆流して、Ｓへと変化する割合が高くなる。
【0197】
このように、反応容器（反応塔塔）内の液体の流れによって、反応成績は異なると考えられる。
【0198】
以上説明した本実施の形態では、液相中に気泡が分散した混相流の解析において、液相に関し、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算された。
【0199】
これにより、液相が気泡に同伴することによる、液相において生じる流れの作用が考慮するためのパラメータβ，ＢＬ１２を取得できる。したがって、反応をメゾスケールで捉えて高精度にシミュレーションできるようになる。
【符号の説明】
【0200】
１００ コンピュータ、１０１ コンピュータ本体、１０２ モニタ、１０３ キーボード、１０４ マウス、１０６ メモリ、１０７ 固定ディスク、１０９ 通信インターフェース、１１１ ＦＤ駆動装置、１１３ ＣＤ−ＲＯＭ駆動装置。