Steady‐State Modeling of Equilibrium Distillation

Abstract

In this chapter, an algorithm for the solution of the mathematical model featuring steady‐state multicomponent distillation columns is analyzed and applied in the case study of the separation of hydrocarbon mixture. The development of the model has assumed each stage outlet streams in thermodynamic equilibrium in the phases liquid and vapor. The modeling of liquid was considerate and the non‐ideality behavior was described by activity. The non‐ideality of gas phase was calculated by Peng‐Robinson equation of state. The model consists of a set of nonlinear algebraic equations. The algorithm and numerical procedure to solve a set of equations are presented in a sequential, general and very simple form. A methodology to produce the good initial guess was defined based on rude simplifications of the system. In the study case, the initial estimates generated by the method are very good, being only about 20% far from the simulation results and considering a tolerance of 10−10, the convergence was obtained with 28 iterations.

1. Introduction

The analysis of a plant, by simulation, within the development of new processes may frequently show beforehand whether it is technically and economically viable. The simulation process in already operating plants may optimize the operational conditions for better quality products, decrease energy consumption and other losses in the process [1].

The design of a multicomponent distillation column by phenomenological models is quite complex due to the large number of parameters and variables involved [2] and also usually, it is required to solve the set of nonlinear equations and differential equations. Mathematical modeling is a powerful and useful tool in the design of this type of equipment, it assists the control and optimization column and therefore, the project and operating costs can be significantly reduced.

The use of distillation as separation method is disseminated by the modern chemical industry. One can find it in almost all industrial chemical processes where liquid separation is required. Common commercial binary distillations are as follows: water/ethylene glycol, benzene/toluene, o‐xylene/m‐xylene, isopentane/n‐pentane, ethylbenzene/styrene, water/acetic acid, ethanol/water, among others [2]. There are a lot of hydrocarbons mixtures that can be cited as commercial multicomponent distillation examples.

With all the foregoing, it is clear the need to provide the theoretical study on a simplified model for evaluating the possible separation processes using distillation columns. The simplest mathematical model for a distillation column is obtained, considering that all stages outlet streams (liquid and vapor) are in thermodynamic equilibrium. What represents, with no much accuracy, what happens in an actual process; however, this study is very important to get an idea of the theoretically best result that can be achieved in the process in question.

Thus, the equations of a distillation column model are obtained from mass and energy balances, mass balances by component and iso‐fugacity equations. The equations that represent such model are highly nonlinear, particularly those describing the phase equilibria and energy balances.

The solution of a set of nonlinear equations is quite difficult and generally requires that good initial guesses are provided in the way the method presents convergence [3]. Thus, the solution of the obtained model is divided into several steps, where in each step is calculated a set of model unknowns (mole fraction, temperature, flow rate, etc.), such that it is not necessary the solution of nonlinear equations sets, but only sets of linear equations, evaluation of explicit algebraic expressions and root finding of single variable. Also will be presented a methodology to generate good initial guesses and an example will be studied using the methodology presented.

The simplest methods used for solving the modeling of distillation columns are the graphical ones like McCabe‐Thiele [4] and Ponchon‐Savarit [5, 6]. The equation obtained by modeling the steady state of equilibrium distillation columns forms a set of highly nonlinear equations (MESH equations, obtained from mass and energy balances, phase equilibrium relations and mole fractions summations) that are normally solved all at the same time by Newton‐Raphson and like‐one method. Another, much used, type of method is that one which uses “tearing equations,” that is, large sparse systems of algebraic equations are split into smaller systems and solved in a sequential form [7–11]; it makes the solution process simpler but, as a consequence, can occur some instability. So, the idea of the method being proposed aims to join the efficiency of Newton‐Raphson‐like methods with the simplicity of methods based on tearing equations, what is possible only with a very good initial estimates generating methodology.

2. Model assumptions

The following assumptions were made when formulating the model of the distillation process:

• Steady state;

• No reaction occurs in the column;

• The vapor and liquid phases are homogeneous in all stages;

• The vapor and liquid leaving any stage are in phase equilibrium;

• Heat transfer only on condenser and reboiler, unless otherwise specified;

• The model does not include effects due to column internals (e.g., pressure drops and flooding/weeping).

3. Modeling

The modeling of a steady‐state distillation column is based mainly on mass and energy balances; in this way, it is needed to understand the equipment layout to obtain such mathematical equations [8, 12, 13]. This model is based on the equations of the column called MESH (material balance equations, phase equilibrium equations, mole fractions summation equations and heat, which means energy balance equations). Aiming to make the model as general as possible, it will be considered that can exist feed stream in any stage and the output streams (bottoms and distillate) can be in liquid phase, vapor phase or both phases [14]. Figure 1 shows a schematic representation of a distillation column and Figure 2 shows a schematic representation of input and output streams in a stage.

Figure 1.

Schematic representation of a distillation column.

Figure 2.

Schematic representation of input and output streams in a generic stage.

