Open access peer-reviewed chapter

Steady‐State Modeling of Equilibrium Distillation

Written By

Vilmar Steffen and Edson Antonio da Silva

Submitted: April 14th, 2016 Reviewed: November 10th, 2016 Published: June 28th, 2017

DOI: 10.5772/66833

Chapter metrics overview

2,858 Chapter Downloads

View Full Metrics


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.


  • distillation
  • steady state
  • phase equilibria
  • modeling
  • simulation

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 [711]; 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 Fjis the flow rate of the feed stream to stage j, Qjis the heat load from stage j(for convention, it was assumed the heat leaving the stage), Ljis the liquid flow rate outputting stage jand inputting stage j+1, Vjis the vapor flow rate outputting stage jand inputting stage j1, Ujis the liquid side flow rate outputting stage jand Wjis the vapor side flow rate outputting stage j. The side streams Uand Ware used to represent the output streams like distillate and bottoms, that is, Dand B. These streams can occur as liquid (LBand LD) or vapor (VBand VD), in which D=LD+VDand B=LB+VB. We are going to consider nstages, 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 U0and L0that would output and stream V1that would input this stage do not exist; similarly, how there is no stage n+1, so streams Vn+1and Wn+1that would output and stream Lnthat would input this stage do not exist. So, it is needed to fix Wn+1=Vn+1=Ln=U0=L0=V1=0.

So, one can obtain the mass balance in stage jas


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


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


The mass balance of component iin stage jis given by


And the global mass balance of component iby


Or putting separated the bottoms and distillate


where xijis the fraction of component iin the liquid phase of stage j, yijis the fraction of component iin the vapor phase of stage jand zijis the fraction of component iin the feed stream of stage j.

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


where HjFis the enthalpy of the feed stream to stage j, HjVis the enthalpy of the vapor stream from stage jand HjLis 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 djthe ratio between the side stream and the stream outputting a stage and inputting the near stages, that is,


The side stream Rjis given by the sum of liquid side stream Ujand vapor side stream Wj, so it is also necessary to define the vapor fraction of side stream ωj, that is,


How V1=0, Eq. (12) is simplified to Eq. (10) by d1=1/rDand how Ln=0, Eq. (12) is simplified to Eq. (11) by dn=1/rB.

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 typeNo side streamOnly liquid side streamOnly vapor side streamVapor and liquid side streams
Generic (stage j)dj=0dj0dj0dj0
Condenser (j = 1)Total refluxTotal condenserPartial condenserPartial condenser
Reboiler (j = n)Total refluxPartial reboilerTotal reboilerPartial reboiler

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 [1519]


where fiVis the fugacity of component iin vapor phase and fiLis the fugacity of component iin liquid phase, which leads to the mathematical relation


where the way for calculating Kidepends on the chosen vapor liquid equilibrium (VLE) formulation, which can be Phi/Phi or Gamma/Phi [1519], respectively,


where ϕ^iLis the fugacity coefficient of component iin liquid phase, ϕ^iVis the fugacity coefficient of component iin vapor phase, ϕisatis the fugacity coefficient of pure component iin saturation state, γiis the activity coefficient of component iin liquid phase, Pis the pressure, Psatis the saturated vapor pressure, Ris the ideal gas constant, ViLis the component ivolume in liquid phase and Tis 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 [1519]. 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 [1519], respectively, by


where mis the number of components, Hiis the enthalpy of pure component i, HEis the excess enthalpy that can be calculated by a model that represents the excess Gibbs free energy (EOS or activity coefficient) and HRis the residual enthalpy calculated using an EOS.

For the calculation of pure component ienthalpy, 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 [1519]. 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 CPis the heat capacity at constant pressure of liquid (L) and ideal gas (IG), Tbpis the bubble point temperature, Tdpis the dew point temperature and ΔHVAPis 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 (TSAT) of pure components at column operation pressure (Pcol) 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 (Tave) pondered on component fractions in feed streams


