4. Amundson_1958

This code follows the approach that was published by Amundson and Pontinen [Amundson1958].

Amundson1958

Amundson, NR and Pontinen, AJ. Multicomponent Distillation Calculations on a Large Digital Computer. Ind. Eng. Chem. 1958;50:730–736.

This approach takes advantage of the following mathematical trick.

4.1. Mathematical Trick

A system of equations such as

\[Ax=b\]

where \(A\) is a matrix of constants, \(x\) is a vector of unknowns, and \(b\) is a vector of constants, can be solved very efficiently if the matrix \(A\) can be made into the following form:

\[\begin{split}A = \begin{bmatrix} d_0 & u_0 & 0 & 0 & 0 \\ l_0 & d_1 & u_1 & 0 & 0 \\ 0 & l_1 & d_2 & u_2 & 0 \\ 0 & 0 & l_2 & d_3 & u_3 \\ 0 & 0 & 0 & l_3 & d_4 \end{bmatrix}\end{split}\]

This special type of matrix is called a tridiagonal_matrix.

It can be solved, for example, if we define the following vectors

\[\begin{split}\begin{align} l &= [l_0, l_1, l_2, l_3] \\ d &= [d_0, d_1, d_2, d_3, d_4] \\ u &= [u_0, u_1, u_2, u_3] \\ b &= [b_0, b_1, b_2, b_3] \end{align}\end{split}\]

or, with python code, like

l = [l_0, l_1, l_2, l_3]
d = [d_0, d_1, d_2, d_3, d_4]
u = [u_0, u_1, u_2, u_3]
b = [b_0, b_1, b_2, b_3]

the \(5\times5\) matrix can be solved efficiently using the following code:

from scipy.sparse import linalg, diags
A = diags(
    diagonals=[l, d, u],
    offsets=[-1, 0, 1],
    shape=(5, 5),
)
x = linalg.spsolve(A, b)

This mathematical trick is used throughout the solution algorithm.

4.2. The algorithm

The complete algorithm is shown in the diagram below.

digraph { input -> initial_guess -> calc_K -> mass_balance -> bubble_pt -> converg_1; converg_1 -> calc_K[label="No",color=red]; converg_1 -> energy_balance[label="Yes",color=green]; energy_balance -> converg_2; converg_2 -> finish[label="Yes",color=green]; converg_2 -> mass_balance[label="No",color=red]; input [shape=polygon,side=4,label="1. Input equilibrium and enthalpy data.\n Input design variables"]; initial_guess [shape=polygon,side=4,label="2. Pick L, V, T\n on every stage"]; calc_K [shape=polygon,side=4,label="3. Calculate K for each component\n on each stage"]; mass_balance [shape=polygon,side=4,label="4. Solve component mass\n balances in matrix form",color=lightblue,style=filled]; bubble_pt [shape=polygon,side=4,label="5. Calculate T on each stage\n using bubble point calculation"]; converg_1 [shape=diamond,label="6. Determine whether\n T is converged\n for all stages"]; energy_balance [shape=polygon,side=4,label="7. Use energy balance to calculate\n L, V on every stage",color=lightblue,style=filled]; converg_2 [shape=diamond,label="8. Determine whether\n L, V are converged\n for all stages"]; finish [shape=ellipse,label="9. Finished"]; }

4.3. Tutorial

4.3.1. Step 1, Input

In Step 1, the property data and design variables are input. The property data is described in Equilibrium Properties. A model with this property data and input design variables is created by the following code:

>>> from distillation.amundson_1958 import Model

>>> model = Model(
        components=['n-Butane', 'n-Pentane'],
        F=1000.,              # feed flow rate [kmol/h]
        P=101325.*2,          # pressure (constant) [Pa]
        z_feed=[0.45, 0.55],  # mole fraction n-Butane = 0.45
        RR=1.,                # reflux ratio [L/D]
        D=400.,               # distillate flow rate [kmol/h]
        N=3,                  # number of equilibrium contacts
        feed_stage=2,         # stage at which feed goes in
)

Currently, the input design flow rates must be specified in kmol/h and the pressure must be specified in Pa. The method called in this process is distillation.amundson_1958.main.Model.__init__() (see Class Method Reference).

4.3.2. Step 2, Make Initial Guess

Before beginning the solution procedure, we need to have an initial guess for the liquid flow rates, vapor flow rates, and temperatures on each stage \(j\). That is, we need to generate initial guesses for \(L_j\), \(V_j\), and \(T_j\).

