多重効用蒸発 −連立非線形方程式−

 化学工学の特徴が最もよくあらわれているのが単位操作における図式解法である。蒸留,抽出,吸収等の単位操作の設計計算は,物質収支と平衡関係を組み合わせた多変数非線形連立方程式を解く問題に帰着される。化学工学の実務では偏微分,複素数等の高度な数学は不要であるが,未知数の多い連立非線形方程式を解く必要が多い。幾何と代数の同一性は受験勉強の数学でおなじみのことであるが,単位操作の図式解法は計算機の無い時代に非線形連立方程式をどうしても解く必要に迫られたケミカルエンジニアの先輩達の偉大な考案であった。図式解法はプロセス設計の実際をわかりやすく示すが,計算機の時代にふさわしく,数値計算の点から図式解法を見直す必要もある。
 ここでは多重効用蒸発プロセスを取り上げ,物質収支と熱収支からなる連立方程式を解く演習をおこなう。

(パルプ製造における黒液の濃縮用の多重効用蒸発装置(北越製紙))

多重効用蒸発の熱および物質収支

 3重効用蒸発缶(順流)を用いてF= 303K, F=0.05重量分率,のある水溶液 (= 2.0 + 0.5c )[kg/s]を xW( = 0.2 + 0.01)[重量分率]まで濃縮する。

簡単のため各缶の伝熱面積および総括伝熱係数は等しく各々 [m2],=2300W/(m2・K) であり, 蒸発潜熱,熱容量も温度依存性を無視して各々 =2.3×106J/kg, p=4.2×103J /(kg・K) とする。B3=325K,s= 403Kとして伝熱面積 [m2],温度B1, TB2,各缶の蒸発量1[kg/s],2, 3,供給蒸気量S を求めよ。

【解】<全物質収支> - = 1 + 2 + 3
<溶質収支>    F= W
上式より, 1 + 2 + 3 = ( 1 -F/W )     (1)
<蒸発潜熱に関する熱収支>
第1缶:  rVS = rV1 + Fcp(B1-F)        (2)
第2缶:  rV1 = rV2 + ( -1)p(B2-B1)    (3)
第3缶:  rV2 = rV3 + (-1-2)p(B3-B2)  (4)
<各缶の伝熱速度>
第1缶:  UA(s -B1) = rVS (5)
第2缶:  UA(B1-B2) = rV1 (6)
第3缶:  UA(B2-B3) = rV2 (7)

以上の7つの非線形連立方程式を解く。

 未知数を次のように置き換える。

B1 B2 1 2 3 S
X(1) X(2) X(3) X(4) X(5) X(6) X(7)

実習レポートの様式 (A4方眼紙使用)

 提出するレポートには3重効用蒸発缶の略図を描いてその図中に各流量,濃度,温度を単位つきでしめしなさい

連立非線形方程式の数値解法

?

 n変数の連立非線形方程式 f (x)=0  すなわち 

を解くことは,f (x)→0となるx =(x1,x2,・・・・xn)T を求めることである。
 一変数の非線形方程式 f(x)=0 の解法であるNewton-Raphson法:

n 元連立非線形方程式にたいし拡張できる。.この疑似Newton法による非線形連立方程式解法の計算手順は以下のようである。
1 (プログラム350-370) 適当な初期ベクトル x(0)を与える。k=0,1,2・・・・に対し以下を反復計算する。
2 (プログラム380-390) f (x (k ))の値を計算する。各ベクトルの値は0でない。
3 (プログラム30-440,サブルーチン3000) ∂fi /∂xj を前進差分で近似して求める。

これをJacobi(ヤコビ)行列と呼ぶ。
4 (プログラム450-550,サブルーチン1000,2000) 次の連立1次方程式をGaussの消去法で解き,解Δx (k )を求める。
J (x (k ))Δx (k ) = -f (x (k ))
5 (プログラム590-650) x(k+1) = x(k) +Δx(k) とおく。
6 (プログラム660-670) Δx(k) が収束判定値εより小さいなら x(k+1)を方程式f (x )=0 の解として計算を終了する。そうでないなら x(k+1)を新しい近似値として2に戻る。

連立非線形方程式解法プログラム

100 '==連立非線形方程式解法プログラム=====(プログラムの大きさ3.2K)
110 ' Pseudo Newton's Method
120 ' by スノウチ ハルオ " パソコンニヨルスウチケイサン ",p.144,サイエンスシャ
210 READ NN
220 DIM ZA(NN, NN), ZB(NN), ZX(NN), ZY(NN), NG(NN)
230 DIM F(NN), ZF1(NN), ZX1(NN), X(NN)
240 '-----------------
280 FOR NI = 1 TO NN: READ ZX1(NI): NEXT
300 NL1 = 100
310 ZE = .00001 (これは1変数非線型方程式の解法のNewton-Raphson
320 ZD = .001         法をn変数に拡張したものである
330 '----------------     3080行でJacobi行列を前進差分近似で求めている
340 FOR NL = 1 TO NL1     したがって,初期値がある程度近くないと収束しな
350 FOR NI = 1 TO NN     いことに注意。またNewton法同様expに弱い。)
360 X(NI) = ZX1(NI)   近似解を代入
370 NEXT NI

