SciELO - Scientific Electronic Library Online

vol.33 número2Use of wiener nonlinear MPC to control a CSTR with multiple steady stateOptimal solvent cycle design in supercritical fluid processes índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados



  • No hay articulos citadosCitado por SciELO

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google


Latin American applied research

versión impresa ISSN 0327-0793

Lat. Am. appl. res. v.33 n.2 Bahía Blanca abr./jun. 2003


Dynamic simulation of high-index models of batch distillation processes

E. F. Costa Jr.1, R. C. Vieira1, A. R. Secchi2 and E. C. Biscaia Jr.1

1 Programa de Engenharia Química - COPPE/UFRJ - Rio de Janeiro - RJ - Brazil
Caixa Postal 68.502 - CEP 21.945-970 - Rio de Janeiro - Brazil
2 Departamento de Engenharia Química - UFRGS - Porto Alegre - RS - Brazil

Abstract - The dynamic modeling of batch distillation columns frequently leads to a mixed system of differential and algebraic equations (DAEs) with differential index greater than one, and this particular feature has many implications on the resolution strategy adopted. As the number of stages and components can be arbitrarily high, those mathematical models can be large scale systems, and the analysis of the system prior to numerical resolution can be cumbersome. Additionally, the consistent initialization step can pose a nontrivial numerical task. For the numerical resolution of this model it is employed the computational code under development at PEQ/COPPE/UFRJ and DEQUI/UFRGS. This code employs Automatic Differentiation (AD) tools to perform index reduction and consistent initialization with minimum interference of the user. The resulting consistency system is solved and the numerical integration of the final index-one DAE is accomplished by means of the integration code DASSLC.

Keywords - Batch Distillation. Index Reduction. Differential-Algebraic Equations (DAEs). Automatic Differentiation.


   Distillation processes are among the most important technologies for the separation of multicomponent mixtures. Distillation can be done either in a continuous operation regime (where the mixture is continuously supplied and the products are continuously led off) or in a discontinuous regime - the so-called batch distillation process - where the apparatus is loaded with the mixture at start-up time and the different fractions are taken out one by one.
   Batch operation is advantageous if: (i) the required operating capacity of a proposed facility is too small to permit continuous operation at a practical rate - pumps, boilers, piping, instrumentation and other auxiliary equipment generally have a minimum industrial operating capacity; (ii) the operating requirements of a facility fluctuates widely in characteristics of feed material as well as processing rate - batch equipment usually has considerably more operating flexibility than continuous equipment. This is the reason for the predominance of batch equipment for multiple purpose solvent recovery or pilot plant applications (Henley and Seader, 1981).
   The description of a typical batch-distillation column is conveniently divided into two parts: (1) the start-up period and (2) the production period. The production period is that part of the process during which a product is withdrawn from the column. Necessary adjustments to bring the column to an operating condition that produces a distillate with a desired purity are made in the start-up period (Holland and Liapis, 1983).
   The modeling of a batch distillation process gives rise to systems of coupled differential and algebraic equations, mostly presenting high index. In this contribution, the numerical resolution of such mathematical models is addressed. In section §2, the mathematical model of the start-up phase of a batch distillation column is presented. In §3 and §4, aspects related to its numerical resolution are discussed, and the computational code under development at PEQ/COPPE, in collaboration with DEQUI/UFRGS, is briefly described. In §5, numerical results are presented to prove the adequacy of the proposed approach.


   In this contribution, it is simulated a batch distillation column with np stages and nc components. A simplified scheme of the equipment is shown in Fig. 1, where the stage number 0 represents the condenser and the stage np+1, the reboiler.
   The main assumptions of the model are: (i) the liquid and the vapour leaving a stage are in thermodynamic equilibrium; (ii) the vapour hold-up is negligible; (iii) there is no heat loss throughout the column; (iv) the liquid and vapour phases behave as ideal mixtures; (v) the pressure on the column is constant.
   Mass and energy balances for each stage are represented by the Eqn. 1 and 2. The Eqn. 3 and 4 represent the total number of moles of stage i and the chemical equilibrium relation. Equations 5 and 6 follow from the assumption of ideal behavior of liquid and vapour phases. The constraint on the molar fraction summation is enforced via Eqn. 7, and Eqn. 8 comes from the assumption of constant number of moles at each stage.