First we calculate the feed temperature (which is a saturated liquid) using a bubble point calculation (see Bubble Point). We set the temperatures \(T_j\) to be the same as the feed temperature. This step is performed automatically when the model is initialized with input parameters. We can check to see the feed temperature calculated for our model by the following

>>> model.T_feed
306.37018410667076

Then, we generate the initial guesses for these values by assuming constant molal overflow (the Lewis method). Computationally, we can do this using the following

>>> model.generate_initial_guess()

We can then check the arrays of liquid flow rates (\(L\)) and vapor flow rates (\(V\)) obtained by the calculations as below

>>> model.L
array([ 400.,  400., 1400.,  600.])
>>> model.V
array([   0., 800., 800., 800.])

or, looking at the liquid flow rates by stage more specifically, .. code-block:: python

>>> for stage, L_j in enumerate(model.L):
...    print('L_j', stage, L_j)
...
L_j 0 400.0
L_j 1 400.0
L_j 2 1400.0
L_j 3 600.0

Here, we see that the liquid flow rates in the enriching section are \((L/D)\times D=1\times 400=400=L\). Since the feed stage is stage 2, we notice that the liquid leaving this stage is equal to \(\overline{L}=L + F=400 + 1000\). Finally, we notice that the liquid leaving the bottom stage, the partial reboiler, is the same as the bottoms flow rate that would be calculated from an overall balance, \(B = F - D = 1000 - 400 = 600\).

>>> for stage, V_j in enumerate(model.V):
...    print('V_j', stage, V_j)
...
V_j 0 0.0
V_j 1 800.0
V_j 2 800.0
V_j 3 800.0

We note that there is no vapor leaving stage 0 because the column is equipped with a total condensor. We note that V[1] = model.RR*model.D + model.D, which arises from the mass balance at the top of the column. Similarly, we also realize that V[3] = model.L[2] - model.B.

4.3.3. Step 3, Update \(K\)-values

In Step 3, we calculate \(K_{i,j}\), the \(K\)-value for each component \(i\) on stage \(j\). We can do this now for our model by the following:

>>> model.update_K_values()

We can think of this step as doing something similar to the following:

for i in model.components:
    for j in range(model.N + 1):
        model.K[i][j] = K_func(i, T[j], model.P_feed)

where the K_func function returns the \(K\)-value for component \(i\) on stage \(j\). The for loops iteratively store the \(K\)-values in each component of the amtrix.

In our example, we can checkout our new \(K\)-values as follows:

>>> import pprint  # for pretty printing
>>> pprint.pprint(model.K)
{'n-Butane': array([1.61320297, 1.61320297, 1.61320297, 1.61320297]),
 'n-Pentane': array([0.49828848, 0.49828848, 0.49828848, 0.49828848])}

As expected, the values for butane are higher for the light component, butane. The values are also the same for all stages because the temperature of all stages are currently fixed at the feed temperature.

Behind the scenes, we then store the current stage temperatures in an attribute T_old which we will need later in Step 6, T_j Convergence Check to determine whether the temperature has converged.

The code for this step is depicted in distillation.amundson_1958.main.Model.update_K_values().

4.3.4. Step 4, Solve Component Mass Balances

The two stages in which tridiagonal matrices are used for computational efficiency are Step 4, Solve Component Mass Balances and Step 7, Solve Energy Balances, which are depicted in colored blue boxes.

In this step, the matrices are tridiagonal for the complete mass balances for each component (as opposed to all components together). The variables that are calculated are \(l_{ij}\) for each component \(i\) and each stage \(j\). This new variable is introduced where \(l_{ij}=x_iL_j\), and can be interpreted as a component flow rate. By formulating the mass balance equations into function of \(l_{ij}\) the system of equations can be converted in to a tridiagonal matrix. Computationally, we can perform this step using our model by the following:

>>> for i in model.components:
...    model.solve_component_mass_bal(i)
...
>>>

First, the tridiagonal matrix is generated. Then, the efficient approach described in Mathematical Trick is used. to solve the equations.

We can see the results of the calculation by probing the model attributes as follows

>>> import pprint # pretty printing
>>> pprint.pprint(model.l)
{'n-Butane': array([288.88918652, 179.07801532, 507.65007098, 161.11081348]),
 'n-Pentane': array([ 74.88288594, 150.28018773, 790.77762525, 475.11711406])}

