McCabe–Thiele法
McCabe–Thiele法により精留塔の理論段数を求める操作は，作図が少しずれただけでとんでもないことになり，面倒です。 そこで，「面倒なことはPythonにやらせよう」ということでPythonに頑張ってもらいました。
機能
基本的にはMcCabe–Thiele法により精留塔の段数を求めるだけですが，色々余計な機能が付いています。
- 多くの化合物を指定できる (名前を文字列で与えればよい)
- 全圧と原料の状態(q値)を変更できる
- 各段での気液組成と温度分布を計算する
- 低沸点成分と高沸点成分を逆に指定した場合は，それを検出して適切な組み合わせに修正する
- 還流比と最小還流比の大小を確認する
動作環境
Python 3が必要です。2では動作しません。
ライブラリ
matplotlibとnumpyは必須です。SciPyもあった方がいいですが，なくても動作します。
iOS上のPython環境「Pythonista 3」でも動作するようになっています。
ソースコード
############################## Parameters ##############################
# The name of the lower-boiling point component #
low_bp = 'methanol'
# The name of the higher-boiling point component #
high_bp = 'water'
########################################################################
# Total pressure, measured in Pa #
P = 101325
assert P > 0, 'Total pressure must be greater than 0 (Pa).' #
########################################################################
# The reflux ratio #
R = 1.5
assert R > 0, 'The reflux ratio must be greater than 0.' #
########################################################################
# The parameter q #
# Superheated vapor feed : q < 0 #
# Saturated vapor feed : q = 0 #
# Partially vaporized feed: 0 < q < 1 #
# Saturated liquid feed : q = 1 #
# Subcooled liquid feed : q > 1 #
q = 1
########################################################################
# The distillate composition #
# (measured in mole fraction of the lower-boiling feed component
# in the liquid phase; the same shall apply hereinafter) #
x_D = .95
assert 0 < x_D < 1, \
'The distillate composition must be greater than 0 and less than 1.' #
########################################################################
# The bottoms composition #
x_W = .05
assert 0 < x_W < 1, \
'The bottoms composition must be greater than 0 and less than 1.' #
########################################################################
# The feed composition #
z_f = .40
assert 0 < z_f < 1, \
'The feed composition must be greater than 0 and less than 1.' #
########################################################################
y_D = x_D
y_W = x_W
antoine_constants = {
# Antoine equation:
# log(p) = A - B / (T - C)
# where p is the vapor pressure in kilopascal
# and T is temperature in kelvin
# Data:
# Substance: (A, B, C)
# (ref: http://memoirs.lib-e.yamaguchi-u.ac.jp/621/01.pdf)
## Alkane
'2-methylbutane' : (5.93330, 1029.602, 38.856),
'pentane' : (5.99028, 1071.187, 40.384),
'2-methylpentane' : (5.99479, 1152.210, 44.579),
'3-methylpentane' : (5.99139, 1162.069, 44.870),
'hexane' : (6.01098, 1176.102, 48.251),
'heptane' : (6.02701, 1267.592, 56.354),
'2,3-dimethylpentane' : (5.98293, 1240.404, 51.056),
'octane' : (6.04394, 1351.938, 64.030),
'2,2,4-trimethylpentane' : (5.92751, 1252.340, 53.060),
## cyclic compound
'cyclohexane' : (6.00569, 1223.273, 48.061),
'benzene' : (6.01905, 1204.637, 53.081),
'toluene' : (6.08436, 1347.620, 53.363),
## ether
'diethyl ether' : (6.04920, 1061.391, 45.090),
'methul t-butyl ether' : (6.070343, 1158.912, 43.200),
'ethyl t-butyl ether' : (6.073724, 1206.874, 49.190),
't-amyl methyl ether' : (6.067822, 1256.258, 50.100),
'diisopropyl ether' : (5.97081, 1137.408, 54.634),
'dibutyl ether' : (5.92274, 1298.256, 82.006),
## ketone
'acetone' : (6.25017, 1214.208, 43.148),
'methyl ethyl ketone' : (6.18397, 1258.940, 51.425),
'diethyl ketone' : (6.14570, 1307.941, 59.182),
'methyl propyl ketone' : (6.13931, 1309.629, 58.585),
'methyl isopropyl ketone': (6.09024, 1265.595, 57.631),
'methyl isobutyl ketone' : (5.81291, 1176.833, 80.225),
## alcohol
'methanol' : (7.24693, 1605.615, 31.317),
'ethanol' : (7.24222, 1595.811, 46.702),
'1-propanol' : (6.87065, 1438.587, 74.598),
'2-propanol' : (6.86634, 1360.183, 75.557),
'1-butanol' : (6.54068, 1335.028, 96.496),
'2-butanol' : (6.35079, 1169.924, 103.413),
'water' : (7.06252, 1650.270, 46.804),
}
import matplotlib.pyplot as plt
import numpy as np
import sys
try:
from scipy import optimize
def inv(func, y, x0=0, eps=None, **kwargs):
return optimize.fsolve(lambda x: func(x) - y, x0)[0]
except ImportError:
def inv(func, y, x0=0, x1=None, eps=2*sys.float_info.epsilon):
assert not x0 == x1, 'x1 cannot be equal to x0.'
f = lambda x: func(x) - y
if not x1:
x1 = x0 + 1
elif x0 > x1:
x0, x1 = x1, x0
# scant method
while abs(f(x1)) > eps:
s = (x0 * f(x1) - x1 * f(x0)) / (func(x1) - func(x0))
x0, x1 = x1, s
if abs(x0 - x1) < eps:
break
return x1
def Antoine(func):
def wrapped(*args, **kwargs):
res = func(*args, **kwargs)
return (10 ** res) * 1000
return wrapped
@Antoine
def vp(consts, T):
return consts[0] - consts[1] / (T - consts[2])
def vp_1(T):
return vp(antoine_constants[low_bp], T)
def vp_2(T):
return vp(antoine_constants[high_bp], T)
def temperature(x):
def total_P(T):
# Raoult's law
p_2 = vp_2(T) * (1 - x)
p_1 = vp_1(T) * x
# Dalton's law
return p_2 + p_1
return inv(total_P, y=P, x0=373.15, eps=1e-5)
@np.vectorize
def vle(x):
T = temperature(x)
return vp_1(T) * x / P
def draw(z_f, x_W, x_D, R, q, P):
plt.close()
t = np.linspace(0, 1, 200)
func_r = lambda x: R / (R + 1) * (x - x_D) + y_D # rectifying section
if q == 1:
Q_x = z_f
Q_y = func_r(z_f)
min_s = (y_D - vle(z_f)) / (x_D - z_f)
else:
func_q = lambda x: -q / (1 - q) * x + z_f / (1 - q)
func_diff = lambda x: func_r(x) - func_q(x)
f0 = func_diff(0)
Q_x = -f0 / (func_diff(1) - f0)
Q_y = func_q(Q_x)
p_x = inv(lambda x: vle(x) - func_q(x), 0, z_f-.1, x1=z_f+.1, eps=1e-5)
p_y = func_q(p_x)
min_s = (y_D - p_y) / (x_D - p_x)
func_s = lambda x: (y_W - Q_y) / (x_W - Q_x) * (x - x_W) + y_W # stripping section
func = lambda x: func_r(x) if x >= Q_x else func_s(x)
min_R = min_s / (1 - min_s)
plt.plot([0, 1], [0, 1], 'k') # diagonal line
plt.plot(t, vle(t), 'b') # VLE line
plt.plot([x_D, x_D], [0, y_D], 'k--') # distillate composition line
plt.plot([x_W, x_W], [0, y_W], 'k--') # bottoms composition line
plt.plot([z_f, z_f], [0, z_f], 'k--') # feed composition line
label_opt = {
'verticalalignment': 'top',
'horizontalalignment': 'center'
}
plt.text(x_D, -.0105, f'{x_D:.2f}', **label_opt)
plt.text(x_W, -.0105, f'{x_W:.2f}', **label_opt)
plt.plot([Q_x, x_D], [Q_y, y_D], 'g') # rectifying operating line
plt.plot([x_W, Q_x], [y_W, Q_y], 'g') # stripping section operating line
plt.plot([z_f, Q_x], [z_f, Q_y], 'b') # q-line
if R > min_R:
y1 = y_D
x1 = x_D
cnt = 1
while (x1 - x_W) > 0.01:
x2 = inv(vle, y=y1, x0=0, x1=x1, eps=1e-5)
y2 = max(func(x2), x2)
T = round(temperature(x2), 1)
if y1 >= Q_y:
if q == 0:
if y2 <= Q_y <= y1:
feed = cnt, T
elif q == 1:
if x2 <= Q_x <= x1:
feed = cnt, T
else:
F_x = inv(func_q, y=y1, x0=x2, x1=x1, eps=1e-5)
if (y1 < vle(F_x)) or (y2 <= func_q(x2) <= y1):
feed = cnt, T
plt.text((x1 - x2) / 2 + x2, y1 + .005,
f'{T:.1f} K',
horizontalalignment='center'
)
yield f'{cnt}: x = {round(x2, 2):.2f}, y = {round(y1, 2):.2f} at {T:.1f} K'
plt.plot([x1, x2], [y1, y1], 'c')
plt.plot([x2, x2], [y1, y2], 'c')
x1, y1 = x2, y2
cnt += 1
yield f'Feed tray: Step #{feed[0]} ({feed[1]:.1f} K)'
else:
yield f'R <= {round(min_R, 2):.2f}'
plt.xlim(xmin=0.0, xmax=1.0)
plt.ylim(ymin=0.0, ymax=1.05)
plt.title(f'{low_bp}-{high_bp} (q = {q}, R = {R}, P = {P:.1f} Pa)')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
def validate_compound():
global low_bp, high_bp
if low_bp == high_bp:
print('Same compound')
return False
for compound in [low_bp, high_bp]:
if compound not in antoine_constants:
print(f'No data: {compound}')
return False
bp_opt = [P, 273.15, 373.15, 1e-5]
bp1 = inv(vp_1, *bp_opt)
bp2 = inv(vp_2, *bp_opt)
if bp1 > bp2:
print(f'b.p.: {low_bp} ({bp1:.1f} K) > {high_bp} ({bp2:.1f} K)')
print(' They will be swapped.')
low_bp, high_bp = high_bp, low_bp
return True
def main():
if not validate_compound():
return
for R in [1.5, 3.0, 0.5]:
print(f'R = {R}')
for data in draw(z_f, x_W, x_D, R, q, P):
print(data)
if __name__ == '__main__':
main()
実行結果
R=1.5
R = 1.5
1: x = 0.82, y = 0.95 at 341.3 K
2: x = 0.64, y = 0.87 at 345.9 K
3: x = 0.45, y = 0.76 at 351.4 K
4: x = 0.33, y = 0.65 at 355.8 K
5: x = 0.22, y = 0.51 at 360.6 K
6: x = 0.12, y = 0.33 at 365.8 K
7: x = 0.05, y = 0.16 at 369.7 K
Feed tray: Step #4 (355.8 K)
R=3.0
R = 3.0
1: x = 0.82, y = 0.95 at 341.3 K
2: x = 0.60, y = 0.86 at 346.9 K
3: x = 0.37, y = 0.69 at 354.5 K
4: x = 0.21, y = 0.49 at 361.2 K
5: x = 0.09, y = 0.27 at 367.2 K
6: x = 0.03, y = 0.11 at 370.9 K
Feed tray: Step #3 (354.5 K)
R=0.5
R = 0.5
R <= 0.73
解説
2～40行目
パラメータの設定です。コメントとassertで書いてある通りに設定します。
42～43行目
y_D = x_D
y_W = x_W
x_D，x_Wに対応する，対角線上のyの値を求めています。
45～89行目
Antoine式で使用するパラメータです。Antoine式については，コメントを参照してください。
パラメータは，
Correlation of Vapor-Liquid Equilibria Using Wilson Equation with Parameters Estimated fromSolubility Parameters and Molar Volumes
を参照しました。
91～93行目
import matplotlib.pyplot as plt
import numpy as np
import sys
必要なライブラリを読み込んでいます。
95～114行目
try:
from scipy import optimize
def inv(func, y, x0=0, eps=None, **kwargs):
return optimize.fsolve(lambda x: func(x) - y, x0)[0]
except ImportError:
def inv(func, y, x0=0, x1=None, eps=2*sys.float_info.epsilon):
assert not x0 == x1, 'x1 cannot be equal to x0.'
f = lambda x: func(x) - y
if not x1:
x1 = x0 + 1
elif x0 > x1:
x0, x1 = x1, x0
# scant method
while abs(f(x1)) > eps:
s = (x0 * f(x1) - x1 * f(x0)) / (func(x1) - func(x0))
x0, x1 = x1, s
if abs(x0 - x1) < eps:
break
return x1
逆関数の値を求める汎用的な関数invの定義です。
SciPyが利用できる環境では，scipy.optimize.fsolveを利用します。
SciPyが利用できない(i.e. ImportErrorが投げられる)環境では，逆関数の計算のロジックを自作します。
Newton法にしたかったのですが，微分の計算が重すぎて実用的な時間内に計算が終了しなかったので割線法で実装しました。
McCabe–Thiele法で出てくる関数はすべて単調増加なので(q=0，1の場合のq-線のみが例外)，ゼロ除算の問題はありません。
少しでもオーバーヘッドを減らすため，110行目の代入演算子の右オペランドの分母をfではなくfuncにしています。
117～121行目
def Antoine(func):
def wrapped(*args, **kwargs):
res = func(*args, **kwargs)
return (10 ** res) * 1000
return wrapped
Antoineデコレータの定義です。Antoine式では圧力(kPa; これはパラメータの問題)の対数が求まるので，これをPaに換算するデコレータです。
現在の実装では1回しか使われていませんが，かつて複数回使われていた時代に，同じロジックを複数回書くことを避けるために(DRY)デコレータを作成した次第です。
125～126行目
def vp(consts, T):
return consts[0] - consts[1] / (T - consts[2])
パラメータと温度から蒸気圧(の対数を取った値)を求める，Antoine式をそのまま書いたものです。
Antoineデコレータを付けることで蒸気圧をそのまま返す関数になります。
129～134行目
def vp_1(T):
return vp(antoine_constants[low_bp], T)
def vp_2(T):
return vp(antoine_constants[high_bp], T)
各成分の蒸気圧を求める関数です。実装はvpに丸投げしています。
137～145行目
def temperature(x):
def total_P(T):
# Raoult's law
p_2 = vp_2(T) * (1 - x)
p_1 = vp_1(T) * x
# Dalton's law
return p_2 + p_1
return inv(total_P, y=P, x0=373.15, eps=1e-5)
特定の液組成で気液平衡になるような温度を求める関数です。
内部のtotal_Pは，ある温度(と液組成)における全圧を求める関数です。
この部分は化学の知識が必要になるので少し詳しく説明します。
まず，ある温度が与えられたとき，Antoine式により純物質の蒸気圧を求めることができます。 さらに液組成が与えられればRaoultの法則より，蒸気中の各成分の分圧が求まります。 ここで，Daltonの法則から，混合気体の全圧は各成分の分圧の和に等しいため，蒸気の全圧を求めることができます。
temperatureでは，total_Pの逆関数を利用して，蒸気の全圧がPに一致する温度を求めています。
149～151行目
def vle(x):
T = temperature(x)
return vp_1(T) * x / P
与えられた液組成で気液平衡に達したときに気組成を返す関数です。
関数名はVapor–Liquid Equilibrium(気液平衡)に拠ります。
まず，特定の液組成で気液平衡になる温度Tを求めます。
その温度での低沸点成分の蒸気圧(Raoultの法則により，液組成をかけたもの)と全圧との比が求める気組成になります。
154～234行目
長いので分けます。
名前がdrawであるので，データは文字列として出力は行わず，ジェネレータとして呼び出し元に通知しています。
155～156行目
plt.close()
t = np.linspace(0, 1, 200)
どうでもいい下準備です。
matplotlib.pyplotをクリアし，液組成を数値計算用に200分割したリストtを用意します。
158～174行目
func_r = lambda x: R / (R + 1) * (x - x_D) + y_D # rectifying section
if q == 1:
Q_x = z_f
Q_y = func_r(z_f)
min_s = (y_D - vle(z_f)) / (x_D - z_f)
else:
func_q = lambda x: -q / (1 - q) * x + z_f / (1 - q)
func_diff = lambda x: func_r(x) - func_q(x)
f0 = func_diff(0)
Q_x = -f0 / (func_diff(1) - f0)
Q_y = func_q(Q_x)
p_x = inv(lambda x: vle(x) - func_q(x), 0, z_f-.1, x1=z_f+.1, eps=1e-5)
p_y = func_q(p_x)
min_s = (y_D - p_y) / (x_D - p_x)
func_s = lambda x: (y_W - Q_y) / (x_W - Q_x) * (x - x_W) + y_W # stripping section
func = lambda x: func_r(x) if x >= Q_x else func_s(x)
min_R = min_s / (1 - min_s)
まず，158行目で濃縮部操作線(rectifying section)を表す関数func_rを定義しています。
159行目から162行目では，q-値が1の場合のq-線と濃縮部操作線の交点の座標と，濃縮部操作線の最小の傾き(min_s)を求めています。
163～171行目では，q-値が1以外の場合のq-線と濃縮部操作線の交点の座標と，濃縮部操作線の最小の傾き(min_s)を求めています。
172～174行目では，回収部操作線の関数，2つの操作線を合わせた関数，最小還流比(min_R)を定義しています。
176～181行目
plt.plot([0, 1], [0, 1], 'k') # diagonal line
plt.plot(t, vle(t), 'b') # VLE line
plt.plot([x_D, x_D], [0, y_D], 'k--') # distillate composition line
plt.plot([x_W, x_W], [0, y_W], 'k--') # bottoms composition line
plt.plot([z_f, z_f], [0, z_f], 'k--') # feed composition line
この部分では，コメントにある通りの5本の線を描画しています。
183～188行目
label_opt = {
'verticalalignment': 'top',
'horizontalalignment': 'center'
}
plt.text(x_D, -.0105, f'{x_D:.2f}', **label_opt)
plt.text(x_W, -.0105, f'{x_W:.2f}', **label_opt)
xDとxWの濃度のラベルを書き込んでいます。-0.0105はPythonistaで調整した位置です。
190～192行目
plt.plot([Q_x, x_D], [Q_y, y_D], 'g') # rectifying operating line
plt.plot([x_W, Q_x], [y_W, Q_y], 'g') # stripping section operating line
plt.plot([z_f, Q_x], [z_f, Q_y], 'b') # q-line
濃縮部操作線，回収部操作線，q-線を描いています。
194～227行目
if R > min_R:
y1 = y_D
x1 = x_D
cnt = 1
while (x1 - x_W) > 0.01:
x2 = inv(vle, y=y1, x0=0, x1=x1, eps=1e-5)
y2 = max(func(x2), x2)
T = round(temperature(x2), 1)
if y1 >= Q_y:
if q == 0:
if y2 <= Q_y <= y1:
feed = cnt, T
elif q == 1:
if x2 <= Q_x <= x1:
feed = cnt, T
else:
F_x = inv(func_q, y=y1, x0=x2, x1=x1, eps=1e-5)
if (y1 < vle(F_x)) or (y2 <= func_q(x2) <= y1):
feed = cnt, T
plt.text((x1 - x2) / 2 + x2, y1 + .005,
f'{T:.1f} K',
horizontalalignment='center'
)
yield f'{cnt}: x = {round(x2, 2):.2f}, y = {round(y1, 2):.2f} at {T:.1f} K'
plt.plot([x1, x2], [y1, y1], 'c')
plt.plot([x2, x2], [y1, y2], 'c')
x1, y1 = x2, y2
cnt += 1
yield f'Feed tray: Step #{feed[0]} ({feed[1]:.1f} K)'
else:
yield f'R <= {round(min_R, 2):.2f}'
McCabe–Thiele法のメインの部分です。還流比が最小還流比より大きい場合のみ実行します。 還流比が最小還流比以下の場合には最小還流比を示したうえで還流比がそれよりも小さいことを通知します。
ステップの描画は，液組成と目標の差が0.01より大きい間繰り返されます。
199行目で気組成が気液平衡になっているような液組成を求め，200行目でその液組成での操作線上の気組成を求めます。
このとき，線が対角線の下まで伸びると見た目が悪いので，maxで対角線より下にならないようにしています。
201行目ではこのときの温度を求めています。
203行目から213行目では，原料を入れる段であるかの判定を行っています。 複数の段が候補となる場合でも，上から順に計算して言っているので最終的に変数に残されるのは最も下の候補の段となります。
215行目から219行目では，グラフに段の温度を書き込み，その段の情報を呼び出し元に通知しています。
221行目から224行目では，段の線を描きこんだ後，気液各組成を更新し，段のカウントをインクリメントします。
229～234行目
plt.xlim(xmin=0.0, xmax=1.0)
plt.ylim(ymin=0.0, ymax=1.05)
plt.title(f'{low_bp}-{high_bp} (q = {q}, R = {R}, P = {P:.1f} Pa)')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
グラフの範囲，タイトル，軸ラベルを設定してグラフを表示します。
237～255行目
def validate_compound():
global low_bp, high_bp
if low_bp == high_bp:
print('Same compound')
return False
for compound in [low_bp, high_bp]:
if compound not in antoine_constants:
print(f'No data: {compound}')
return False
bp_opt = [P, 273.15, 373.15, 1e-5]
bp1 = inv(vp_1, *bp_opt)
bp2 = inv(vp_2, *bp_opt)
if bp1 > bp2:
print(f'b.p.: {low_bp} ({bp1:.1f} K) > {high_bp} ({bp2:.1f} K)')
print(' They will be swapped.')
low_bp, high_bp = high_bp, low_bp
return True
各成分の正当性を検証します。戻り値は，計算を継続していいかを示すBoolean値です。
まず，低沸点成分と高沸点成分が同一である場合は問答無用で弾きます。
また，いずれかの成分のAntoine式のパラメータが登録されていな場合，継続不能と判断します。
さらに，247行目から253行目では，沸点の高低を検証します。 各成分の沸点を計算し，低沸点成分と高沸点成分が逆である場合には警告を表示したうえで成分を入れ替えます。
258～265行目
def main():
if not validate_compound():
return
for R in [1.5, 3.0, 0.5]:
print(f'R = {R}')
for data in draw(z_f, x_W, x_D, R, q, P):
print(data)
メインの関数です。成分の検証を行った後，還流比を変えながらMcCabe–Thiele法で作図します。
267～268行目
ボイラープレートコード的な部分です。
Tweet