プロセス解析プログラム
【課題】気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようにする。
【解決手段】本実施の形態のシミュレーション装置は、液相中に気泡が分散した混相流の解析において利用される。当該解析では、液相が、気泡に同伴する領域とそれ以外の領域とに分類される。そして、各分類について、体積率およびモル濃度が計算される。
【解決手段】本実施の形態のシミュレーション装置は、液相中に気泡が分散した混相流の解析において利用される。当該解析では、液相が、気泡に同伴する領域とそれ以外の領域とに分類される。そして、各分類について、体積率およびモル濃度が計算される。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、気液混合相におけるプロセス解析プログラムに関する。
【背景技術】
【0002】
化学プラントの酸化、水素化、発酵、廃水処理などのさまざまなプロセスにおいて、気泡塔と呼ばれる装置が用いられている。この気泡塔は、液相中に気泡が分散した「気泡流」を利用した設備であり、他の気液反応装置に比較してより大きな気液界面積が利用でき、また構造を単純化できるため装置の管理が容易になるといった利点を有している。
【0003】
このような気泡塔の解析では、気泡流構造が複雑であるので、実験的手法を用いて大規模な解析を行なうことは困難である。そのため、数値計算(シミュレーション)を用いて、気泡塔の解析を行なうことが注目されている。このようなシミュレーションにより気泡塔を解析する際には、対象の気泡塔についての計算モデルを高度化することが重要となっている。
【0004】
このような状況下において、たとえば非特許文献1〜3や特許文献1には、二流体モデルを用いた種々の計算モデルが提案されている。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特開2010−40013号公報
【非特許文献】
【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】島田他、「気泡塔内気泡流数値予測に必要となる相関式の感度解析」、社団法人 化学工学会、化学工学論文集、第29巻、第6号、pp.778−785、2003年
【非特許文献3】島田他、「化学反応・ガス吸収・熱輸送を伴う気泡塔内気泡流の数値解法」、社団法人 化学工学会、化学工学論文集、第31巻、第6号、pp.377−387、2005年
【非特許文献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】
具体的には、たとえば上記したような気泡塔が、直径0.4m以上、高さ4m以上のものとなった場合に、シミュレーション結果の実際の反応系からの乖離が顕著に現われていた。
【0009】
本発明は係る実情に鑑み考え出されたものであり、その目的は、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようにすることである。
【課題を解決するための手段】
【0010】
本発明に従ったプロセス解析プログラムは、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算するステップは、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算するステップと、質量保存式および状態保存式に基づいて、液相の体積率を計算するステップと、質量保存式および状態保存式に基づいて、気相の体積率を計算するステップと、化学種保存式に基づいてモル濃度を計算するステップとを含み、液相の体積率を計算するステップは、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算するステップは、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。
【0011】
好ましくは、第1液相領域の体積率は、気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる。
【0012】
好ましくは、パラメータは、0.01〜1.0の範囲の値に設定される。
本発明に従ったプロセス解析装置は、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算する手段は、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算する手段と、質量保存式および状態保存式に基づいて、液相の体積率を計算する手段と、質量保存式および状態保存式に基づいて、気相の体積率を計算する手段と、化学種保存式に基づいてモル濃度を計算する手段とを含み、液相の体積率を計算する手段は、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算する手段は、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。
【発明の効果】
【0013】
本発明によれば、液相中に気泡が分散した混相流を解析する際に、液相において、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算される。
【0014】
これにより、気泡によって液体が流れを変えられる作用を考慮でき、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようになる。
【図面の簡単な説明】
【0015】
【図1】本発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【図2】反応における各成分の濃度の時間変化の一例を示す図である。
【図3】反応における各成分の濃度の時間変化の他の例を示す図である。
【図4】反応における各成分の濃度の時間変化のさらに他の例を示す図である。
【図5】反応の中間生成物の、生成量(溶液中の濃度)の時間変化を示す図である。
【図6】気泡が液体を同伴して流れる状態の概念図である。
【図7】液相の第1液相領域と第2液相領域を説明するための図である。
【図8】本発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図9】本発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。
【図10】気泡流入開始後1秒におけるpHの分布を示す図である。
【図11】気泡流入開始後1秒における速度ベクトルの分布を示す図である。
【図12】pHの分布を示す図である。
【図13】速度ベクトルの分布を示す図である。
【発明を実施するための形態】
【0016】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【0017】
<解析対象>
この発明の実施の形態に従うプロセス解析は、気泡塔などの液相中に気泡が分散した混相流を生じる系を解析対象とする。
【0018】
図1は、この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【0019】
図1を参照して、解析対象の一例として、反応器10内に原料1(液相)を供給するとともに、原料1中に原料2(気相)が吹き込まれる系SYSを考える。反応器10の底部には、原料1を供給するための供給ライン12が接続されている。この供給ライン12の経路上には、原料1を供給/遮断するための開閉弁14が設けられている。なお、図1に示す系SYSでは、バッチ処理またはセミバッチ処理が行なわれるので、原料1は、当該バッチ処理毎に供給される。
【0020】
また、反応器10には、気相の原料2を供給するための供給ライン16が接続されており、この供給ライン16は、反応器10内で水平方向に延びる吹き出しライン18と接続される。吹き出しライン18には、原料1を液相L内に吹き込むための孔が垂直上方向に向けて多数設けられている。これにより、気相の原料2が供給されると、原料1の液相中に原料2の気泡が分散し、この気泡が上昇することにより流れが生じる。また、気泡中の反応成分が液相中に溶解することにより、酸化、水素化、発酵、廃水処理などの化学反応が生じる。
【0021】
反応器10の底部には、化学反応によって生成された製品を取り出すための取り出しライン20が接続されており、反応器10内で生じた生成品は、この取り出しライン20を通じて次工程へ運ばれる。また、この反応器10の頂点部には、反応器10内の圧力および液面を調整するため、開閉弁26を設けた通気ライン24が接続されている。
【0022】
<本発明の概要>
図1に示したような気泡塔内では、たとえば以下のa)〜e)に示すような現象が起こっていると考えられる。
【0023】
a) 底部から液体が供給される
b) 塔内部に空気が供給される
c) 液体(液相)中を、空気が気泡として分散しながら流れる
d) 気泡から液体(液相)中にガスが溶解する
e) 溶解したガスと液体が反応する
本明細書では、これだけの複雑な過程を含む反応系について、まず、ラボの規模で反応系を検討し、そして、ミニプラント、実機へと順にスケールアップさせることとした。つまり、ラボやミニプラントで得た知見を実機の合理化設計へと適切に展開するために、上記現象を様々なスケールで捉え、スケールごとに重要となる因子を議論することにより、スケールアップのあり方を調整することとした。
【0024】
そこで、本明細書では、次の表1に挙げられた、ミクロ、メゾ、マクロの各スケールにおいて検討を行なう。なお、表1では、生じると考えられる主な現象によって、スケールが定義されている。ミクロスケールは、反応機構レベルで現象を捉えるスケールである。マクロスケールは、反応塔内における液混合全体的に現象を捉えるスケールである。メゾは、ミクロとマクロの中間のスケールであり、ガスの溶解や、その際に各気泡が液を同伴して移動するなど、相間の界面での現象を捉えるスケールである。
【0025】
【表1】
【0026】
(ミクロスケール)
同じ反応機構を持っていた場合であっても、気泡塔内の液体の流れによって反応成績が変化する。その原因について、以下説明する。
【0027】
反応容器における流体の流れとして、供給した液体が瞬時にして容器全体に混合する流れを「完全混合」と、供給した液体が流出部へと一次元的に移動する「押し出し流れ」とがある。
【0028】
ラボで使用される攪拌槽反応器における流れは、完全混合とみなすことができる。
一方、押し出し流れは、反応容器内の液体が全く混合していない状態と捉えることができる。たとえば、管内の流れなどは、押し出し流れの一例と言える。
【0029】
そして、通常の反応容器内の流れは、完全混合と押し出し流れの間に位置づけられる。
ここで、次の(EQ)式に示すような、出発物質Aが中間体Rを経て、最終生成物Sへと変化する「A→R→S」という逐次反応を対象として、Rの選択率を検討する。なお、k1,k2は、それぞれの反応の反応速度定数である。
【0030】
【数1】
【0031】
供給されたAは、反応容器の流出部に進むにつれてRへと変化する。Rの一部は、反応が進んでSへと変化する。
【0032】
反応容器において混合があると、供給したAが反応せずに流出部へと流れやすくなる一方で、Aが反応して生成したRが入口に逆流して、Sへと変化する割合が増える。
【0033】
反応液に限らず、流体の流れの状態を表す指標として混合段数(混合槽数)という値が知られている。この混合段数とは、筒状の胴体を有する容器(以下、適宜「流通器」と称する)の中を一方向に流れる流体が、流通器を流れる間に混合される度合いを示す値であり、上記流通器を、直列に連結された互いに同一体積の複数の完全混合槽にモデル化した場合の、当該完全混合槽の数である。混合段数については、文献A(橋本健治、「反応工学」、改訂版、培風館、1993年3月20日、p.179−187、195―197)および文献B(Octave Levenspiel, “Chemical reaction engineering”, Fifth printing, USA, John Wiley and Sons, INC., May, 1967 p.261-265)にその説明が記載されている。そして、この混合段数が1(最小)の場合が、押し出し流れの状態であり、混合段数が無限大の場合が、完全混合の状態である。
【0034】
(EQ)式の反応における、A成分、O成分、R成分、および、S成分の各濃度についての、従来の技術に従ったシミュレーション結果について、説明する。
【0035】
図2は、k1とk2の値が同じである場合の、各成分の濃度の時間変化を示す図である。なお、図2(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図2(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。また、各図では、A成分の濃度が実線で、O成分の濃度が破線で、R成分の濃度が一点鎖線で、そして、S成分の濃度が点線で、記載されている。
【0036】
図2(B)の不均一の場合には、反応することなく残留するA成分およびR成分がほとんど見られないのに対し、図2(A)の均一の場合には、或る程度の量のA成分およびR成分が反応することなく残留している。
【0037】
図3は、k1がk2よりも大きい場合の、各成分の濃度の時間変化を示す図である。なお、図3(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図3(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。図3では、各成分の濃度変化が、図2と同様の線の種類で記載されている。
【0038】
また、図4は、k1がk2よりも小さい場合の、各成分の濃度の時間変化を示す図である。なお、図4(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図4(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。図4では、各成分の濃度変化が、図2と同様の線の種類で記載されている。
【0039】
図3と図4に示されたいずれの場合においても、反応容器内が均一である場合と不均一である場合とでは、各成分の濃度変化の挙動に差異を生じている。
【0040】
つまり、図2〜図4から理解されるように、反応系における各成分の濃度の時間変化は、k1とk2の間の大小関係に関わらず、反応容器内が均一であるか不均一であるかによって相違することとなる。
【0041】
図5は、(EQ)式の中間生成物であるR成分の、生成量(溶液中の濃度)の時間変化を示す図である。なお、図5(A)は、図2(A),図2(B),図3(A),図3(B),図4(A),図4(B)のそれぞれに示されたR成分についての濃度変化を示す。そして、図5(B)は、図5(A)に示されたそれぞれの状況下でのR成分の濃度変化の、反応初期における変化を拡大して(つまり、図5(A)の各グラフを横軸方向に拡大して)、示すものである。
【0042】
図5(A)および図5(B)に示されるように、中間生成物の生成量は、反応条件によって、大きく異なるといえる。ここでいう反応条件とは、k1とk2の間の大小関係と、反応容器における流体の流れ(完全混合/押し出し流れ)との組合せである。また、図5(A)および図5(B)において、R成分の生成量の挙動の差異は、反応初期だけでなく、反応中期においても見られている。
【0043】
なお、k1とk2の大小関係のいずれであっても(k1=k2、k1>k2、または、k1<k2)、R成分の濃度は、反応初期には、完全混合よりも押し出し流れの方が高くなっているが、反応中期以降は、完全混合の方が高くなる傾向にある。つまり、同じ反応率において完全混合よりも押し出し流れの方が、Sを生成する反応についての選択率が高いと言える。
【0044】
以上のことから、反応系における生成物の生成量を予測するためには、k1とk2の間の大小関係とともに、反応容器における流体の流れについても検討することは重要であるといえる。
【0045】
(メゾスケール)
上記したように、反応における選択率は反応容器内の流体の流れに影響を受け、反応機構が同じであれば、完全混合よりも選択率の低い流れは存在しないものと考えられる。しかしながら、実機における反応成績には、完全混合よりも選択率が劣る場合が観測されている。
【0046】
ここで、上記したような考察は、反応容器内では液体が単相で流れているという想定に基づいたものであった。これに対して、実際の気泡塔内では、気泡と液体が混在した混相流の状態になっている。このことから、上記したような、完全混合よりも選択率の低い系が存在したという結果は、気泡が液体を同伴して流れたことによる特別な作用に起因するものではないかと考えられる。
【0047】
図6は、気泡が液体を同伴して流れる状態の概念図である。図6では、液相L内を、気泡Bが、液体ELを同伴して、上方へ流れている状態が示されている。
【0048】
従来の計算モデルでは、計算セルサイズΔxに比べて、小さな気泡を分散相とし、また、その大きさに応じてN群に分類することにより、各相の体積率αを、(1)式で表していた。
【0049】
【数2】
【0050】
(1)式では、気泡流に含まれる相を1つの連続気相(添え字G)、1つの連続液相(添え字L)、そして、N群の気泡(Bm:Bm=1,2,・・・,N)の合計(N+2)種のグループに分類がなされていた。
【0051】
本実施の形態では、液相において、気泡の影響を受ける可能性を考慮して、体積率の検討において領域の分類を図る。具体的には、図7に示されるように、液相において、気泡の影響を受ける液体ELを第1液相領域とし、気泡の影響を受けにくいと考えられるバルクの液相Lを第2液相領域とする。本実施の形態では、(1)式においてαLで記述されていた液相の体積率は、第1液相領域の体積率αL1と第2液相領域の体積率αL2とを用いて、次の式のように記述される。
【0052】
【数3】
【0053】
そして、第1液相領域の体積率αL1は、それを同伴する気泡の体積率αBに依存すると考えられる。したがって、各液体ELの体積率αL1は、各気泡の体積率αBを用いて、次の式のように記述される。
【0054】
【数4】
【0055】
ここで、βは、気泡の影響範囲を表すパラメータである。
そして、液相の体積率は、次の式で示すように、第1液相領域の体積率αL1の総量と第2液相領域の体積率αL2の総量の和として記述される。
【0056】
【数5】
【0057】
上記したように、本実施の形態では、計算モデルにおける体積率について、気泡の影響の有無を考慮し、液相の領域を第1液相領域と第2液相領域に分割して検討している。そして、第1液相領域の体積率αL1を、気泡の体積率αBとパラメータβによって、つまり、系のスケールに依存しない値によって、言及している。
【0058】
これにより、本実施の形態では、系のスケールに依存することなく、液相における気泡の影響を考慮できるため、より実際の反応における挙動に沿ったシミュレーションが可能となる。
【0059】
なお、本明細書では、気泡Bと区別するために、反応系内のバルクの気相を適宜気相Gと表記する。
【0060】
<ハードウェア構成>
本実施の形態に従うプロセス解析装置は、代表的にコンピュータによって実現される。
【0061】
図8は、この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータ100の概略構成図である。
【0062】
図8を参照して、コンピュータ100は、FD(Flexible Disk)駆動装置111およびCD−ROM(Compact Disk-Read Only Memory)駆動装置113を搭載したコンピュータ本体101と、モニタ102と、キーボード103と、マウス104とを含む。
【0063】
コンピュータ本体101は、FD駆動装置111およびCD−ROM駆動装置113に加えて、相互にバスで接続された、演算装置であるCPU(Central Processing Unit)105と、メモリ106と、記憶装置である固定ディスク107と、通信インターフェース109とを含む。
【0064】
本実施の形態に従うプロセス解析処理は、CPU105がメモリ106などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、FD112やCD−ROM114などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、FD駆動装置111やCD−ROM駆動装置113などにより記録媒体から読取られて、または通信インターフェース109にて受信されて、固定ディスク107に格納される。さらに、このようなプログラムは、固定ディスク107からメモリ106に読出されて、CPU105により実行される。
【0065】
CPU105は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ106は、CPU105のプログラム実行に応じて各種の情報を記憶する。
【0066】
モニタ102は、CPU105が出力する情報を表示するための表示部であって、一例としてLCD(Liquid Crystal Display)やCRT(Cathode Ray Tube)などから構成される。すなわち、モニタ102には、プロセス制御の状態などが表示される。
【0067】
マウス104は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード103は、入力されるキーに応じたユーザから指令を受付ける。
【0068】
通信インターフェース109は、コンピュータ100と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。
【0069】
なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。
【0070】
また、コンピュータ100を本実施の形態に従うように機能させるプログラムが記録される記憶媒体は、上記したFDやCD−ROMに限定されない。記憶媒体の他の例としては、DVD−ROM(Digital Versatile Disk - Read Only Memory)、USB(Universal Serial Bus)メモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、MO(Magnetic Optical Disc)、MD(Mini Disc)、IC(Integrated Circuit)カード(メモリカードを除く)、光カード、マスクROM、EPROM、EEPROM(Electronically Erasable Programmable Read-Only Memory)などの、不揮発的にプログラムを格納する媒体が挙げられる。
【0071】
<計算モデル>
以下、本実施の形態に従うプロセス解析に用いられる計算モデルについて説明する。
【0072】
(1.基礎式)
本実施の形態に従うプロセス解析に用いる計算モデルでは、気泡流に含まれる相を1つの連続気相(添え字G)、液相領域(添え字L)、N群の気泡(Bm,Bm=1,2,・・・,N)の合計(N+2)種のグループに分類する。
【0073】
本実施の形態に従うプロセス解析では、各相((N+2)種の各グループ)について、質量保存式、運動量保存式、状態方程式、および、化学種保存式を考えることで、状態値の時間的変化を計算する。
【0074】
以下の計算モデルでは、計算セルサイズΔxに比べて、気液界面の長さスケールが大きな場合には、その界面で隔てられた二相を連続気相と連続液相(上記した第2液相領域)とみなす。一方、計算セルサイズΔxに比べて、小さな気泡は分散相とし、さらにその大きさに応じてN群に分類する。各相の体積率αは、体積率の定義から上記(1)式を満たす。
【0075】
気液両相は相変化のない等温非圧縮性ニュートン流体と仮定し、気泡内部の分子および乱流拡散を無視すると、セル体積平均に基づいて、(2)式〜(6)式に示す質量保存式および運動量保存式が成立する。
【0076】
【数6】
【0077】
ここで、下付添字のGは連続気相を、Lは連続液相を、Bmはグループm(m=1,2,…,N)の気泡を意味する。また、tは時間、Vは速度、ρは密度、ΓCBは合体・分裂による単位体積・単位時間当りの質量生成量、ΓAEは吸収・蒸発による単位体積・単位時間当りの質量生成量を表す。なお、質量生成量Γの下付添字は輸送に関与する二相を意味し、たとえばLBmは気泡グループ Bmから連続液相Lへの輸送を表す。また、DBm/DtおよびDC/Dtは、以下の実質微分を意味する。
【0078】
【数7】
【0079】
本実施の形態に従うプロセス解析で用いる計算モデルについても、2つの連続相に対して一流体近似を採用している。そのため、2つの連続相に対する運動量保存式としては、連続相混合体についての(6)式のみを用いる。(6)式における連続相混合体の体積率αCは、αC=αG+αLとして与え、混合体の密度ρCを次の(7)式として与える。
【0080】
【数8】
【0081】
また、本実施の形態に従うプロセス解析で用いる計算モデルでは、化学種の量を評価するパラメータとして、モル濃度・モル分率・分密度・分圧・質量分率・絶対湿度があげられる。ここでは、質量分率Yとモル濃度Cを用いて、状態方程式を構成する場合を考える。グループmの気泡のみが存在する場合、Yの定義から次の(8)式が成立する。
【0082】
【数9】
【0083】
ここで、下付添字のiは、化学種を表す。また、ρBmi(P,TBm)は、純物質時の状態式である(8)式から、平均密度に関して次の(9)式が得られる。
【0084】
【数10】
【0085】
連続液相と連続気相の平均密度も同様の手順で得られる。本実施の形態では、各相(Bm,L,G)を一括して表す下付添字kを用い、以下の(10)式のように表記する。
【0086】
【数11】
【0087】
(10)式では、添字kを用いて表記する便宜上、TcをTL,TG(一流体近似からTc=TL=TG)に置き換えている。一方、モル濃度Ckiと平均密度との関係は、次の(11)式を採用している。
【0088】
【数12】
【0089】
ここで、Mは分子量をあらわす。変数YとCは、(10)式を用いる方法と(11)式を用いる方法は等価である。これらの二式を比較した場合、変分が単純となる(11)式を用いる方が簡便である。また、多くの場合、反応速度はCに直接依存するため、多くの等温系液相流れの計算でCが使用されている。しかしながら、圧縮性流体を扱う場合はCに温度・圧力依存性を加味する必要がある。一方、(10)式を用いる場合は、ρki(P,Tk)を物性値推算法・データベースなどを用いて設定することで、密度の温度・圧力依存性を加味できる。このため、圧縮性・非圧縮性の区別なくYの輸送計算ができる利点があり、ρki(P,Tk)が明確なガス燃焼流れなどの計算に良く使用されている。工業的なガス吸収操作では、ガスを溶質として液相側で反応させる場合が極めて多い。本実施の形態では、気相(Bm,G)の化学種を質量分率Y、液相(L)の化学種をモル濃度Cで評価し、化学種保存式として、次の(12)式〜(14)式を用いる。
【0090】
【数13】
【0091】
なお、(13)式において、Kは気泡(BM)または液相に溶解している気体(C)を表す。これにより、(13)式では、気体成分について、気泡(BM)の状態および液相内の溶質の状態の運動が考慮されている。
【0092】
ここで、wは反応に伴う単位体積当たりのモル生成率、Dは拡散係数、BLBmiは気泡グループと連続液相における単位体積・単位時間当たりの化学種iの生成量、ΦBmiとΦGiはΓGBmに伴う化学種iの生成量を表す。(12)式〜(14)式の右辺第一項は生成項、(12)式の第二項と(14)式の第三項はΓに伴う生成項、(13)式と(14)式の第二項は拡散項、(13)式の第三項は温度・圧力依存項である。
【0093】
なお、本実施の形態では、図7を参照して上記したように、液相を第1液相領域と第2液相領域に分割して検討を行なうため、(13)式が改良される。改良にあたり、検討を容易にするため、次の仮定を設ける。
【0094】
仮定1:連続気相から液相への物質移動はない
仮定2:各モル濃度は圧力の影響を受けない
仮定3:気相からの物質移動は第1液相領域への移動のみが生じ、その後、第1液相領域と第2液相領域の間で交換が生じる
仮定1によれば次の(15)式が導かれ、仮定2によれば次の(16)式が導かれる。
【0095】
【数14】
【0096】
また、モル濃度に関する(13)式を第1液相領域と第2液相領域に分割して検討すると、次の(17)式および(18)式が導かれる。
【0097】
【数15】
【0098】
ここで、Ljiは液相の領域jにいる成分iを意味している。また、BL12iは第1液相領域と第2液相領域の間での物質交換量であり、BL1Bmiは気相からの物質移動についての物質交換量である。気相からの物質移動が無い成分については、の場合、BL1Bmi=0となる。したがって、この場合、(17)式と(18)式は次の(19)式と(20)式のように書ける。
【0099】
【数16】
【0100】
(17)式〜(20)式を解くには、本実施の形態において新たに採用したパラメータβとBL12を指定する必要がある。
【0101】
βは、気泡の形状、上昇速度、密集度、乱れの状況など、多くの流れの作用を受ける複雑なパラメータであるが、既往の研究(例えば、R. Clift , J. R. Grace , M. E. Weber“Bubbles,Drops and Particles”, Dover Publications, 2005)からb = 0.01〜1.0の範囲が考えられる。
【0102】
BL12についても、気泡周りに関する詳細な知見から与える必要がある。これらの情報は、以下の研究内容から取得できる。
【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】
(2.構成式)
界面運動量輸送MLBmFには、次の(21)式で表される構成式が用いられる。
【0106】
【数17】
【0107】
また、粘性および乱流拡散項Fμには、次の(22)式で表される構成式が用いられる。
【0108】
【数18】
【0109】
ここで、上付添字Tは転置であり、μcは連続相混合体の粘度である。
また、表面張力FSには、次の(23)式が用いられる。
【0110】
【数19】
【0111】
また、化学種輸送量Bおよびこれに伴う質量輸送ΓAEには、次の(24)式〜(27)式が用いられる。
【0112】
【数20】
【0113】
(3.解法)
上述した式を用いて、対象プロセスにおける状態値の時間的変化を計算するための方法について説明する。
【0114】
(3−1.化学種保存式の解法)
まず、化学種を求めるための解法を提示する。
【0115】
(12)式〜(14)式および(15)式〜(20)式を参照して説明した化学種保存式には、対流項、反応項、相間輸送項、拡散項、温度・圧力依存項、体積輸送に伴う項が含まれている。各項の特性時間が大きく異なる場合、同じ時間積分法で各項を計算する方法は非効率である。そこで、各項を相間輸送項(上付数字INT)、反応項(RT)、温度・圧力依存項(PT)、それ以降の項(OT)に分類し別々に取り扱うことで、計算の効率向上を図る。
【0116】
YまたはCを一括してΨとすると、その時間変化は、次の(28)式で評価できる。
【0117】
【数21】
【0118】
ここで、上付添字(n+1)は時刻(n+1)Δtを表す。
(28)式は、後述する圧力・速度計算時に必要となる。圧力・速度計算後、Ψの時間更新を行なう。時間更新時には、計算不安定化を回避するため、陰解法を用いて相間輸送項による計算を行なった後、次の(29)式を用いて相間輸送項以外の項で時間更新する。
【0119】
【数22】
【0120】
ここで、上付添字(*)は、相間輸送項のみを用いて更新する値を意味する。
以下に、相間輸送項、反応項、温度・圧力依存項、それ以外の項の評価法を整理する。
【0121】
(3−1−1.相間輸送項)
気泡塔内気泡流では、気泡から液相へのガス吸収が最も支配的な相間輸送項である。そこで、簡便のため、以下では相間の化学種iの輸送を気泡から液相へのガス吸収に限定する。(12)式〜(14)式等を参照して説明したような相間輸送項に後述のモル流束に関する相関式を代入すると、以下の(30)式および(31)式のように整理できる。
【0122】
【数23】
【0123】
ここで、ξ,Ψ,ζ,ηは代入操作によって生じるYBmi,CLiの各係数をまとめたものである。時間変化項を評価する際には、(30)式〜(31)式の右辺におけるYBmi,CLiを時刻nで評価した、次の(32)式および(33)式を用いる。
【0124】
【数24】
【0125】
一方、時間更新時には本項による計算不安定化を回避するため、(30)式〜(31)式の右辺を陰的に扱った、次の(34)式および(35)式を用いる。
【0126】
【数25】
【0127】
(34)式と(35)式を行列表示に置換え、さらに、YBmi (*),CLi (*)について解くと、次の(36)式が得られる。
【0128】
【数26】
【0129】
(3−1−2.反応項)
ここでは、液相内での反応のみを取扱う。
【0130】
化学反応速度式は非線形連立方程式で、特性時間が運動量輸送式離散化的に用いる時間刻み幅Δtより短い場合が多い。計算不安定化と計算誤差を回避するため、反応項の離散化には半陰解法を用いる。さらに、ΔtをJ分割したΔtJ(=Δt/J)を定義し、l=1,…,Jとして、次の(37)式による計算をCLi (l+1)に関して繰返す。
【0131】
【数27】
【0132】
ここで、上付添字の(n)は時刻nΔtを、(l)は時刻nΔt+lΔtJを表しており、l=0は時刻nΔtを、l=Jは時刻(n+1)Δtを意味する。
【0133】
(37)式の計算過程においてCLiが平衡状態に達したとみなされる場合は、J回以下でも繰返し計算を終了できる。反応項による時間変化は、次の(38)式で評価する。
【0134】
【数28】
【0135】
(3−1−3.温度・圧力依存項)
連続液相のモル濃度に対して温度・圧力に関する変分を取ると、次の(39)式が得られる。
【0136】
【数29】
【0137】
(δCLi/δTC)(n),(δCLj/δP)(n)は物性値推算法・データベースなどから取得する。(DCTC/Dt)(n+1)および(DCP/Dt)(n+1)の評価法に関しては、後述する。Yを使用する場合は、本項は存在しない。
【0138】
(3−1−4.上記以外の項)
上記以外の項として、拡散項および質量輸送に伴う項がある。これらの項では、次の(40)式を用い、各変数を陽的に扱う。
【0139】
【数30】
【0140】
(3−2.密度の化学種・温度・圧力依存性を考慮した解法)
上記の説明により、化学種の時間変化および時間更新計算法が構築できたので、これを質量・運動量・エネルギ保存式に組込み、全体の計算手順を整理する。
【0141】
密度の化学種・温度・圧力依存性を考慮して、(2)式〜(4)式で示した質量保存式を変形する。
【0142】
(15)式で表した状態方程式から、密度の変分に関して、次の(41)式が成り立つ。
【0143】
【数31】
【0144】
そして、(41)式により、密度のラグランジュ微分は、次の(42)式のように書くことができる。
【0145】
【数32】
【0146】
ここでは、Tcと同様に、Dc/Dtの、DL/Dt,DG/Dtに置換えられている。そして、モル濃度を用いる場合は、(11)式から、次の(43)式のように記述することができる。
【0147】
【数33】
【0148】
また、(2)式で表される気泡の質量保存式と上記(42)式により、次の(44)式が得られる。
【0149】
【数34】
【0150】
そして、(44)式を質量保存式として、解法を構成する。
次に、各項を評価する時刻を、(45)式に基づいて定める。
【0151】
【数35】
【0152】
次に、運動量およびエネルギ保存式の各項を評価する時刻を、(46)式および(47)式のように定める。
【0153】
【数36】
【0154】
そして、連続液相と連続気相の基礎式に関しても同様の操作を行なう。ただし、連続液相ではモル濃度を用いているため、(45)式に対応する式は、次の(48)式となる。
【0155】
【数37】
【0156】
Tomiyama et al. (1991)は、温度が速い過渡変化を伴う場合、温度および密度を評価する時刻が数値解に有意な差を生じさせることを指摘している。(45)式〜(48)式においても密度変化の影響が各式に可能な限り反映されるように、各項を評価する時刻を選択している。
【0157】
圧力のラグランジュ微分項である(DkP/Dt)(n+1)は、次の(49)式で評価する。
【0158】
【数38】
【0159】
Shimada and Tomiyama (2005)では、大きな静水圧差下の気泡流計算を通して、(49)式の右辺第1項を無視する近似の妥当性を確認している。(28)式と(47)式,(48)式により、時刻nの値のみを用いて(DkYki/Dt)(n+1),(DkP/Dt)(n+1)が求められる。このため、質量保存式(気泡の場合(44)式)の右辺は既知量となり、エネルギ保存式および化学種保存式を用いずに質量・運動量保存式のみを用いて新しい時刻の速度Vk(n+1)・圧力P(n+1)を計算できる。(1)式を考慮して(45)式に相当する式をすべての相に関して辺々を加え合せると、以下の混合体に対する質量保存式である(50)式を得る。
【0160】
【数39】
【0161】
次に、(21)式を運動量保存式(気泡の場合(46)式)に代入し、さらに、VBm (n+1)およびVG (n+1)について整理すると、次の(51)式および(52)式が得られる。
【0162】
【数40】
【0163】
ここで、Am,Bm,Cm,Dは、時刻nで評価される項をまとめたものである。そして、(39)式と(40)式を行列表示すると、次の(53)式が得られる。
【0164】
【数41】
【0165】
そして、(53)式をVBm(n+1),VC(n+1)について解き、整理すると、次の(54)式が得られる。
【0166】
【数42】
【0167】
ここで、ρの上付添字(〜)およびGの上付添字(〜)は、(41)式の右辺に逆行列を乗じた項を置換えたものである。(50)式および(54)式からVk(n+1),P(n+1)を求める手順には、種々の非圧縮性流れの数値解法を利用できる。ここでは、SMAC法(Simplified Marker and Cell method)(Amsden and. Hallow, 1970)に基づく求解手順を整理しておく。
【0168】
まず、圧力修正量δPを次の(55)式で定義する。
【0169】
【数43】
【0170】
次に、(54)式における圧力を既知の圧力P(n)で置換えた次の(56)式で、速度の推定値Vkの上付添字(〜)を定義する。
【0171】
【数44】
【0172】
(50)式に、(54)式〜(56)式を代入すると、以下のδPに関するポアソン方程式である(57)式が得られる。
【0173】
【数45】
【0174】
化学種と同様に温度を時間更新する際も、計算安定性確保のため界面エネルギ輸送項を陰的に評価する必要がある。温度計算の詳細については、文献(Shimada and Tomiyama, 2005)などで開示されている技術を採用することができる。
【0175】
<各モデルのシミュレーション方法>
図9は、この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。なお、図9に示す処理は、解析対策を空間的に分割して得られるそれぞれのセルについて実行される。
【0176】
図9を参照して、まず、プロセス解析装置は、ステップS100で、各相の密度ρ,体積率α,速度V,圧力P,モル濃度C,質量分率Yを初期化する。すなわち、シミュレーション対象となるプロセスの初期状態における値を各変数にセットする。その後、処理は、ステップS102へ進められる。
【0177】
ステップS102では、(10)式および(11)式に従って、各化学種についての(δρ/δP)nを算出する。その後、処理はステップS104に進められる。
【0178】
ステップS104では、(28)式、(47)式および(48)式に従って、(DkYki/Dt)(n+1),(DkP/Dt)(n+1)を計算する。その後、処理はステップS106へ進められる。
【0179】
ステップS106では、(56)式に従って、Vkの上付添字(〜)を計算する。その後、処理はステップS108へ進められる。
【0180】
ステップS108では、(57)式に従って、圧力修正量δPを計算する。その後、処理はステップS110へ進められる。
【0181】
ステップS110では、(54)式および(55)式に従って、速度Vk(n+1)および圧力P(n+1)を計算する。その後、処理はステップS112へ進められる。
【0182】
ステップS112〜ステップS116では、気泡の体積率αBmn+1と、気相の体積率αGn+1と、液相の体積率αLn+1が、それぞれ計算される。なお、これらの体積率は、(1)式、(45)式、および、(48)式に従って求められる。そして、処理はステップS118へ進められる。
【0183】
ステップS118では、図6および図7を参照して説明した第1液相領域ELと第2液相領域のそれぞれの体積率αLjn+1が、以下の(58)式および(59)式に従って計算される。そして、処理はステップS120へ進められる。
【0184】
【数46】
【0185】
ステップS120では、(17)式〜(20)式に従って第1液相領域と第2液相領域それぞれのモル濃度CLjin+1が計算される。そして、ステップS122へ処理が進められる。
【0186】
ステップS122では、各化学種の質量分率Ykn+1が計算される。そして、ステップS124へ処理が進められる。
【0187】
ステップS124では、(10)式および(11)式に従って、密度ρkn+1が計算される。その後、処理はステップS126へ進められる。
【0188】
ステップS126では、すべての計算処理が終了したか否かが判断される。すべての計算処理が終了していなければ(ステップS126においてNO)、時刻nに1が加算され(ステップS128)、次の時刻における各パラメータを算出するために処理はステップS102に戻される。そして、ステップS102以降の処理が再度実行される。一方、すべての計算処理が終了していれば(ステップS128においてYES)、処理は終了する。
【0189】
このように、解析対象を空間的に分離して得られるそれぞれのセルについて、状態値が時系列に計算される。
【0190】
<検証実験>
本願発明者らは、模擬装置を用いて、本実施の形態におけるパラメータβ,BL12を取得するためのシミュレーションについての検証を行なった。
【0191】
検証は、水のSO2ガス吸収操作を対象としている。
図10は、気泡流入開始後1秒におけるpHの分布を示す図である。なお、図10(A)は実験結果であり、図10(B)は、本実施の形態に従うプロセス解析方法による計算結果である。そして、図10(A)および図10(B)では、pHの値がハッチングの違いで表現されている。具体的には、図中、pH4.4からpH5.2までのpHの変化が、P1〜P6の6段階で示されている。これらは、P1がpH4.4に対応し、P6がpH5.2に対応する。また、図11は、気泡流入開始後1秒における速度ベクトルの分布を示す図である。なお、図11(A)は実験結果であり、図11(B)は、本実施の形態に従うプロセス解析方法による計算結果である。
【0192】
図10および図11を参照して、気泡流入開始後1秒では、pHと速度ベクトルのいずれについても、反射の影響により、気泡周辺でのパラメータの計算結果が実験値に沿ったものとはなっていない。
【0193】
図12および図13に、気泡流入開始後2秒における実験結果および計算結果を示す。なお、図12は、pHの分布を示す図である。図12(A)は実験結果であり、図12(B)は、本実施の形態に従うプロセス解析方法による計算結果である。図12(A)および図12(B)においても、pHの変化は、図10(A)および図10(B)と同様に、6段階のハッチングで示されている。また、図13は、速度ベクトルの分布を示す図である。なお、図13(A)は実験結果であり、図13(B)は、本実施の形態に従うプロセス解析方法による計算結果である。
【0194】
図12および図13を参照して、気泡流入開始後2秒では、実験結果において、気泡周辺に速度の高い領域が見られ、この傾向は、計算結果でも見られている。このことから、計算結果は、実験結果に沿ったものと言える。
【0195】
<実施の形態の効果>
たとえば、混相流における「A→R→S」の逐次反応を対象としてRの選択率を考えた場合、反応容器に供給されたAは、当該反応容器の流出部に進むにつれてRへと変化する。そして、Rの一部は反応が進んでSへと変化する。
【0196】
このような反応系において、反応容器内で混合があると、供給したAが反応せず(Rへと変化することなく)に流出部へと流れやすくなる。これと同時に、Aが反応することによって生じたRのうち、反応容器の入口に逆流して、Sへと変化する割合が高くなる。
【0197】
このように、反応容器(反応塔塔)内の液体の流れによって、反応成績は異なると考えられる。
【0198】
以上説明した本実施の形態では、液相中に気泡が分散した混相流の解析において、液相に関し、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算された。
【0199】
これにより、液相が気泡に同伴することによる、液相において生じる流れの作用が考慮するためのパラメータβ,BL12を取得できる。したがって、反応をメゾスケールで捉えて高精度にシミュレーションできるようになる。
【符号の説明】
【0200】
100 コンピュータ、101 コンピュータ本体、102 モニタ、103 キーボード、104 マウス、106 メモリ、107 固定ディスク、109 通信インターフェース、111 FD駆動装置、113 CD−ROM駆動装置。
【技術分野】
【0001】
本発明は、気液混合相におけるプロセス解析プログラムに関する。
【背景技術】
【0002】
化学プラントの酸化、水素化、発酵、廃水処理などのさまざまなプロセスにおいて、気泡塔と呼ばれる装置が用いられている。この気泡塔は、液相中に気泡が分散した「気泡流」を利用した設備であり、他の気液反応装置に比較してより大きな気液界面積が利用でき、また構造を単純化できるため装置の管理が容易になるといった利点を有している。
【0003】
このような気泡塔の解析では、気泡流構造が複雑であるので、実験的手法を用いて大規模な解析を行なうことは困難である。そのため、数値計算(シミュレーション)を用いて、気泡塔の解析を行なうことが注目されている。このようなシミュレーションにより気泡塔を解析する際には、対象の気泡塔についての計算モデルを高度化することが重要となっている。
【0004】
このような状況下において、たとえば非特許文献1〜3や特許文献1には、二流体モデルを用いた種々の計算モデルが提案されている。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特開2010−40013号公報
【非特許文献】
【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】島田他、「気泡塔内気泡流数値予測に必要となる相関式の感度解析」、社団法人 化学工学会、化学工学論文集、第29巻、第6号、pp.778−785、2003年
【非特許文献3】島田他、「化学反応・ガス吸収・熱輸送を伴う気泡塔内気泡流の数値解法」、社団法人 化学工学会、化学工学論文集、第31巻、第6号、pp.377−387、2005年
【非特許文献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】
具体的には、たとえば上記したような気泡塔が、直径0.4m以上、高さ4m以上のものとなった場合に、シミュレーション結果の実際の反応系からの乖離が顕著に現われていた。
【0009】
本発明は係る実情に鑑み考え出されたものであり、その目的は、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようにすることである。
【課題を解決するための手段】
【0010】
本発明に従ったプロセス解析プログラムは、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算するステップは、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算するステップと、質量保存式および状態保存式に基づいて、液相の体積率を計算するステップと、質量保存式および状態保存式に基づいて、気相の体積率を計算するステップと、化学種保存式に基づいてモル濃度を計算するステップとを含み、液相の体積率を計算するステップは、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算するステップは、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。
【0011】
好ましくは、第1液相領域の体積率は、気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる。
【0012】
好ましくは、パラメータは、0.01〜1.0の範囲の値に設定される。
本発明に従ったプロセス解析装置は、液相と気相を収容する反応容器内における、液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、状態値を時系列で計算する手段は、質量保存式および状態保存式に基づいて、液相中の気泡の体積率を計算する手段と、質量保存式および状態保存式に基づいて、液相の体積率を計算する手段と、質量保存式および状態保存式に基づいて、気相の体積率を計算する手段と、化学種保存式に基づいてモル濃度を計算する手段とを含み、液相の体積率を計算する手段は、液相内で気泡の影響を受ける第1液相領域の体積率と、液相内の第1液相領域以外の領域である第2液相領域の体積率とを計算し、モル濃度を計算する手段は、第1液相領域のモル濃度と第2液相領域のモル濃度とを計算する。
【発明の効果】
【0013】
本発明によれば、液相中に気泡が分散した混相流を解析する際に、液相において、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算される。
【0014】
これにより、気泡によって液体が流れを変えられる作用を考慮でき、気液混合相のシミュレーションにおいて、系のサイズによらず、液相中に気泡が分散した混相流での現象をより高精度にシミュレーションできるようになる。
【図面の簡単な説明】
【0015】
【図1】本発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【図2】反応における各成分の濃度の時間変化の一例を示す図である。
【図3】反応における各成分の濃度の時間変化の他の例を示す図である。
【図4】反応における各成分の濃度の時間変化のさらに他の例を示す図である。
【図5】反応の中間生成物の、生成量(溶液中の濃度)の時間変化を示す図である。
【図6】気泡が液体を同伴して流れる状態の概念図である。
【図7】液相の第1液相領域と第2液相領域を説明するための図である。
【図8】本発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図9】本発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。
【図10】気泡流入開始後1秒におけるpHの分布を示す図である。
【図11】気泡流入開始後1秒における速度ベクトルの分布を示す図である。
【図12】pHの分布を示す図である。
【図13】速度ベクトルの分布を示す図である。
【発明を実施するための形態】
【0016】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【0017】
<解析対象>
この発明の実施の形態に従うプロセス解析は、気泡塔などの液相中に気泡が分散した混相流を生じる系を解析対象とする。
【0018】
図1は、この発明の実施の形態に従うプロセス解析が解析対象とする系の一例を示す模式図である。
【0019】
図1を参照して、解析対象の一例として、反応器10内に原料1(液相)を供給するとともに、原料1中に原料2(気相)が吹き込まれる系SYSを考える。反応器10の底部には、原料1を供給するための供給ライン12が接続されている。この供給ライン12の経路上には、原料1を供給/遮断するための開閉弁14が設けられている。なお、図1に示す系SYSでは、バッチ処理またはセミバッチ処理が行なわれるので、原料1は、当該バッチ処理毎に供給される。
【0020】
また、反応器10には、気相の原料2を供給するための供給ライン16が接続されており、この供給ライン16は、反応器10内で水平方向に延びる吹き出しライン18と接続される。吹き出しライン18には、原料1を液相L内に吹き込むための孔が垂直上方向に向けて多数設けられている。これにより、気相の原料2が供給されると、原料1の液相中に原料2の気泡が分散し、この気泡が上昇することにより流れが生じる。また、気泡中の反応成分が液相中に溶解することにより、酸化、水素化、発酵、廃水処理などの化学反応が生じる。
【0021】
反応器10の底部には、化学反応によって生成された製品を取り出すための取り出しライン20が接続されており、反応器10内で生じた生成品は、この取り出しライン20を通じて次工程へ運ばれる。また、この反応器10の頂点部には、反応器10内の圧力および液面を調整するため、開閉弁26を設けた通気ライン24が接続されている。
【0022】
<本発明の概要>
図1に示したような気泡塔内では、たとえば以下のa)〜e)に示すような現象が起こっていると考えられる。
【0023】
a) 底部から液体が供給される
b) 塔内部に空気が供給される
c) 液体(液相)中を、空気が気泡として分散しながら流れる
d) 気泡から液体(液相)中にガスが溶解する
e) 溶解したガスと液体が反応する
本明細書では、これだけの複雑な過程を含む反応系について、まず、ラボの規模で反応系を検討し、そして、ミニプラント、実機へと順にスケールアップさせることとした。つまり、ラボやミニプラントで得た知見を実機の合理化設計へと適切に展開するために、上記現象を様々なスケールで捉え、スケールごとに重要となる因子を議論することにより、スケールアップのあり方を調整することとした。
【0024】
そこで、本明細書では、次の表1に挙げられた、ミクロ、メゾ、マクロの各スケールにおいて検討を行なう。なお、表1では、生じると考えられる主な現象によって、スケールが定義されている。ミクロスケールは、反応機構レベルで現象を捉えるスケールである。マクロスケールは、反応塔内における液混合全体的に現象を捉えるスケールである。メゾは、ミクロとマクロの中間のスケールであり、ガスの溶解や、その際に各気泡が液を同伴して移動するなど、相間の界面での現象を捉えるスケールである。
【0025】
【表1】
【0026】
(ミクロスケール)
同じ反応機構を持っていた場合であっても、気泡塔内の液体の流れによって反応成績が変化する。その原因について、以下説明する。
【0027】
反応容器における流体の流れとして、供給した液体が瞬時にして容器全体に混合する流れを「完全混合」と、供給した液体が流出部へと一次元的に移動する「押し出し流れ」とがある。
【0028】
ラボで使用される攪拌槽反応器における流れは、完全混合とみなすことができる。
一方、押し出し流れは、反応容器内の液体が全く混合していない状態と捉えることができる。たとえば、管内の流れなどは、押し出し流れの一例と言える。
【0029】
そして、通常の反応容器内の流れは、完全混合と押し出し流れの間に位置づけられる。
ここで、次の(EQ)式に示すような、出発物質Aが中間体Rを経て、最終生成物Sへと変化する「A→R→S」という逐次反応を対象として、Rの選択率を検討する。なお、k1,k2は、それぞれの反応の反応速度定数である。
【0030】
【数1】
【0031】
供給されたAは、反応容器の流出部に進むにつれてRへと変化する。Rの一部は、反応が進んでSへと変化する。
【0032】
反応容器において混合があると、供給したAが反応せずに流出部へと流れやすくなる一方で、Aが反応して生成したRが入口に逆流して、Sへと変化する割合が増える。
【0033】
反応液に限らず、流体の流れの状態を表す指標として混合段数(混合槽数)という値が知られている。この混合段数とは、筒状の胴体を有する容器(以下、適宜「流通器」と称する)の中を一方向に流れる流体が、流通器を流れる間に混合される度合いを示す値であり、上記流通器を、直列に連結された互いに同一体積の複数の完全混合槽にモデル化した場合の、当該完全混合槽の数である。混合段数については、文献A(橋本健治、「反応工学」、改訂版、培風館、1993年3月20日、p.179−187、195―197)および文献B(Octave Levenspiel, “Chemical reaction engineering”, Fifth printing, USA, John Wiley and Sons, INC., May, 1967 p.261-265)にその説明が記載されている。そして、この混合段数が1(最小)の場合が、押し出し流れの状態であり、混合段数が無限大の場合が、完全混合の状態である。
【0034】
(EQ)式の反応における、A成分、O成分、R成分、および、S成分の各濃度についての、従来の技術に従ったシミュレーション結果について、説明する。
【0035】
図2は、k1とk2の値が同じである場合の、各成分の濃度の時間変化を示す図である。なお、図2(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図2(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。また、各図では、A成分の濃度が実線で、O成分の濃度が破線で、R成分の濃度が一点鎖線で、そして、S成分の濃度が点線で、記載されている。
【0036】
図2(B)の不均一の場合には、反応することなく残留するA成分およびR成分がほとんど見られないのに対し、図2(A)の均一の場合には、或る程度の量のA成分およびR成分が反応することなく残留している。
【0037】
図3は、k1がk2よりも大きい場合の、各成分の濃度の時間変化を示す図である。なお、図3(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図3(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。図3では、各成分の濃度変化が、図2と同様の線の種類で記載されている。
【0038】
また、図4は、k1がk2よりも小さい場合の、各成分の濃度の時間変化を示す図である。なお、図4(A)は、反応容器内が均一の場合(つまり、完全混合の場合)を示し、図4(B)は、反応容器内の不均一の場合(つまり、押し出し流れの場合)を示す。図4では、各成分の濃度変化が、図2と同様の線の種類で記載されている。
【0039】
図3と図4に示されたいずれの場合においても、反応容器内が均一である場合と不均一である場合とでは、各成分の濃度変化の挙動に差異を生じている。
【0040】
つまり、図2〜図4から理解されるように、反応系における各成分の濃度の時間変化は、k1とk2の間の大小関係に関わらず、反応容器内が均一であるか不均一であるかによって相違することとなる。
【0041】
図5は、(EQ)式の中間生成物であるR成分の、生成量(溶液中の濃度)の時間変化を示す図である。なお、図5(A)は、図2(A),図2(B),図3(A),図3(B),図4(A),図4(B)のそれぞれに示されたR成分についての濃度変化を示す。そして、図5(B)は、図5(A)に示されたそれぞれの状況下でのR成分の濃度変化の、反応初期における変化を拡大して(つまり、図5(A)の各グラフを横軸方向に拡大して)、示すものである。
【0042】
図5(A)および図5(B)に示されるように、中間生成物の生成量は、反応条件によって、大きく異なるといえる。ここでいう反応条件とは、k1とk2の間の大小関係と、反応容器における流体の流れ(完全混合/押し出し流れ)との組合せである。また、図5(A)および図5(B)において、R成分の生成量の挙動の差異は、反応初期だけでなく、反応中期においても見られている。
【0043】
なお、k1とk2の大小関係のいずれであっても(k1=k2、k1>k2、または、k1<k2)、R成分の濃度は、反応初期には、完全混合よりも押し出し流れの方が高くなっているが、反応中期以降は、完全混合の方が高くなる傾向にある。つまり、同じ反応率において完全混合よりも押し出し流れの方が、Sを生成する反応についての選択率が高いと言える。
【0044】
以上のことから、反応系における生成物の生成量を予測するためには、k1とk2の間の大小関係とともに、反応容器における流体の流れについても検討することは重要であるといえる。
【0045】
(メゾスケール)
上記したように、反応における選択率は反応容器内の流体の流れに影響を受け、反応機構が同じであれば、完全混合よりも選択率の低い流れは存在しないものと考えられる。しかしながら、実機における反応成績には、完全混合よりも選択率が劣る場合が観測されている。
【0046】
ここで、上記したような考察は、反応容器内では液体が単相で流れているという想定に基づいたものであった。これに対して、実際の気泡塔内では、気泡と液体が混在した混相流の状態になっている。このことから、上記したような、完全混合よりも選択率の低い系が存在したという結果は、気泡が液体を同伴して流れたことによる特別な作用に起因するものではないかと考えられる。
【0047】
図6は、気泡が液体を同伴して流れる状態の概念図である。図6では、液相L内を、気泡Bが、液体ELを同伴して、上方へ流れている状態が示されている。
【0048】
従来の計算モデルでは、計算セルサイズΔxに比べて、小さな気泡を分散相とし、また、その大きさに応じてN群に分類することにより、各相の体積率αを、(1)式で表していた。
【0049】
【数2】
【0050】
(1)式では、気泡流に含まれる相を1つの連続気相(添え字G)、1つの連続液相(添え字L)、そして、N群の気泡(Bm:Bm=1,2,・・・,N)の合計(N+2)種のグループに分類がなされていた。
【0051】
本実施の形態では、液相において、気泡の影響を受ける可能性を考慮して、体積率の検討において領域の分類を図る。具体的には、図7に示されるように、液相において、気泡の影響を受ける液体ELを第1液相領域とし、気泡の影響を受けにくいと考えられるバルクの液相Lを第2液相領域とする。本実施の形態では、(1)式においてαLで記述されていた液相の体積率は、第1液相領域の体積率αL1と第2液相領域の体積率αL2とを用いて、次の式のように記述される。
【0052】
【数3】
【0053】
そして、第1液相領域の体積率αL1は、それを同伴する気泡の体積率αBに依存すると考えられる。したがって、各液体ELの体積率αL1は、各気泡の体積率αBを用いて、次の式のように記述される。
【0054】
【数4】
【0055】
ここで、βは、気泡の影響範囲を表すパラメータである。
そして、液相の体積率は、次の式で示すように、第1液相領域の体積率αL1の総量と第2液相領域の体積率αL2の総量の和として記述される。
【0056】
【数5】
【0057】
上記したように、本実施の形態では、計算モデルにおける体積率について、気泡の影響の有無を考慮し、液相の領域を第1液相領域と第2液相領域に分割して検討している。そして、第1液相領域の体積率αL1を、気泡の体積率αBとパラメータβによって、つまり、系のスケールに依存しない値によって、言及している。
【0058】
これにより、本実施の形態では、系のスケールに依存することなく、液相における気泡の影響を考慮できるため、より実際の反応における挙動に沿ったシミュレーションが可能となる。
【0059】
なお、本明細書では、気泡Bと区別するために、反応系内のバルクの気相を適宜気相Gと表記する。
【0060】
<ハードウェア構成>
本実施の形態に従うプロセス解析装置は、代表的にコンピュータによって実現される。
【0061】
図8は、この発明の実施の形態に従うプロセス解析装置を実現するための代表的なハードウェア構成であるコンピュータ100の概略構成図である。
【0062】
図8を参照して、コンピュータ100は、FD(Flexible Disk)駆動装置111およびCD−ROM(Compact Disk-Read Only Memory)駆動装置113を搭載したコンピュータ本体101と、モニタ102と、キーボード103と、マウス104とを含む。
【0063】
コンピュータ本体101は、FD駆動装置111およびCD−ROM駆動装置113に加えて、相互にバスで接続された、演算装置であるCPU(Central Processing Unit)105と、メモリ106と、記憶装置である固定ディスク107と、通信インターフェース109とを含む。
【0064】
本実施の形態に従うプロセス解析処理は、CPU105がメモリ106などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、FD112やCD−ROM114などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、FD駆動装置111やCD−ROM駆動装置113などにより記録媒体から読取られて、または通信インターフェース109にて受信されて、固定ディスク107に格納される。さらに、このようなプログラムは、固定ディスク107からメモリ106に読出されて、CPU105により実行される。
【0065】
CPU105は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ106は、CPU105のプログラム実行に応じて各種の情報を記憶する。
【0066】
モニタ102は、CPU105が出力する情報を表示するための表示部であって、一例としてLCD(Liquid Crystal Display)やCRT(Cathode Ray Tube)などから構成される。すなわち、モニタ102には、プロセス制御の状態などが表示される。
【0067】
マウス104は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード103は、入力されるキーに応じたユーザから指令を受付ける。
【0068】
通信インターフェース109は、コンピュータ100と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。
【0069】
なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。
【0070】
また、コンピュータ100を本実施の形態に従うように機能させるプログラムが記録される記憶媒体は、上記したFDやCD−ROMに限定されない。記憶媒体の他の例としては、DVD−ROM(Digital Versatile Disk - Read Only Memory)、USB(Universal Serial Bus)メモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、MO(Magnetic Optical Disc)、MD(Mini Disc)、IC(Integrated Circuit)カード(メモリカードを除く)、光カード、マスクROM、EPROM、EEPROM(Electronically Erasable Programmable Read-Only Memory)などの、不揮発的にプログラムを格納する媒体が挙げられる。
【0071】
<計算モデル>
以下、本実施の形態に従うプロセス解析に用いられる計算モデルについて説明する。
【0072】
(1.基礎式)
本実施の形態に従うプロセス解析に用いる計算モデルでは、気泡流に含まれる相を1つの連続気相(添え字G)、液相領域(添え字L)、N群の気泡(Bm,Bm=1,2,・・・,N)の合計(N+2)種のグループに分類する。
【0073】
本実施の形態に従うプロセス解析では、各相((N+2)種の各グループ)について、質量保存式、運動量保存式、状態方程式、および、化学種保存式を考えることで、状態値の時間的変化を計算する。
【0074】
以下の計算モデルでは、計算セルサイズΔxに比べて、気液界面の長さスケールが大きな場合には、その界面で隔てられた二相を連続気相と連続液相(上記した第2液相領域)とみなす。一方、計算セルサイズΔxに比べて、小さな気泡は分散相とし、さらにその大きさに応じてN群に分類する。各相の体積率αは、体積率の定義から上記(1)式を満たす。
【0075】
気液両相は相変化のない等温非圧縮性ニュートン流体と仮定し、気泡内部の分子および乱流拡散を無視すると、セル体積平均に基づいて、(2)式〜(6)式に示す質量保存式および運動量保存式が成立する。
【0076】
【数6】
【0077】
ここで、下付添字のGは連続気相を、Lは連続液相を、Bmはグループm(m=1,2,…,N)の気泡を意味する。また、tは時間、Vは速度、ρは密度、ΓCBは合体・分裂による単位体積・単位時間当りの質量生成量、ΓAEは吸収・蒸発による単位体積・単位時間当りの質量生成量を表す。なお、質量生成量Γの下付添字は輸送に関与する二相を意味し、たとえばLBmは気泡グループ Bmから連続液相Lへの輸送を表す。また、DBm/DtおよびDC/Dtは、以下の実質微分を意味する。
【0078】
【数7】
【0079】
本実施の形態に従うプロセス解析で用いる計算モデルについても、2つの連続相に対して一流体近似を採用している。そのため、2つの連続相に対する運動量保存式としては、連続相混合体についての(6)式のみを用いる。(6)式における連続相混合体の体積率αCは、αC=αG+αLとして与え、混合体の密度ρCを次の(7)式として与える。
【0080】
【数8】
【0081】
また、本実施の形態に従うプロセス解析で用いる計算モデルでは、化学種の量を評価するパラメータとして、モル濃度・モル分率・分密度・分圧・質量分率・絶対湿度があげられる。ここでは、質量分率Yとモル濃度Cを用いて、状態方程式を構成する場合を考える。グループmの気泡のみが存在する場合、Yの定義から次の(8)式が成立する。
【0082】
【数9】
【0083】
ここで、下付添字のiは、化学種を表す。また、ρBmi(P,TBm)は、純物質時の状態式である(8)式から、平均密度に関して次の(9)式が得られる。
【0084】
【数10】
【0085】
連続液相と連続気相の平均密度も同様の手順で得られる。本実施の形態では、各相(Bm,L,G)を一括して表す下付添字kを用い、以下の(10)式のように表記する。
【0086】
【数11】
【0087】
(10)式では、添字kを用いて表記する便宜上、TcをTL,TG(一流体近似からTc=TL=TG)に置き換えている。一方、モル濃度Ckiと平均密度との関係は、次の(11)式を採用している。
【0088】
【数12】
【0089】
ここで、Mは分子量をあらわす。変数YとCは、(10)式を用いる方法と(11)式を用いる方法は等価である。これらの二式を比較した場合、変分が単純となる(11)式を用いる方が簡便である。また、多くの場合、反応速度はCに直接依存するため、多くの等温系液相流れの計算でCが使用されている。しかしながら、圧縮性流体を扱う場合はCに温度・圧力依存性を加味する必要がある。一方、(10)式を用いる場合は、ρki(P,Tk)を物性値推算法・データベースなどを用いて設定することで、密度の温度・圧力依存性を加味できる。このため、圧縮性・非圧縮性の区別なくYの輸送計算ができる利点があり、ρki(P,Tk)が明確なガス燃焼流れなどの計算に良く使用されている。工業的なガス吸収操作では、ガスを溶質として液相側で反応させる場合が極めて多い。本実施の形態では、気相(Bm,G)の化学種を質量分率Y、液相(L)の化学種をモル濃度Cで評価し、化学種保存式として、次の(12)式〜(14)式を用いる。
【0090】
【数13】
【0091】
なお、(13)式において、Kは気泡(BM)または液相に溶解している気体(C)を表す。これにより、(13)式では、気体成分について、気泡(BM)の状態および液相内の溶質の状態の運動が考慮されている。
【0092】
ここで、wは反応に伴う単位体積当たりのモル生成率、Dは拡散係数、BLBmiは気泡グループと連続液相における単位体積・単位時間当たりの化学種iの生成量、ΦBmiとΦGiはΓGBmに伴う化学種iの生成量を表す。(12)式〜(14)式の右辺第一項は生成項、(12)式の第二項と(14)式の第三項はΓに伴う生成項、(13)式と(14)式の第二項は拡散項、(13)式の第三項は温度・圧力依存項である。
【0093】
なお、本実施の形態では、図7を参照して上記したように、液相を第1液相領域と第2液相領域に分割して検討を行なうため、(13)式が改良される。改良にあたり、検討を容易にするため、次の仮定を設ける。
【0094】
仮定1:連続気相から液相への物質移動はない
仮定2:各モル濃度は圧力の影響を受けない
仮定3:気相からの物質移動は第1液相領域への移動のみが生じ、その後、第1液相領域と第2液相領域の間で交換が生じる
仮定1によれば次の(15)式が導かれ、仮定2によれば次の(16)式が導かれる。
【0095】
【数14】
【0096】
また、モル濃度に関する(13)式を第1液相領域と第2液相領域に分割して検討すると、次の(17)式および(18)式が導かれる。
【0097】
【数15】
【0098】
ここで、Ljiは液相の領域jにいる成分iを意味している。また、BL12iは第1液相領域と第2液相領域の間での物質交換量であり、BL1Bmiは気相からの物質移動についての物質交換量である。気相からの物質移動が無い成分については、の場合、BL1Bmi=0となる。したがって、この場合、(17)式と(18)式は次の(19)式と(20)式のように書ける。
【0099】
【数16】
【0100】
(17)式〜(20)式を解くには、本実施の形態において新たに採用したパラメータβとBL12を指定する必要がある。
【0101】
βは、気泡の形状、上昇速度、密集度、乱れの状況など、多くの流れの作用を受ける複雑なパラメータであるが、既往の研究(例えば、R. Clift , J. R. Grace , M. E. Weber“Bubbles,Drops and Particles”, Dover Publications, 2005)からb = 0.01〜1.0の範囲が考えられる。
【0102】
BL12についても、気泡周りに関する詳細な知見から与える必要がある。これらの情報は、以下の研究内容から取得できる。
【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】
(2.構成式)
界面運動量輸送MLBmFには、次の(21)式で表される構成式が用いられる。
【0106】
【数17】
【0107】
また、粘性および乱流拡散項Fμには、次の(22)式で表される構成式が用いられる。
【0108】
【数18】
【0109】
ここで、上付添字Tは転置であり、μcは連続相混合体の粘度である。
また、表面張力FSには、次の(23)式が用いられる。
【0110】
【数19】
【0111】
また、化学種輸送量Bおよびこれに伴う質量輸送ΓAEには、次の(24)式〜(27)式が用いられる。
【0112】
【数20】
【0113】
(3.解法)
上述した式を用いて、対象プロセスにおける状態値の時間的変化を計算するための方法について説明する。
【0114】
(3−1.化学種保存式の解法)
まず、化学種を求めるための解法を提示する。
【0115】
(12)式〜(14)式および(15)式〜(20)式を参照して説明した化学種保存式には、対流項、反応項、相間輸送項、拡散項、温度・圧力依存項、体積輸送に伴う項が含まれている。各項の特性時間が大きく異なる場合、同じ時間積分法で各項を計算する方法は非効率である。そこで、各項を相間輸送項(上付数字INT)、反応項(RT)、温度・圧力依存項(PT)、それ以降の項(OT)に分類し別々に取り扱うことで、計算の効率向上を図る。
【0116】
YまたはCを一括してΨとすると、その時間変化は、次の(28)式で評価できる。
【0117】
【数21】
【0118】
ここで、上付添字(n+1)は時刻(n+1)Δtを表す。
(28)式は、後述する圧力・速度計算時に必要となる。圧力・速度計算後、Ψの時間更新を行なう。時間更新時には、計算不安定化を回避するため、陰解法を用いて相間輸送項による計算を行なった後、次の(29)式を用いて相間輸送項以外の項で時間更新する。
【0119】
【数22】
【0120】
ここで、上付添字(*)は、相間輸送項のみを用いて更新する値を意味する。
以下に、相間輸送項、反応項、温度・圧力依存項、それ以外の項の評価法を整理する。
【0121】
(3−1−1.相間輸送項)
気泡塔内気泡流では、気泡から液相へのガス吸収が最も支配的な相間輸送項である。そこで、簡便のため、以下では相間の化学種iの輸送を気泡から液相へのガス吸収に限定する。(12)式〜(14)式等を参照して説明したような相間輸送項に後述のモル流束に関する相関式を代入すると、以下の(30)式および(31)式のように整理できる。
【0122】
【数23】
【0123】
ここで、ξ,Ψ,ζ,ηは代入操作によって生じるYBmi,CLiの各係数をまとめたものである。時間変化項を評価する際には、(30)式〜(31)式の右辺におけるYBmi,CLiを時刻nで評価した、次の(32)式および(33)式を用いる。
【0124】
【数24】
【0125】
一方、時間更新時には本項による計算不安定化を回避するため、(30)式〜(31)式の右辺を陰的に扱った、次の(34)式および(35)式を用いる。
【0126】
【数25】
【0127】
(34)式と(35)式を行列表示に置換え、さらに、YBmi (*),CLi (*)について解くと、次の(36)式が得られる。
【0128】
【数26】
【0129】
(3−1−2.反応項)
ここでは、液相内での反応のみを取扱う。
【0130】
化学反応速度式は非線形連立方程式で、特性時間が運動量輸送式離散化的に用いる時間刻み幅Δtより短い場合が多い。計算不安定化と計算誤差を回避するため、反応項の離散化には半陰解法を用いる。さらに、ΔtをJ分割したΔtJ(=Δt/J)を定義し、l=1,…,Jとして、次の(37)式による計算をCLi (l+1)に関して繰返す。
【0131】
【数27】
【0132】
ここで、上付添字の(n)は時刻nΔtを、(l)は時刻nΔt+lΔtJを表しており、l=0は時刻nΔtを、l=Jは時刻(n+1)Δtを意味する。
【0133】
(37)式の計算過程においてCLiが平衡状態に達したとみなされる場合は、J回以下でも繰返し計算を終了できる。反応項による時間変化は、次の(38)式で評価する。
【0134】
【数28】
【0135】
(3−1−3.温度・圧力依存項)
連続液相のモル濃度に対して温度・圧力に関する変分を取ると、次の(39)式が得られる。
【0136】
【数29】
【0137】
(δCLi/δTC)(n),(δCLj/δP)(n)は物性値推算法・データベースなどから取得する。(DCTC/Dt)(n+1)および(DCP/Dt)(n+1)の評価法に関しては、後述する。Yを使用する場合は、本項は存在しない。
【0138】
(3−1−4.上記以外の項)
上記以外の項として、拡散項および質量輸送に伴う項がある。これらの項では、次の(40)式を用い、各変数を陽的に扱う。
【0139】
【数30】
【0140】
(3−2.密度の化学種・温度・圧力依存性を考慮した解法)
上記の説明により、化学種の時間変化および時間更新計算法が構築できたので、これを質量・運動量・エネルギ保存式に組込み、全体の計算手順を整理する。
【0141】
密度の化学種・温度・圧力依存性を考慮して、(2)式〜(4)式で示した質量保存式を変形する。
【0142】
(15)式で表した状態方程式から、密度の変分に関して、次の(41)式が成り立つ。
【0143】
【数31】
【0144】
そして、(41)式により、密度のラグランジュ微分は、次の(42)式のように書くことができる。
【0145】
【数32】
【0146】
ここでは、Tcと同様に、Dc/Dtの、DL/Dt,DG/Dtに置換えられている。そして、モル濃度を用いる場合は、(11)式から、次の(43)式のように記述することができる。
【0147】
【数33】
【0148】
また、(2)式で表される気泡の質量保存式と上記(42)式により、次の(44)式が得られる。
【0149】
【数34】
【0150】
そして、(44)式を質量保存式として、解法を構成する。
次に、各項を評価する時刻を、(45)式に基づいて定める。
【0151】
【数35】
【0152】
次に、運動量およびエネルギ保存式の各項を評価する時刻を、(46)式および(47)式のように定める。
【0153】
【数36】
【0154】
そして、連続液相と連続気相の基礎式に関しても同様の操作を行なう。ただし、連続液相ではモル濃度を用いているため、(45)式に対応する式は、次の(48)式となる。
【0155】
【数37】
【0156】
Tomiyama et al. (1991)は、温度が速い過渡変化を伴う場合、温度および密度を評価する時刻が数値解に有意な差を生じさせることを指摘している。(45)式〜(48)式においても密度変化の影響が各式に可能な限り反映されるように、各項を評価する時刻を選択している。
【0157】
圧力のラグランジュ微分項である(DkP/Dt)(n+1)は、次の(49)式で評価する。
【0158】
【数38】
【0159】
Shimada and Tomiyama (2005)では、大きな静水圧差下の気泡流計算を通して、(49)式の右辺第1項を無視する近似の妥当性を確認している。(28)式と(47)式,(48)式により、時刻nの値のみを用いて(DkYki/Dt)(n+1),(DkP/Dt)(n+1)が求められる。このため、質量保存式(気泡の場合(44)式)の右辺は既知量となり、エネルギ保存式および化学種保存式を用いずに質量・運動量保存式のみを用いて新しい時刻の速度Vk(n+1)・圧力P(n+1)を計算できる。(1)式を考慮して(45)式に相当する式をすべての相に関して辺々を加え合せると、以下の混合体に対する質量保存式である(50)式を得る。
【0160】
【数39】
【0161】
次に、(21)式を運動量保存式(気泡の場合(46)式)に代入し、さらに、VBm (n+1)およびVG (n+1)について整理すると、次の(51)式および(52)式が得られる。
【0162】
【数40】
【0163】
ここで、Am,Bm,Cm,Dは、時刻nで評価される項をまとめたものである。そして、(39)式と(40)式を行列表示すると、次の(53)式が得られる。
【0164】
【数41】
【0165】
そして、(53)式をVBm(n+1),VC(n+1)について解き、整理すると、次の(54)式が得られる。
【0166】
【数42】
【0167】
ここで、ρの上付添字(〜)およびGの上付添字(〜)は、(41)式の右辺に逆行列を乗じた項を置換えたものである。(50)式および(54)式からVk(n+1),P(n+1)を求める手順には、種々の非圧縮性流れの数値解法を利用できる。ここでは、SMAC法(Simplified Marker and Cell method)(Amsden and. Hallow, 1970)に基づく求解手順を整理しておく。
【0168】
まず、圧力修正量δPを次の(55)式で定義する。
【0169】
【数43】
【0170】
次に、(54)式における圧力を既知の圧力P(n)で置換えた次の(56)式で、速度の推定値Vkの上付添字(〜)を定義する。
【0171】
【数44】
【0172】
(50)式に、(54)式〜(56)式を代入すると、以下のδPに関するポアソン方程式である(57)式が得られる。
【0173】
【数45】
【0174】
化学種と同様に温度を時間更新する際も、計算安定性確保のため界面エネルギ輸送項を陰的に評価する必要がある。温度計算の詳細については、文献(Shimada and Tomiyama, 2005)などで開示されている技術を採用することができる。
【0175】
<各モデルのシミュレーション方法>
図9は、この発明の実施の形態に従う計算モデルを用いたシミュレーションの手順を示すフローチャートである。なお、図9に示す処理は、解析対策を空間的に分割して得られるそれぞれのセルについて実行される。
【0176】
図9を参照して、まず、プロセス解析装置は、ステップS100で、各相の密度ρ,体積率α,速度V,圧力P,モル濃度C,質量分率Yを初期化する。すなわち、シミュレーション対象となるプロセスの初期状態における値を各変数にセットする。その後、処理は、ステップS102へ進められる。
【0177】
ステップS102では、(10)式および(11)式に従って、各化学種についての(δρ/δP)nを算出する。その後、処理はステップS104に進められる。
【0178】
ステップS104では、(28)式、(47)式および(48)式に従って、(DkYki/Dt)(n+1),(DkP/Dt)(n+1)を計算する。その後、処理はステップS106へ進められる。
【0179】
ステップS106では、(56)式に従って、Vkの上付添字(〜)を計算する。その後、処理はステップS108へ進められる。
【0180】
ステップS108では、(57)式に従って、圧力修正量δPを計算する。その後、処理はステップS110へ進められる。
【0181】
ステップS110では、(54)式および(55)式に従って、速度Vk(n+1)および圧力P(n+1)を計算する。その後、処理はステップS112へ進められる。
【0182】
ステップS112〜ステップS116では、気泡の体積率αBmn+1と、気相の体積率αGn+1と、液相の体積率αLn+1が、それぞれ計算される。なお、これらの体積率は、(1)式、(45)式、および、(48)式に従って求められる。そして、処理はステップS118へ進められる。
【0183】
ステップS118では、図6および図7を参照して説明した第1液相領域ELと第2液相領域のそれぞれの体積率αLjn+1が、以下の(58)式および(59)式に従って計算される。そして、処理はステップS120へ進められる。
【0184】
【数46】
【0185】
ステップS120では、(17)式〜(20)式に従って第1液相領域と第2液相領域それぞれのモル濃度CLjin+1が計算される。そして、ステップS122へ処理が進められる。
【0186】
ステップS122では、各化学種の質量分率Ykn+1が計算される。そして、ステップS124へ処理が進められる。
【0187】
ステップS124では、(10)式および(11)式に従って、密度ρkn+1が計算される。その後、処理はステップS126へ進められる。
【0188】
ステップS126では、すべての計算処理が終了したか否かが判断される。すべての計算処理が終了していなければ(ステップS126においてNO)、時刻nに1が加算され(ステップS128)、次の時刻における各パラメータを算出するために処理はステップS102に戻される。そして、ステップS102以降の処理が再度実行される。一方、すべての計算処理が終了していれば(ステップS128においてYES)、処理は終了する。
【0189】
このように、解析対象を空間的に分離して得られるそれぞれのセルについて、状態値が時系列に計算される。
【0190】
<検証実験>
本願発明者らは、模擬装置を用いて、本実施の形態におけるパラメータβ,BL12を取得するためのシミュレーションについての検証を行なった。
【0191】
検証は、水のSO2ガス吸収操作を対象としている。
図10は、気泡流入開始後1秒におけるpHの分布を示す図である。なお、図10(A)は実験結果であり、図10(B)は、本実施の形態に従うプロセス解析方法による計算結果である。そして、図10(A)および図10(B)では、pHの値がハッチングの違いで表現されている。具体的には、図中、pH4.4からpH5.2までのpHの変化が、P1〜P6の6段階で示されている。これらは、P1がpH4.4に対応し、P6がpH5.2に対応する。また、図11は、気泡流入開始後1秒における速度ベクトルの分布を示す図である。なお、図11(A)は実験結果であり、図11(B)は、本実施の形態に従うプロセス解析方法による計算結果である。
【0192】
図10および図11を参照して、気泡流入開始後1秒では、pHと速度ベクトルのいずれについても、反射の影響により、気泡周辺でのパラメータの計算結果が実験値に沿ったものとはなっていない。
【0193】
図12および図13に、気泡流入開始後2秒における実験結果および計算結果を示す。なお、図12は、pHの分布を示す図である。図12(A)は実験結果であり、図12(B)は、本実施の形態に従うプロセス解析方法による計算結果である。図12(A)および図12(B)においても、pHの変化は、図10(A)および図10(B)と同様に、6段階のハッチングで示されている。また、図13は、速度ベクトルの分布を示す図である。なお、図13(A)は実験結果であり、図13(B)は、本実施の形態に従うプロセス解析方法による計算結果である。
【0194】
図12および図13を参照して、気泡流入開始後2秒では、実験結果において、気泡周辺に速度の高い領域が見られ、この傾向は、計算結果でも見られている。このことから、計算結果は、実験結果に沿ったものと言える。
【0195】
<実施の形態の効果>
たとえば、混相流における「A→R→S」の逐次反応を対象としてRの選択率を考えた場合、反応容器に供給されたAは、当該反応容器の流出部に進むにつれてRへと変化する。そして、Rの一部は反応が進んでSへと変化する。
【0196】
このような反応系において、反応容器内で混合があると、供給したAが反応せず(Rへと変化することなく)に流出部へと流れやすくなる。これと同時に、Aが反応することによって生じたRのうち、反応容器の入口に逆流して、Sへと変化する割合が高くなる。
【0197】
このように、反応容器(反応塔塔)内の液体の流れによって、反応成績は異なると考えられる。
【0198】
以上説明した本実施の形態では、液相中に気泡が分散した混相流の解析において、液相に関し、気泡に同伴する領域とそれ以外の領域とに分類されて、体積率およびモル濃度が計算された。
【0199】
これにより、液相が気泡に同伴することによる、液相において生じる流れの作用が考慮するためのパラメータβ,BL12を取得できる。したがって、反応をメゾスケールで捉えて高精度にシミュレーションできるようになる。
【符号の説明】
【0200】
100 コンピュータ、101 コンピュータ本体、102 モニタ、103 キーボード、104 マウス、106 メモリ、107 固定ディスク、109 通信インターフェース、111 FD駆動装置、113 CD−ROM駆動装置。
【特許請求の範囲】
【請求項1】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算するステップは、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算するステップと、
前記化学種保存式に基づいてモル濃度を計算するステップとを含み、
前記液相の体積率を計算するステップは、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、
前記モル濃度を計算するステップは、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析プログラム。
【請求項2】
前記第1液相領域の体積率は、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる、請求項1に記載のプロセス解析プログラム。
【請求項3】
前記パラメータは、0.01〜1.0の範囲の値に設定される、請求項2に記載のプロセス解析プログラム。
【請求項4】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、
空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算する手段は、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算する手段と、
前記化学種保存式に基づいてモル濃度を計算する手段とを含み、
前記液相の体積率を計算する手段は、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、
前記モル濃度を計算する手段は、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析装置。
【請求項1】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析プログラムであって、空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算するステップをコンピュータに実行させ、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算するステップは、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算するステップと、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算するステップと、
前記化学種保存式に基づいてモル濃度を計算するステップとを含み、
前記液相の体積率を計算するステップは、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、
前記モル濃度を計算するステップは、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析プログラム。
【請求項2】
前記第1液相領域の体積率は、前記気泡の体積率と当該気泡の影響範囲を表すパラメータとの積によって求められる、請求項1に記載のプロセス解析プログラム。
【請求項3】
前記パラメータは、0.01〜1.0の範囲の値に設定される、請求項2に記載のプロセス解析プログラム。
【請求項4】
液相と気相を収容する反応容器内における、前記液相中に気泡が分散した混相流を解析するためのプロセス解析装置であって、
空間的に分割して得られる各セルについて、計算モデルに従って状態値を時系列で計算する手段を備え、
前記計算モデルは、各相についての質量保存式、運動量保存式、状態方程式、および、化学種保存式を含み、
前記状態値を時系列で計算する手段は、
前記質量保存式および前記状態保存式に基づいて、前記液相中の気泡の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記液相の体積率を計算する手段と、
前記質量保存式および前記状態保存式に基づいて、前記気相の体積率を計算する手段と、
前記化学種保存式に基づいてモル濃度を計算する手段とを含み、
前記液相の体積率を計算する手段は、前記液相内で前記気泡に同伴される第1液相領域の体積率と、前記液相内の前記第1液相領域以外の領域である第2液相領域の体積率とを計算し、
前記モル濃度を計算する手段は、前記第1液相領域のモル濃度と前記第2液相領域のモル濃度とを計算する、プロセス解析装置。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【公開番号】特開2013−45423(P2013−45423A)
【公開日】平成25年3月4日(2013.3.4)
【国際特許分類】
【出願番号】特願2011−185059(P2011−185059)
【出願日】平成23年8月26日(2011.8.26)
【出願人】(000002093)住友化学株式会社 (8,981)
【出願人】(504258527)国立大学法人 鹿児島大学 (284)
【Fターム(参考)】
【公開日】平成25年3月4日(2013.3.4)
【国際特許分類】
【出願日】平成23年8月26日(2011.8.26)
【出願人】(000002093)住友化学株式会社 (8,981)
【出願人】(504258527)国立大学法人 鹿児島大学 (284)
【Fターム(参考)】
[ Back to top ]