Where $F_j$is the flow rate of the feed stream to stage $j$, $Q_j$is the heat load from stage $j$(for convention, it was assumed the heat leaving the stage), $L_j$is the liquid flow rate outputting stage $j$and inputting stage $j+1$, $V_j$is the vapor flow rate outputting stage $j$and inputting stage $j-1$, $U_j$is the liquid side flow rate outputting stage $j$and $W_j$is the vapor side flow rate outputting stage $j$. The side streams $U$and $W$are used to represent the output streams like distillate and bottoms, that is, $D$and $B$. These streams can occur as liquid ($L_B$and $L_D$) or vapor ($V_B$and $V_D$), in which $D=L_D+V_D$and $B=L_B+V_B$. We are going to consider $n$stages, with the stage numbering starting on condenser ($j=1$) and continuing until reboiler is reached ($j=n$). The generalization from Figure 2 to represent all stage is shown in Figure 3.

Figure 3.

Schematic representation of input and output streams in all stages of a distillation column.

Considering the schematic representation of Figure 3, it is possible to consider side stream in any stage, not only in the first but also in the last stages. A little care is needed in utilizing Figure 3 as a base for obtaining the mass and energy balances, that is, how there is no stage $0$, so the streams $U_0$and $L_0$that would output and stream $V_1$that would input this stage do not exist; similarly, how there is no stage $n+1$, so streams $V_{n+1}$and $W_{n+1}$that would output and stream $L_n$that would input this stage do not exist. So, it is needed to fix $W_{n+1}=V_{n+1}=L_n=U_0=L_0=V_1=0$.

So, one can obtain the mass balance in stage $j$as

$F_j+V_{j+1}+L_{j-1}-V_j-W_j-L_j-U_j=0.$E1

The sum of mass balance in all stage gives us the global mass balance

${\displaystyle \sum _{j=1}^n{F}_j}-{\displaystyle \sum _{j=1}^n\left(W_j+U_j\right)}=0.$E2

Or putting separated the bottoms and distillate (side streams in condenser and reboiler)

${\displaystyle \sum _{j=1}^n{F}_j}-B-D-{\displaystyle \sum _{j=2}^{n-1}\left(W_j+U_j\right)}=0$E3

The mass balance of component $i$in stage $j$is given by

$z_{ij}{F}_j+y_{ij+1}{V}_{j+1}+x_{ij-1}{L}_{j-1}-y_{ij}\left(V_j+W_j\right)-x_{ij}\left(L_j+U_j\right)=0$E4

And the global mass balance of component $i$by

${\displaystyle \sum _{j=1}^n{z}_{ij}{F}_j}-{\displaystyle \sum _{j=1}^n\left(y_{ij}{W}_j+x_{ij}{U}_j\right)}=0.$E5

Or putting separated the bottoms and distillate

${\displaystyle \sum _{j=1}^n{z}_{ij}{F}_j}-x_{in}{L}_B-y_{in}{V}_B-x_{i1}{L}_D+y_{i1}{V}_D-{\displaystyle \sum _{j=2}^{n-1}\left(x_{ij}{W}_j+x_{ij}{U}_j\right)}=0$E6

where $x_{ij}$is the fraction of component $i$in the liquid phase of stage $j$, $y_{ij}$is the fraction of component $i$in the vapor phase of stage $j$and $z_{ij}$is the fraction of component $i$in the feed stream of stage $j$.

The energy balance in stage $j$, ignoring changes in kinetic and potential energies, is given by

$F_j{H}_j^F+V_{j+1}{H}_{j+1}^V+L_{j-1}{H}_{j-1}^L-\left(V_j+W_j\right)H_j^V-\left(L_j+U_j\right)H_j^L-Q_j=0$E7

where $H_j^F$is the enthalpy of the feed stream to stage $j$, $H_j^V$is the enthalpy of the vapor stream from stage $j$and $H_j^L$is the enthalpy of the liquid stream from stage $j$.

The global energy balance can be obtained by the sum of energy balance of all stages

${\displaystyle \sum _{j=1}^n{F}_j}{H}_j^F-{\displaystyle \sum _{j=1}^n\left(W_j{H}_j^V+U_j{H}_j^L+Q_j\right)}=0$E8

Or in a equivalent form

${\displaystyle \sum _{j=1}^n{F}_j}{H}_j^F-L_B{H}_n^L-V_B{H}_n^V-L_D{H}_1^L-V_D{H}_1^V-Q_C-Q_R-{\displaystyle \sum _{j=2}^{n-1}\left(W_j{H}_j^V+U_j{H}_j^L+Q_j\right)}=0$E9

Normally for the operation of a distillation column, are specified the reflux and the reboil ratios, respectively, defined by

$r_D=\frac{L_1}{D}$E10

$r_B=\frac{V_n}{B}$E11

To generalize this part will be defined $d_j$the ratio between the side stream and the stream outputting a stage and inputting the near stages, that is,

$d_j=\frac{R_j}{L_j+V_j}=\frac{U_j+W_j}{L_j+V_j}$E12

The side stream $R_j$is given by the sum of liquid side stream $U_j$and vapor side stream $W_j$, so it is also necessary to define the vapor fraction of side stream ${\omega }_j$, that is,

$W_j={\omega }_j{R}_j$E13

$U_j=\left(1-{\omega }_j\right)R_j$E14

