蒸留塔のモデル化と設計
バイオ・リファイナリー（再生可能資源であるバイオマスを原料にバイオ燃料や樹脂などを製造するプラントや技術)のシミュレーションソフト"BioSTEAM"で、蒸留塔のモデル化と設計について説明しています。
オリジナルのページはDistillation modeling & designです。
ソースコードは以下の実行環境で確認しています。
- Visual Studio Code バージョン: 1.104.2
- 拡張機能:Jupyter バージョン 2025.8.0
- Python 3.12.10
- biosteam 2.52.13
- graphviz-14.0.2
蒸留塔のモデル化と設計
蒸留塔は化学工業で最も幅広く用いられている装置で、石油化学、医薬品、バイオマス利用工業で応用されています。混合物を分離したい場合に、多くの成分で蒸留は有効な選択肢となります。このため、新しい分離技術(例えば膜分離や、吸着)は、しばしば蒸留を基盤としたプロセスと比較されます。ここでは、炭化水素の分離や小さな有機分子の精製のために、蒸留塔を設計しシミュレーションしてみます。
簡易蒸留モデル
一般的に、蒸留は複雑でモデル化が難しい装置です。段数を M、成分数を N とすると、物質収支・エネルギー収支・相平衡に対して、\( M\cdot (3N+1) \) 個の非線形かつ強く連成した方程式を解くことになります。蒸留塔の最適化はさらに挑戦的な課題となります。そこで、最初は蒸留塔の設計やモデル化に、単純化した仮定を用いる“簡易的”手法を利用します。良く使われる簡易的な手法として、２成分系の蒸留に用いられる「マッケーブ・シール」法と「フェンスキー・アンダーウッド・ギランド」法という２つの方法があります。これらの手法は場合によっては驚くほど高い精度を示すことがあり、高精度が求められない迅速な解析にもしばしば利用されます。
マッケーブ・シール法による２成分系の蒸留
マッケーブ・ティーレ法では、１モルの液体が蒸発すると１モルの蒸気が凝縮する、すなわち、モル流量が一定である、と仮定しています。この仮定のもとでは、蒸留塔の各段は、沸点／露点計算と単純な物質収支を用いて順番に解くことができます。この仮定は、蒸発熱が近い成分同士の場合には比較的よく成り立ちます。
２成分系
とは、水に混合しているエタノール、のように溶媒と取り出したい成分の２つの蒸留特性に絞ってモデル化することを指していて、BioSTEAM では、３成分以上の系に対しては、注目する２つの成分のうち、軽い(沸点の低い)成分より揮発性の高い化学種は留出液側へ、重い方の成分より重い化学種は塔底液側へ分配されると仮定することで、マッケーブ・ティーレ法を拡張しています。これは注目している２つの成分以外の成分が微量の場合は、この仮定が妥当であることが多いとされています。
エタノールを含む発酵液からエタノールを分離するための蒸留塔をモデル化してみましょう。
- 2成分系の蒸留塔
- 結果
- マッケーブ・シール図
import biosteam as bst
import numpy as np
bst.nbtutorial()
# 最初に使用する化学成分を指定します。
bst.settings.set_thermo(['Water', 'Ethanol', 'Glycerol', 'Yeast'], db='BioSTEAM')
# 供給ストリームを定義します。
feed = bst.Stream(
'feed',
Water=1.08e+03,
Ethanol=586,
Glycerol=10,
Yeast=50,
units='kg/hr',
)
feed.vle(V=0, P=101325)
# 蒸留塔を定義します。
McCabeThiele = bst.BinaryDistillation(
'McCabeThiele', ins=feed,
outs=('distillate', 'bottoms_product'),
LHK=('Ethanol', 'Water'), # 注目している２つの成分のうち、沸点の低い成分(軽質成分)と高い成分を指定
y_top=0.79, # 軽質成分の、蒸留後の濃度を指定
x_bot=0.001, # 軽質成分の塔底での濃度を指定
k=1.25, # 実際の還流比の最小還流比に対する比率
is_divided=True, # 精留塔と分離塔が分離されているかどうか
)
McCabeThiele.simulate()
McCabeThiele.diagram()
McCabeThiele.show()
BinaryDistillation: McCabeThiele
ins...
[0] feed
phases: ('g', 'l'), T: 356.95 K, P: 101325 Pa
flow (kmol/hr): (l) Water 59.9
Ethanol 12.7
Glycerol 0.109
Yeast 2.21
outs...
[0] distillate
phase: 'g', T: 351.63 K, P: 101325 Pa
flow (kmol/hr): Water 3.37
Ethanol 12.7
[1] bottoms_product
phase: 'l', T: 372.85 K, P: 101325 Pa
flow (kmol/hr): Water 56.6
Ethanol 0.0566
Glycerol 0.109
Yeast 2.21
BinaryDistillationの特徴として、重質なキー成分に指定した水より沸点の高いグリセロールは全て塔底に分離されます。モデルの詳細な説明はこちら BinaryDistillationクラス。
McCabeThiele.results()| Divided Distillation Column | Units | McCabeThiele | |
|---|---|---|---|
| Electricity | Power | kW | 0.139 |
| Cost | USD/hr | 0.0109 | |
| Cooling water | Duty | kJ/hr | -6.37e+05 |
| Flow | kmol/hr | 435 | |
| Cost | USD/hr | 0.212 | |
| Low pressure steam | Duty | kJ/hr | 1.4e+06 |
| Flow | kmol/hr | 36.3 | |
| Cost | USD/hr | 8.63 | |
| Design | Theoretical feed stage | 24 | |
| Theoretical stages | 29 | ||
| Minimum reflux | Ratio | 0.8 | |
| Reflux | Ratio | 1 | |
| Rectifier stage | 41 | ||
| Stripper stages | 10 | ||
| Rectifier height | ft | 73 | |
| Stripper height | ft | 27.3 | |
| Rectifier diameter | ft | 3 | |
| Stripper diameter | ft | 3 | |
| Rectifier wall thickness | in | 0.5 | |
| Stripper wall thickness | in | 0.312 | |
| Rectifier weight | lb | 9.16e+03 | |
| Stripper weight | lb | 3.6e+03 | |
| Purchase cost | Rectifier trays | USD | 2.92e+04 |
| Stripper trays | USD | 1.03e+04 | |
| Rectifier tower | USD | 5.86e+04 | |
| Stripper platform and ladders | USD | 2.13e+04 | |
| Stripper tower | USD | 3.41e+04 | |
| Rectifier platform and ladders | USD | 9.69e+03 | |
| Condenser - Floating head | USD | 5.19e+03 | |
| Reflux drum - Vertical pressure vessel | USD | 7.58e+03 | |
| Reflux drum - Platform and ladders | USD | 1.85e+03 | |
| Pump - Pump | USD | 4.36e+03 | |
| Pump - Motor | USD | 141 | |
| Reboiler - Floating head | USD | 5.2e+03 | |
| Total purchase cost | USD | 1.88e+05 | |
| Utility cost | USD/hr | 8.85 |
結果の見方は機器ユニットの計算結果を見てください。
McCabeThiele.plot_stages()
フェンスキー–アンダーウッド–ギリランド法による多成分系の蒸留
フェンスキー–アンダーウッド–ギリランド法(FUG 法)は、多成分蒸留を扱うため、モル数の変化なし、蒸留段間での平均的な相対揮発度は一定であると仮定しています。この方法は、段ごとに揮発度が“良好に振る舞い”、かつ成分間の潜熱が類似している場合に、妥当な推算結果を与えることができます。ただし、アゼトロープ
(蒸留しても組成が変わらない状態)がある場合は適用できません。では、この簡易(ショートカットな)蒸留設計法であるFUG法を使って、炭化水素混合物を C5 未満成分（−C5）と C6 以上成分（C6+）に分離するモデル化を行ってみましょう。ここでは ShortcutColumnクラスを使っています。
- 多成分系の蒸留
- BinaryDistillationとの比較
- 沸点
import biosteam as bst
import numpy as np
bst.nbtutorial()
hydrocarbons = [
'Ethane', 'Propane', 'n-Butane',
'n-Pentane', 'n-Hexane', 'n-Heptane',
'Cyclohexane', 'Cycloheptane', 'Benzene',
'Toluene', 'Hydrogen'
]
bst.settings.set_thermo(hydrocarbons, pkg='Peng Robinson')
feed = bst.Stream('feed',)
feed.imol[hydrocarbons] = 242.6 * np.array([
7.54105070e-10, 2.56035674e-06, 5.34074413e-03,
1.850257e-01, 3.06487804e-01, 1.27207300e-01,
3.97638023e-02, 6.96044980e-04, 3.35385551e-01,
9.04125972e-05, 1.00013935e-18
])
feed.vle(T=340.3, P=366463)
FUG = bst.ShortcutColumn(
'FUG',
ins=[feed],
outs=['distillate', 'bottoms_product'],
LHK=['Cyclohexane', 'n-Heptane'],
P=feed.P,
Lr=0.98, # Light key recovery
Hr=0.98, # Heavy key recovery
k=1.25 # Reflux to minimum reflux
)
FUG.simulate()
FUG.show('cmol100')ShortcutColumn: FUG
ins...
[0] feed
phases: ('g', 'l'), T: 340.3 K, P: 366463 Pa
composition (%): (l) Ethane 7.54e-08
Propane 0.000256
n-Butane 0.534
n-Pentane 18.5
n-Hexane 30.6
n-Heptane 12.7
Cyclohexane 3.98
Cycloheptane 0.0696
Benzene 33.5
Toluene 0.00904
Hydrogen 1e-16
------------ 243 kmol/hr
outs...
[0] distillate
phase: 'g', T: 401.49 K, P: 366463 Pa
composition (%): n-Pentane 21.4
n-Hexane 35.5
n-Heptane 0.294
Cyclohexane 4.51
Cycloheptane 3.31e-08
Benzene 38.3
Toluene 2.47e-07
------------ 210 kmol/hr
[1] bottoms_product
phase: 'l', T: 402.88 K, P: 366463 Pa
composition (%): Ethane 5.56e-07
Propane 0.00189
n-Butane 3.94
n-Pentane 1.43e-09
n-Hexane 0.00763
n-Heptane 91.9
Cyclohexane 0.586
Cycloheptane 0.513
Benzene 2.97
Toluene 0.0667
Hydrogen 7.37e-16
------------ 32.9 kmol/hr
C5 未満成分（−C5）と C6 以上成分（C6+）に分離するとなっていますが、結果を見るとdistillate(塔頂)側にもC6のn-ヘキサンやベンゼンが混入していたり、軽質なn-ブタンやプロパンがbottoms_product(塔底)に混入してしまっています。これが現実に近いのかどうかは分かりませんが、キー成分より軽質、重質な成分を塔頂、塔底に分離してしまうBinaryDistillationと比較してみましょう。
D1 = bst.BinaryDistillation(
'D1',
ins=[feed],
outs=['distillate', 'bottoms_product'],
LHK=['Cyclohexane', 'n-Heptane'],
P=feed.P,
Lr=0.98, # Light key recovery
Hr=0.98, # Heavy key recovery
k=1.25 # Reflux to minimum reflux
)
D1.simulate()
D1.show('cmol100')BinaryDistillation: D1
ins...
[0] feed
phases: ('g', 'l'), T: 340.3 K, P: 366463 Pa
composition (%): (l) Ethane 7.54e-08
Propane 0.000256
n-Butane 0.534
n-Pentane 18.5
n-Hexane 30.6
n-Heptane 12.7
Cyclohexane 3.98
Cycloheptane 0.0696
Benzene 33.5
Toluene 0.00904
Hydrogen 1e-16
------------ 243 kmol/hr
outs...
[0] distillate
phase: 'g', T: 405.23 K, P: 366463 Pa
composition (%): Ethane 8.63e-08
Propane 0.000293
n-Butane 0.611
n-Pentane 21.2
n-Hexane 35.1
n-Heptane 0.291
Cyclohexane 4.46
Benzene 38.4
Hydrogen 1.14e-16
----------- 212 kmol/hr
[1] bottoms_product
phase: 'l', T: 417.02 K, P: 366463 Pa
composition (%): n-Heptane 98.7
Cyclohexane 0.63
Cycloheptane 0.551
Toluene 0.0716
------------ 30.6 kmol/hr
ここで、BioSTEAMのデータベースで沸点を調べておくと
import thermosteam as tmo
bp_list = []
for name in hydrocarbons:
chem = tmo.Chemical(name)
Tb_K = chem.Tb
Tb_C = Tb_K - 273.15
bp_list.append((name, Tb_K, Tb_C))
# 沸点（K）でソート
bp_list_sorted = sorted(bp_list, key=lambda x: x[1])
# 表示
print(f"{'Chemical':<15} {'Boiling Point [K]':>20}")
for name, Tb_K, Tb_C in bp_list_sorted:
print(f"{name:<15} {Tb_C:>20.2f}")Chemical Boiling Point [K] Hydrogen -252.78 Ethane -88.58 Propane -42.11 n-Butane -0.49 n-Pentane 36.06 n-Hexane 68.72 Benzene 80.07 Cyclohexane 80.71 n-Heptane 98.40 Toluene 110.60 Cycloheptane 118.80ですので、おおよそ、Cyclohexaneとn-Heptaneが境界となって、それより軽質なものと重質なものに分離されています。このように結果が違いますので、模擬したい状況に応じて使い分けが必要と思います。
より厳密な蒸留塔モデル
厳密平衡蒸留塔
モデルは、物質収支(Material)、平衡関係(Equilibrium)、各物質の割合の総和(Summation)、エンタルピー(Enthalpy)の頭文字を取ったMESH
を完全な方程式を解きます。そのために、蒸留段数、供給段/サイドカットの位置(段)、サイドカットの分割比を決める必要があります。その上で、さらに2 つの自由度が残るため、2 つの変数を指定する必要があります。通常は供給量に対する塔底製品の流量、コンデンサー部の還流量を与えます。その他の設定としては、リボイラーのボイラー比、コンデンサー、リボイラーの温度があります。
簡易(ショートカット)から厳密蒸留モデルへの変換
簡易(ショートカット)なモデルの計算結果を使って厳密な蒸留塔モデルに変換するto_rigorous_column
メソッドを使ってみましょう。
to_rigorous_column
メソッド- FUG 法から
to_rigorous_column
メソッド
MESH = McCabeThiele.to_rigorous_column()
# 以下と等価:
# MESH = MESHDistillation(
# ins=[feed],
# outs=[distillate, bottoms_product],
# N_stages=29,
# feed_stages=[24],
# LHK=('Ethanol', 'Water'),
# stage_specifications={0: ('Reflux', 1.0), -1: ('Flow', 0.78624)},
# P=[101325 + i*690.0 for i in range(29)]
# )
MESH.simulate()
MESH.show('cmol100')MESHDistillation
ins...
[0] feed
phases: ('g', 'l'), T: 356.95 K, P: 101325 Pa
composition (%): (l) Water 79.9
Ethanol 17
Glycerol 0.145
Yeast 2.95
-------- 75 kmol/hr
outs...
[0] distillate
phase: 'g', T: 351.57 K, P: 101325 Pa
composition (%): Water 21.1
Ethanol 78.9
Glycerol 1.95e-17
-------- 15.6 kmol/hr
[1] bottoms_product
phase: 'l', T: 375.75 K, P: 120645 Pa
composition (%): Water 95.3
Ethanol 0.751
Glycerol 0.183
Yeast 3.72
-------- 59.4 kmol/hr
より厳密に計算すると、塔頂にもグリセロールが混入するようです。実際には重質な成分でも、沸点以下の温度で蒸気圧がゼロではなく、さらに、有限段数/有限還流比の現実の蒸留塔ではこのような現象が起こるそうです。
次に、FUG 法の例も同様に厳密モデルに変換してみます。
MESH = FUG.to_rigorous_column()
# 以下と等価:
# MESH = MESHDistillation(
# ins=[feed],
# outs=[distillate, bottoms_product],
# N_stages=42,
# feed_stages=[18],
# LHK=('Cyclohexane', 'n-Heptane'),
# stage_specifications={0: ('Reflux', 1.629), -1: ('Flow', 0.13563)},
# P=[366463.0 + i*690.0 for i in range(42)]
# )
MESH.simulate()
MESH.show('cmol100')MESHDistillation
ins...
[0] feed
phases: ('g', 'l'), T: 340.3 K, P: 366463 Pa
composition (%): (l) Ethane 7.54e-08
Propane 0.000256
n-Butane 0.534
n-Pentane 18.5
n-Hexane 30.6
n-Heptane 12.7
Cyclohexane 3.98
Cycloheptane 0.0696
Benzene 33.5
Toluene 0.00904
Hydrogen 1e-16
------------ 243 kmol/hr
outs...
[0] distillate
phase: 'g', T: 401.2 K, P: 366463 Pa
composition (%): Ethane 8.72e-08
Propane 0.000296
n-Butane 0.618
n-Pentane 21.4
n-Hexane 35.5
n-Heptane 0.306
Cyclohexane 4.36
Cycloheptane 5.09e-08
Benzene 37.8
Toluene 3.77e-07
Hydrogen 1.16e-16
------------ 210 kmol/hr
[1] bottoms_product
phase: 'l', T: 554.64 K, P: 394753 Pa
composition (%): Ethane 1e-59
Propane 4.39e-37
n-Butane 1.06e-21
n-Pentane 1.95e-11
n-Hexane 0.00348
n-Heptane 91.8
Cyclohexane 1.51
Cycloheptane 0.513
Benzene 6.07
Toluene 0.0667
Hydrogen 1.15e-317
------------ 32.9 kmol/hr
BinaryDistillationとは逆に、全ての成分が塔頂からも塔底からも検出される結果になりました。
シミュレーションアルゴリズムへの導入
蒸留塔の解法となると、使える手法はすべて総動員となります。この複雑な問題を解くために、さまざまな方程式分割法と数値計算法のソフトウェアを用います。BioSTEAM では、バブルポイント法、サムレート法、インサイドアウト法、そして同時修正法（方程式指向法）を使用できます。
- バブルポイント法 バブルポイント法は通常の蒸留ではうまく機能するが、吸収／脱離では破綻することがある。
- サムレート法 サムレート法は吸着／脱離では非常に良好に機能するが、通常の蒸留では性能がかなり低い。
- インサイドアウト法 インサイドアウト法は、一般的に蒸留に対して最も堅牢な手法の一つ。BioSTEAM では、吸収／脱離用と通常の蒸留用の 2 種類のバージョンが準備されている。
- 同時修正法（方程式指向法） 同時修正法（トラストリージョン最適化を用いる）は、初期推定値が解に十分近ければ、幅広い問題に対して良好に収束できる。
BioSTEAM は、デフォルトでまずインサイドアウト法を試みます。簡易フラッシュ計算で塔を初期化した後、分配係数の分布に基づいて使用するバージョンを判断します。インサイドアウト法が収束しない場合は、他のアルゴリズムを試行します。
では、各手法の、炭化水素系蒸留塔の収束プロファイルを見てみましょう。
- 各手法の収束プロファイル
import matplotlib.pyplot as plt
feed = bst.Stream(
ID='feed', phase='l',
T=329.54, P=372615,
Ethane=1.829e-07, Propane=0.0006211,
**{'n-Butane': 1.296, 'n-Pentane': 44.89,
'n-Hexane': 74.35, 'n-Heptane': 30.86},
Cyclohexane=9.647, Cycloheptane=0.1689, Benzene=81.36,
Toluene=0.02193, Hydrogen=2.426e-16,
units='kmol/hr'
)
N_stages= 42
P_condenser = 351625
dP = 723.78
MESH = bst.MESHDistillation(
'MESH',
N_stages=N_stages, ins=[feed], feed_stages=[29],
outs=['distillate', 'bottoms'],
stage_specifications={
0: ('Reflux', 0.857),
-1: ('Flow', 0.76) # Bottoms product flow rate as a fraction of column feed
},
LHK=('Cyclohexane', 'n-Heptane'),
P=[P_condenser + i*dP for i in range(N_stages)],
use_cache=True,
)
x0 = MESH.hot_start()
residual_profiles = []
algorithms =(
'inside out',
'phenomena', # Equivalent to Wang-Hanke's bubble point algorithm
'simultaneous correction',
)
maxiter = 10000
maxtime = 20
residual_profiles = [
MESH.convergence_analysis(
maxiter, maxtime, x0=x0, algorithm=alg,
verbose=False,
plot=False,
) for alg in algorithms
]
colors = bst.utils.GG_colors
profile_colors = [colors.blue, colors.red, colors.purple, colors.green]
yticks = [-60, -50, -40, -30, -20, -10, 0, 10, 20]
fig, ax = plt.subplots(1, 1)
tmax = 0
for alg, color, residual_profile in zip(algorithms, profile_colors, residual_profiles):
tmax = max(tmax, residual_profile.time[-1])
plt.scatter(
residual_profile.time,
residual_profile.log_residual,
color=color.RGBn,
label=alg or 'default',
s=5,
)
plt.ylim([yticks[0], yticks[-1]])
plt.xlim([0, tmax])
plt.yticks(yticks)
plt.xlabel('Time [s]')
plt.ylabel('ln(residual)')
plt.legend()
plt.show()
このテストケースでは、多数の化学種と段数を含むため、インサイドアウト法が最も早く収束しました。
次のステップとなる題材
混合物中の分子間力はアゼトロープを生じることがあり、これが気液分離をより困難にします。典型的な例として、発酵液からバイオエタノールを分離する操作が挙げられます。エタノールのアゼトロープの存在により、単一の蒸留塔では分離が不可能となる場合があります。しかし、熱力学の知識を活用すれば、アゼトロープを破る、あるいはシフトさせることで、蒸留による実現可能な分離プロセスを設計することができます。今後の例では、圧力変換蒸留や抽出蒸留がどのようにアゼトロープをシフトまたは破壊し、アゼトロープ混合物の分離を可能にするかを示します。