Since \(\sum_i l_{i,j} = L_j\), it is worthwhile to compare our calculations to what was calculated for \(L_j\) with the Lewis method (i.e. assuming Constant Molal Overflow)

>>> for stage in model.stages:
...     print(stage, model.L[stage], model.l['n-Butane'][stage]+model.l['n-Pentane'][stage])
...
0 400.0 363.7720724598715
1 400.0 329.3582030449366
2 1400.0 1298.4276962288486
3 600.0 636.2279275401286
>>>

From this analysis, we realize that our initial guess for the liquid and vapor flow rates was not that bad.

4.3.5. Step 5

In Step 2, Make Initial Guess, we performed a Bubble Point calculation to determine the temperature of the feed. We do the same thing here, except we do it multiple times (i.e., for each stage). However, we have not calculated the liquid phase mole fractions explicitly yet. We can use the liquid component mass balances calculated in Step 4, Solve Component Mass Balances to calculate the mole fractions as

\[x_{ij} = \frac{l_{ij}}{\sum_k l_{kj}}\]

and then perform the Bubble Point calculations with these mole fractions.

Computationally, we can perform this step with our model using the following

>>> model.T
array([306.37018411, 306.37018411, 306.37018411, 306.37018411])
>>> model.update_T_values()
>>> model.T
array([295.34137166 302.94861315 308.72858323 314.97039634])

And we can see how the temperatures change from the initial guess (the feed temperature) to the results of the first calculation.

4.3.6. Step 6, \(T_j\) Convergence Check

In this step, we determine whether the temperatures on all stages \(j\), \(T_j\), have converged. Mathematically, we require the following

\[\sqrt{\left(T_{j,\mathrm{new}} - T_{j,\mathrm{old}}\right)^2} < \epsilon\]

The temperature tolerance \(\epsilon\) is attribute distillation.amundson_1958.main.Model.temperature_tol; It is currently set to 0.01.

We can see from Step 5 that the temperature has clearly changed by more than 0.01~K between the first iteration. When we perform this step, we find that the temperature has not converged.

>>> model.T_is_converged()
False

Following the algorithmic diagram, we go back to Step 3 and perform a few more iterations.

Its convenient to perform these iterations with a while loop, as follows

while not model.T_is_converged():
...     model.update_K_values()
...     for i in model.components:
...         model.solve_component_mass_bal(i)
...     model.update_T_values()
...     print(model.T)
...
[295.50131711 303.41982848 309.35773209 315.96224522]
[295.6166139  303.62102821 309.52271272 316.15309479]
[295.65094925 303.67756389 309.56508571 316.19923696]
[295.66006297 303.69238402 309.57598554 316.21097278]
[295.66242518 303.69621481 309.57879188 316.21398756]

The temperature is converged after 5 more iterations, as we can see from the output above. Since the temperature has converged, we can now proceed to the next step.

4.3.7. Step 7, Solve Energy Balances

In Step 7, we solve the energy balances. Here, all the balances can again be combined into one banded matrix.

For a general stage \(j\), the energy balance is

\[L_j h_j + V_jH_j = V_{j+1}H_{j+1} + L_{j-1}h_{j-1} + F_j h_{j,\mathrm{feed}} + Q_j\]

The liquid flow rates can be substituted in using the following relationship

\[L_k = V_{k+1} - D + \sum_{m=0}^k F_m\]

So that the system of equations only becomes functions of \(V_j\) and \(V_{j+1}\).

Todo

Finish description of energy balance

A helper function for creation and solution of these matrices is provided in the model. In other words, the energy balances can be solved by the following code:

>>> model.L
array([ 400.,  400., 1400.,  600.])
>>> model.V
array([  0., 800., 800., 800.])
>>> model.solve_energy_balances()
>>> model.L
array([ 400.        ,  312.3605968 , 1300.10206406,  600.        ])
>>> model.V
array([  0.        , 730.64500782, 712.3605968 , 700.10206406])
>>>

Where we have added extra steps so that we can check the values of model.L and model.V before and after solving the energy balances.

4.3.8. Step 8, \(V_j\), \(L_j\) Convergence Check

In this step, we determine if the simulation has converged. If the following holds true for all stages \(j\)

\[\sqrt{\left(\frac{X_{j,\mathrm{new}} - X_{j,\mathrm{old}}}{X_{j,\mathrm{new}}}\right)^2} < \epsilon\]

for each variable \(X=V\) and \(X=L\). The flow rate tolerance, \(\epsilon\), is found in the attribute distillation.amundson_1958.main.Model.flow_rate_tol. The code for this step is found in distillation.amundson_1958.main.Model.flow_rates_converged()