380 '----------------
390 GOSUB 4000
400 FOR NI = 1 TO NN
410 ZF1(NI) = F(NI)
420 NEXT NI
430 '----------------
440 GOSUB 3000
450 '----------------連立一次方程式をGaussの消去法で解く
460 '----------------
470 FOR NI = 1 TO NN
480 NG(NI) = NI
490 NEXT NI
500 GOSUB 1000
510 '----------------
520 FOR NI = 1 TO NN
530 NJ = NG(NI): ZB(NI) = -ZF1(NJ)
540 NEXT NI
550 GOSUB 2000
560 '----------------
570 '----------------
580 PRINT "ハンプク カイスウ "; NL
590 NC = 0
600 FOR NI = 1 TO NN
610 ZW = X(NI)
620 ZX1(NI) = ZX1(NI) + ZW
630 IF ABS(ZW) < ZE OR ABS(ZW) < ZE * ABS(ZX1(NI)) THEN NC = NC + 1
640 PRINT "X("; NI; ")="; ZX1(NI)
650 NEXT NI
660 '----------------
670 IF NC = NN THEN 730
680 NEXT NL
690 '----------------
700 BEEP: PRINT "ハンプク="; NL1; "デ シュウソクシナイ"
710 GOTO 770
720 '----------------
730 PRINT "*****ホウテイシキ ノ カイ*****"
740 FOR NI = 1 TO NN
750 PRINT " X("; NI; ")= "; ZX1(NI)
760 NEXT NI
770 END
1000 '---------------
1010 '---------------
1020 FOR NK = 1 TO NN - 1
1030 '-------------
1040 ZW = ABS(ZA(NK, NK)): NK1 = NK
1050 FOR NI = NK + 1 TO NN
1060 ZW1 = ABS(ZA(NI, NK))
1070 IF ZW1 < ZW THEN 1090
1080 ZW = ZW1: NK1 = NI
1090 NEXT NI
1100 '-------------
1110 IF NK = NK1 THEN 1160
1120 FOR NJ = 1 TO NN
1130 ZW = ZA(NK, NJ): ZA(NK, NJ) = ZA(NK1, NJ): ZA(NK1, NJ) = ZW
1140 NEXT NJ
1150 NI1 = NG(NK): NG(NK) = NG(NK1): NG(NK1) = NI1
1160 '-------------
1170 ZW = ZA(NK, NK)
1180 FOR NI = NK + 1 TO NN
1190 ZA(NI, NK) = ZA(NI, NK) / ZW
1200 FOR NJ = NK + 1 TO NN
1210 ZA(NI, NJ) = ZA(NI, NJ) - ZA(NI, NK) * ZA(NK, NJ)
1220 NEXT NJ
1230 NEXT NI
1240 NEXT NK
1250 RETURN
2000 '-------SOLVE
2010 '---------------
2020 ZY(1) = ZB(1)
2030 FOR NI = 2 TO NN
2040 ZW = ZB(NI)
2050 FOR NK = 1 TO NI - 1
2060 ZW = ZW - ZA(NI, NK) * ZY(NK)
2070 NEXT NK
2080 ZY(NI) = ZW
2090 NEXT NI
2100 '----------------
2110 X(NN) = ZY(NN) / ZA(NN, NN)
2120 FOR NI = NN - 1 TO 1 STEP -1
2130 ZW = ZY(NI)
2140 FOR NK = NI + 1 TO NN
2150 ZW = ZW - ZA(NI, NK) * X(NK)
2160 NEXT NK
2170 X(NI) = ZW / ZA(NI, NI)
2180 NEXT NI
2190 RETURN
3000 '----------------Jacobi行列の計算
3010 '----------------
3020 FOR NJ = 1 TO NN
3030 IF NJ <> 1 THEN X(NJ - 1) = ZX1(NJ - 1)
3040 X(NJ) = X(NJ) + ZD
3050 '--------------
3060 GOSUB 4000
3070 FOR NI = 1 TO NN
3080 ZA(NI, NJ) = (F(NI) - ZF1(NI)) / ZD
3090 NEXT NI
3100 NEXT NJ
3110 RETURN

4000 '**** 方程式 **** 関数値ベクトルf(x)の計算
4005 ' 例題:3重効用蒸発缶
4010 '(Z,N ノツクヘンスウハツカワナイ )
4020 '--ミチスウ ノ カズ X()
4030 DATA 7
4040 '初期値A, TB1,TB2, V1, V2, V3, Vs
4051 DATA ??, ???, ???, ???, ???, ???, ??? (初期値の選定には「常識」
4060 '--レンリツ ホウテイシキ                を はたらかせて)
4070 U = 2300: R = 2300000!: TF = 303: CP = 4200!: TS = 403: TB3 = 325
4080 FF = ???: XXF =0.05: XXW =0.?? (FF:Fの値, XXW:xwの値)
4100 F(1) = X(4) + X(5) + X(6) - (1 - XXF / XXW) * FF
4110 F(2) = R * X(7) - R * X(4) - FF * CP * (X(2) - TF)
4120 F(3) = R * X(4) - R * X(5) - (FF - X(4)) * CP * (X(3) - X(2))
4130 F(4) = R * X(5) - R * X(6) - (FF - X(4) - X(5)) * CP * (TB3 - X(3))
4140 F(5) = R * X(7) - U * X(1) * (TS - X(2))
4150 F(6) = R * X(4) - U * X(1) * (X(2) - X(3))
4170 F(7) = R * X(5) - U * X(1) * (X(3) - TB3)
4220 RETURN

?


▲目次