How $V_1=0$, Eq. (12) is simplified to Eq. (10) by $d_1=1/r_D$and how $L_n=0$, Eq. (12) is simplified to Eq. (11) by $d_n=1/r_B$.

Talking about side streams, there are four possibilities: no side stream, only liquid side stream, only vapor side stream and both liquid and vapor side stream. More details are shown in Table 1 and Figure 4.

Figure 4.

Four types of side streams.

Stage type	No side stream	Only liquid side stream	Only vapor side stream	Vapor and liquid side streams
Generic (stage j)	$d_j=0$	$d_j\ne 0$	$d_j\ne 0$	$d_j\ne 0$
	$W_j=0$	${\omega }_j=0$	${\omega }_j=1$	$0<{\omega }_j<1$
	$U_j=0$	$R_j=d_j\left(L_j+V_j\right)$	$R_j=d_j\left(L_j+V_j\right)$	$R_j=d_j\left(L_j+V_j\right)$
		$W_j=0$	$W_j=R_j$	$W_j={\omega }_j{R}_j$
		$U_j=R_j$	$U_j=0$	$U_j=\left(1-{\omega }_j\right)R_j$
Condenser (j = 1)	Total reflux	Total condenser	Partial condenser	Partial condenser
	$d_1=0$	$d_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_D$}\right.$	$d_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_D$}\right.$	$d_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_D$}\right.$
	$W_1=V_D=0$	${\omega }_1={\omega }_D=0$	${\omega }_1={\omega }_D=1$	$0<{\omega }_1={\omega }_D<1$
	$U_1=L_D=0$	$R_1=d_1{L}_1=D$	$R_1=d_1{L}_1=D$	$R_1=d_1{L}_1=D$
		$W_1=V_D=0$	$W_1=V_D=R_1$	$W_1=V_D={\omega }_D{R}_1$
		$U_1=L_D=R_1$	$U_1=L_D=0$	$U_1=L_D=\left(1-{\omega }_D\right)R_1$
Reboiler (j = n)	Total reflux	Partial reboiler	Total reboiler	Partial reboiler
	$d_n=0$	$d_n=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_B$}\right.$	$d_n=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_B$}\right.$	$d_n=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r_B$}\right.$
	$W_n=V_B=0$	${\omega }_n={\omega }_B=0$	${\omega }_n={\omega }_B=1$	$0<{\omega }_n={\omega }_B<1$
	$U_n=L_B=0$	$R_n=d_n{V}_n=B$	$R_n=d_n{V}_n=B$	$R_n=d_n{V}_n=B$
		$W_n=V_B=0$	$W_n=V_B=R_n$	$W_n=V_B={\omega }_B{R}_n$
		$U_n=L_B=R_n$	$U_n=L_B=0$	$U_n=L_B=\left(1-{\omega }_B\right)R_n$

Table 1.

Four types of side streams.

4. Thermodynamics

Thermodynamics plays a key role on the modeling of phase equilibrium. For the phases of liquid and vapor in thermodynamic equilibrium, the fraction of component in each phase is connected by the iso‐fugacity relation [15–19]

$f_i^L=f_i^V$E15

where $f_i^V$is the fugacity of component $i$in vapor phase and $f_i^L$is the fugacity of component $i$in liquid phase, which leads to the mathematical relation

$y_i=K_i{x}_i$E16

where the way for calculating $K_i$depends on the chosen vapor liquid equilibrium (VLE) formulation, which can be Phi/Phi or Gamma/Phi [15–19], respectively,

$K_i=\frac{{\hat{\varphi }}_i^L}{{\widehat{\varphi }}_i^V}$E17

$K_i=\frac{{\gamma }_i{\varphi }_i^{sat}{P}_i^{sat}}{{\widehat{\varphi }}_i^V}\exp \left({\displaystyle {\int }_{P_i^{sat}}^P\frac{V_i^L}{RT}}dP\right)$E18

where ${\widehat{\varphi }}_i^L$is the fugacity coefficient of component $i$in liquid phase, ${\widehat{\varphi }}_i^V$is the fugacity coefficient of component $i$in vapor phase, ${\varphi }_i^{sat}$is the fugacity coefficient of pure component $i$in saturation state, ${\gamma }_i$is the activity coefficient of component $i$in liquid phase, $P$is the pressure, $P^{sat}$is the saturated vapor pressure, $R$is the ideal gas constant, $V_i^L$is the component $i$volume in liquid phase and $T$is the temperature.

For the calculation of fugacity coefficient is required to use an equation of state (EOS) and for the calculation of activity coefficient is required to use a model that represents the excess Gibbs free energy. The saturated vapor pressure is calculated using one of the equation that describes the relation between vapor pressure and temperature for pure components, like Antoine, Wagner, Riedel, Harlecher‐Braun, among others [15–19]. These equations are based on Clapeyron equation and their constant is obtained by experimental data fitting.

Another very important calculation that is needed to resort thermodynamics is the enthalpy. The liquid and vapor phases enthalpies are calculated [15–19], respectively, by

$H^L={\displaystyle \sum _{i=1}^m{x}_i{H}_i^L}+H^E$E19

