Steady‐State Modeling of Equilibrium Distillation

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.


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][8][9][10][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.
• 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).

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.
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 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.
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 So, one can obtain the mass balance in stage j as The sum of mass balance in all stage gives us the global mass balance  ðW j þ U j Þ ¼ 0: Or putting separated the bottoms and distillate (side streams in condenser and reboiler) The mass balance of component i in stage j is given by And the global mass balance of component i by Or putting separated the bottoms and distillate The energy balance in stage j, ignoring changes in kinetic and potential energies, is given by where H F j is the enthalpy of the feed stream to stage j, H V j is the enthalpy of the vapor stream from stage j and H L j 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 Or in a equivalent form Normally for the operation of a distillation column, are specified the reflux and the reboil ratios, respectively, defined by 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, 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 ω j , that is, Stage type No side stream Only liquid side stream Only vapor side stream Vapor and liquid side streams Condenser (j=1) Total reflux Total condenser Partial condenser Partial condenser Reboiler (j=n) Total reflux Partial reboiler Total reboiler Partial reboiler 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

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][16][17][18][19] f where f V i is the fugacity of component i in vapor phase and f L i is the fugacity of component i in liquid phase, which leads to the mathematical relation 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][16][17][18][19], respectively, 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][16][17][18][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][16][17][18][19], respectively, by 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 (ΔH f ), by convention, for an element in its most stable form standard state is zero [15][16][17][18][19]. The calculation of the enthalpy of a component is based on the elemental reference state, shown in Figure 5.
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 ΔH VAP is the vaporization enthalpy.

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

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 bubblepoint 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 ð0Þ ) 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 And a minimal temperature equal to T ave minus the average absolute deviation (AAD) That is, 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 vapor; • Saturated liquid and vapor, q is equal to liquid fraction; • Subcooled liquid, q > 1; • Superheated vapor, q < 0.
Using the mass balance in a stage and the global mass balance, one can find the expression 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 is almost constant along the column. And, it can be used to generate initial estimates for side flow rates 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 So V n and L 1 can be estimated using, respectively, the reboil and reflux ratio 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, …, 2 The liquid flow rates are estimated using the mass balance by stage for j ¼ 2, …, n−1 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.

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.
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, …, n (d n ¼ 1=r B and d 1 ¼ 1=r D ); • The vapor fraction of side stream, ω j for j ¼ 1, …, n; • Column pressure, P col ; • The heat transfer for intermediate stages, Q j for j ¼ 2, …, n−1(they are normally adiabatic stages, that is, Q j ¼ 0); • Feed flow rates, F j for j ¼ 1, …, n; • Feed temperature, T Feed j for j ¼ 1, …, n;  • Feed quality, q j for j ¼ 1, …, n; • Feed fractions, z i j for j ¼ 1, …, n and i ¼ 1, …, 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.
Using the mass balance by stage and the phase equilibrium relation One obtain 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 a 1 u 1 0 0 0 l 2 a 2 u 2 0 0 0 ⋱ ⋱ ⋱ 0 0 0 l n−1 a n−1 u n−1 0 0 0 l n a n 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 Figure 7. Algorithm for the simulation of steady-state distillation columns.
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, The solution of this step is made with aid of a iterative root-finding method, like secant 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.
The value of vapor flow rate from stage n is obtained by using the mass balance in this stage And the values of vapor flow rate from stages 2, …, 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 where, 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 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 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

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 3 ), n-butane(n-C 4 ), isopentane (i-C 5 ) and n-pentane (n-C 5 ). The operational conditions are presented in Table 3.
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 nonideality, 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.
For a tolerance of 1 · 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 9 shows a comparison between initial guess and simulation results for stage temperatures. Initial guess profile is extremely close to results. 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. Side streams Condenser (stage 1)  Table 4. Physical and critical properties for the components involved in the case study.    The liquid mole fraction profiles calculated are shown in Figure 11. Some important numeric results of the simulation are shown in Table 5.
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.
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.
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.
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).

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 Figure 15. Heat load to reboiler in function of the reflux ratio. 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.