空気分離は、私たちの身の回りにある空気から、酸素・窒素・アルゴンなどを取り出す技術です。化学工場では窒素、医療用や溶断では酸素、溶接ではアルゴンが必要になってきます。
空気分離で最も生産量が多いのは、「深冷分離法」です。深冷分離法は、窒素、酸素、アルゴンの沸点の違いを利用して-170や-190℃といった環境の蒸留塔で分ける方式です。
このページでは、プログラムを自作して、深冷分離法の運転条件を詳細に調べましたので、是非ご覧ください。
目次
空気分離とは？空気から酸素・窒素・アルゴンを取り出す技術
深冷分離法：空気を低温で液化し、沸点の違いで分ける方式
空気分離とは、空気から窒素、酸素などを取り出す技術の総称です。空気分離の方法には、いくつか種類がありますが、最も生産量が多いものは「深冷分離法」です。
英語ではAir Separation Unitと呼ばれ、ASUと略されます。
以降では、「ASU=深冷分離法」として記事を書きます。
ASUは白い箱が目立つので、すぐ分かりますね
深冷分離法は、窒素、酸素、アルゴンの沸点の違いを利用して蒸留塔で分けています。窒素の沸点はかなり低い（大気圧で-196℃）ので、ASUはかなりの低温で蒸留を行います。そのため、「深冷」分離法と呼ばれます。複数の蒸留設備を組み合わせて高純度の製品を得るため、大規模な連続製造や酸素・窒素・アルゴンの同時回収に適しています。
また、装置内部の温度が非常に低くなっているため、太陽からの入熱を可能な限り下げるため、ASUの蒸留塔は白い箱で囲われています。また、箱の中はパーライトが詰め込まれているそうです（私は見たことはありませんが）。
原料となる空気の成分と、分離して取り出せるガス
原料となる空気には、体積比（モル比）で約78％が窒素、約21％が酸素、約0.93％がアルゴンが含まれており、その他、水分や二酸化炭素なども含まれます。水分や二酸化炭素は-190℃といった低温にすると当然固体になってしまい、バルブや配管を詰まらせてしまうため、水分と二酸化炭素は、ASUの蒸留塔の前に吸着装置によって除去されます。
| 成分 | 乾燥空気中の体積割合の目安 | 主な用途 |
|---|---|---|
| 窒素 | 約78％ | 酸化防止、パージ |
| 酸素 | 約21％ | 製鋼、医療 |
| アルゴン | 約0.93％ | 溶接、半導体製造 |
ASUでは原料となる空気はタダですが、大量の電気を消費します。入口となるタービン式圧縮機は400kWを超える大型のモーターをもつものが使われます。
深冷分離装置のプロセス・仕組み
ASUの概略フローについて
ASUは蒸留によって窒素、酸素、アルゴンを分離させています。ASUの構成は資料によって微妙に異なっていますが、「高圧塔」「低圧塔」「粗アルゴン塔」「純アルゴン塔」の４本の蒸留塔で解説している資料が多いです。
「高圧塔」「低圧塔」「粗アルゴン塔」「純アルゴン塔」を利用した基本的なASUのフローは下図のような形です。このフロー図の数値は書籍「分離プロセスの最適化とスケールアップの進め方」で紹介されているASUの製品濃度に沿うよう、私が蒸留プログラムを自作して計算させたものです。
各塔の濃度プロファイルやプログラムは後で紹介します。
| ASU | ① 原料空気 | ② 廃棄窒素 | ③ 純アルゴン | ④ 純酸素 | ⑤ 純窒素 |
|---|---|---|---|---|---|
| 流量 (mol/hr) | 10000 | 4000 | 0.5 | 1999.5 | 4000 |
| 窒素 | 78% | 95.004% | 0.070% | 0.000% | 99.996% |
| 酸素 | 21% | 3.308% | 0.034% | 98.425% | 0.02% |
| アルゴン | 1% | 1.688% | 99.896% | 1.575% | 0.01% |
計算誤差で合計が100%にならない部分があります
「分離プロセスの最適化とスケールアップの進め方」にあるASUのプロファイル
- 高圧塔は550kPa、低圧塔は130kPa
- 高圧塔 塔頂 O2：0.1ppm、塔底 O2：15%
- 低圧塔 塔頂 O2：0.1ppm、塔底 O2：99.8%
- 低圧塔でAr濃度7～13%となる中部から粗アルゴン塔へサイドカット
- 粗アルゴン塔：塔頂 Ar：98%
空気分離では原料コストは各社変わらないので、製造コストを極限にまで下げることが差別化になります。そのため、製造メーカーはいかに高純度のガスを低コストで作られるかを競っています。例えば、蒸留塔のリボイラーやコンデンサーは、蒸留塔ごとにうまくペアを作って熱回収の効率を極限にまで高めています。
概略フローは図の簡略化のため、低圧塔のリボイラーだけ高圧塔のコンデンサーとペアさせています。
しかし、ペアを作るだけでは-190℃といった超低温状態を作り出したり、維持することはできません。そのため、ASUでは膨張タービンを持っており、原料空気の一部を通すことでジュール・トムソン膨張を利用して、熱をASU系外に放出しています。だいたい0.5MPaG⇒0.05MPaGまで膨張させることで、空気は-170℃⇒-195℃程度まで温度が下がります。
高圧塔
ASUの高圧塔は、原料となる空気（0.5MPaG程度）が初めに入る蒸留塔です。塔内は高圧（0.5MPaG程度）のため、塔内温度はおよそ-175～-178℃とASUのなかでは比較的高温となっています。
| 高圧塔 | 塔底 | 塔頂 |
|---|---|---|
| 温度 | -174.99℃ | -178.97℃ |
| 窒素濃度（mol） | 62.003% | 99.996% |
| 酸素濃度（mol） | 36.271% | 0.002% |
| アルゴン濃度（mol） | 1.726% | 0.001% |
高圧塔の塔頂からは99.99%以上の純度の高い窒素が得られ、液体窒素の状態で製品タンクにためられます。半導体製造設備では2000tonや5000ton級の液体窒素貯槽があります。
2000tonを超える液体窒素の貯槽はとてつもなく存在感があります。
高圧塔では運転条件を変えても塔内組成が大きく変わることはありません。これは、窒素が酸素、アルゴンより原料空気中の量が多いことと、沸点が大きくことなることによります。
高圧塔の計算条件はこちら
| ファクター | 設定値 |
|---|---|
| 圧力 | 0.4MPaG |
| 段数 | 21 |
| FEED段 | 10 |
| FEED流量 | 9500mol/hr |
| TOP抜き出し量 | 4000mol/hr |
| 還流比 | 3.0 |
低圧塔
低圧塔の下部からは99%程度の酸素が得られます。しかし、アルゴンと酸素の沸点の差がわずかしかないため、アルゴン濃度を十分に下げることは難しいです。
上部からは、濃縮しきれない分が廃ガスとして系外から取り払われます。廃ガスの一部はASU装置入口の脱水・脱CO2の吸着器を再生させるために利用されます。
| 低圧塔 | 塔底 | 塔頂 |
|---|---|---|
| 温度 | -182.55℃ | -195.03℃ |
| 窒素濃度（mol） | 0.000% | 95.004% |
| 酸素濃度（mol） | 98.425% | 3.308% |
| アルゴン濃度（mol） | 1.575% | 1.688% |
吸着器は運転中１つと再生中１つのスイング運転をしています
また、低圧塔の中部からは一部がサイドカットされ、粗アルゴン塔へFEEDされます。サイドカットのポイントは①窒素濃度は限りなく下げる②アルゴン濃度は可能な限り高い濃度で抜き出すことが重要です。窒素がわずかでも残っていると高濃度アルゴンを手に入れることが難しくなります（窒素の方がアルゴンより軽いため）また、アルゴン濃度が低すぎても濃縮するのが大変になります。
低圧塔の計算条件
| ファクター | 設定値 |
|---|---|
| 圧力 | 0.05MPaG |
| 段数 | 25 |
| FEED段 | 2 |
| SIDE CUT段 | 15 |
| FEED流量 | 6099.5mol/hr |
| TOP抜き出し量 | 4000mol/hr |
| SIDE CUT抜き出し量 | 100mol/hr |
| 還流比 | 5.0 |
粗アルゴン塔、純アルゴン塔
低圧塔からサイドカットされたものは粗アルゴン塔、されには純アルゴン塔へ入れられます。２つのアルゴン塔のポイントは還流比が高いことです。
酸素とアルゴンの沸点は-183.0℃、-185.7℃と3℃足らずしかありません。そのため、シミュレーションでは還流比$50$という高い還流比をかけることでようやく99.9%程度のアルゴンを手に入れることができています。
| 粗アルゴン塔 | 塔底 | 塔頂 |
|---|---|---|
| 温度 | -182.89℃ | -185.26℃ |
| 窒素濃度（mol） | 0.000% | 0.035% |
| 酸素濃度（mol） | 88.005% | 5.390% |
| アルゴン濃度（mol） | 11.995% | 94.575% |
| 純アルゴン塔 | 塔底 | 塔頂 |
|---|---|---|
| 温度 | -185.11℃ | -185.40℃ |
| 窒素濃度（mol） | 0.000% | 0.070% |
| 酸素濃度（mol） | 10.746% | 0.034% |
| アルゴン濃度（mol） | 89.254% | 99.896% |
粗アルゴン塔の計算条件
| ファクター | 設定値 |
|---|---|
| 圧力 | 0.05MPa |
| 段数 | 21 |
| FEED段 | 19 |
| FEED流量 | 100mol/hr |
| TOP抜き出し量 | 1mol/hr |
| 還流比 | 50.0 |
純アルゴン塔
| ファクター | 設定値 |
|---|---|
| 圧力 | 0.05MPa |
| 段数 | 21 |
| FEED段 | 19 |
| FEED流量 | 1mol/hr |
| TOP抜き出し量 | 0.5mol/hr |
| 還流比 | 50.0 |
利用したプログラム
本ホームページでは、FORTRANを利用してASUの塔内環境を再現しています。
コードはこちら
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,K,L,II,JJ,KK,LL,NCOMP
INTEGER NPLATE_1,NPLATE_FEED_1
INTEGER NPLATE_2,NPLATE_FEED_2
INTEGER NPLATE_3,NPLATE_FEED_3
INTEGER NPLATE_4,NPLATE_FEED_4
INTEGER NPLATE_SCUT_1,NPLATE_SFEED_1
INTEGER NPLATE_SCUT_2,NPLATE_SFEED_2
INTEGER NPLATE_SCUT_3,NPLATE_SFEED_3
INTEGER NPLATE_SCUT_4,NPLATE_SFEED_4
CHARACTER(18) CHARA(10)
INTEGER MATERIAL(10)
REAL(8) P0
REAL(8) EX_RATIO
REAL(8) FEED_TOTAL
REAL(8) FEED_EX
REAL(8) FEED_STREAM_1,TOP_1,REF_1,BTM_1,SCUT_1,SFEED_1
REAL(8) FEED_STREAM_2,TOP_2,REF_2,BTM_2,SCUT_2,SFEED_2
REAL(8) FEED_STREAM_3,TOP_3,REF_3,BTM_3,SCUT_3,SFEED_3
REAL(8) FEED_STREAM_4,TOP_4,REF_4,BTM_4,SCUT_4,SFEED_4
REAL(8),ALLOCATABLE,DIMENSION(:) :: FEED_1,LIQ_1,VAP_1
REAL(8),ALLOCATABLE,DIMENSION(:) :: TEMP_1
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: xLIQ_1
REAL(8),ALLOCATABLE,DIMENSION(:) :: xLIQ_FEED_1,xLIQ_SFEED_1
REAL(8),ALLOCATABLE,DIMENSION(:) :: FEED_2,LIQ_2,VAP_2
REAL(8),ALLOCATABLE,DIMENSION(:) :: TEMP_2
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: xLIQ_2
REAL(8),ALLOCATABLE,DIMENSION(:) :: xLIQ_FEED_2,xLIQ_SFEED_2
REAL(8),ALLOCATABLE,DIMENSION(:) :: FEED_3,LIQ_3,VAP_3
REAL(8),ALLOCATABLE,DIMENSION(:) :: TEMP_3
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: xLIQ_3
REAL(8),ALLOCATABLE,DIMENSION(:) :: xLIQ_FEED_3,xLIQ_SFEED_3
REAL(8),ALLOCATABLE,DIMENSION(:) :: FEED_4,LIQ_4,VAP_4
REAL(8),ALLOCATABLE,DIMENSION(:) :: TEMP_4
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: xLIQ_4
REAL(8),ALLOCATABLE,DIMENSION(:) :: xLIQ_FEED_4,xLIQ_SFEED_4
P0 = 1013.25d0
NCOMP = 3
EX_RATIO = 0.050d0
FEED_TOTAL = 10000.0d0
FEED_EX = FEED_TOTAL * EX_RATIO
NPLATE_1 = 21
NPLATE_FEED_1 = 10
NPLATE_SCUT_1 = 10
NPLATE_SFEED_1 = 10
NPLATE_2 = 25
NPLATE_FEED_2 = 2
NPLATE_SCUT_2 = 15
NPLATE_SFEED_2 = 15
NPLATE_3 = 21
NPLATE_FEED_3 = 19
NPLATE_SCUT_3 = 19
NPLATE_SFEED_3 = 19
NPLATE_4 = 21
NPLATE_FEED_4 = 19
NPLATE_SCUT_4 = 19
NPLATE_SFEED_4 = 19
ALLOCATE (xLIQ_1(NCOMP,NPLATE_1))
ALLOCATE (xLIQ_FEED_1(NCOMP),xLIQ_SFEED_1(NCOMP))
ALLOCATE (FEED_1(NPLATE_1),LIQ_1(NPLATE_1),VAP_1(NPLATE_1))
ALLOCATE (TEMP_1(NPLATE_1))
ALLOCATE (xLIQ_2(NCOMP,NPLATE_2))
ALLOCATE (xLIQ_FEED_2(NCOMP),xLIQ_SFEED_2(NCOMP))
ALLOCATE (FEED_2(NPLATE_2),LIQ_2(NPLATE_2),VAP_2(NPLATE_2))
ALLOCATE (TEMP_2(NPLATE_2))
ALLOCATE (xLIQ_3(NCOMP,NPLATE_3))
ALLOCATE (xLIQ_FEED_3(NCOMP),xLIQ_SFEED_3(NCOMP))
ALLOCATE (FEED_3(NPLATE_3),LIQ_3(NPLATE_3),VAP_3(NPLATE_3))
ALLOCATE (TEMP_3(NPLATE_3))
ALLOCATE (xLIQ_4(NCOMP,NPLATE_4))
ALLOCATE (xLIQ_FEED_4(NCOMP),xLIQ_SFEED_4(NCOMP))
ALLOCATE (FEED_4(NPLATE_4),LIQ_4(NPLATE_4),VAP_4(NPLATE_4))
ALLOCATE (TEMP_4(NPLATE_4))
FEED_STREAM_1 = FEED_TOTAL - FEED_EX
TOP_1 = 0.40d0 * FEED_TOTAL
REF_1 = 3.0d0
SCUT_1 = 0.0d0 * FEED_TOTAL
SFEED_1 = 0.0d0 * FEED_TOTAL
BTM_1 = FEED_STREAM_1 - TOP_1 - SCUT_1 + SFEED_1
xLIQ_FEED_1(1) = 0.78d0
xLIQ_FEED_1(2) = 0.21d0
xLIQ_FEED_1(3) = 0.010d0
xLIQ_SFEED_1(1) = 0.9d0
xLIQ_SFEED_1(2) = 0.1d0
xLIQ_SFEED_1(3) = 0.0d0
FEED_STREAM_2 = BTM_1 + FEED_EX
TOP_2 = 0.4d0 * FEED_TOTAL
REF_2 = 5.0d0
SCUT_2 = 0.01d0 * FEED_TOTAL
SFEED_2 = 0.0d0 * FEED_TOTAL
BTM_2 = FEED_STREAM_2 - TOP_2 - SCUT_2 + SFEED_2
FEED_STREAM_3 = SCUT_2
TOP_3 = 0.0001d0 * FEED_TOTAL
REF_3 = 50.0d0
SCUT_3 = 0.0d0 * FEED_TOTAL
SFEED_3 = 0.0d0 * FEED_TOTAL
BTM_3 = FEED_STREAM_3 - TOP_3 - SCUT_3 + SFEED_3
FEED_STREAM_4 = TOP_3
TOP_4 = FEED_STREAM_4 * 0.90d0
REF_4 = 50.0d0
SCUT_4 = 0.0d0 * FEED_TOTAL
SFEED_4 = 0.0d0 * FEED_TOTAL
BTM_4 = FEED_STREAM_4 - TOP_4 - SCUT_4 + SFEED_4
CALL COLUMN(NCOMP,5.0d0*P0,
1 NPLATE_1,NPLATE_FEED_1,NPLATE_SCUT_1,NPLATE_SFEED_1,
1 xLIQ_FEED_1,FEED_STREAM_1,TOP_1,BTM_1,REF_1,SCUT_1,SFEED_1,
1 xLIQ_1,xLIQ_SFEED_1,VAP_1,LIQ_1,TEMP_1,
1 MATERIAL,CHARA,1)
c CALL WRITE_RESULT(NCOMP,NPLATE_1,MATERIAL,CHARA,
c 1 NPLATE_FEED_1,NPLATE_SCUT_1,NPLATE_SFEED_1,
c 1 FEED_STREAM_1,
c 1 VAP_1,LIQ_1,SCUT_1,SFEED_1,
c 1 TOP_1,TEMP_1,xLIQ_1,xLIQ_FEED_1,xLIQ_SFEED_1)
xLIQ_FEED_2(1) = (xLIQ_1(1,NPLATE_1) * BTM_1
1 + 0.780d0 * FEED_EX) / FEED_STREAM_2
xLIQ_FEED_2(2) = (xLIQ_1(2,NPLATE_1) * BTM_1
1 + 0.210d0 * FEED_EX) / FEED_STREAM_2
xLIQ_FEED_2(3) = (xLIQ_1(3,NPLATE_1) * BTM_1
1 + 0.010d0 * FEED_EX) / FEED_STREAM_2
xLIQ_SFEED_2(1)= 0.78d0
xLIQ_SFEED_2(2)= 0.21d0
xLIQ_SFEED_2(3)= 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_2,NPLATE_FEED_2,NPLATE_SCUT_2,NPLATE_SFEED_2,
1 xLIQ_FEED_2,FEED_STREAM_2,TOP_2,BTM_2,REF_2,SCUT_2,SFEED_2,
1 xLIQ_2,xLIQ_SFEED_2,VAP_2,LIQ_2,TEMP_2,
1 MATERIAL,CHARA,2)
c CALL WRITE_RESULT(NCOMP,NPLATE_2,MATERIAL,CHARA,
c 1 NPLATE_FEED_2,NPLATE_SCUT_2,NPLATE_SFEED_2,
c 1 FEED_STREAM_2,
c 1 VAP_2,LIQ_2,SCUT_2,SFEED_2,
c 1 TOP_2,TEMP_2,xLIQ_2,xLIQ_FEED_2,xLIQ_SFEED_2)
xLIQ_FEED_3(1) = xLIQ_2(1,NPLATE_SCUT_2)
xLIQ_FEED_3(2) = xLIQ_2(2,NPLATE_SCUT_2)
xLIQ_FEED_3(3) = xLIQ_2(3,NPLATE_SCUT_2)
xLIQ_SFEED_3(1) = 0.78d0
xLIQ_SFEED_3(2) = 0.21d0
xLIQ_SFEED_3(3) = 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_3,NPLATE_FEED_3,NPLATE_SCUT_3,NPLATE_SFEED_3,
1 xLIQ_FEED_3,FEED_STREAM_3,TOP_3,BTM_3,REF_3,SCUT_3,SFEED_3,
1 xLIQ_3,xLIQ_SFEED_3,VAP_3,LIQ_3,TEMP_3,
1 MATERIAL,CHARA,3)
c CALL WRITE_RESULT(NCOMP,NPLATE_3,MATERIAL,CHARA,
c 1 NPLATE_FEED_3,NPLATE_SCUT_3,NPLATE_SFEED_3,
c 1 FEED_STREAM_3,
c 1 VAP_3,LIQ_3,SCUT_3,SFEED_3,
c 1 TOP_3,TEMP_3,xLIQ_3,xLIQ_FEED_3,xLIQ_SFEED_3)
FEED_STREAM_4 = TOP_3
TOP_4 = FEED_STREAM_4 * 0.5d0
REF_4 = 50.0d0
SCUT_4 = 0.0d0
SFEED_4 = 0.0d0
BTM_4 = FEED_STREAM_4 - TOP_4 - SCUT_4 + SFEED_4
xLIQ_FEED_4(1) = xLIQ_3(1,1)
xLIQ_FEED_4(2) = xLIQ_3(2,1)
xLIQ_FEED_4(3) = xLIQ_3(3,1)
xLIQ_SFEED_4(1) = 0.78d0
xLIQ_SFEED_4(2) = 0.21d0
xLIQ_SFEED_4(3) = 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_4,NPLATE_FEED_4,NPLATE_SCUT_4,NPLATE_SFEED_4,
1 xLIQ_FEED_4,FEED_STREAM_4,TOP_4,BTM_4,REF_4,SCUT_4,SFEED_4,
1 xLIQ_4,xLIQ_SFEED_4,VAP_4,LIQ_4,TEMP_4,
1 MATERIAL,CHARA,4)
DO I = 1,2
FEED_STREAM_2 = BTM_1 + FEED_EX + BTM_3 + BTM_4
TOP_2 = 0.4d0 * FEED_TOTAL
REF_2 = 5.0d0
SCUT_2 = 0.01d0 * FEED_TOTAL
SFEED_2 = 0.0d0 * FEED_TOTAL
BTM_2 = FEED_STREAM_2 - TOP_2 - SCUT_2 + SFEED_2
xLIQ_FEED_2(1) = (xLIQ_1(1,NPLATE_1) * BTM_1
1 + 0.780d0 * FEED_EX
1 + xLIQ_3(1,NPLATE_3) * BTM_3
1 + xLIQ_4(1,NPLATE_4) * BTM_4)
1 / FEED_STREAM_2
xLIQ_FEED_2(2) = (xLIQ_1(2,NPLATE_1) * BTM_1
1 + 0.210d0 * FEED_EX
1 + xLIQ_3(2,NPLATE_3) * BTM_3
1 + xLIQ_4(2,NPLATE_4) * BTM_4)
1 / FEED_STREAM_2
xLIQ_FEED_2(3) = (xLIQ_1(3,NPLATE_1) * BTM_1
1 + 0.010d0 * FEED_EX
1 + xLIQ_3(3,NPLATE_3) * BTM_3
1 + xLIQ_4(3,NPLATE_4) * BTM_4)
1 / FEED_STREAM_2
xLIQ_SFEED_2(1)= 0.78d0
xLIQ_SFEED_2(2)= 0.21d0
xLIQ_SFEED_2(3)= 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_2,NPLATE_FEED_2,NPLATE_SCUT_2,NPLATE_SFEED_2,
1 xLIQ_FEED_2,FEED_STREAM_2,TOP_2,BTM_2,REF_2,SCUT_2,SFEED_2,
1 xLIQ_2,xLIQ_SFEED_2,VAP_2,LIQ_2,TEMP_2,
1 MATERIAL,CHARA,2)
xLIQ_FEED_3(1) = xLIQ_2(1,NPLATE_SCUT_2)
xLIQ_FEED_3(2) = xLIQ_2(2,NPLATE_SCUT_2)
xLIQ_FEED_3(3) = xLIQ_2(3,NPLATE_SCUT_2)
xLIQ_SFEED_3(1) = 0.78d0
xLIQ_SFEED_3(2) = 0.21d0
xLIQ_SFEED_3(3) = 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_3,NPLATE_FEED_3,NPLATE_SCUT_3,NPLATE_SFEED_3,
1 xLIQ_FEED_3,FEED_STREAM_3,TOP_3,BTM_3,REF_3,SCUT_3,SFEED_3,
1 xLIQ_3,xLIQ_SFEED_3,VAP_3,LIQ_3,TEMP_3,
1 MATERIAL,CHARA,3)
FEED_STREAM_4 = TOP_3
TOP_4 = FEED_STREAM_4 * 0.5d0
REF_4 = 50.0d0
SCUT_4 = 0.0d0
SFEED_4 = 0.0d0
BTM_4 = FEED_STREAM_4 - TOP_4 - SCUT_4 + SFEED_4
xLIQ_FEED_4(1) = xLIQ_3(1,1)
xLIQ_FEED_4(2) = xLIQ_3(2,1)
xLIQ_FEED_4(3) = xLIQ_3(3,1)
xLIQ_SFEED_4(1) = 0.78d0
xLIQ_SFEED_4(2) = 0.21d0
xLIQ_SFEED_4(3) = 0.01d0
CALL COLUMN(NCOMP,1.05d0*P0,
1 NPLATE_4,NPLATE_FEED_4,NPLATE_SCUT_4,NPLATE_SFEED_4,
1 xLIQ_FEED_4,FEED_STREAM_4,TOP_4,BTM_4,REF_4,SCUT_4,SFEED_4,
1 xLIQ_4,xLIQ_SFEED_4,VAP_4,LIQ_4,TEMP_4,
1 MATERIAL,CHARA,4)
ENDDO
CALL WRITE_RESULT(NCOMP,NPLATE_1,MATERIAL,CHARA,
1 NPLATE_FEED_1,NPLATE_SCUT_1,NPLATE_SFEED_1,
1 FEED_STREAM_1,
1 VAP_1,LIQ_1,SCUT_1,SFEED_1,
1 TOP_1,TEMP_1,xLIQ_1,xLIQ_FEED_1,xLIQ_SFEED_1,1)
CALL WRITE_RESULT(NCOMP,NPLATE_2,MATERIAL,CHARA,
1 NPLATE_FEED_2,NPLATE_SCUT_2,NPLATE_SFEED_2,
1 FEED_STREAM_2,
1 VAP_2,LIQ_2,SCUT_2,SFEED_2,
1 TOP_2,TEMP_2,xLIQ_2,xLIQ_FEED_2,xLIQ_SFEED_2,2)
CALL WRITE_RESULT(NCOMP,NPLATE_3,MATERIAL,CHARA,
1 NPLATE_FEED_3,NPLATE_SCUT_3,NPLATE_SFEED_3,
1 FEED_STREAM_3,
1 VAP_3,LIQ_3,SCUT_3,SFEED_3,
1 TOP_3,TEMP_3,xLIQ_3,xLIQ_FEED_3,xLIQ_SFEED_3,3)
CALL WRITE_RESULT(NCOMP,NPLATE_4,MATERIAL,CHARA,
1 NPLATE_FEED_4,NPLATE_SCUT_4,NPLATE_SFEED_4,
1 FEED_STREAM_4,
1 VAP_4,LIQ_4,SCUT_4,SFEED_4,
1 TOP_4,TEMP_4,xLIQ_4,xLIQ_FEED_4,xLIQ_SFEED_4,4)
END PROGRAM
SUBROUTINE COLUMN(NCOMP,P0,
1 NPLATE,NPLATE_FEED,NPLATE_SCUT,NPLATE_SFEED,
1 xLIQ_FEED,FEED_STREAM,TOP,BTM,REF,SCUT,SFEED,
1 xLIQ,xLIQ_SFEED,VAP,LIQ,TEMP,
1 MATERIAL,CHARA,IFLAG)
IMPLICIT NONE
INTEGER IW,I,J,II,JJ,NCOMP,N_JACOB,ITE
INTEGER IFLAG,INFO
INTEGER NPLATE,NPLATE_FEED,NPLATE_SFEED,NPLATE_SCUT
REAL(8) FEED_STREAM,TOP,BTM,REF,SCUT,SFEED
REAL(8) P0,allError,delta,ErrorTolerance
REAL(8) temp_work
REAL(8) work(3)
REAL(8) FAC_NOMAL
REAL(8) FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE)
REAL(8) S_CUT(NPLATE),S_FEED(NPLATE)
REAL(8) xLIQ_FEED(NCOMP),xLIQ_SFEED(NCOMP),xLIQ(NCOMP,NPLATE)
REAL(8) M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
REAL(8) M_SCUT_LIQ(NCOMP,NPLATE)
REAL(8) M_SFEED_LIQ(NCOMP,NPLATE)
REAL(8) M_FEED(NCOMP,NPLATE),M_TOP(NCOMP)
REAL(8) TEMP(NPLATE)
REAL(8) ANTOINE(4,NCOMP),WILLSON(NCOMP,NCOMP)
REAL(8) CP(3,NCOMP),TB(NCOMP)
REAL(8) D_LIQ(NPLATE),D_VAP(NPLATE)
INTEGER MATERIAL(10)
CHARACTER(18) CHARA(10)
REAL(8) E(NPLATE),M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
REAL(8) D_M_VAP(NCOMP,NPLATE),D_M_LIQ(NCOMP,NPLATE)
REAL(8) D_TEMP(NPLATE)
REAL(8) D_xLIQ(NCOMP,NPLATE)
REAL(8) ALPHA(NCOMP,NPLATE),D_ALPHA(NCOMP,NPLATE)
REAL(8) JACOB((2*NCOMP+1)*NPLATE,(2*NCOMP+1)*NPLATE)
REAL(8) FX((2*NCOMP+1)*NPLATE)
INTEGER IPIV((2*NCOMP+1)*NPLATE)
REAL(8) working(NPLATE)
REAL(8) D_E(NPLATE),D_S(NPLATE)
REAL(8) D_M(NCOMP,NPLATE),D_Q(NCOMP,NPLATE)
IW = 6
delta = 0.010d0
ErrorTolerance = 0.10d0**(6.0d0)
FAC_NOMAL = FEED_STREAM
MATERIAL(1) = 1
MATERIAL(2) = 2
MATERIAL(3) = 3
CHARA(1) = "N2"
CHARA(2) = "O2"
CHARA(3) = "Ar"
FEED_STREAM = 100.0d0
TOP = TOP / FAC_NOMAL * 100.0d0
BTM = BTM / FAC_NOMAL * 100.0d0
SCUT = SCUT / FAC_NOMAL * 100.0d0
SFEED = SFEED / FAC_NOMAL * 100.0d0
CALL SET_PROPERTY(NCOMP,MATERIAL,ANTOINE,WILLSON,CP,TB)
CALL SET_xLIQ_START(NCOMP,NPLATE,
1 xLIQ_FEED,xLIQ_SFEED,xLIQ,
1 TB,NPLATE_FEED,FEED_STREAM,TOP,BTM)
CALL SET_INITIAL(NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 SCUT,SFEED,NPLATE_SCUT,NPLATE_SFEED,
1 TB,xLIQ_FEED,xLIQ_SFEED,
1 xLIQ,VAP,LIQ,FEED,TOP,BTM,S_CUT,S_FEED,
1 REF,M_VAP,M_LIQ,M_SCUT_LIQ,M_SFEED_LIQ,M_FEED,M_TOP,TEMP,
1 IFLAG)
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0,TEMP,xLIQ,ALPHA)
N_JACOB = (2*NCOMP+1)*NPLATE
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,'(15x,a12,i3,a12)') "START TOWER",IFLAG, " CALCULATION"
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
CCC LOOP
CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IF(IFLAG.EQ.2) THEN
ITE = 2
ENDIF
DO ITE = 1,10000
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0
1 ,TEMP,xLIQ,ALPHA)
CALL CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM,
1 M_VAP,M_LIQ,M_FEED,M_TOP,M_SCUT_LIQ,M_SFEED_LIQ,
1 VAP,LIQ,S_CUT,S_FEED,ALPHA,CP,E,M,Q)
DO II = 1,(2*NCOMP+1)*NPLATE
DO JJ = 1,(2*NCOMP+1)*NPLATE
JACOB(II,JJ) = 0.0d0
ENDDO
ENDDO
DO II = 1,(2*NCOMP+1)*NPLATE
FX(II) = 0.0d0
ENDDO
DO JJ = 1,(2*NCOMP+1)*NPLATE
CALL CALC_DEL_VEC(JJ,NCOMP,NPLATE,delta
1 ,TEMP,M_VAP,M_LIQ,D_TEMP,D_M_VAP,D_M_LIQ)
CALL CALC_LIQ_VAP_xLIQ(NCOMP,NPLATE
1 ,D_VAP,D_LIQ,S_CUT,S_FEED,D_M_VAP,D_M_LIQ,D_xLIQ)
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0
1 ,D_TEMP,D_xLIQ,D_ALPHA)
CALL CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM
1 ,D_M_VAP,D_M_LIQ,M_FEED,M_TOP,M_SCUT_LIQ,M_SFEED_LIQ
1 ,D_VAP,D_LIQ,S_CUT,S_FEED
1 ,D_ALPHA,CP,D_E,D_M,D_Q)
CALL CALC_JACOB_COL(JJ,NCOMP,NPLATE
1 ,E,M,Q,D_E,D_M,D_Q,delta,JACOB,FX)
ENDDO
CALL DGESV(N_JACOB,1,JACOB,N_JACOB,IPIV,FX,N_JACOB,INFO)
CALL CALC_NEW_VEC(NCOMP,NPLATE
1 ,TEMP,M_VAP,M_LIQ,M_SCUT_LIQ,M_TOP
1 ,VAP,LIQ,S_CUT,FX,xLIQ,TOP)
c IF(mod(ITE,10).LE.5) THEN
c IF(IFLAG.EQ.2) THEN
c xLIQ(2,NPLATE) = 1.0d0
c xLIQ(3,NPLATE) = 0.0d0
c xLIQ(2,NPLATE-1) = 1.0d0
c xLIQ(3,NPLATE-1) = 0.0d0
c
c ENDIF
c ENDIF
DO J = 1,NPLATE
working(J) = abs(E(J))
DO I = 1,NCOMP
working(J) = working(J) + abs(M(I,J)) + abs(Q(I,J))
ENDDO
ENDDO
allError = 0.0d0
DO J = 1,NPLATE
allError = allError + working(J)
ENDDO
IF(allError .LE. REAL(NPLATE) * ErrorTolerance) THEN
GOTO 100
ELSE IF(ITE<10) THEN
GOTO 300
ELSE IF(ITE .LE. 2000 .AND. MOD(ITE,10) .EQ.0) THEN
GOTO 300
ELSE
GOTO 300
ENDIF
200 WRITE(IW,'(15x,i10,10x,e15.5)')
1 ITE,allError/real(NPLATE)
300 CONTINUE
CCCCCCCCCCCCCCCCCCCC
CCCC
CCCC ITE LOOP END
ENDDO
CCCC
CCCCCCCCCCCCCCCCCCCC
100 WRITE(IW,'(10x,a5,i10,10x,e15.5,a23)')
1 "ITE",ITE,allError/real(NPLATE)
1 ,"< 1.0 E-05 (Converged)"
FEED_STREAM = FAC_NOMAL
TOP = TOP * FAC_NOMAL / 100.0d0
BTM = BTM * FAC_NOMAL / 100.0d0
SCUT = SCUT * FAC_NOMAL / 100.0d0
SFEED = SFEED * FAC_NOMAL / 100.0d0
DO J = 1, NPLATE
LIQ(J) = LIQ(J) * FAC_NOMAL / 100.0d0
VAP(J) = VAP(J) * FAC_NOMAL / 100.0d0
ENDDO
END
SUBROUTINE CALC_ALPHA
1 (NCOMP,NPLATE,ANTOINE,WILLSON,P0,TEMP,xLIQ,ALPHA)
IMPLICIT NONE
INTEGER IW,I,J,K,II,JJ,NCOMP,N_JACOB,ITE
INTEGER IFLAG,INFO
INTEGER NPLATE,NPLATE_FEED,NPLATE_SFEED,NPLATE_SCUT
REAL(8) FEED_STREAM,TOP,BTM,REF,SCUT,SFEED
REAL(8) workingFirst,workingBunbo,workingThird
REAL(8) P0
REAL(8) xLIQ(NCOMP,NPLATE)
REAL(8) ALPHA(NCOMP,NPLATE)
REAL(8) ANTOINE(4,NCOMP)
REAL(8) WILLSON(NCOMP,NCOMP)
REAL(8) P(NCOMP)
REAL(8) GAM(NCOMP,NPLATE)
REAL(8) TEMP(NPLATE)
DO J = 1, NPLATE
TEMP(J) = TEMP(J) + 273.150d0
DO I = 1, NCOMP
P(I) = ANTOINE(1,I) + ANTOINE(2,I)/TEMP(J)
1 + ANTOINE(3,I)*TEMP(J) + ANTOINE(4,I)*log(TEMP(J))
P(I) = exp(P(I)) * 1013.25d0
ENDDO
TEMP(J) = TEMP(J) - 273.150d0
DO I = 1 , NCOMP
ALPHA(I,J) = P(I)/P0
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE SET_PROPERTY(NCOMP,MATERIAL,ANTOINE,WILLSON,CP,TB)
IMPLICIT NONE
INTEGER IW,I,J,K,II,JJ,NCOMP,N_JACOB,ITE
REAL(8) DB_ANTOINE(4,10),DB_WILLSON(10,10)
REAL(8) DB_CP(3,10),DB_TB(10)
REAL(8) ANTOINE(4,NCOMP),WILLSON(NCOMP,NCOMP)
REAL(8) CP(3,NCOMP),TB(NCOMP)
INTEGER MATERIAL(10)
DB_ANTOINE(1,1) = 44.20590d0
DB_ANTOINE(2,1) = -1058.65d0
DB_ANTOINE(3,1) = 0.041088d0
DB_ANTOINE(4,1) = -7.750d0
DB_ANTOINE(1,2) = 41.8710d0
DB_ANTOINE(2,2) = -1224.02d0
DB_ANTOINE(3,2) = 0.0304130d0
DB_ANTOINE(4,2) = -6.895d0
DB_ANTOINE(1,3) = 41.8128d0
DB_ANTOINE(2,3) = -1170.9d0
DB_ANTOINE(3,3) = 0.032033d0
DB_ANTOINE(4,3) = -6.98d0
DB_WILLSON(1,1) = 1.00000d0
DB_WILLSON(1,2) = 0.99990d0
DB_WILLSON(1,3) = 0.38208d0
DB_WILLSON(2,1) = 1.00000d0
DB_WILLSON(2,2) = 1.00000d0
DB_WILLSON(2,3) = 0.05840d0
DB_WILLSON(3,1) = 1.09178d0
DB_WILLSON(3,2) = 0.99900d0
DB_WILLSON(3,3) = 1.00000d0
DB_CP(1,1) = 55.8d0
DB_CP(2,1) = 0.030d0
DB_CP(3,1) = 0.045d0
c
DB_CP(1,2) = 68.2d0
DB_CP(2,2) = 0.029d0
DB_CP(3,2) = 0.053d0
c
DB_CP(1,3) = 65.2d0
DB_CP(2,3) = 0.021d0
DB_CP(3,3) = 0.032d0
DB_TB(1) = -150.0d0
DB_TB(2) = -140.0d0
DB_TB(3) = -100.0d0
DO I = 1,NCOMP
ANTOINE(1,I) = DB_ANTOINE(1,MATERIAL(I))
ANTOINE(2,I) = DB_ANTOINE(2,MATERIAL(I))
ANTOINE(3,I) = DB_ANTOINE(3,MATERIAL(I))
ANTOINE(4,I) = DB_ANTOINE(4,MATERIAL(I))
ENDDO
DO I = 1,NCOMP
DO J = 1,NCOMP
WILLSON(I,J) = DB_WILLSON(MATERIAL(I),MATERIAL(J))
ENDDO
ENDDO
DO I = 1,NCOMP
CP(1,I) = DB_CP(1,MATERIAL(I))
CP(2,I) = DB_CP(2,MATERIAL(I))
CP(3,I) = DB_CP(3,MATERIAL(I))
ENDDO
DO I = 1,NCOMP
TB(I) = DB_TB(MATERIAL(I))
TB(I) = DB_TB(MATERIAL(I))
TB(I) = DB_TB(MATERIAL(I))
ENDDO
ENDSUBROUTINE
SUBROUTINE SET_INITIAL(NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 SCUT,SFEED,NPLATE_SCUT,NPLATE_SFEED,
1 TB,xLIQ_FEED,xLIQ_SFEED,
1 xLIQ,VAP,LIQ,FEED,TOP,BTM,S_CUT,S_FEED,
1 REF,M_VAP,M_LIQ,M_SCUT_LIQ,M_SFEED_LIQ,M_FEED,M_TOP,TEMP,
1 IFLAG)
IMPLICIT NONE
INTEGER IW,I,J,II,JJ,NCOMP,N_JACOB,ITE
INTEGER IFLAG,INFO
INTEGER NPLATE,NPLATE_FEED,NPLATE_SFEED,NPLATE_SCUT
REAL(8) FEED_STREAM,TOP,BTM,REF,SCUT,SFEED
REAL(8) P0,allError,delta,ErrorTolerance
REAL(8) S_CUT_TOTAL,S_FEED_TOTAL,T_MAX,T_MIN
REAL(8) FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE)
REAL(8) S_CUT(NPLATE),S_FEED(NPLATE)
REAL(8) xLIQ_FEED(NCOMP),xLIQ_SFEED(NCOMP),xLIQ(NCOMP,NPLATE)
REAL(8) M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
REAL(8) M_SCUT_LIQ(NCOMP,NPLATE)
REAL(8) M_SFEED_LIQ(NCOMP,NPLATE)
REAL(8) M_FEED(NCOMP,NPLATE),M_TOP(NCOMP)
REAL(8) TEMP(NPLATE)
REAL(8) TB(NCOMP)
DO I = 1, NPLATE
FEED(I) = 0.0d0
S_CUT(I) = 0.0d0
S_FEED(I) = 0.0d0
ENDDO
FEED(NPLATE_FEED)=FEED_STREAM
S_CUT(NPLATE_SCUT)=SCUT
S_FEED(NPLATE_SFEED)=SFEED
VAP(1) = 0.0d0
DO I = 2, NPLATE
VAP(I) = (REF + 1.0d0) * TOP
ENDDO
DO I = 1, NPLATE-1
IF(I.EQ.1) THEN
LIQ(1) = REF * TOP
ELSE
LIQ(I) = LIQ(I-1) + FEED(I) - S_CUT(I) + S_FEED(I)
ENDIF
ENDDO
c WRITE(*,'(a30)') "INITIAL"
c DO I=1,NPLATE
c WRITE(*,'(i5,4f10.5)') I,LIQ(I),FEED(I),S_CUT(I),S_FEED(I)
c ENDDO
DO I = 1,NPLATE
S_CUT(I) = 0.0d0
S_FEED(I) = 0.0d0
ENDDO
S_CUT(NPLATE_SCUT) = SCUT
S_FEED(NPLATE_SFEED) = SFEED
S_CUT_TOTAL = 0.0d0
S_FEED_TOTAL = 0.0d0
DO I = 1, NPLATE
S_CUT_TOTAL = S_CUT_TOTAL + S_CUT(I)
S_FEED_TOTAL = S_FEED_TOTAL + S_FEED(I)
ENDDO
LIQ(NPLATE) = FEED_STREAM-TOP-S_CUT_TOTAL+S_FEED_TOTAL
T_MIN = -185.0d0
T_MAX = -180.0d0
DO I = 2, NCOMP
IF(TB(I).LT.T_MIN) THEN
T_MIN = TB(I)
ENDIF
IF(TB(I).GT.T_MAX) THEN
T_MAX = TB(I)
ENDIF
ENDDO
DO I = 1, NPLATE
TEMP(I) = T_MIN + (T_MAX - T_MIN) / NPLATE * real(I)
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
M_VAP (I,J) = xLIQ_FEED(I)*VAP(J)
M_LIQ (I,J) = xLIQ_FEED(I)*LIQ(J)
M_SCUT_LIQ(I,J) = xLIQ_FEED(I)*S_CUT(J)
M_SFEED_LIQ(I,J) = xLIQ_SFEED(I)*S_FEED(J)
M_FEED(I,J) = xLIQ_FEED(I)*FEED(J)
ENDDO
ENDDO
IF(IFLAG.EQ.2) THEN
DO J = 1,NPLATE
DO I = 1,NCOMP
M_VAP (2,J) = real(J)/real(NPLATE) *VAP(J)
M_LIQ (2,J) = real(J)/real(NPLATE) *LIQ(J)
M_VAP (3,J) = real(NPLATE-J+1)/real(NPLATE) *VAP(J)
M_LIQ (3,J) = real(NPLATE-J+1)/real(NPLATE) *LIQ(J)
ENDDO
ENDDO
c DO J = 1, NPLATE
c WRITE(*,'(i5,2f10.5)') J,M_VAP(2,J),M_VAP(3,J)
c ENDDO
c STOP
ENDIF
DO I = 1,NCOMP
M_TOP(I) = TOP*xLIQ(I,1)
ENDDO
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM,
1 M_VAP,M_LIQ,M_FEED,M_TOP,M_SCUT_LIQ,M_SFEED_LIQ,
1 VAP,LIQ,S_CUT,S_FEED,ALPHA,CP,E,M,Q)
IMPLICIT NONE
INTEGER I,J
INTEGER NPLATE,NCOMP
DOUBLE PRECISION H_VAP(NPLATE),H_LIQ(NPLATE),H_FEED(NPLATE)
DOUBLE PRECISION H_SCUT(NPLATE),H_SFEED(NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_SCUT_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_SFEED_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_TOP(NCOMP),M_FEED(NCOMP,NPLATE)
DOUBLE PRECISION ALPHA(NCOMP,NPLATE)
DOUBLE PRECISION LIQ(NPLATE),VAP(NPLATE)
DOUBLE PRECISION S_CUT(NPLATE),S_FEED(NPLATE)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION CP(3,NCOMP)
DOUBLE PRECISION E(NPLATE),M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
DOUBLE PRECISION REF,TOP,BTM
DO J = 1,NPLATE
LIQ(J) = 0.0d0
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
ENDDO
ENDDO
DO J = 1, NPLATE
H_VAP(J) = 0.0d0
H_LIQ(J) = 0.0d0
H_FEED(J) = 0.0d0
H_SCUT(J) = 0.0d0
H_SFEED(J) = 0.0d0
DO I = 1, NCOMP
H_VAP(J) = H_VAP(J) + M_VAP(I,J)*(CP(1,I)+TEMP(J)*CP(2,I))
H_LIQ(J) = H_LIQ(J) + M_LIQ(I,J) * CP(3,I)
H_FEED(J) = H_FEED(J) + M_FEED(I,J)* CP(3,I)
H_SCUT(J) = H_SCUT(J) + M_SCUT_LIQ(I,J)* CP(3,I)
H_SFEED(J) = H_SFEED(J) + M_SFEED_LIQ(I,J)* CP(3,I)
ENDDO
ENDDO
DO I=1,NCOMP
M(I,1) = M_VAP(I,1) + M_LIQ(I,1) + M_TOP(I)
1 - (M_VAP(I,2))
DO J = 2,NPLATE-1
M(I,J) = M_VAP(I,J) + M_LIQ(I,J)
1 - (M_VAP(I,J+1) + M_LIQ(I,J-1)
1 + M_FEED(I,J) - M_SCUT_LIQ(I,J)+ M_SFEED_LIQ(I,J))
ENDDO
M(I,NPLATE) = M_VAP(I,NPLATE) + M_LIQ(I,NPLATE)
1 - M_LIQ(I,NPLATE-1)
ENDDO
E(1) = LIQ(1) - REF * ( TOP)
E(NPLATE) = LIQ(NPLATE) - BTM
DO I = 2, NPLATE-1
E(I) = H_VAP(I) + H_LIQ(I)
1 - (H_VAP(I+1) + H_LIQ(I-1)
1 -H_SCUT(I) + H_SFEED(I) + H_FEED(I))
ENDDO
E(NPLATE) = LIQ(NPLATE) - BTM
Q(1,1) = LIQ(1)
DO I = 1,NCOMP
Q(1,1) = Q(1,1) - ALPHA(I,1)*M_LIQ(I,1)
ENDDO
DO I = 2,NCOMP
Q(I,1) = M_LIQ(I,1)*VAP(2) - LIQ(1)*M_VAP(I,2)
ENDDO
DO J = 2,NPLATE
DO I = 1,NCOMP
Q(I,J) = ALPHA(I,J)*VAP(J)/LIQ(J)*M_LIQ(I,J) - M_VAP(I,J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_LIQ_VAP_xLIQ(NCOMP,NPLATE
1 ,VAP,LIQ,S_CUT,S_FEED,M_VAP,M_LIQ,xLIQ)
IMPLICIT NONE
INTEGER I,J
INTEGER NPLATE,NCOMP
DOUBLE PRECISION VAP(NPLATE),LIQ(NPLATE)
DOUBLE PRECISION S_CUT(NPLATE),S_FEED(NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
c DO J = 1,NPLATE
c LIQ(J) = LIQ(J) - S_CUT(J)
c ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_JACOB_COL(JJ,NCOMP,NPLATE,E,M,Q,D_E,D_M,D_Q,
1 delta,JACOB,FX)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION E(NPLATE)
DOUBLE PRECISION M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
DOUBLE PRECISION D_E(NPLATE)
DOUBLE PRECISION D_M(NCOMP,NPLATE),D_Q(NCOMP,NPLATE)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DOUBLE PRECISION JACOB((2*NCOMP+1)*NPLATE,(2*NCOMP+1)*NPLATE)
DOUBLE PRECISION FX((2*NCOMP+1)*NPLATE)
DOUBLE PRECISION delta
c
IF(JJ.LE.NPLATE) THEN
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
ELSE
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE - NPLATE * NCOMP
JJJ = (JJJ+NCOMP-1)/NCOMP
c
ENDIF
IF(JJ.LE.NPLATE) THEN
DO KK = 1,NPLATE
JACOB(KK,JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ ) = -E(JJ)
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
DO KK = 1,NPLATE
JACOB(KK, JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL ,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ) = -M(II,JJJ)
c
c
c
ELSE
DO KK = 1,NPLATE
JACOB(KK, JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL
1 ,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL
1 ,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ) = -Q(II,JJJ)
ENDIF
ENDSUBROUTINE
SUBROUTINE CALC_NEW_VEC(NCOMP,NPLATE
1 ,TEMP,M_VAP,M_LIQ,M_SCUT_LIQ,M_TOP
1 ,VAP,LIQ,S_CUT,FX,xLIQ,TOP)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_SCUT_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_TOP(NCOMP)
DOUBLE PRECISION S_CUT(NPLATE)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION VAP(NPLATE),LIQ(NPLATE)
DOUBLE PRECISION FX((2*NCOMP+1)*NPLATE)
DOUBLE PRECISION TOP
DO J = 1,NPLATE
TEMP(J) = TEMP(J)
1 + FX(J)/100.0d0
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
M_VAP(I,J) = M_VAP(I,J)
1 + FX(NPLATE +NCOMP*(J-1)+I)/100.0d0
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
M_LIQ(I,J) = M_LIQ(I,J)
1 + FX((NCOMP+1)*NPLATE+NCOMP*(J-1)+I)/100.0d0
ENDDO
ENDDO
DO I = 1, NCOMP
M_TOP(I) = TOP * xLIQ(I,1)
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
IF(M_VAP(I,J).LE.0.0d0) THEN
M_VAP(I,J) = 0.0d0
ENDIF
IF(M_LIQ(I,J).LE.0.0d0) THEN
M_LIQ(I,J) = 0.0d0
ENDIF
M_VAP(I,1) = 0.0d0
ENDDO
ENDDO
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
M_SCUT_LIQ(I,J) = xLIQ(I,J)*S_CUT(J)
ENDDO
ENDDO
c DO I = 1,NPLATE
c WRITE(*,'(i5,f10.5,5x,3f10.5)') I, S_CUT(I),
c 1 (M_SCUT_LIQ(J,I),J=1,3)
c ENDDO
c STOP
ENDSUBROUTINE
SUBROUTINE CALC_DEL_VEC(JJ,NCOMP,NPLATE,delta
1 ,TEMP,M_VAP,M_LIQ,D_TEMP,D_M_VAP,D_M_LIQ)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION D_M_VAP(NCOMP,NPLATE),D_M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION D_TEMP(NPLATE)
DOUBLE PRECISION delta
DO I = 1,NCOMP
DO J = 1,NPLATE
D_M_VAP(I,J) = M_VAP(I,J)
D_M_LIQ(I,J) = M_LIQ(I,J)
ENDDO
ENDDO
DO J = 1,NPLATE
D_TEMP(J) = TEMP(J)
ENDDO
IF(JJ.LE.NPLATE) THEN
D_TEMP(JJ) = D_TEMP(JJ) + delta
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
D_M_VAP(II,JJJ) = D_M_VAP(II,JJJ) + delta
ELSE
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE - NPLATE * NCOMP
JJJ = (JJJ+NCOMP-1)/NCOMP
D_M_LIQ(II,JJJ) = D_M_LIQ(II,JJJ) + delta
ENDIF
ENDSUBROUTINE
SUBROUTINE SET_xLIQ_START(NCOMP,NPLATE,
1 xLIQ_FEED,xLIQ_SFEED,xLIQ,
1 TB,NPLATE_FEED,FEED_STREAM,TOP,BTM)
IMPLICIT NONE
INTEGER I,J,K,LL,JJJ
INTEGER NPLATE,NCOMP,NPLATE_FEED
INTEGER N_TB(NCOMP)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DOUBLE PRECISION xLIQ_FEED(NCOMP),xLIQ_SFEED(NCOMP)
DOUBLE PRECISION TB(NCOMP),TB_work(NCOMP)
DOUBLE PRECISION FEED_STREAM,TOP,BTM
DOUBLE PRECISION T_SAVE
DO I = 1, NPLATE
DO J = 1, NCOMP
xLIQ(J,I) = xLIQ_FEED(J)
ENDDO
ENDDO
DO I = 1 , NCOMP
N_TB(I) = I
TB_work(I) = TB(I)
ENDDO
DO I = 1, NCOMP-1
DO J = I+1, NCOMP
IF(TB_work(I).GT.TB_work(J)) THEN
T_SAVE = TB_work(I)
TB_work(I) = TB_work(J)
TB_work(J) = T_SAVE
K = N_TB(I)
N_TB(I) = N_TB(J)
N_TB(J) = K
ENDIF
ENDDO
ENDDO
DO I = 1, NCOMP/2
J = N_TB(I)
xLIQ(J,1) = FEED_STREAM/TOP*xLIQ_FEED(J)
IF(xLIQ(J,1).GE.1.0d0) THEN
xLIQ(J,1) = 1.0d0
ENDIF
xLIQ(J,NPLATE) = 0.00d0
ENDDO
DO I = NCOMP/2+1,NCOMP
J = N_TB(I)
xLIQ(J,1) = 0.0d0
xLIQ(J,NPLATE) = FEED_STREAM/BTM*xLIQ_FEED(J)
IF(xLIQ(J,NPLATE).GE.1.0d0) THEN
xLIQ(J,NPLATE) = 1.0d0
ENDIF
ENDDO
DO J = 2,NPLATE_FEED
DO I = 1,NCOMP
xLIQ(I,J) = (xLIQ(I,NPLATE_FEED)-xLIQ(I,1))
1 / real(NPLATE_FEED-1)*real(J) + xLIQ(I,1)
ENDDO
ENDDO
DO J = NPLATE_FEED,NPLATE-1
DO I = 1,NCOMP
xLIQ(I,J) = (xLIQ(I,NPLATE)-xLIQ(I,NPLATE_FEED))
1 / real(NPLATE-NPLATE_FEED)*real(J-NPLATE_FEED) + xLIQ_FEED(I)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE WRITE_RESULT(NCOMP,NPLATE,MATERIAL,CHARA,
1 NPLATE_FEED,NPLATE_SCUT,NPLATE_SFEED,FEED_STREAM,
1 VAP,LIQ,SCUT,SFEED,
1 TOP,TEMP,xLIQ,xLIQ_FEED,xLIQ_SFEED,IFLAG)
IMPLICIT NONE
INTEGER IW,I,J,NCOMP,IFLAG
INTEGER NPLATE,NPLATE_FEED,NPLATE_SFEED,NPLATE_SCUT
REAL(8) FEED_STREAM,TOP,BTM,REF,SCUT,SFEED
REAL(8) working1,working2,working3,working4,working5
REAL(8) FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE)
REAL(8) xLIQ_FEED(NCOMP),xLIQ(NCOMP,NPLATE)
REAL(8) xLIQ_SFEED(NCOMP)
REAL(8) ALPHA(NCOMP,NPLATE)
REAL(8) TEMP(NPLATE)
REAL(8) E(NPLATE),S(NPLATE)
REAL(8) M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
REAL(8) M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
REAL(8) working(8)
REAL(8) M_SCUT_LIQ(NCOMP,NPLATE)
REAL(8) M_SFEED_LIQ(NCOMP,NPLATE)
REAL(8) M_FEED(NCOMP,NPLATE),M_TOP(NCOMP)
INTEGER MATERIAL(10)
CHARACTER(18) CHARA(10)
IW = 6
c WRITE(IW,*) ""
c WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
c WRITE(IW,'(15x,a35)') "END CALCULATION"
c WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
c WRITE(IW,*) ""
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,'(15x,a12,i3,a12)') "RESULT TOWER",IFLAG, " CALCULATION"
WRITE(IW,*) ""
WRITE(IW,'(10x,a40)') "CONCENTRATION of COMPONENT (mol/hr)"
WRITE(IW,*) ""
WRITE(IW,'(15x,6a15)') "FEED","SFEED","TOP","BTM","SCUT","Err"
BTM = FEED_STREAM - TOP - SCUT + SFEED
DO I = 1,NCOMP
WRITE(IW,'(13x,a5,6f15.5)')
1 CHARA(MATERIAL(I)),FEED_STREAM*xLIQ_FEED(I)
1 ,SFEED*xLIQ_SFEED(I),TOP*xLIQ(I,1),BTM*xLIQ(I,NPLATE)
1 ,SCUT*xLIQ(I,NPLATE_SCUT)
1 ,FEED_STREAM*xLIQ_FEED(I)
1 +SFEED*xLIQ_SFEED(I)-TOP*xLIQ(I,1)-BTM*xLIQ(I,NPLATE)
1 -SCUT*xLIQ(I,NPLATE_SCUT)
ENDDO
working1 = 0.0d0
working2 = 0.0d0
working3 = 0.0d0
working4 = 0.0d0
working5 = 0.0d0
DO I = 1,NCOMP
working1 = working1 + FEED_STREAM*xLIQ_FEED(I)
working2 = working2 + SFEED*xLIQ_SFEED(I)
working3 = working3 + TOP*xLIQ(I,1)
working4 = working4 + BTM*xLIQ(I,NPLATE)
working5 = working5 + SCUT*xLIQ(I,NPLATE_SCUT)
ENDDO
WRITE(IW,*) ""
WRITE(IW,'(13x,a5,5f15.5)') "TOTAL",
1 working1,working2,working3,working4,working5
WRITE(IW,*) ""
WRITE(IW,*) ""
WRITE(IW,'(14x,a20)') "CONDITION IN COLUMN"
WRITE(IW,*) ""
WRITE(IW,'(10x,a10,a20,25x,3a20)')
1 "PLATE","xLIQ","TEMP","VAP","LIQ"
WRITE(IW,'(30x,5(a9,x))')
1 (CHARA(MATERIAL(J)),J=1,NCOMP)
WRITE(IW,'(106x,f20.5)') TOP
c
DO I = 1,NPLATE
IF(I.EQ.1) THEN
WRITE(IW,'(10x,i10,a,5x,3f10.5,10x,f20.5,a20,f20.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),"--------",LIQ(I)
ELSEIF(I.EQ.NPLATE_FEED) THEN
WRITE(IW,'(10x,a7,i3,a,5x,3f10.5,10x,3f20.5)')
1 "FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSEIF(I.EQ.NPLATE_SCUT) THEN
WRITE(IW,'(9x,a8,i3,a,5x,3f10.5,10x,3f20.5)')
1 "S_CUT<<<",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSEIF(I.EQ.NPLATE_SFEED) THEN
WRITE(IW,'(8x,a9,i3,a,5x,3f10.5,10x,3f20.5)')
1 "S_FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSE
WRITE(IW,'(10x,i10,a,5x,3f10.5,10x,3f20.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ENDIF
ENDDO
ENDSUBROUTINE
SUBROUTINE MATRIX_OUT(A,M)
IMPLICIT NONE
INTEGER I,J,II,JJ,M,MN
INTEGER NDim,NDimI,NDimJ
DOUBLE PRECISION A(M,M)
CHARACTER*10 LINE,inputname
CHARACTER*10 orbital(50)
C
LINE = '-------'
C
IF(M.LE.20) THEN
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,M)
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,M)
DO I = 1, M
WRITE(*,'(i3,x,a2,20(f12.5,x))') I,"|",(A(I,J),J=1,M)
ENDDO
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,M)
WRITE(*,*) ""
WRITE(*,*) ""
ELSE
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,20)
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,20)
DO I = 1, M
WRITE(*,'(i3,x,a2,20(f6.2,x))')
1 I,"|",(A(I,J),J=1,20)
ENDDO
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=1,20)
WRITE(*,*) ""
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=21,M)
DO I = 1,M
WRITE(*,'(i3,x,a2,20(f6.2,x))')
1 I,"|",(A(I,J),J=21,M)
ENDDO
WRITE(*,'(22a7)') LINE,LINE,(LINE,J=21,M)
WRITE(*,*) ""
WRITE(*,*) ""
ENDIF
END
まとめ
このページでは、ASU（空気分離装置）を数値シミュレーションし、蒸留塔内部の温度や濃度、流量プロファイルがどのようなものかを調べてみました。
プラントで使われている窒素などはまず間違いなく、ASU装置を経て出てきています。たとえボンベから持ってきていても、ボンベに入る前はASUの中にいたはずです。
ASUを持っている工場は少数かと思いますが、もしASUの運転を体験したい場合は、このホームページのプログラムを動かしてみてください。（lapack環境が必要ですが。。）