$H^V={\displaystyle \sum _{i=1}^m{y}_i{H}_i^V}+H^R$E20

where $m$is the number of components, $H_i$is the enthalpy of pure component $i$, $H^E$is the excess enthalpy that can be calculated by a model that represents the excess Gibbs free energy (EOS or activity coefficient) and $H^R$is the residual enthalpy calculated using an EOS.

For the calculation of pure component $i$enthalpy, it is necessary to specify a reference state. How we are considering that there is no reaction in the column and it is possible to put the enthalpies of all components equal to zero in a reference state (same conditions of temperature, pressure and phase). But, the intention here is to make the model as general as possible it will be chosen elemental reference state that can be used also in a reactive system. In this reference state (298K and 1 atm), the standard enthalpy of formation ($\text{Δ}H{°}_f$), by convention, for an element in its most stable form standard state is zero [15–19]. The calculation of the enthalpy of a component is based on the elemental reference state, shown in Figure 5.

Figure 5.

Calculation of the enthalpy of a component based on the elemental reference state.

Where $C_P$is the heat capacity at constant pressure of liquid (L) and ideal gas (IG), $T_{bp}$is the bubble point temperature, $T_{dp}$is the dew point temperature and $\text{Δ}{H}^{VAP}$is the vaporization enthalpy.

5. Materials and methods

In this section, will be presented the methods for generating the initial estimates and for solution of equations showed above.

5.1. Initial estimates

The resolution of equations that models the steady‐state equilibrium distillation column involves a set of highly nonlinear equation, mainly on phase equilibria and energy balances. The algorithms, for solution of this type of problem, request good initial estimates in order that it can be possible to reach a solution. Moreover, here we are working in an algorithm that avoids the use of methods for solving nonlinear equations systems, what makes the quality of initial estimate even more important.

There are a lot of methods for solving models of steady‐state equilibrium distillation column, considering various levels of layout complexity, number of components involved and accuracy of properties calculation. Simpler models do not need initial estimates, but for more complex models, a good initial estimate is fundamental.

McCabe‐Thiele is a graphical method for combining the equilibrium curve with mass balance, assuming that there are two sections in the distillation column (between reboiler and feed stage and between feed stage and condenser) where molar vapor and liquid flow rates are constant, in addition to the assumption that there is no heat loss, eliminates the need of energy balances [4], something like the non‐heat effect presented below. Ponchon‐Savarit is a graphical method that includes energy balances, utilizing for this an enthalpy‐concentration diagram [5, 6]. How, Ponchon‐Savarit method utilizes energy balances and it is more accurate than McCabe‐Thiele method. These methods do not need initial estimates but, unfortunately, are applicable only for distillation of binary mixtures.

For complex systems, it is suggested that the procedure for solving this models should be based on the solution of a system of nonlinear equations using an appropriated method for solving systems of nonlinear equations like Newton‐Raphson. This system of equations is composed by MESH equations or combinations of them. But, the solution convergence of this type of problem is totally dependent on the quality of initial guess.

There are a lot of methods that use a technique called “tearing equations” that split large and sparse systems of algebraic equations into smaller system [20]. They are relatively simple, but are restricted to ideal and nearly ideal mixtures. The methods of Lewis‐Matheson [9], Thiele‐Geddes [10] and theta [8] are based on equation tearing for solving simple distillation columns with one feed and two product stream. The bubble‐point method receives this name because it tears the MESH equations in a way that a new set of stage temperatures is computed from bubble‐point equations [11]. Similarly, the sum‐rates method calculates, at each new iteration, the values of liquid streams by the summation of components flow rates in liquid phase [7].

The MESH equations wrote in this work will be rearranged for using tearing‐equation method like bubble‐point method, as can be seen in the next section. So, the idea here is to propose a method for generating good initial guesses aiming to avoid instability, normally presented for simpler method. That is, we are trying to join the efficiency of Newton‐Raphson‐like methods with the simplicity of methods based on tearing equations. It is already demonstrated that, with good initial guesses, it is possible to use tearing equation for very complex models like steady state for reactive distillation columns [13].

The algorithm at issue needs initial estimates for temperature, liquid and vapor streams and side streams. For initial estimates of temperature will be considered a model based on saturation temperature ($T^{SAT}$) of pure components at column operation pressure ($P_{col}$) and liquid, vapor and side streams will be estimated considering a non‐heat effect.

For the initial estimates of stage temperatures ($T^{\left(0\right)}$) will be considered a linear profile, in which, approximately in the middle of the column, we have an average temperature ($T_{ave}$) pondered on component fractions in feed streams

$T_{ave}=\frac{{\displaystyle \sum _{i=1}^m\left(T_i^{SAT}{\displaystyle \sum _{j=1}^n{z}_{ij}{F}_j}\right)}}{{\displaystyle \sum _{j=1}^n{F}_j}}$E21

And a minimal temperature equal to $T_{ave}$minus the average absolute deviation (AAD)

$T_{\min }=T_{ave}-AAD=T_{ave}-\frac{{\displaystyle \sum _{i=1}^m\left(\left|T_i^{SAT}-T_{ave}\right|{\displaystyle \sum _{j=1}^n{z}_{ij}{F}_j}\right)}}{{\displaystyle \sum _{j=1}^n{F}_j}}$E22

