らい・ぶらり

精留塔設計：McCabe-Thiele法をExcelマクロで再現し塔内濃度分布を算出

McCabe-Thieleによる階段作図をExcelで再現するマクロの作成

本ページは下記資料の内容を理解していることを前提に話を進めています。精留塔とは気液平衡とは何か？については下記資料をご覧ください。
また、本ページは「McCabe-Thiele法による階段作図」に関する技術を講義で受け、単位を取得した人向けに書きました。

[blogcard url=”http://shimaphoto03.com/science/distillation/”] [blogcard url=”http://shimaphoto03.com/science/mccabe/”]

はしがき

工学部化学系の学生は「分離工学」あるいは「化学工学」などといった講義の名前で「精留塔の理論段数を求める方法」として「McCabe-Thiele法による階段作図」を学び、単位をとるわけですが、実際の化学工場の設計となると、この階段作図はほとんど行われません。

理由は２点あります。

①コンピュータ技術の発展により、多元連立方程式を解くことが容易になったので、作図による設計は非効率になった。

②理論段数を求めるという行為は「新規精留塔の設計」に当たるのだが、近年は精留塔を新しく建設する案件が減ってきた。

①についてですが、
精留塔の理論段数を求める問題は「気液平衡関係と物質収支に関する多元連立方程式を解く」問題に帰着します。しかし、多元連立方程式を紙とペンだけで解くのは非常に大変なので「McCabe-Thiele法による階段作図」という、図による解法が1925年に開発され、それがそのまま21世紀の講義の内容に受け継がれています。

その結果として現在、大学の定期試験で「McCabe-Thiele法による階段作図により理論段数を求めよ」という問題が解けると単位が手に入ってしまいます。
しかもAAが手に入ってしまって、精留塔を「知った気」でいると、技術者として働き始めてから苦労してしまいます。

なぜなら、現実の精留塔は「気液平衡関係と物質収支」だけでモデル化できるわけがなく、「熱収支」「気液間物質移動速度論（段効率）」「第三成分の影響」「気液平衡の非理想性（ラウールモデルからの乖離）」などが関わってきており、こういった現象の理解と定式化、解法の技術が必要なのです。

現代は何十元もの多元連立方程式を数秒で解いてしまうコンピュータが一家に一台レベルになってしまっているので、「現象の理解と定式化と解法の技術」がそのまま競争力に直結する時代なのでは、と思います。

②についてですが、
「McCabe-Thiele法による階段作図」が開発された1925年ごろは、化学工業産業の黎明期であり、新規精留塔建設の案件が多かったと思いますが、現代は、新しく精留塔を建設するというよりは「今ある精留塔を使って新しいものを作りたい。できるか？できないか？できないならどうすればできるようになるか？」という課題の方が多いかと思います。

そういった背景からも、「McCabe-Thiele法による階段作図により理論段数を求めよ」という問題が解けると単位が手に入ってしまう現在の大学教育はレベルが低く、その低いレベルで学位をとってしまった「新卒」は「レベルが低い」と指差されてしまっても致し方ないのでは、と思います。

McCabe-Thieleの階段作図のマクロ化

はしがきが長くなってしまいましたが、ここからが本題です。

現代のコンピュータ技術の発展により、何十元もの多元連立方程式を解くことが容易になった。とお話ししました。
本ページはそのコンピュータ技術の発展の恩恵を、直接設計につなげるための「１ステップ目」として、McCabe-Thieleの階段作図をコンピュータ上でできるようにするというのが趣旨となります。

今回、マクロで解いている方程式の導出（＝気液平衡関係と物質収支の導き方）はこちらにて解説しております。
[blogcard url=”http://shimaphoto03.com/science/mccabe/”]

マクロを使用するExcelシートスタイル

マクロを適用するExcelシートのスタイルは以下の通りです。罫線を引いてある部分に数値を入れます。それ以外の部分は自動で計算されます。

プログラム

上記のスタイルで下記のマクロを実行するとMcCabe-Thieleの階段作図と同じことができるようになります。

※予めグラフの範囲の指定をしておいてください。

Sub Distillation()

F = Cells(3, 2)
D = Cells(4, 2)
W = F – D
Cells(5, 2).FormulaR1C1 = “=R3C2-R4C2”
q = Cells(7, 2)
R = Cells(8, 2)
zF = Cells(9, 2)
L = R * D
Cells(6, 2).FormulaR1C1 = “=R8C2*R4C2”
alpha = Cells(3, 8)

‘nF:原料供給段
nF = Cells(10, 2)

‘n:理論段数
n = Cells(11, 2)

‘nD:留出部段数
nD = nF – 1
Cells(12, 2) = nD

‘nw:回収部段数
nW = n – nF + 1
Cells(13, 2) = nW

‘nrb=リボイラ段
nrb = n + 1
Cells(14, 2) = nrb

Cells(16, 2) = 1
Cells(16, 3) = 1

For i = 1 To nrb
Cells(16 + i, 1) = i
Cells(16 + i, 2) = 0.5
xi = Cells(16 + i, 2)

‘yi = (alpha * x) / (1 + (alpha – 1) * x)
Cells(16 + i, 3).FormulaR1C1 = “=(R3C8 * RC[-1]) / (1 + (R3C8 – 1) * RC[-1])”

Next i

xD = Cells(16 + 1, 2)
xW = Cells(16 + nrb, 2)

Cells(16 + nrb + 1, 2) = 0
Cells(16 + nrb + 1, 3) = 0

‘F1 = D * xD + W * xW – F * zF
Cells(16 + 1, 4).FormulaR1C1 = “=R4C2*R17C2+R5C2*R” & 16 + nrb & “C2-R3C2*R9C2”

For i = 2 To nD
‘F2 = (L * x) / (L + D) + (D * xD) / (L + D) – y
Cells(16 + i, 4).FormulaR1C1 = “=(R6C2*RC[-2])/(R6C2+R4C2)+(R4C2*R17C2)/(R6C2+R4C2)-R[1]C[-1]”
Next i

For i = nF To n
x = Cells(16 + i, 2)
y = Cells(16 + i + 1, 3)
‘F3 = (L + q * F) * x / (L + q * F – W) – (W * xW) / (L + q * F – W) – y
Cells(16 + i, 4).FormulaR1C1 = “=(R6C2+R7C2*R3C2)*RC[-2]/(R6C2+R7C2*R3C2-R5C2)-(R5C2*R” & 16 + nrb & “C2)/(R6C2+R7C2*R3C2-R5C2)-R[1]C[-1]”
Next i

‘F4 = xD-y2
Cells(16 + nrb, 4).FormulaR1C1 = “=R17C2-R” & 16 + 2 & “C3”

Cells(16 + nrb + 1, 4).FormulaR1C1 = “=SUMSQ(R17C4:R” & 16 + nrb & “C4)”

Cells(16 + 1, 2).Select
Selection.Font.Bold = True
Cells(16 + nrb, 2).Select
Selection.Font.Bold = True

SolverReset
SolverOk SetCell:=Cells(17 + nrb, 4).Address, MaxMinVal:=2, ValueOf:=0, ByChange:=Range(Cells(16 + 1, 2), Cells(16 + nrb, 2)).Address, _
Engine:=1, EngineDesc:=”GRG Nonlinear”
SolverOk SetCell:=Cells(17 + nrb, 4).Address, MaxMinVal:=2, ValueOf:=0, ByChange:=Range(Cells(16 + 1, 2), Cells(16 + nrb, 2)).Address, _
Engine:=1, EngineDesc:=”GRG Nonlinear”
SolverSolve

‘operating line
Cells(3, 5).FormulaR1C1 = “=R6C2/(R6C2+R4C2)”
aD = Cells(3, 5)
Cells(4, 5).FormulaR1C1 = “=(R4C2*R17C2)/(R6C2+R4C2)”
bD = Cells(4, 5)
Cells(5, 5).FormulaR1C1 = “=(R6C2+R7C2*R3C2)/(R6C2+R7C2*R3C2-R5C2)”
aW = Cells(5, 5)
Cells(6, 5).FormulaR1C1 = “=-(R5C2*R” & 16 + nrb & “C2)/(R6C2+R7C2*R3C2-R5C2)”
bW = Cells(6, 5)
Cells(7, 5).FormulaR1C1 = “=-R7C2/(1-R7C2)”
aq = Cells(7, 5)
Cells(8, 5).FormulaR1C1 = “=R9C2/(1-R7C2)”
bq = Cells(8, 5)

For i = 0 To 1000
Cells(i + 3, 13) = i / 1000

If i / 1000 < Cells(16 + 1, 2) Then Cells(i + 3, 14).FormulaR1C1 = “=R6C2*RC[-1]/(R6C2+R4C2)+(R4C2*R17C2)/(R6C2+R4C2)” End If Cells(i + 3, 16) = i / 1000 If i / 1000 > Cells(16 + nrb, 2) Then
Cells(i + 3, 17).FormulaR1C1 = “=(R6C2+R7C2*R3C2)*RC[-1]/(R6C2+R7C2*R3C2-R5C2)-(R5C2*R” & 16 + nrb & “C2)/(R6C2+R7C2*R3C2-R5C2)”
End If

Cells(i + 3, 19) = i / 1000

Cells(i + 3, 20).FormulaR1C1 = “=-R7C2*RC[-1]/(1-R7C2)+R9C2/(1-R7C2)”

Next i

End Sub

実行後の結果

例題

今回作成したマクロが技術者の仕事においてどんな風に役に立ちそうか、経験は浅いですが、浅いなりにも考えてみました。

『最近研究所で開発された新製品「銘柄Y」の工業化プロジェクトを進めている。
このプロジェクトでは、4年前に撤退してしまった「銘柄X」を精製するために使用していた段数8段の精留塔を使用することを検討している。

「銘柄Y」の反応工程より送られてくる製品は不純物よりも低沸点である。流量は1kmol/h、純度は60mol%であり、不純物との相対揮発度は2.5であり、気液平衡関係はラウールモデルに従うこととしても差し支えないことが、昨年までの本プロジェクトに関する報告書に掲載されていた。なお、精留塔供給前に熱交換器にて沸点まで予熱させる予定である。
また、供給量の需要から試算すると、留出流量は0.5kmol/h必要であり、純度は95mol%まで高めたいようである。

この精留塔には、原料供給のための配管設備が3段目、5段目、7段目に用意されている。どの段に接続すればよいか検討せよ。
また、他の段に新しく配管を用意する予算が用意されているようであるが、その工事が本当に必要であるかどうか検討せよ。』

解答例

F = 1kmol/h, z_F = 0.6,　q = 1（沸点の液供給）、D = 0.5 kmol/h(留出流量)、α＝2.5の条件において、x_D = 0.95以上を求めている。という問題に書き下すことができる。

還流比が与えられていないので、還流比を比較的低めの「2」と仮定してシミュレーションを行う。
原料供給段を3段目～8段目に変化させ、シミュレーションを行った結果を示す。(gif動画です。)

3段目供給

5段目供給

7段目供給

以上、シミュレーション結果より、7段目に原料供給段を設置すればよいことがわかる。各段で気液平衡が成立していれば、純度は97.1mol%まで高めることができる。さらに要求されるスペック95mol%より多少の余裕がある。
段効率など精留塔個々の事情はあるが、少なくとも配管新設の工事は必要ないと言える。