Figure 1. Scheme of the batch distillation column.









   In the previous equations, j= 1, ..., nc and i= 0, ..., np+1 except for Eqn. 8, for which i= 0, ..., np. Additionally, i represents stage, j represents component, yi,j represents molar fraction in vapour phase, hi and Hi are enthalpies in liquid and vapour phases, Vi and Li are the vapour and liquid fluxes leaving the ith stage, ni,j is the amount of the jth component in the ith stage, Ti is temperature and Qi is the heat removed from the ith stage. In the simulated column, all heat duties Qi are known except for the condenser heat duty, Q0. Hence, the state vector is [Li Vi Ni ni,,j yi,,j Ti Hi hi Q0]T.
The equilibrium constant and the enthalpies of each component, needed at Eqn. 4 to 6, are function of T according to the following relations:




   Numerical values for the constants employed in Eqn. 9 to 11 can be found in Holland and Liapis (1983, pp. 211-212), and are reproduced in Table 1. Additional constants needed for simulation of the model are the heat duty of the reboiler Qnp+1 = - 3.5 106 BTU/h and the heat removed from each stage, Qi = 0 (i= 1, ..., np).


   The mathematical model of the start-up of a batch distillation column is typically a high-index DAE system with nonlinear state equations and equilibrium constraints. The differential equations generally come from energy and mass balances in the equilibrium stages, and the algebraic equations are flux constraints and equilibrium relationships. Additionally, the resulting DAEs can easily become a large-scale system as the number of stages and components get high.
   Experience has proved that working directly with the differential-algebraic equations (DAEs) is advantageous for many reasons: (i) the systems are simpler and the variables have a physical significance; (ii) time consuming manipulations are avoided and fundamental information, that could be lost through differentiation, is preserved; (iii) constitutive relations, such as phase equilibrium relations, kinetic equations and equilibrium isotherms, are easily included and/or changed whenever necessary (Vieira and Biscaia Jr., 2001).

