Last active
May 24, 2022 08:29
-
-
Save marcosfelt/fd7625441f1501f107d7cb3f7a09217d to your computer and use it in GitHub Desktop.
Dynamic simulation of a distillation column in python
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from gekko import GEKKO | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from typing import List, Union, Optional | |
class BinaryDistillationColumn: | |
""" | |
Notes | |
------ | |
Stages numbered from the top with condenser as index 0 and reboiler as index n_trays+1 | |
""" | |
def __init__( | |
self, | |
m: Optional[GEKKO], | |
# feed_rate: Union[List[float], float], | |
# feed_composition: float, | |
# dynamic: bool = False, | |
n_trays: int = 30, | |
feed_stage: int = 17, | |
feed_pressure: float = 1000, | |
tray_holdup: float = 0.25, | |
condenser_holdup: float = 0.5, | |
reboiler_holdup: float = 0.1, | |
feed_temperature: float = 343, | |
q_reb: Optional[float] = None, | |
temp_cond: Optional[float] = None, | |
# q_cond: Optional[float] = None, | |
pressure_drop_stripping: float = 1.2, | |
pressure_drop_rectifying: float = 2.0, | |
pressure_top: float = 1000.0, | |
murphee_efficiency_rectifying: float = 0.35, | |
murphee_efficiency_stripping: float = 0.62, | |
): | |
self.N = n_trays | |
self.feed_stage = feed_stage | |
self.Ptop = pressure_top | |
self.dPrect = pressure_drop_rectifying | |
self.dPstrip = pressure_drop_stripping | |
self.murphee_efficiency_rectifying = murphee_efficiency_rectifying | |
self.murphee_efficiency_stripping = murphee_efficiency_stripping | |
### Define constants ### | |
# Reflux Ratio | |
self.rr = m.Param(value=0.7) | |
# Feed flowrate (mol/min) | |
self.Feed = m.Param(value=2) | |
# Mole fraction of feed | |
self.x_Feed = m.Param(value=0.5) | |
# Total molar holdup on each tray | |
self.n = [m.Const(value=tray_holdup) for _ in range(self.N + 2)] | |
self.n[self.N + 1].value = reboiler_holdup | |
self.n[0].value = condenser_holdup | |
# mole fraction of component A | |
self.x = [] | |
for i in range(self.N + 2): | |
self.x.append(m.Var(0.3)) | |
# Temperatures | |
self.T = [ | |
m.Var( | |
value=273.15 + 80, | |
name=f"Temperature[{i}]", | |
lb=273.15 + 65, | |
ub=273.15 + 100, | |
) | |
for i in range(self.N + 2) | |
] | |
### Define intermediates ### | |
# Pressures | |
self.P = [ | |
m.Intermediate(self.Ptop - i * self.dPrect) | |
if i < self.feed_stage | |
else m.Intermediate(self.Ptop - i * self.dPstrip) | |
for i in range(self.N + 2) | |
] | |
# Distillate flowrate (mol/min) | |
self.D = m.Intermediate(0.5 * self.Feed) | |
# Liquid flowrate in rectification section (mol/min) | |
self.L = m.Intermediate(self.rr * self.D) | |
# Vapor Flowrate in column (mol/min) | |
self.V = m.Intermediate(self.L + self.D) | |
# Liquid flowrate in stripping section (mol/min) | |
self.FL = m.Intermediate(self.Feed + self.L) | |
# Bottoms | |
self.B = m.Intermediate(self.Feed - self.D) | |
def p_star(T, c): | |
# From NIST | |
# Methanol = 1 | |
if c == 0: | |
A, B, C = 5.15853, 1569.613, -34.846 | |
# Water = 2 | |
elif c == 1: | |
A, B, C = 4.6543, 1435.264, -64.848 | |
return 1000 * 10 ** ((A - (B / (T + C)))) | |
self.p_star_0 = [p_star(temp, 0) for temp in self.T] | |
self.p_star_1 = [p_star(temp, 1) for temp in self.T] | |
# vapor mole fraction of Component A | |
# For trays calculated using murphee effciiencey below | |
# For condenser and reboiler, from the equilibrium assumption and mole balances | |
# 1) vol = (yA/xA) / (yB/xB) | |
# 2) xA + xB = 1 | |
# 3) yA + yB = 1 | |
self.y = [m.Var(value=0.3) for _ in range(self.N + 2)] | |
# Condenser and reboiler | |
for i in [0, self.N + 1]: | |
# Relative volatility = (yA/xA)/(yB/xB) = KA/KB = alpha(A,B) | |
vol = m.Var(value=3.5) | |
m.Equation(vol == self.p_star_0[i] / self.p_star_1[i]) | |
m.Equation(self.y[i] == self.x[i] * vol / (1 + (vol - 1) * self.x[i])) | |
# Component balances | |
m.Equations(self.create_component_expressions()) | |
# Pressure balances | |
m.Equations(self.create_pressure_balances()) | |
# Murphee efficiencies | |
m.Equation(self.create_murphee_efficiencies()) | |
def create_component_expressions(self): | |
equations = [] | |
for i in range(self.N + 2): | |
# Right side of ODE | |
if i == 0: | |
# condenser | |
rhs = self.V * (self.y[1] - self.x[0]) | |
elif i == self.N + 1: | |
# reboiler | |
rhs = ( | |
self.FL * self.x[self.N] | |
- self.B * self.x[self.N + 1] | |
- self.V * self.y[self.N + 1] | |
) | |
elif i == self.feed_stage: | |
# feed stage | |
rhs = ( | |
self.Feed * self.x_Feed | |
+ self.L * self.x[i - 1] | |
- self.FL * self.x[i] | |
- self.V * (self.y[i] - self.y[i + 1]) | |
) | |
else: | |
# non-feed trays | |
if i > self.feed_stage: | |
l = self.FL | |
else: | |
l = self.L | |
rhs = l * (self.x[i - 1] - self.x[i]) - self.V * ( | |
self.y[i] - self.y[i + 1] | |
) | |
# Left hand side of ODE | |
lhs = self.n[i] * self.x[i].dt() | |
equations.append(lhs == rhs) | |
return equations | |
def create_pressure_balances(self): | |
equations = [] | |
for i in range(self.N + 2): | |
rhs = self.x[i] * self.p_star_0[i] + (1 - self.x[i]) * self.p_star_1[i] | |
lhs = self.P[i] | |
equations.append(lhs == rhs) | |
return equations | |
def create_murphee_efficiencies(self): | |
equations = [] | |
for i in range(1, self.N + 1): | |
if i < self.feed_stage: | |
me = self.murphee_efficiency_rectifying | |
else: | |
me = self.murphee_efficiency_stripping | |
rhs = ( | |
me * self.p_star_0[i] / self.P[i] * self.x[i] + (1 - me) * self.y[i + 1] | |
) | |
lhs = self.y[i] | |
equations.append(lhs == rhs) | |
return equations | |
def plot(self, time): | |
fig = plt.figure(figsize=(10, 10)) | |
# Flows | |
ax = fig.add_subplot(2, 2, 1) | |
ax.plot(time, self.Feed.value, label="Feed") | |
ax.plot(time, self.D.value, label="Distillate") | |
ax.plot(time, self.L.value, label="Reflux") | |
ax.plot(time, self.B.value, label="Bottoms") | |
ax.legend(loc="best") | |
ax.set_ylabel("Molar flow rate (mol/min)") | |
# Pressures | |
ax = fig.add_subplot(2, 2, 2) | |
ax.plot(time, self.P[0].value, label="Reboiler") | |
for tray in [5, 10, 15, 20, 25, 31]: | |
ax.plot(time, self.P[tray].value, label=f"Tray {tray-1}") | |
ax.plot(time, self.P[self.N + 1].value, label="Condenser") | |
ax.set_ylabel("Pressure [mbar]") | |
ax.legend(loc="best") | |
# Compositions | |
ax = fig.add_subplot(2, 2, 3) | |
ax.plot(time, self.x[0].value, "r--", label="Condenser") | |
for tray in [5, 10, 15, 20, 25, 31]: | |
ax.plot(time, self.x[tray].value, label=f"Tray {tray-1}") | |
ax.plot(time, self.x[self.N + 1].value, "k-", label="Reboiler") | |
ax.set_ylabel("Composition") | |
ax.legend(loc="best") | |
# Temperatures | |
ax = fig.add_subplot(2, 2, 4) | |
ax.plot(time, self.T[0].value, label="Reboiler") | |
for tray in [5, 10, 15, 20, 25, 31]: | |
ax.plot(time, self.T[tray].value, label=f"Tray {tray-1}") | |
ax.plot(time, self.T[self.N + 1].value, label="Condenser") | |
ax.set_ylabel("Temperatures [mbar]") | |
ax.legend(loc="best") | |
fig.supxlabel("Time (min)") | |
return fig, ax | |
if __name__ == "__main__": | |
def p_star(T, c): | |
# From NIST | |
# Methanol = 1 | |
if c == 0: | |
A, B, C = 5.15853, 1569.613, -34.846 | |
# Water = 2 | |
elif c == 1: | |
A, B, C = 4.6543, 1435.264, -64.848 | |
return 1000 * 10 ** ((A - (B / (T + C)))) | |
T_degC = 80 | |
print(f"Temperature: {T_degC+273.1}K") | |
print("Vapor pressure 1:", p_star(273.1 + T_degC, 0)) | |
print("Vapor pressure 2:", p_star(273.1 + T_degC, 1)) | |
print("Relative volatility:", p_star(273.1 + T_degC, 0) / p_star(273.1 + T_degC, 1)) | |
m = GEKKO(remote=False) | |
column = BinaryDistillationColumn(m) | |
# steady state solution | |
m.solve() | |
print(column.x) # with RR=0.7 | |
# switch to dynamic simulation | |
m.options.imode = 4 | |
nt = 61 | |
m.time = np.linspace(0, 60, 61) | |
# step change in reflux ratio | |
rr_step = np.ones(nt) * 0.7 | |
rr_step[10:] = 3.0 | |
column.rr.value = rr_step | |
# xfeed_step = np.ones(nt) * column.x_Feed.value | |
# xfeed_step[10:] = xfeed_step[0] * 1.2 | |
# column.x_Feed.value = xfeed_step | |
# Solve dynamics | |
m.solve() | |
# Make plots | |
fig, ax = column.plot(m.time) | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Distillation model
This is a model of a simple binary distillation column.
To install:
To run