That is,

$T_j^{(0)}=T_{\min }+\frac{2\left(j-1\right)}{n}\left(T_{ave}-T_{\min }\right)$E23

For non‐heat effects model is made the assumption of constant vaporization enthalpy, what can guarantee that the change in the liquid and vapor flow rates in each section of the column is due to only the feed quality, feed flow rate and side flow rates [21], as illustrated in Figure 6. Where the feed quality is defined as

• Saturated liquid;

  $q=1$E24

• Saturated vapor;

  $q=0$E25

• Saturated liquid and vapor, $q$is equal to liquid fraction;

  $0<q<1$E26

• Subcooled liquid, $q>1$;

  $q=\left(\Delta H^{VAP}+{\displaystyle {\int }_{T_{bp}}^T{C}_P^{L}dT}\right)/\Delta H^{VAP}$E27

• Superheated vapor, $q<0$.

  $q=-\left({\displaystyle {\int }_{T_{dp}}^T{C}_P^{V}dT}\right)/\Delta H^{VAP}$E28

Figure 6.

Non‐heat effect behavior in a stage for: (a) saturated liquid; (b) saturated vapor; (c) saturated liquid and vapor; (d) subcooled liquid and (e) superheated vapor.

Using the mass balance in a stage and the global mass balance, one can find the expression

$R_j=\frac{\frac{d_j}{d_j+1}\left(F_j+V_{j+1}+L_{j-1}\right)}{{\displaystyle \sum _{j=1}^n\frac{d_j}{d_j+1}\left(F_j+V_{j+1}+L_{j-1}\right)}}{\displaystyle \sum _{j=1}^n{F}_j}$E29

At this point, we want to know just a rude value of side flow rates. So, let us make the rude assumption that feed and side streams are 50% in each phase and assume that the non‐heat effects are valid. In this case, the expression

$\frac{1}{d_j+1}\left(F_j+V_{j+1}+L_{j-1}\right)$E30

is almost constant along the column. And, it can be used to generate initial estimates for side flow rates

$R_j^{\left(0\right)}=\frac{d_j}{{\displaystyle \sum _{j=1}^n{d}_j}}{\displaystyle \sum _{j=1}^n{F}_j}$E31

If there are side streams only in the reboiler and in the condenser or the side streams in the intermediate stages are known, the non‐heat effects model can be used to generate initial estimates directly to bottoms and distillate flow rates

$B^{\left(0\right)}=\frac{r_D{F}_1+\left(r_D+1\right)F_n+{\displaystyle \sum _{j=2}^{n-1}\left[\left(q_j+r_D\right)F_j-\left(r_D+1\right)U_j-W_j\right]}}{r_D+r_B+1}$E32

$D^{\left(0\right)}=\frac{\left(r_B+1\right)F_1+r_D{F}_n+{\displaystyle \sum _{j=2}^{n-1}\left[\left(r_B+1-q_j\right)F_j+\left(r_D+1\right)U_j+W_j\right]}}{r_D+r_B+1}$E33

So $V_n$and $L_1$can be estimated using, respectively, the reboil and reflux ratio

$V_n^{\left(0\right)}=r_B{B}^{\left(0\right)}$E34

$L_1^{\left(0\right)}=r_D{D}^{\left(0\right)}$E35

Using the non‐heat effects model, it also can be obtained the initial estimates for vapor flow rates, by mass balances in vapor phase, for $j=n-1,\mathrm{\dots },2$

$V_j^{\left(0\right)}=\frac{\left(1-q_j\right)F_j+V_{j+1}^{\left(0\right)}}{d_j+1}$E36

The liquid flow rates are estimated using the mass balance by stage for $j=2,\mathrm{\dots },n-1$

$L_j^{\left(0\right)}=\frac{F_j+V_{j+1}^{\left(0\right)}+L_{j-1}^{\left(0\right)}}{d_j+1}-V_j^{\left(0\right)}$E37

With the initial estimates for side, liquid and vapor stream and temperature, it can be started the iterative process of steady‐state model solution. The only unknowns lasting, to generate initial estimates, are the components fractions in each stream, but, in the algorithm, we are going to work; these initial estimates, showed how to generate, are sufficient for the calculation of components fractions. So, for these unknowns, it is not necessary to generate initial estimates.

5.2. Algorithm

When one talks about the solution of algebraic equation, it is necessary to keep in mind that for the calculation of a number unknowns is necessary the same number of equations. In a complex system of equations, it is not easy to be sure that we have the right number of equations and unknowns. To make it easier to do, this balance is presented in Table 2.