Obviously, our flow rates have not converged after just one iteration, as we have seen that several values have changed by \(\\approx 100\) kmol/h. We will need to come up with a looping algorithm to get the simulation to converge, as described in the next section.

4.4. Putting it all together

So far, we can combine our code from the tutorial into the following

>>> from distillation.amundson_1958 import Model

>>> model = Model(
        components=['n-Butane', 'n-Pentane'],
        F=1000.,              # feed flow rate [kmol/h]
        P=101325.*2,          # pressure (constant) [Pa]
        z_feed=[0.45, 0.55],  # mole fraction n-Butane = 0.45
        RR=1.,                # reflux ratio [L/D]
        D=400.,               # distillate flow rate [kmol/h]
        N=3,                  # number of equilibrium contacts
        feed_stage=2,         # stage at which feed goes in
)
>>> model.generate_initial_guess()
>>> while not model.T_is_converged():
...     model.update_K_values()
...     for i in model.components:
...         model.solve_component_mass_bal(i)
...     model.update_T_values()
...
>>> model.solve_energy_balances()

But we need to figure out how to do the last iteration to automatically try to converge. It turns out, we can do this with a nested while loop, as shown below

>>> while not model.flow_rates_converged():
...     for i in self.components:
...         model.solve_component_mass_bal(i)
...     model.update_T_values()
...     while not model.T_is_converged():
...         model.update_K_values()
...         for i in model.components:
...             model.solve_component_mass_bal(i)
...         model.update_T_values()
...     model.solve_energy_balances()
...

This is what happens when the user invokes model.run() For our example, it turns out that it takes about 4 loops through the step 4-8 loop, and then we will be converged. The final results for the flow rates can be calculated as

>>> model.L
array([ 400.        ,  353.80944637, 1340.92773179,  600.        ])
>>> model.V
array([  0.        , 733.53885641, 753.80944637, 740.92773179])

which is still quite similar to what we calculated by the Lewis method. We can also take a look at the stage temperatures and mole fractions as below

>>> model.T
array([295.60405689, 303.58875055, 309.37291267, 315.8254746 ])
>>> list(model.x_ij_expr('n-Butane', i) for i in model.stages)
[0.7648435805973592, 0.5587863474554224, 0.38180690272177853, 0.24010427960176078]
>>> list(model.y_ij_expr('n-Butane', i) for i in model.stages)
[0.902879383384724, 0.8342045394354831, 0.6681268674620511, 0.4965317369291433]

As expected, the distillate is greatly enriched in the light component (n-Butane).

4.5. Class Method Reference

class distillation.amundson_1958.main.Model(components: list = None, F: float = 0.0, P: float = 101325.0, z_feed: list = None, RR: float = 1, D: float = 0, N: int = 1, feed_stage: int = 0, T_feed_guess: float = 300.0)[source]
H_j_rule(stage)[source]

Enthalpy of vapor on stage j. Calculated for ideal mixture

\[H_j = \sum_i y_{ij}H^*_i(T_j)\]

where the asterisk indicates the pure component enthalpy

Todo

convert y mole fractions to dynamic expression

Returns

\(H_j\) [J/kmol]

H_pure_rule(c, T)[source]

Rule for vapor enthalpy of pure component

Q_condenser_rule()[source]

Condenser requirement can be determined from balances around total condenser

T_is_converged()[source]

In this step, we determine whether the temperatures on all stages \(j\), \(T_j\), have converged. Mathematically, we require the following

\[\sqrt{\left(T_{j,\mathrm{new}} - T_{j,\mathrm{old}}\right)^2} < \epsilon\]

The temperature tolerance \(\epsilon\) is attribute distillation.amundson_1958.main.Model.temperature_tol; It is currently set to 0.01.

Returns

True if T is converged, else False

flow_rates_converged()[source]

Determine if flow rates are converged

Use the mathematical criterion in Model.is_below_relative_error()

generate_initial_guess()[source]

Before beginning the solution procedure, we need to have an initial guess for the liquid flow rates, vapor flow rates, and temperatures on each stage \(j\). That is, we need to generate initial guesses for \(L_j\), \(V_j\), and \(T_j\).

First we calculate the feed temperature (which is a saturated liquid) using a bubble point calculation (see Bubble Point). We set the temperatures \(T_j\) to be the same as the feed temperature. This step is performed automatically when the model is initialized with input parameters. We can check to see the feed temperature calculated for our model by the following