Table 1. Constants employed to compute equilibrium conditions.

   When solving DAEs one must be concerned about the index of the system and about the consistency of the initial conditions. The index of a DAE system is the number of times that all or part of the system must be differentiated with respect to time in order to convert it into an explicit set of first order ODEs. Index-1 systems may be solved with modified ODE methods, while high-index systems (systems with index 2 or greater) require specific methods. Generally, high-index systems are rewritten in an index-1 formulation and solved with standard integration codes like DASSL (Petzold, 1989) or DASSLC (Secchi, 1992), written in FORTRAN and C, respectively.
   During the index reduction, some extra algebraic equations are obtained, which generally correspond to derivatives of the original algebraic equations. Those hidden algebraic equations along with the original DAEs compose the extended system. Consistent initial values must satisfy not only the DAE system but also the underlying extended system (Brenan et al., 1989).
   Several rigorous methods for initialization have been proposed in the literature, and they all depend to some extent on the characterization of the DAE system. In other words, in order to perform the numerical integration of a DAE system with standard integration codes, it is necessary to: (i) perform the characterization of the DAE system, in order to determine its index and, if required, the set of equations to be differentiated; (ii) perform the differentiation; (iii) check for feasibility of the initial states arbitrarily set; (iv) solve the consistency system for a consistent initial state.
   The mathematical modeling of the batch distillation column gives rise to a differential algebraic (DAE) system: Eqn. 1 and 2 are differential equations, while the remaining equations are algebraic. The differential variables of this system are ni,j and hi and the algebraic variables are Ni, yi,j, Hi, Ti, Vi, Q0 and Li. Hence, the dimension of the system is (np+2)(2nc+6)-1. It can be noticed that the algebraic variable Q0 does not appear in any algebraic equation. Additionally, the 2(np+1) algebraic variables Vi e Li only appear at (np+1) equations (represented as Eqn. 8). Hence, the DAE system presents high index.
   Consistent initialization and index reduction techniques are to be considered prior to numerical resolution of the batch distillation model with standard integration tools, as DASSL or DASSLC. In the Appendix, a procedure for index reduction is presented, which requires several steps of differentiation and algebraic manipulation. The main purpose of the developed code is to aid the researchers, so that they could avoid performing such cumbersome, tedious and error prone algebraic manipulations.
   It should be pointed out that the procedure in the Appendix is not the one actually implemented on the code, but instead, the one that a researcher could carry out when dealing with this high-index DAE model. The main difference between the approaches is that a person can substitute an equation into another during the index reduction procedure. The code, on the other hand, includes all the generated equations and variables on the extended system, increasing its dimension.


   In the present contribution the characterization/initialization module described by Costa Jr. et al. (2001) is used. Based on Pantelides' algorithm (Pantelides, 1988), it uses a graph-theoretical approach to identify the minimum set of equations to be differentiated in order to generate an index-one formulation of the original DAEs. Through AD code ADOL-C (Griewank et al., 1996), the differentiation is performed automatically at an affordable computational cost, and, most important of all, the differentiation technique employed does not incur truncation error. The required structural pattern of the jacobian matrix is computed via numerical perturbation of the original equations and of the exact extended system.
   The user must supply the DAE system F(z,z')=0 according to DASSLC standards. Then, the code applies Pantelides' algorithm, constructs the extended system and informs the number r of degrees of freedom of the system. The user is prompted to give r arbitrary initial conditions, on which the code will perform a structural check for feasibility. It must be stressed that structural feasibility is a necessary but not sufficient condition for feasibility. That comes from the fact that structural computation cannot determine if a structurally regular matrix is not in fact a singular matrix.
   In the present code, when an initial arbitrary set is found to be structurally feasible, the code will carry out the resolution of the extended system by a Newton-type method. The consistency system was constructed via automatic differentiation, and hence it is exact. The initial conditions obtained from its resolution can be provided directly to the DASSLC numerical integration code. If the initial conditions are unfeasible, the user is warned and prompted to supply another arbitrary set.


   Using the present code, it is possible to integrate high-index DAE systems directly using DASSLC code. All index reductions and manipulation, needed in order to achieve consistent initialization, are performed automatically. The user must only give the original DAEs according to DASSLC standards, determine which are the arbitrarily set initial variables (based on the calculated number of degrees of freedom) and provide initial values for those variables.
   Characterization results for an example with 5 components and 12 stages are shown in Table 2. The feasible arbitrary initial set used in this example corresponds to all the variables ni,j ( i= 0, ..., 13 and j= 1, ..., 5), performing 70 state variables. Table 3 presents the numerical values actually set for ni,j , and Table 4 presents the initial conditions calculated by the code for the remaining variables.

Table 2. Characterization results.

   After the differentiation of the equations indicated in Table 2, the extended system presents 335 equations and 405 variables, 70 of which must be assigned arbitrarily. It should be pointed out that the original system has 223 equations and 307 unknowns, among state variables (223) and some of their derivatives (84). The 405-307=98 new variables represent derivatives of former algebraic variables (, and ). In a more general case, they could also comprise some higher order derivatives of differential variables.

Table 3. Initial conditions arbitrarily set.

Table 4. Computed initial conditions.

   It should be emphasized that other sets could be determined arbitrarily. There are (approximately 1080) possible sets to be assigned, and only some of them are feasible. It is for this reason that the code asks the user to inform the desired initial set to be tested for feasibility. In order to illustrate this idea, some sets other than the one used in this contribution have been verified for structural feasibility. The results of this analysis are presented in Table 5.

Table 5. Some sets verified for structural feasibility.

   The simulation results are shown in the next figures. Figures 2 to 4 show the mole fraction of the components in the condenser, in plate 6 and in the reboiler. These results can be compared with the ones presented by Holland and Liapis (1983) for the same model, and it can be noticed that reported results have been reproduced. The direct resolution of the DAE model arising from the mathematical modeling of the process eliminates the algebraic manipulation step to the researcher. With the present code, even higher index problems can be solved with minor intervention of the user.

Figure 2. Number of moles of each component in the condenser.

Figure 3. Number of moles of each component in the intermediate plate.

Figure 4. Number of moles of each component in the reboiler.

   In Fig. 5, it is shown the liquid flowrate leaving the condenser (L0). It follows from the constant mole number per stage assumption that all Li variables have equal values.

Figure 5. Liquid flowrate leaving the condenser L0.

   In Fig. 6, it is presented the heat duty of the condenser. It can be seen that the steady state value of Q0 is equal to the heat duty of the reboiler, which is coherent with the consideration of negligible heat loss throughout the column. Figures 7 and 8 show the behavior of the temperature and of the enthalpy of liquid and vapour phases in the condenser, in the reboiler and in the sixth stage of the column.

Figure 6. Condenser heat duty.

Figure 7. Temperature profiles in a few stages.

Figure 8. Enthalpy of liquid and vapour phases.


   The optimized operation and the control of a transient distillation process require the numerical resolution of the mathematical model representative of the process, a DAE system. On the presented formulation, the model presents high index. Characterization of the system, index reduction and consistent initialization are issues to be dealt with prior to the numerical resolution of the model. The drawbacks of this formulation may be overcome using the developed characterization/initialization module. This module has been coupled with the DASSLC code for numerical integration of index 1 DAEs. As a consequence, high-index DAE systems can be simulated with minor intervention of the user.
In this contribution, the start-up of a batch distillation process has been analyzed. Numerical results agree with reported profiles for concentration in the reboiler and in the condenser.
   The user must supply the DAE system according to DASSLC standards, choose r variables to be assigned arbitrarily at t=0 and provide numerical values for those variables. The code automatically checks for feasibility of the chosen set, computes a complete initial set of initial conditions and allows the numerical integration of the resulting DAE system by existing numerical codes. The automatic differentiation tool ADOL-C has been utilized to perform all the required differentiation.
   It should be pointed out that the direct resolution of high-index systems makes it feasible to perform the numerical resolution of DAE models of arbitrary index and dimension. Consequently, it encourages the development of more detailed mathematical models.

APPENDIX - Reducing the Index of the DAE Model

   In order to reduce the index of the DAE system comprised of Eqn. 1 to 8, the following procedure could be used. Equation 4 can be rewritten as:

Replacing yi,j by the right hand side of this expression at Eqn. 7, it follows that:


Differentiating Eqn. 12, 5 and 3 for the first time, the following expressions are obtained, respectively.




Substituting Eqn. 13 into 14, it is obtained:


Substituting Eqn.15 into 16, it is obtained:


Substituting Eqn.1 into 17, it is obtained:


Substituting Eqn.18 into 2, it is obtained:


Substituting Eqn.8 into 19, it is obtained:

i= 0, ..., np+1

   The algebraic variable Q0 and the np+1 molar fluxes Li are present at these np+2 equations (represented as Eqn. 20). In the above equations, variables with index i-1 and i+1 do not exist when i = 0 and i = np+1, respectively. Via this procedure, it has been built an index 1 system, composed of Eqn. 1-6 coupled with Eqn. 8 and 20. This system can be integrated with a standard code as DASSL or DASSLC. As only one differentiation of each equation has been necessary to build an index 1 system, the original system presents index 2.


Brenan, K., Campbell, S. , Petzold, L., Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, New York, Elsevier Sc. Publ. Co. (1989).         [ Links ]
Costa Jr., E.F., Vieira, R.C., Secchi, A.R., Biscaia, E.C., "Automatic Structural Characterization of DAE Systems", Proceedings of 11th ESCAPE, Denmark (2001)         [ Links ]
Griewank, A., Juedes, D., Mitev, H., Utke, J., Vogel, O., Walther, A., "ADOL-C: A Package for the Automatic Differentiation of Algorithms Written in C/C++", ACM TOMS 22, 131-167 (1996).         [ Links ]
Henley, E.J., Seader, J.D., Equilibrium-Stage Separation Operations in Chemical Engineering, New York , John Wiley & Sons Inc. (1981).         [ Links ]
Holland, C., Liapis, A.I., Computer Methods for Solving Dynamic Separation Problems, McGraw Hill Publishers (1983).         [ Links ]
Pantelides, C.C., "The Consistent Initialization of Differential-Algebraic Systems", SIAM Journal of Scientific and Statistical Computing, 9, 213-231 (1988).         [ Links ]
Petzold, L.R., "DASSL code, version 1989", Computing and Mathematics Research Division, Lawrence Livermore National Laboratory, L316, PO Box 808, Livermore, CA 94559 (1989).         [ Links ]
Secchi, A.R., "DASSLC: User's Manual, a Differential-Algebraic System Solver", Technical Report,, UFRGS, Porto Alegre, RS/Brazil (1992).         [ Links ]
Vieira, R.C., Biscaia Jr., E.C., "Direct Methods for Consistent Initialization of DAE Systems", Comp. Chem. Eng., 25, 1299-1311 (2001).
        [ Links ]

Received: September 16, 2001.
Accepted for publication: August 29, 2002.
Recommended by Guest Editors: J. Cerdá, S. Díaz and A. Bandoni.