And a minimal temperature equal to Taveminus 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 liquid;


  • Saturated vapor;


  • Saturated liquid and vapor, qis equal to liquid fraction;


  • Subcooled liquid, q>1;


  • Superheated vapor, q<0.


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


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 Vnand L1can 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=n1,,2


The liquid flow rates are estimated using the mass balance by stage for j=2,,n1


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 typeNumber of unknownsStageEquation(s) used to calculate
Temperaturesn1to nThe restriction that the sum of vapor mole fraction is equal to unity, for stages 1to n
Vapor flow ratesn12to nEnergy balance on stages to n1and mass balance on stage n
Liquid flow ratesn11to n1Mass balance for stage 1to n1
Side flow rates (liquid and vapor)2n1to nRatio between side streams and streams outputting stages for nstages. And vapor fraction of side stream for nstages
Liquid mole fractionsnm1to nMass balance by component in nstages for mcomponents
Vapor mole fractionsnm1to nIso‐fugacity relation for mcomponents in nstages
Heat transfer21and nEnergy balance in 1stages and 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, djfor j=1,,n(dn=1/rBand d1=1/rD);

  • The vapor fraction of side stream, ωjfor j=1,,n;

  • Column pressure, Pcol;

  • The heat transfer for intermediate stages, Qjfor j=2,,n1(they are normally adiabatic stages, that is, Qj=0);

  • Feed flow rates, Fjfor j=1,,n;

  • Feed temperature, TjFeedfor j=1,,n;

  • Feed quality, qjfor j=1,,n;

  • Feed fractions, zijfor j=1,,nand 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.

Figure 7.

Algorithm for the simulation of steady‐state distillation columns.

Using the mass balance by stage and the phase equilibrium relation


One obtain




Remembering that the L0, Ln, V1and Vn+1streams do not exist, it is obtained a tri‐diagonal linear system of equations for each component


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


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 1and nare 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 nis obtained by using the mass balance in this stage


And the values of vapor flow rate from stages 2,,n1are calculated using energy balance with the aid of mass balances to eliminate the liquid flow rates terms. The calculation is made from n1back to 2




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 n1(remember, there is no Ln) 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


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 (C3), n‐butane(n‐C4), isopentane (i‐C5) and n‐pentane (n‐C5). The operational conditions are presented in Table 3.

PressureAll stages13.8 bar
Type1 condenser
10 adiabatic stages (Qj=0)
1 reboiler
FeedLocationStage 6
ConditionSaturated liquid
Flow rate100 kmol/h
Temperature359.3 K
Mole fractionsC3 (0.4)
n‐C4 (0.4)
i‐C5 (0.1)
n‐C5 (0.1)
Side streamsCondenser (stage 1)d1=1/rD
Stages 2–11dj=0
Reboiler (stage 12)d12=1/rB

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 aand 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.

Properties for calculation of fugacity coefficientsTc (K)369.8425.2460.4469.7
Pc (bar)42.538.033.933.7
Tbp (K)231.1272.7301.0309.2
VL (m3/kmol)0.07580.10040.11640.1152
CPIG(kJ/(kmol K)) T (K)C1-4.2249.487-9.525-3.626
CPL(kJ/(kmol K)) T (K)C187.19110.4115.6133.7
PsatT (K) and P (bar)A9.10589.05809.01369.2173

Table 4.

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

For a tolerance of 1×1010, 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.