Unknown type	Number of unknowns	Stage	Equation(s) used to calculate
Temperatures	$n$	$1$to $n$	The restriction that the sum of vapor mole fraction is equal to unity, for stages $1$to $n$
Vapor flow rates	$n-1$	$2$to $n$	Energy balance on stages to $n-1$and mass balance on stage $n$
Liquid flow rates	$n-1$	$1$to $n-1$	Mass balance for stage $1$to $n-1$
Side flow rates (liquid and vapor)	$2n$	$1$to $n$	Ratio between side streams and streams outputting stages for $n$stages. And vapor fraction of side stream for $n$stages
Liquid mole fractions	$nm$	$1$to $n$	Mass balance by component in $n$stages for $m$components
Vapor mole fractions	$nm$	$1$to $n$	Iso‐fugacity relation for $m$components in $n$stages
Heat transfer	$2$	$1$and $n$	Energy balance in $1$stages and $n$
Total	$\left(2m+5\right)n$	$n$	$\left(2m+5\right)n$

Table 2.

Balance of equations and unknowns.

The first step of algorithm is the imputation data. The data needed to be imputed are as follows:

• Number of components, $m$;

• Number of stages, $n$;

• The side streams ratios, $d_j$for $j=1,\mathrm{\dots },n$($d_n=1/r_B$and $d_1=1/r_D$);

• The vapor fraction of side stream, ${\omega }_j$for $j=1,\mathrm{\dots },n$;

• Column pressure, $P_{col}$;

• The heat transfer for intermediate stages, $Q_j$for $j=2,\mathrm{\dots },n-1$(they are normally adiabatic stages, that is, $Q_j=0$);

• Feed flow rates, $F_j$for $j=1,\mathrm{\dots },n$;

• Feed temperature, $T_j^{Feed}$for $j=1,\mathrm{\dots },n$;

• Feed quality, $q_j$for $j=1,\mathrm{\dots },n$;

• Feed fractions, $z_{ij}$for $j=1,\mathrm{\dots },n$and $i=1,\mathrm{\dots },m$;

• Physical, critical and other properties of components, for evaluation of enthalpies and phase equilibrium.

At this point, a more watchful reader must be thinking: it was used the restriction of vapor mole fraction sum, but was not used the same restriction for the liquid mole fraction, that is, right? Actually, this restriction is implicitly used, because the mass balance of a stage is the sum of mass balances by component in that stage and the sum of liquid mole fractions is forced to be equal to unity by a normalization step used in the algorithm presented ahead.

The sequence of calculation of the algorithm is presented in Figure 7.

Figure 7.

Algorithm for the simulation of steady‐state distillation columns.

Using the mass balance by stage and the phase equilibrium relation

$y_{ij}=K_{ij}{x}_{ij}$E38

One obtain

$l_j{x}_{ij-1}+m_j{x}_{ij}+u_j{x}_{ij}=b_j$E39

where

$l_j=L_{j-1}$E40

$a_j=-\left[\left(L_j+U_j\right)+\left(V_j+W_j\right)K_{ij}\right]$E41

$u_j=V_{j+1}{K}_{ij+1}$E42

$b_j=-F_j{z}_{ij}$E43

Remembering that the $L_0$, $L_n$, $V_1$and $V_{n+1}$streams do not exist, it is obtained a tri‐diagonal linear system of equations for each component

$\left[\begin{array}{ccccc}{a}_1& u_1& 0& 0& 0\\ l_2& a_2& u_2& 0& 0\\ 0& \ddots & \ddots & \ddots & 0\\ 0& 0& l_{n-1}& a_{n-1}& u_{n-1}\\ 0& 0& 0& l_n& a_n\end{array}\right]\cdot \left[\begin{array}{c}{x}_{i1}\\ x_{i2}\\ ⋮\\ x_{in}\\ x_{in}\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\\ b_n\end{array}\right]$E44

Because the initial guess imprecision (mainly on temperatures) and the mole fractions of each component are calculated separately, especially in the first iterations, they may have values without any physical meaning (such as, negative or a sum different from unity,). So that, the convergence process of algorithm may be accelerated by normalization, undertaken by the following equation

$x_{ij}=\frac{\left|x_{ij}\right|}{{\displaystyle \sum _{i=1}^m\left|x_{ij}\right|}}$E45

The temperatures are the only set of unknowns that cannot be calculated by a linear method. The temperatures are calculated stage by stage using the restriction that the sum of mole fractions in vapor phase must be equal to unit with aid of phase equilibrium relation, that is,

$F\left(T_j\right)=1-{\displaystyle \sum _{i=1}^m{K}_{ij}{x}_{ij}}=0$E46

The solution of this step is made with aid of a iterative root‐finding method, like secant

$T_j^{k+1}=T_j^k-F\left(T_j^k\right)\left[\frac{T_j^k-T_j^{k-1}}{F\left(T_j^k\right)-F\left(T_j^{k-1}\right)}\right]$E47

With the values of temperature and liquid mole fractions, it can be calculated the vapor mole fractions simply by using the phase equilibrium relation, that is, Eq. (38).

The enthalpies of the all streams are calculated as previously shown. These enthalpies are used in the energy balances. The energy balances in stages $1$and $n$are used to calculate the heat duties on these stages. In the others stages, the heat transfer must be zero (adiabatic stages) or must be fixed in the input data.

$Q_C=Q_1=F_1{H}_1^F+V_2{H}_2^V-W_1{H}_1^V-\left(L_1+U_1\right)H_1^L$E48