>>> model.T_feed
306.37018410667076

Then, we generate the initial guesses for these values by assuming constant molal overflow (the Lewis method). Computationally, we can do this using the following

>>> model.generate_initial_guess()

We can then check the arrays of liquid flow rates (\(L\)) and vapor flow rates (\(V\)) obtained by the calculations as below

>>> model.L
array([ 400.,  400., 1400.,  600.])
>>> model.V
array([   0., 800., 800., 800.])

or, looking at the liquid flow rates by stage more specifically, .. code-block:: python

>>> for stage, L_j in enumerate(model.L):
...    print('L_j', stage, L_j)
...
L_j 0 400.0
L_j 1 400.0
L_j 2 1400.0
L_j 3 600.0

Here, we see that the liquid flow rates in the enriching section are \((L/D)\times D=1\times 400=400=L\). Since the feed stage is stage 2, we notice that the liquid leaving this stage is equal to \(\overline{L}=L + F=400 + 1000\). Finally, we notice that the liquid leaving the bottom stage, the partial reboiler, is the same as the bottoms flow rate that would be calculated from an overall balance, \(B = F - D = 1000 - 400 = 600\).

>>> for stage, V_j in enumerate(model.V):
...    print('V_j', stage, V_j)
...
V_j 0 0.0
V_j 1 800.0
V_j 2 800.0
V_j 3 800.0

We note that there is no vapor leaving stage 0 because the column is equipped with a total condensor. We note that V[1] = model.RR*model.D + model.D, which arises from the mass balance at the top of the column. Similarly, we also realize that V[3] = model.L[2] - model.B.

h_feed_rule(stage)[source]

Enthalpy of liquid in feed mixture Calculated for ideal mixture

\[h = \sum_i x_{ij}h^*_i(T_j)\]

where the asterisk indicates the pure component enthalpy

Returns

\(h\) [J/kmol]

h_j_rule(stage)[source]

Enthalpy of liquid on stage j. Calculated for ideal mixture

\[h_j = \sum_i x_{ij}h^*_i(T_j)\]

where the asterisk indicates the pure component enthalpy

Returns

\(h_j\) [J/kmol]

h_pure_rule(c, T)[source]

rule for liquid enthalpy of pure component

is_below_relative_error(new, old)[source]

Determine relative error between two vectors

\[\sqrt{\left(\frac{X_{\mathrm{new}} - X_{\mathrm{old}}}{X_{\mathrm{new}}}\right)^2} < \epsilon\]

The flow rate tolerance, \(\epsilon\), is found in the attribute distillation.amundson_1958.main.Model.flow_rate_tol

Parameters

  • new –

  • old –

Return type

bool

solve_component_mass_bal(component)[source]

Solve component mass balances

solve_energy_balances()[source]

Solve energy balances

update_K_values()[source]

In Step 3, we calculate \(K_{i,j}\), the \(K\)-value for each component \(i\) on stage \(j\). We can do this now for our model by the following:

>>> model.update_K_values()

We can think of this step as doing something similar to the following:

for i in model.components:
    for j in range(model.N + 1):
        model.K[i][j] = K_func(i, T[j], model.P_feed)

where the K_func function returns the \(K\)-value for component \(i\) on stage \(j\). The for loops iteratively store the \(K\)-values in each component of the amtrix.

In our example, we can checkout our new \(K\)-values as follows:

>>> import pprint  # for pretty printing
>>> pprint.pprint(model.K)
{'n-Butane': array([1.61320297, 1.61320297, 1.61320297, 1.61320297]),
 'n-Pentane': array([0.49828848, 0.49828848, 0.49828848, 0.49828848])}

As expected, the values for butane are higher for the light component, butane. The values are also the same for all stages because the temperature of all stages are currently fixed at the feed temperature.

Behind the scenes, we then store the current stage temperatures in an attribute T_old which we will need later in Step 6, T_j Convergence Check to determine whether the temperature has converged.

The code for this step is depicted in distillation.amundson_1958.main.Model.update_K_values().

update_T_values()[source]

Update temperatures in all stages by performing bubble point calculation

Todo

vectorize with matrix multiplication

x_ij_expr(i, j)[source]
Parameters

  • i – component name

  • j – stage number

Returns

mole fraction on stage

y_ij_expr(i, j)[source]
Parameters

  • i – component name

  • j – stage number

Returns

gas-phase mole fraction on stage