This workKing [23]
Flow rate (kmol/h)D40.00040.0
Heat duty (kJ/h)QC4.1133 × 106
QR-4.0774 × 106
Temperature (K)Condenser315.91316.15
Distillate mole fractionsC30.929660.980
Bottoms mole fractionsC30.046880.013

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.


  1. 1. Marquini MF, Mariani DC, Meireles AJA, Santos OAA, Jorge LMM. Simulation and analysis of an industrial system of columns for ethanol distillation. Acta Scientiarum. Technology. 2007;29(1):23–28. doi:10.4025/actascitechnol.v29i1.81
  2. 2. Seader JD, Henley EJ, Roper DK. Separation process principles: chemical and biochemical operations. 3rd ed. New York: John Wiley & Sons; 2011. 821 p.
  3. 3. Corazza FC, Oliveira JV, Corazza ML. Application of a subdivision algorithm for solving nonlinear algebraic systems. Acta Scientiarum. Technology. 2008;30(1):27–38. doi:10.4025/actascitechnol.v30i1.3198
  4. 4. McCabe WL, Thiele EW. Graphical Design of Fractionating Columns. Industrial & Engineering Chemistry. 1925;17(6):605–611. doi:10.1021/ie50186a023
  5. 5. Ponchon M. Application of Graphs to the Study of Commercial Fractional Distillation. La Technique Moderne. 1921;13:20, 55.
  6. 6. Savarit R. Arts et Metiers. 1922:65, 142, 178, 241, 266, 307.
  7. 7. Burningham DW, Otto FD. Which computer design for absorbers. Hydrocarbon Processing. 1967;46(10):163–170.
  8. 8. Holland CD. Fundamentals of Multicomponent Distillation. New York: McGraw‐Hill; 1980. 626 p.
  9. 9. Lewis WK, Matheson GL. Studies in distillation design of rectifying columns for natural and refinery gasoline. Industrial & Engineering Chemistry. 1932;24:496–498.
  10. 10. Thiele EW, Geddes RL. Computation of distillation apparatus for hydrocarbon mixtures. Industrial & Engineering Chemistry. 1933;25(3):289–295. doi:10.1021/ie50279a011
  11. 11. Wang JC, Henke GE. Tridiagonal matrix for distillation. Hydrocarbon Process. 1966;45(8):155–163.
  12. 12. Seader JD, Henley EJ. Separations Process Principles. New York: John Wiley & Sons; 1998. 886 p.
  13. 13. Steffen V, Silva EA. Steady‐state modeling of reactive distillation columns. Acta Scientiarum. Technology. 2012;34(1):61–69. doi:10.4025/actascitechnol.v34i1.9535
  14. 14. Alfradique MF, Castier M. Modeling and simulation of reactive distillation columns using computer algebra. Computers and Chemical Engineering. 2005;29(9):1875–1884. doi:10.1016/j.compchemeng.2005.04.002
  15. 15. Sandler SI. Chemical, Biochemical and Engineering Thermodynamics. 4th ed. New York: John Wiley & Sons, Inc.; 2006. 945 p.
  16. 16. Poling BE, Prausnitz JM, O'Connell JP. The Properties of Gases and Liquids. 5th ed. New York: McGraw‐Hill; 2004.
  17. 17. Smith JM, Van Ness HC, Abbott MM. Introduction to Chemical Engineering Thermodynamics. 7th ed. New York: MacGraw Hill; 2005. 817 p.
  18. 18. Koretsky MD. Engineering and Chemical Thermodynamics. 2nd ed. New York: John Wiley & Sons, Inc.; 2013. 704 p.
  19. 19. Elliott JR, Lira CT. Introductory Chemical Engineering Thermodynamics. 2nd ed. New York: Prentice Hall; 2012. 912 p.
  20. 20. Krong G. The Piecewise Solution of Large‐Scale Systems Hardcover. London: MacDonald & Co; 1963. 1963 p.
  21. 21. Chen F, Huss RS, Malone MF, Doherty MF. Simulation of kinetic effects in reactive distillation. Computers and Chemical Engineering. 200;24(11):2457–2472. doi:10.1016/S0098‐1354(00)00609‐8
  22. 22. Yaws CL. Yaws’ Handbook of Thermodynamics and Physical Properties of Chemical Compounds. New York: Knovel; 2003. 823 p.
  23. 23. King CJ. Separation Processes. 2nd ed. New York: McGrawHill; 1980. 850 p.

Written By

Vilmar Steffen and Edson Antonio da Silva

Submitted: April 14th, 2016 Reviewed: November 10th, 2016 Published: June 28th, 2017