$Q_R=Q_n=F_n{H}_n^F+L_{n-1}{H}_{n-1}^L-\left(V_n+W_n\right)H_n^V-U_n{H}_n^L$E49

The value of vapor flow rate from stage $n$is obtained by using the mass balance in this stage

$V_n=F_n+L_{n-1}-U_n-W_n$E50

And the values of vapor flow rate from stages $2,\mathrm{\dots },n-1$are calculated using energy balance with the aid of mass balances to eliminate the liquid flow rates terms. The calculation is made from $n-1$back to $2$

$V_j=\frac{{\beta }_j{V}_{j+1}+{\gamma }_j}{{\alpha }_j}$E51

where,

${\gamma }_j=\left(H_j^F-H_j^L\right)F_j+\left(H_{j-1}^L-H_j^L\right){\displaystyle \sum _{k=1}^{j-1}\left(F_k-U_k-W_k\right)}+\left(H_j^L-H_j^V\right)W_j-Q_j$E52

${\beta }_j=H_{j+1}^V-H_j^L$E53

${\alpha }_j=H_j^V-H_{j-1}^L$E54

The energy balance normally is used to calculate the temperatures, but it would be necessary a solution of a system of highly nonlinear equations. Instead of this, the energy balances are used to calculate the vapor flow rates, with aid of mass balances and some algebraic rearrangement, by sequential evaluations (one vapor flow rate at a time).

The side flow rates are calculated by substituting the ratio between the side stream and the stream outputting a stage in the mass balance

$R_j=\left(\frac{d_j}{d_j+1}\right)\left(F_j+V_{j+1}+L_{j-1}\right)$E55

$W_j={\omega }_j{R}_j$E56

$U_j=\left(1-{\omega }_j\right)R_j$E57

The liquid flow rates are calculated sequentially from stage 1 to stage $n-1$(remember, there is no $L_n$) by using the mass balance in the respective stages

$L_j=F_j+V_{j+1}-V_j+L_{j-1}-R_j$E58

The iterative process finishes when the relative variation of some main unknowns of the model is very low. The smallness of the variation depends on the chosen tolerance

${\displaystyle \sum _{j=1}^n{\left[\frac{T_j^{(k)}-T_j^{(k-1)}}{T_j^{(k)}}\right]}^2}+{\displaystyle \sum _{j=1}^{n-1}{\left[\frac{L_j^{(k)}-L_j^{(k-1)}}{L_j^{(k)}}\right]}^2}+{\displaystyle \sum _{j=2}^n{\left[\frac{V_j^{(k)}-V_j^{(k-1)}}{V_j^{(k)}}\right]}^2}\le tolerance$E59

6. Results and discussion: case study

For the evaluation of the initial estimates and algorithm in question, it will be tested with an example of hydrocarbons separation. Where a feed stream contains four hydrocarbons: propane (C₃), n‐butane(n‐C₄), isopentane (i‐C₅) and n‐pentane (n‐C₅). The operational conditions are presented in Table 3.

Variables	Specifications
Pressure	All stages	13.8 bar
Ratio	Reflux	$r_D=$5.0
	Reboil	$r_B=$3.2531
Stages	Number	12
	Type	1 condenser; 10 adiabatic stages ($Q_j=0$); 1 reboiler
Feed	Location	Stage 6
	Condition	Saturated liquid
	Flow rate	100 kmol/h
	Temperature	359.3 K
	Mole fractions	C3 (0.4)
		n‐C4 (0.4)
		i‐C5 (0.1)
		n‐C5 (0.1)
Side streams	Condenser (stage 1)	$d_1=1/r_D$
		${\omega }_1={\omega }_D=0$
		$D=L_D$
	Stages 2–11	$d_j=0$
		$U_j=W_j=0$
	Reboiler (stage 12)	$d_{12}=1/r_B$
		${\omega }_{12}={\omega }_B=0$
		$B=L_B$

Table 3.

Operational conditions of distillation column.

For representing the nonideal behavior of vapor phase, it was used Peng‐Robinson equation of state, where were considered as mixing rule binary iteration parameters with geometric mean for parameter $a$and arithmetic mean for parameter $b$. And for representing the liquid non‐ideality, it was utilized UNIFAC method. The UNIFAC parameters, physical and critical properties for the components involved in the case of study, were taken from [16, 22], in which physical and critical properties are presented in Table 4.

$C_P=C_1+C_{2}T+C_3{T}^2+C_4{T}^3$E60

$\ln \left(P^{sat}\right)=A-\frac{B}{T+C}$E61

Properties		Propane	n‐Butane	isoPentane	n‐Pentane
Properties for calculation of fugacity coefficients	Tc (K)	369.8	425.2	460.4	469.7
	Pc (bar)	42.5	38.0	33.9	33.7
	ω	0.153	0.199	0.227	0.251
	Tbp (K)	231.1	272.7	301.0	309.2
	V^L (m^3/kmol)	0.0758	0.1004	0.1164	0.1152
$C_P^{IG}$(kJ/(kmol K)) T (K)	C1	-4.224	9.487	-9.525	-3.626
	C2	3.063	3.313	5.066	4.873
	C3	-1.586	-1.108	-2.729	-2.580
	C4	3.215	-0.2822	5.723	5.305
