Dynamic simulation of a distillation column in python

distill.py

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()

Distillation model

This is a model of a simple binary distillation column.

To install:

pip install gekko numpy matplotlib

To run

python distill.py