$C_P^L$(kJ/(kmol K)) T (K)	C1	87.19	110.4	115.6	133.7
	C2	-0.6110	-0.6567	-0.1138	-0.7024
	C3	4.128	5.016	5.773	5.904
$P^{sat}$T (K) and P (bar)	A	9.1058	9.0580	9.0136	9.2173
	B	1872.46	2154.90	2348.67	2477.07
	C	-25.16	-34.42	-40.05	-39.94

Table 4.

Physical and critical properties for the components involved in the case study.

For a tolerance of $1\times {10}^{-10}$, it was needed 28 iterations to achieve convergence. A graphic of error versus iteration is shown in Figure 8. In the start of the procedure, few steps present a high error reduction, after this a linear convergence rate is obtained. That is not as good as a quadratic convergence of an algorithm that uses a method for solution of nonlinear equations, like Newton‐Raphson, but is justified by the simplicity of the solution algorithm.

Figure 8.

Numerical convergence.

Figure 9 shows a comparison between initial guess and simulation results for stage temperatures. Initial guess profile is extremely close to results.

Figure 9.

Temperature profiles.

Figure 10 shows comparisons between initial guess and simulation results for molar flow rates of the liquid and vapor phases. Initial guess of rates of both phases was extremely close to the values calculated by simulation.

Figure 10.

Liquid and vapor flow rates profiles.

The liquid mole fraction profiles calculated are shown in Figure 11. Some important numeric results of the simulation are shown in Table 5.

Figure 11.

Liquid mole fractions profiles.

Variables	Specification	Value
		This work	King [23]
Flow rate (kmol/h)	D	40.000	40.0
	B	60.000	60.0
Heat duty (kJ/h)	QC	4.1133 × 10^6	–
	QR	-4.0774 × 10^6	–
Temperature (K)	Condenser	315.91	316.15
	Reboiler	376.10	377.15
Distillate mole fractions	C3	0.92966	0.980
	n‐C4	0.06936	0.020
	i‐C5	0.00070	0.000
	n‐C5	0.00029	0.000
Bottoms mole fractions	C3	0.04688	0.013
	n‐C4	0.62044	0.653
	i‐C5	0.16621	0.167
	n‐C5	0.16648	0.167

Table 5.

Some simulation results.

In this case of study, the initial estimates generated are no more than 20% far from the final result, what confirm the goodness of the methodology used for generating the guesses. The results obtained here are very close to that obtained by King [23], some differences can be justified by different levels of accuracy of the thermodynamic modeling (the thermodynamic modeling of the cited reference is simpler).

This case being studied aims to separate the propane from the other three hydrocarbons. For evaluating the influence of feed stage is presented in Figure 12, which presents the fraction of propane in distillate stream in function of the stage where the feed stream occurs. One can see in this figure that the best stage to put the feed stream is in the middle of the column. It is easy to understand, if the feed occurs near of condenser, there will be a great amount of components other than propane in the stages near of condenser, so part of it eventually outputs the column in the distillate stream and if the feed occurs near of reboiler, there will be a great amount of propane in the stages near of reboiler, so part of propane is present on bottoms stream.

Figure 12.

Fraction of propane in distillate stream in function of feed stage.

The number of stage also has a great influence on the fraction of propane in distillate stream. It is obvious that, the more stages there are in the column, a greater mole fraction of propane there will be in distillate. But, it is possible to see in Figure 13 that after a certain number of stages, the increase on that fraction is too small. Figure 12 confirms what was verified in Figure 13, because there is a comparison between column with the same number of stage, one with feed stream in stage six and another in the middle of column. And again the feed in the middle of column presents a better separation.

Figure 13.

Fraction of propane in distillate stream in function of the number of stages.

Another parameter that has a strong influence on the mole fraction of propane in distillate stream is the reflux ratio, see Figure 14. Higher purity levels are achieved by increasing the reflux ratio, with a clear limit. But greater reflux ratio greater is the heat load to the reboiler, see Figure 15. After a certain value of reflux ratio, the increase in propane mole fraction is very small.

Figure 14.

Fraction of propane in distillate stream in function of the reflux ratio.

Figure 15.

Heat load to reboiler in function of the reflux ratio.

With everything that has been exposed, one can see that achieving a better separation level is result mainly by spending more in operating or investment costs (greater reflux ratio, more stages, etc.) or changing the layout configuration of the operating system (feed stage location).

7. Conclusion

An algorithm was provided in this work for the solution of a distillation column operating in steady state; in this algorithm, the high nonlinear equations are solved in a very simple form. Equations in the model were divided into sets and each set was solved separately. The solution procedure uses an algorithm for solution of systems of tri‐diagonal linear equation, explicit calculations and a method for root finding of equations of one unknown variable. A methodology was also provided to produce initial guess which constitutes a critical step in the solution of nonlinear equations system. The modeling allows a variety of cases depending on the types of condenser and reboiler, number and conditions of feed stream, side streams, etc.

The suggested methodology for the production of initial estimates was efficient with values close to those calculated by simulation. This fact accelerates and increases the convergence warranty.