Copyright (c), H.A. Kooijman, R. Taylor

October 1998


This book accompanies the ChemSep program, which was developed to allow students to do separation calculations on ordinary personal computers. This book is not a guide where we show how to use ChemSep (see the ChemSep user manual for that) but it is intended to supply technical background to help the user in his selection of models and correlations. It is hoped that sensible selections can be made by providing information on, descriptions of, and references to the models and correlations that are employed in ChemSep. Although we have tried to be as extensive as possible, it is impossible to describe all models and their underlaying theory, so references are given for further reading. There are probably many more literature models and correlations than are available in ChemSep, but we have tried to be as comprehensive as we could. Sometimes a choice had to be made in which models to implement without having any criteria to discriminate between models. Furthermore, not all models are applicable to a particular regime of operation. We try to adapt ChemSep as much as possible to comply with all model limitations and user requirements. This book serves as a replacement for the ``manual" information files that we used to distribute with ChemSep. Therefore, some parts of this document might still be incomplete or unorganized and any suggestions or remarks are welcome. Of course, any remarks on the ChemSep program are welcome as well.

This book is written in LaTeX, a complete typesetting language, and set in the standard Times-Roman 11 point font. It is also provided with ChemSep in ASCII text form (file CHEMSEP.TXT) for online reference which was generated with a LaTeX to ASCII converter. The conversion is limited, with the result that the ASCII text file contains some unconverted LaTeX formatting. A PostScript file (BOOK.PS) also generated by LaTeX can be downloaded from our ftp site. Harry Kooijman
Ross Taylor


Many people have helped to shape ChemSep. The project was started in 1988 at Delft University (The Netherlands) by Harry Kooijman, Arno Haket and Ross Taylor. The purpose was to make an interactive interface for doing equilibrium stage calculations on the PC platform. It had to be easy enough for use by students with little computer exposure and yet sufficiently comprehensive to solve the various problems encountered in a course on separation processes.

We would like to express our appreciation to Professor Hans Wesselingh (now at the University of Groningen, the Netherlands) who initially promoted the project and made various resources available and encouraged us by letting students use the program for their course work. This has been an indispensable source of feedback that has helped us to improve the program. We also like to thank Peter Verheijen for his enthusiasm and contributions in the early years of the project. Also, various students have worked on projects to check and improve the programs and documentation, which was very helpful. Finally, ChemSep owes its very existence to the Internet which enabled the authors to keep in touch and continue development while living on different continents.


Chapter 1. Solving Nonlinear Equations

In this chapter we discuss the methods employed in ChemSep to solve the separation problems at hand. Side-issues such as how to start an iterative method or how ChemSep solves the resulting linear system of equations also pass the revue.

1.1 Newton's method

ChemSep uses Newton's method to solve the system of (MESH) equations derived from the flash or column problems. Newton's method is a Simultaneous Correction (SC) method that each time corrects all the variables. To use it, the equations to be solved are written in the form
  F(x) = 0
where F is a vector consisting of all the equations to be solved and x is, again, the vector of variables. A Taylor series expansion of the function vector around the point x_o which the functions are evaluated gives (ignoring second and higher order terms):
  F(x) = F(x ) + J(x - x )
            o           o

where J is the Jacobian matrix of partial derivatives of F with respect to the independent variables x:
  J   = (d F  / d x )
   ij       i      j

If x is the actual solution to the system of equations, then F(x)=0 and we can rewrite the above equation as:
  J(x - x ) = - F(x )
         o         o

This linear system of equations may be solved for a new estimate of the vector x. If the new vector, x, obtained in this way does not actually satisfy the set of equations, F, then the procedure can be repeated using the calculated x as a new x. The entire procedure is summarized below.
  1. Set iteration counter, k, to zero, estimate x_o
  2. Solve linearized equations for x_k+1
  3. Check for convergence; if not obtained, increment k and return to step 2.
Solving the linear system does not require a full matrix inversion of the Jacobian and is normaly done with Gaussian elimination or some type of decomposition technique. If the Jacobian has a lot of zero entries (i.e. it is sparse) then the linear system can be much more efficiently solved by using a sparse linear solver. For Jacobians with specific structures special solvers can be employed which are more efficient than a complete elimination or decomposition. One very important property of the Newton's method is that the convergence is scale invariant and independent of the ordering of the equations. This means that the same convergence is obtained if one of the equations is multiplied by some number or if the equations are reordered in a different manner. This is very important, because this means we are free to order the equations to obtain a special Jacobian which might enable the use of a special solver. It also makes the method applicable to a wider range of problems and without requiring the user to scale equations or variables.

An important drawback of the Newton's method can be its sensitivity to the initial guess, x_o, since quadratic convergence is only achieved ``close" to the solution. In order to obtain convergence, Newton's method requires that reasonable initial estimates be provided for all independent variables. It is obviously impractical to expect the user of a SC method to guess this number of quantities. Thus, the designer of a computer code implementing a SC method must provide one or more methods of generating initial estimates of all the unknown variables. Several techniques have been developed to improve the convergence away from the solution and to prevent the method from taking a step in a wrong direction. The simplest and most common technique is to ``damp" the step of each variable to some range or fraction of the Newton step. However, this damping also reduces the methods effectiveness.

Simultaneous correction procedures have shown themselves to be generally fast and reliable, having a locally quadratic convergence rate in the case of Newton's method, and these methods are much less sensitive to difficulties associated with nonideal problems than are tearing methods. They also lend themselves to be easier extended with optimization, parametric sensitivity, or continuation methods.

1.2 Continuation method

A simple implementation of a continuation method is incorporated in ChemSep for more difficult problems. Continuation methods use a parameter to make a path from a known solution for a simplified model to the desired solution of the complete model. For example the Newton homotopy starts with the initial guess as model and follows the path to the solution of the real problem by solving
  0 = (1-t) (F(X ) - F(X)) + t F(X)

where t varies from 0 (where X=X_o) to 1 (where F(X)=0). Better continuation methods can be formulated while using a parameter which has some physical significance. In separations problems the most appropriate choice would be the degree to which mass transfer (between the present phases) prevails. In the equilibrium model the stage efficiency and in the nonequilibrium model the mass transfer rates represent this degree. Thus, they will be multiplied with a parameter t which will vary from 0 (no separation at all) to 1 (actual separation).

Chapter 2. Property Models

This chapter discusses the thermodynamic and physical property models available in ChemSep. The selection of these models can be quite important for the results produced by ChemSep. Most formulae are repeated here but additional reading is available in two main sources: References are in between parentheses, by combining the letter A or B with the page number, for example (A43). The model "types" are grouped by ChemSep menu.

2.1 Thermodynamic Properties

2.1.1 K-value models

2.1.2 Activity coefficient models

Here we discuss the activity coefficient models available in ChemSep. For an in depth discussion of these models see the standard references. For the calculation of activity coefficients and their derivatives (for diffusion calculations) see also Kooijman and Taylor (1991).

2.1.3 Vapour pressure models

2.1.4 Equations of State (EOS)

Three types of equations of state may be selected in ChemSep; Ideal Gas, Virial, and Cubic EOS. The fugacity coefficient of an ideal gas mixture (B3) is unity (since the fugacity represents the deviation from an ideal gas, and we use the natural logarithm of the fugacity as the fugacity coefficient). The pressure relation for an ideal gas is:
  P = (R T / V)
The Virial and cubic EOS are discussed in the sections below.

2.1.5 Virial EOS

2.1.6 Cubic EOS

2.1.7 Enthalpy

2.2 Physical Properties

A number of different polynomials is implemented in ChemSep to evaluate physical properties over a certain temperature range. These temperature correlations are assigned a unique number in the range of 0-255 (see Table 2.1). For each up to 5 parameters (A-E) are available. Table 2.2 shows which pure component properties can be modeled with temperature correlations and their typical correlation number.

Table 2.1. Temperature correlations

Equation Parameter(s) Formula
2 A,B A + B T
3 A-C A + B T + C T^2
4 A-D A + B T + C T^2 + D T^3
10 A-C exp ( A - B / C+T )
100 A-E A + B T + C T^2 + D T^3 + E T^4
101 A-E exp ( A + B / T + C ln T + D T^E )
102 A-D A T^B / 1 + C/T + D/T^2
103 A-D A + B exp ( -C / T^D )
104 A-E A + B / T + C / T^3 + D / T^8 + E / T^9
105 A-D A / B^(1 + (1 - T / C)^D)
106 A-E A (1-T_r)^(B + C T_r + D T_r^2 + E T_r^3 )
107 A-E A + B ( C / T/sinh (C / T)^2 + D ( D / T/cosh (D / T)^2
All types of equations may be used for any of the physical properties but, of course, some formulas were specifically developed for prediction of particular properties. Besides the parameters A-E the temperature limits of the correlation must also be present. If the temperature specified falls out of the temperature range of a correlation (or the temperature limits are missing/incomplete) normally an alternative (default) method will be used automatically.

Table 2.2. Component properties with the typical correlation number

Liquid density 105
Vapour pressure 101
Heat of vaporisation 106
Liquid heat capacity 100
Ideal gas heat capacity 107
Second virial coefficient 104
Liquid viscosity 101
Vapour viscosity 102
Liquid thermal conductivity 100
Vapour thermal conductivity 102
Surface tension 106
Ideal gas heat capacity (Reid Prausnitz and Poling) 4
Antoine 10
Liquid viscosity (Reid, Prausnitz and Sherwood) 2
Physical properties models can be selected manually or the automatic selection can be used (which is the default). Below we discuss the models for calculating physical properties for pure components and mixtures, for vapour or liquid phases. ChemSep uses an automatic selection when no model is selected at all and the selection is left as *'s. Depending on range, phase, conditions, data availability, and required property ChemSep will make a guess of the best model to use. ChemSep does allow you to pick default models, and will use them if the model's range is valid. In case a property cannot be computed with a specific model it will use an estimation method or a fixed estimate (it is a good habit to check predicted physical properties when possible).

Table 2.3. Default physical property correlations

Mixture liquid density Rackett
Component liquid density Polynomial
Vapour density Cubic EOS
Mixture liquid viscosity Molar averaging
Component vapour viscosity Polynomial/Letsou-Stiel
Mixture vapour viscosity Brokaw
Component vapour viscosity Polynomial
Mixture liquid thermal conductivity Molar average
Component liquid thermal conductivity Polynomial
Mixture vapour thermal conductivity Molar average
Component vapour thermal conductivity Polynomial/9B-3
Liquid diffusivity Kooijman-Taylor/Wilke-Chang
Vapour diffusivity Fuller et al.
Mixture surface tension Molar average
Component surface tension Polynomial
Liquid-liquid interfacial tension Jufu et al.
Certain methods require mixture (critical) properties, commonly used mixing rules are:

  T    = sum      x  T   
   c,m        i=1  i  c,i

  V    = sum      x  V   
   c,m        i=1  i  c,i

  Z    = sum      x  Z   
   c,m        i=1  i  c,i

  P    = Z    R T    / V   
   c,m    c,m    c,m    c,m

  M  = sum      x  M 
   m        i=1  i  i

which will be referred as the "normal" mixing rules. Reduced properties will be calculated by:
  T  = T / T 
   r        c

  P  = P / P 
   r        c

  V  = V / V 
   r        c

unless specified otherwise.

2.2.1 Liquid density

Mixture liquid densities (in kmol/m^3) are calculated with: Pure component liquid densities are computed from the Peng-Robinson EOS for temperatures above a components critical temperature, otherwise with one of the folowing methods:

The pure component liquid densities are corrected for pressure effects with the correction of Thomson et al. (1982) as described for the Hankinson and Thompson method for mixtures.

2.2.2 Vapour density

Vapour densities are computed with the equation of state selected for the thermodynamic properties (possible selections are Ideal gas EOS, Virial EOS, and Cubic EOS).

2.2.3 Liquid Heat Capacity

The mixture liquid heat capacity is the molar average of the component liquid heat capacities, which are generally computed from a temperature correlation. Alternatively the liquid heat capacity could be computed from a corresponding states method and the ideal gas capacity. Rowlinson (1969, see A140) proposed a Lee-Kesler heat capacity departure function which was later modified to:

   L       ig                      -1                               1/3   -1               -1
  C     - C    = 1.45 + 0.45 (1-T )   + 0.25 w [ 17.11 + 25.2 (1-T )    T    + 1.742 (1-T )   ]
    p,i      p                   r                                r      r               r

However, in ChemSep the temperature correlation is used for all temperatures to prevent problems arrising from using different liquid heat capacity methods in the same column (which especially trouble nonequilibrium models). Liquid heat capacities could also be computed from the selected thermodynamic models to circumvent this problem.

2.2.4 Vapour Heat Capacity

The mixture vapour heat capacity is the molar average of the component vapour heat capacities, which are computed from the ideal gas heat capacity (RPP) 4 parameter temperature correlation. If no parameters for this correlation are present, the vapour heat capacity temperature correlation is used (if within the temperature range).

2.2.5 Liquid Viscosity

Mixture liquid viscosity are computed from DIPPR procedure 8H from the pure component liquid viscosities from:

        L        c             L
  ln eta   = sum      z  ln eta  
         m        i=1  i        i

where z_i are either the mole fractions (for molar averaging, the default) or alternatively the weight fractions for mass averaging. A better method is from Teja and Rice (1981, A479). However, this method requires interaction parameters. Here a different mixing rule (for T_cijV_cij) is used which improves the model predictions with unity interaction coefficients:

  w  = sum      x  w 
   m       i=1   i  i

  M  = sum      x  M 
   m       i=1   i  i

               c        c
  V   = sum      sum      x  x  V   
   cm       i=1      j=1   i  j  cij

               1/3      1/3  3
  V    = (( V      + V      )  / 8)
   cij       ci       cj

                c        c
  T   = (sum      sum      x  x  T   V    / V  )
   cm        i=1      j=1   i  j  cij cij    cm

  T   V    = psi   (T   V   + T   V   / 2)
   cij cij      ij   ci  ci    cj  cj

where psi_ij is set to unity for all components. The liquid viscosity of the mixture is computed from two reference components
  ln (epsilon  eta ) = ln (epsilon  eta ) + [ ln (epsilon  eta ) - ln (epsilon  eta ) ] (( w  -w  / w  - w  ))
             m    m               1    1                 2    2               1    1        m   1    2    1

with epsilon defined as

  epsilon  = (V      / sqrt(T   M ))
         i     ci            ci  i

and the reference component vioscosities are evaluated at T T_ci / T_cm. Component liquid viscosities are calculated from the liquid viscosity temperature correlation if the temperature is within the valid range. Otherwise the component viscosity is computed with DIPPR procedure 8G, the Letsou-Stiel method (1973, see A471):

                 1/6                2/3
  xi = 2173.424 T       / sqrt(M ) P      
                    c,i         i      c,i

    (0)                               2      -5
  xi    = ( 1.5174 - 2.135 T  + 0.75 T   ) 10  
                            r          r

    (1)                              2      -5
  xi    = ( 4.2552 - 7.674 T  + 3.4 T   ) 10  
                            r         r

     L        (0)       (1)
  eta   = ( xi    + w xi    ) / xi

Alternatively the simple temperature correlation given in Reid et al. (RPS liquid viscosity, see A439) can be used:
  log eta = A + B / T
A high pressure correction by Lucas (A436) is used to correct the influence of the pressure on the liquid viscosity:

  eta = (1 + D (Delta P  / 2.118)  / 1 +C w  Delta P ) eta  
                       r                   i        r     SL

where eta_SL is the viscosity of the saturated liquid at P_vp, and
  Delta P  = (P - P  ) / P  
         r         vp     ci

                        -4           -0.03877
  A = 0.9991 - [4.674 10  /(1.0523 T          - 1.0513)]

                                        2             3
  C = - 0.07921 + 2.1616 T  - 13.4040 T   + 44.1706 T  
                          r            r             r

              4             5             6             7
  - 84.8291 T   + 96.1209 T   - 59.8127 T   + 15.6719 T  
             r             r             r             r

                        2.573 0.2906
  D = [0.3257/(1.0039-T      )      ]-0.2086

2.2.6 Vapour Viscosity

Mixture vapour viscosities are computed using DIPPR procedure 8D-1 from component viscosites as follows:

     L        c           L
  eta   = sum      (x  eta   / sum x  phi  )
      m        i=1   i     i        i    ij

where the interaction parameters phi_ij can be calculated by Wilke's (1950) method:

                                               1/4 2
  phi   = ((1 + sqrt( \eta  / \eta  ) (M  / M )   )  / sqrt( 8 (1 + M  / M ) ) )
     ij                   i       j     i    j                       i    j

or by Brokaw's method:
  phi   = S A sqrt( \eta  / \eta  )
     ij                 i       j

  sm = ( 4 / (1+M /M ) (1+M /M ) )   
                 j  i      i  j

                                                0.45                                  0.45
  A = (sm / sqrt(M /M )) ( 1 + ((M /M  - (M /M )    ) / 2 (1 + M /M )) + ((1 + (M /M )    ) / sqrt(sm) (1 + M /M )) )
                  i  j            i  j     i  j                 i  j             i  j                        i  j

If the Lennard-Jones energy parameter, epsilon (in Kelvin), and the Stockmayers polar parameter, delta, are known, S is calculated from:

                                                                                                  2                                    2
  S = (1 + sqrt( (T/\epsilon ) (T/\epsilon ) ) + (delta  delta /4) / sqrt(1 + T/\epsilon  + \delta   / 4) sqrt(1 + T/\epsilon  + \delta   / 4))
                            i             j            i      j                         i          i                         j          j

otherwise it is approximated by S=1. epsilon and delta can be estimated from:

  epsilon = 65.3 T    Z      
                  c,i     c,i

                  5     2
  delta = 1.744 10 9 (mu  / V    T   )
                             b,i  b,i

Where mu is the dipole moment in Debye. Vapour viscosities are a function of pressure and a correction is normally applied. Mixture properties are computed with the "normal" mixing rules. DIPPR procedure 8E can be used to compute the high pressure viscosity:
  rho  = 1 / V   
     c        c,m

  rho  = rho / rho 
     r            c

                  1/6                2/3
  xi = 2173.4241 T       / sqrt(M ) P      
                     c,m         m      c,m

  A = exp (1.4439 rho ) - exp (-1.111 rho     )
                     r                       r

  B = 1.08 10   A / xi
  eta   = eta + B

where rho is the vapour mixture molar density.

Both Wilke's and Brokaw's method require pure component viscosities. These are normally obtained from the vapour viscosity temperature correlations, as long as the temperature is within the valid temperature range. If not, then the viscosity can be computed with the Chapman-Enskog kinetic theory (see Hirschfelder et al. 1954 and A391-393):

  T  = T / epsilon

               * -b               *                *
  Omega  = a (T )   + c / exp (d T ) + e / exp (f T )

     V           -             2                    2    *
  eta  = 26.69 10 7 M T / sigma  (Omega  + 0.2 delta  / T )

where the collision integral constants are a=1.16145, b=0.14874, c=0.52487, d=0.77320, e=2.16178, and f=2.43787. The viscosity may also be computed with the Yoon and Thodos method (DIPPR procedure 8B):

                   1/6                2/3
  xi  = 2173.4241 T       / sqrt(M ) P      
    i                 c,i         i      c,i

     V            b                                     8
  eta   = (1 + a T   - c exp (d T-r) + e exp (f T ) / 10  xi)
      i            r                             r

where the constants a-f are given in Table 2.4.

Table 2.4. Constants for the Yoon-Thodos method

Hydrogen Helium Others
a=47.65 a=52.57 a=46.1
b=0.657 b=0.656 b=0.618
c=20.0 c=18.9 c=20.4
d=-0.858 d=-1.144 d=-0.449
e=19.0 e=17.9 e=19.4
f=-3.995 f=-5.182 f=-4.058
Another method for calculating the vapor viscosity is the Lucas (A397) method:

          -7          0.618
  eta = 10   [0.807 T       - 0.357 exp (-0.449T ) +
                     r                          r

                                   o   o
  0.340 exp (-4.058 T ) + 0.018] F   F   / xi
                     r            p   q

                      3    -5    4   1/6
  xi = 0.176 (( T  / M  (10   P )  ))   
                 c             c

where F_p^o and F_q^o are polarity and quantum correction factors. The polarity correction depends on the reduced dipole moment:

                             -30 2    -5         2
  mu  = 52.46 ((mu / 3.336 10   )  (10   P ) / T  )
    r                                     c     c

If mu_r is smaller than 0.022 then the correction factor is unity, else if it is smaller than 0.075 it is given by

    o                         1.72
  F   = 1 + 30.55 (0.292 - Z )    
   p                        c

else by

    o                         1.72
  F   = 1 + 30.55 (0.292 - Z )     [0.96 + 0.1 (T  - 0.7)]
   p                        c                    r

The quantum correction is only used for quantum gases He, H_2, and D_2,

    o         0.15                     2 1/M
  F   = 1.22 Q     (1 + 0.00385[(T -12) ]    sign(T -12) )
   q                              r                r

where Q=1.38 (He), Q=0.76 (H_2), Q=0.52 (D_2). There is also a specific correction for high pressures (A421) by Lucas.

  eta = Y F  F  eta 
           p  q

               e       f          d -1
  Y = 1 + (a P   / b P   + (1+c P  )  )
              r       r          r

               o       -3     o
  F  = (1 + (F   - 1) Y   / F  )
   p          p              p

               o        -1               4      o
  F  = (1 + (F   - 1) [Y   - 0.007 (ln Y) ] / F  )
   q          q                                q

where eta^o refers to the low-pressure viscosity (note that the original Lucas method has a different rule for Y if T_r is below unity, however, this introduces a discontinuity which is avoided here). The parameters a through f are evaluated with:

               -3                    -0.3286
  a = (1.245 10   / T ) exp 5.1726 T        
                     r              r

  b = a (1.6553 T  - 1.2723)

  c = (0.4489 / T ) exp 3.0578 T         
                 r              r

  d = (1.7368 / T ) exp 2.2310 T        
                 r              r

  e = 1.3088

  f = 0.9425 exp -0.1853 T       

where, in case T_r is below unity, T_r is taken to be unity. For mixtures the Lucas model uses the following mixing rules:

  T   = sum      y  T  
   cm       i=1   i  ci

  V   = sum      y  V  
   cm       i=1   i  ci

  Z   = sum      y  Z  
   cm       i=1   i  ci

  P   = R T   Z   / V  
   cm      cm  cm    cm

  M  = sum      y  M 
   m       i=1   i  i

     o          c       o
  F    = sum      y  F   
   pm        i=1   i  pi

     o            c       o
  F    = A sum      y  F   
   qm          i=1   i  qi

where A is a correction factor depending on the molecular weights of the components in the mixture. Let H denote the component of highest molecular weight and L of lowest, then if M_H/M_L > 9:

  A = 1 - 0.01 (M  / M )    
                 H    L

else A=1. The mixture vapor viscosity is computed with the Lucas method as for a component which has the mixture properties T_cm, P_cm, M_m, F_pm^o, and F_qm^o. Therefore, the method is not interpolative in the same way as the techniques of Wilke and Brokaw (that is, the method does not necessarily lead to pure component viscosity eta_i when all y_j=0 except y_i=1).

2.2.7 Liquid Thermal Conductivity

The mixture liquid thermal conductivity, lambda^L_m (W/m K), can be computed using the following methods from the component liquid thermal conductivities: A correction is applied when the pressure is larger than 3.5 bar:

                    1.2                                     1.4
  lambda   = (0.63 T     P  / (30 + P ) + 0.98 + 0.0079 P  T     ) lambda
        hp             r  r          r                   r     r

This is DIPPR procedure 9G-1 where the mixture parameters are computed by the "normal" mixing rules. Component liquid thermal conductivities are calculated from one of the following methods:

2.2.8 Vapour Thermal Conductivity

If the system pressure is larger than 1 atmosphere a corection is applied according to DIPPR procedure 9C-1. Mixture parameters are computed using the "normal" mixing rules. Critical and reduced densities are computed from:
  rho  = (1 / V   )
     c         c,m

  rho  = (rho / rho )
     r             c

If the reduced density is below 0.5 then a=2.702, b=0.535, and c=-1; if the reduced density is witin [0.5,2] then a=2.528, b=0.67, and c=-1.069; otherwise a=0.574, b=1.155, and c=2.016. The high pressure thermal conductivity correction is then calculated from:

                      -8                                   1/6       2/3         5
  Delta lambda = (a 10   (exp (b rho ) + c) / (( sqrt(M ) T       / P       ) ) Z     )
                                    r                  m      c,m       c,m       c,m

which must be added to the calculated thermal conductivity for low pressure.

Pure component vapour thermal conductivities are estimated from the following methods:

2.2.9 Liquid Diffusivity

Generalized Maxwell-Stefan binary diffusion coefficients DD_ij are computed from the Kooijman-Taylor (1990) correlation where

    k      o
  DD    = D   , k=i
     ij     ij

    k      o
  DD    = D   , k=j
     ij     ji

    k           o    o
  DD    = sqrt(D    D   ), k != i, k != j
     ij          ik   jk

             c       k    x_k
  DD   = sum      (DD    )   
    ij        i=1     ij

Liquid binary infinitive diffusion coefficients (D^o_ij) are normally computed by the Wilke-Chang method unless selected otherwise. The following models are available:

2.2.10 Vapour Diffusivity

Generalized Maxwell-Stefan binary diffusion coefficients D_ij are equal to the normal binary diffusion coefficients (since the gas is considered an ideal system for which the thermodynamic matrix is the identity matrix). Normally these are computed with the Fuller-Schettler-Giddings method (see A587) but if Fuller volume parameters are missing the Wilke-Lee modification of the Chapman-Enskog kinetic theory is used.

2.2.11 Surface Tension

Mixture: Component surface tensions are only determined for temperatures below the component's critical temperature, otherwise it is assumed that the component does not contribute to the mixture surface tension (i.e. sigma_i=0). The following methods are available:

2.2.12 Liquid-Liquid Interfacial Tension

This property is only required for simulating Liquid-Liquid extractors with the nonequilibrium model. API method 10B3 uses the calculated surface tensions for both liquid phases and the interfacial tension, sigma^,, is computed from

  sigma  = sigma  + sigma  - 1.1 sqrt( \sigma  \sigma  )
                1        2                   1       2

This method generally overpredicts the interfacial tension for aqueous systems. We use a general method from Jufu et al. (1986):

              ,,     ,
  X = - ln ( x    + x   + x   )
                1     2    3r

                                     ,,        ,
  sigma , = (K R T X / A   exp (X) (x    q  + x   q  + x   q ))
                        w0             1  1     2  2    3r  3

with A_w0=2.5 10^5 (m^2/mol), R=8.3144 (J/mol/K), K=0.9414 (-), and q_i is the UNIQUAC surface area parameters of the components i. The components are ordered in such a manner that component 1 and 2 are the dominating components in the two liquid phases. Then the rest of the components are lumped into one mole fraction, x_3. This lumped mole fraction is taken for the phase which has the largest x_3 (the richest). q_3 is the molar averaged q for that phase for all components except 1 and 2.

Symbol List

a, b    Cubic EOS parameters
B         Second virial coefficient (m3/kmol)
c         Number of components
C_p       Mass heat capacity (J/kg.K)
D_ij    Binary diffusion coefficient (m^2/s)

DD_ij  Binary Maxwell-Stafan diffusion coefficient (m^2/s)
D^o_ij  Infinite dilution binary diffusion coefficient (m^2/s)
K_i       K-value of component i, equilibrium ratio (K_i=y_i/x_i)
k_ij    Binary interaction coefficient (for EOS)
M         Molecular mass (kg/kmol)
R         Gas constant = 8134 (J/kmol K)
r         Radius of gyration (Angstrom)
P         Pressure (Pa)
P^*,P_vap   Vapour pressure (Pa)
P_i       Parachor (m^3 kg^1/4 / s^1/2) of component i
PF        Poynting correction
q         UNIQUAC surface area parameter
T         Temperature (K)
T_r       Reduced temperature (T_r=T/T_c)
T_b       Normal boiling temperature (K)
V         Molar volume (m^3/kmol)
V_b       Molar volume at normal boiling point (m^3/kmol)
V_s       Saturated molar volume (m^3/kmol)
w         Weight fraction (of component)
x         Liquid mole fraction (of component)
y         Vapour mole fraction (of component)
Z         Compressibility
Z_R       Racket parameter
alpha    Attractive parameter in EOS
w    Acentric factor
Omega_a, Omega_b  EOS parameters
Omega_v  Collision integral for viscosity
Omega_D  Collision integral for diffusion
gamma    Activity coefficient
delta    Stockmayer parameter
epsilon  Molecular energy parameter (K)
lambda   Thermal conductivity (W/m.K)
rho      Molar density (kmol/m^3)
eta      Viscosity (Pa.s)
phi_i    Fugacity coefficient of component i
phi_s    Association factor for solvent s (Hayduk-Laudie)
phi_ij   Interaction parameter for viscosities
phi_i    Fugacity coefficient of component i
phi^*  Pure fugacity coefficient at saturation
sigma    Surface tension (N/m)
            Collision diameter (Angstrom)
sigma_b  Surface tension at T_b (N/m)
sigma ,  Liquid-liquid interfacial tension (N/m)
mu       Dipole moment (Debye)
xi       Inverse viscosity (defined in text)
L    Liquid
V, G  Vapour, gas
*    Saturated liquid,
b    at normal boiling point
c    critical
i    of component i
j    of component j
m    mixture
r    reduced
s    saturated liquid
EOS    Equation of State
RK     Redlich-Kwong
SRK    Soave Redlich-Kwong
PR     Peng-Robinson


Digulio, Teja, Chem. Eng. J., Vol. 38 (1988) pp. 205.

E.N. Fuller, K. Ensley, J.C. Giddings, ``A New Method for Prediction of Binary Gas-Phase Diffusion Coefficients'', Ind. Eng. Chem., Vol. 58, (1966) pp. 19--27.

E.N. Fuller, P.D. Schettler, J.C. Giddings, ``Diffusion of Halogonated Hydrocarbons in Helium. The Effect of Structure on Collision Cross sections'', J. Phys. Chem., Vol. 73 (1969) pp. 3679--3685.

W. Hayduk, H. Laudie, ``Prediction of Diffusion Coefficients for Nonelectrolytes in Dilute Aqueous Solutions'', AIChE J., Vol. 20, (1974) pp. 611--615.

W. Hayduk, B.S. Minhas, ``Correlations for Prediction of Molecular Diffusivities in Liquids'', Can. J. Chem. Eng., 60, 295-299 (1982); Correction, Vol. 61, (1983) pp. 132.

J.O. Hirschfelder, C.F. Curtis, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York (1954).

F. Jufu, L. Buqiang, W. Zihao, ``Estimation of Fluid-Fluid Interfacial Tensions of Multicomponent Mixtures'', Chem. Eng. Sci., Vol. 41, No. 10 (1986) pp. 2673--2679.

R. Taylor, H.A. Kooijman, ``Composition Derivatives of Activity Coefficient Models (For the estimation of Thermodynamic Factors in Diffusion)'', Chem. Eng. Comm., Vol. 102, (1991) pp. 87--106.

A. Letsou, L.I. Stiel, AIChE J., Vol. 19 (1973) pp. 409.

Lielmezs, Herrick, Chem. Eng. J., Vol. 32 (1986) pp. 165.

J.M. Prausnitz, T. Anderson, E. Grens, C. Eckert, R. Hsieh, J. O'Connell, Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria, Prentice-Hall (1980).

J.S. Rowlinson, Liquids and Liquid Mixtures, 2nd Ed., Butterworth, London (1969).

R.C. Reid, J.M. Prausnitz, T.K. Sherwood, Properties of Gases and Liquids, 3rd Ed., McGraw-Hill, New York (1977).

R.C. Reid, J.M. Prausnitz and B.E. Poling, The Properties of Gases and Liquids, 4th Ed., McGraw-Hill, New York (1988).

M.A. Siddiqi, K. Lucas, ``Correlations for Prediction of Diffusion in Liquids'', Can. J. Chem. Eng., Vol. 64 (1986) pp. 839--843.

A.S. Teja, P. Rice, ``Generalized Corresponding States Method for the Viscosities of Liquid Mixtures'', Ind. Eng. Chem. Fundam., Vol. 20 (1981) pp. 77-81.

S.M. Walas, Phase Equilibria in Chemical Engineering, Butterworth Publishers, London (1985).

C.R. Wilke, J. Chem. Phys., Vol. 18 (1950) pp. 617.

C.R. Wilke, P. Chang, ``Correlation of Diffusion Coefficients in Dilute Solutions'', AIChE J., Vol. 1 (1955) pp. 264-270.

C.R. Wilke, C.Y. Lee, Ind. Eng. Chem., Vol. 47 (1955) pp. 1253.

P.H. Winterfeld, L.E. Scriven, H.T. Davis, AIChE J., Vol. 24 (1978) pp. 1010.

Chapter 3. Flash Calculations

A flash is a one stage operation where a (multiple phase) feed is "flashed" to a certain temperature and/or pressure and the resulting phases are separated. The flash in ChemSep deals only with two different phases leaving, a vapour and a liquid. Liquid-Liquid or multiphase Vapour-Liquid-Liquid flashes are currently not yet supported in ChemSep. For more information see also the general references given at the end of this chapter.

3.1 Equations

The vapour and liquid streams leaving the flash are assumed to be in equilibrium with each other. The equations that model equilibrium flashes are summarized below:

where F is the molar feedrate with component mole fractions z_i. V and L are the leaving vapour and liquid flows with mole fractions y_i and x_i, respectively. Equilibrium ratios K_i and enthalpies H are computed from property models as discussed in chapter 2. Q is defined as the heat added to the feed before the flash. If we count the equations listed, we will find that there are 2c+3 equations, where c is the number of components. As flash variables we have (depending on the type of flash):

Since we have 2c+3 equations, two of the 2c+5 variables above need not be specified. ChemSep allows the following nine flash specifications:

3.2 Solution of the Flash Equations

FLASH uses Newton's method for solving flash problems as well as simpler bubble and dew point calculations. The vector of variables used in the PQ-FLASH is:

  (X)  = (V, y  , y  ... y  , T, x  , x  ... x  , L)
              1    2      c       1    2      c

the vector of functions, (F), is:

  (F)  = (TMB ,CMB , CMB  ... CMB , H, EQM , EQM  ... E , SUM),
                  1     2        c        1     2      c

The structure of the Jacobian matrix [J] is shown below:


Henley, E.J., J.D. Seader, Equilibrium-Stage Separation Operations in Chemical Engineering, Wiley (1981).

King, C.J., Separation Processes, Second Edition, McGraw Hill (1980).

Chapter 4. Equilibrium Columns

This chapter describes the equilibrium stage model for column operations such as distillation, absorption and extraction. The equations that ChemSep solves are discussed as well as other issues.

4.1 Introduction

Multicomponent separation processes like distillation, absorption and extraction have been modelled using the equilibrium stage concept for a century. The equilibrium stage model was first used by Sorel in 1893 to describe the rectification of alcohol. Since that time it has been applied with ever increasing frequency to all manner of separation processes: distillation (including rectification, stripping, simple (single feed, two product columns), complex (multiple feed, multiple product columns), extractive, azeotropic and petroleum refinery distillation), absorption, stripping, liquid-liquid and supercritical extraction.

The equations that model equilibrium stages are called the MESH equations. The MESH equations for the interior stages of a column together with equations for the reboiler and condenser (if they are needed) are solved together with any specification equations to yield, for each stage, the vapour mole fractions; the liquid mole fractions; the stage temperature and the vapour and liquid flowrates.

Since the late 1950's, hardly a year has gone by without the publication of at least one (and usually more than one) new algorithm for solving the equilibrium stage model equations. One of the incentives for the continued activity has always been (and remains) a desire to solve problems with which existing methods have trouble. The evolution of algorithms for solving the MESH equations has been influenced by, among other things: the availability (or lack) of sufficient computer storage and power, the development of mathematical techniques that can be exploited, the complexity of physical property (K-value and enthalpy) correlations and the form of the model equations being solved.

It is not completely clear who first implemented a simultaneous correction method for solving multicomponent distillation and absorption problems. As is so often the case, it would appear that the problem was being tackled by a number of people independently. Simultaneous solution of all the MESH equations was suggested as a method of last resort by Friday and Smith (1964) in a classic paper analysing the reasons why other algorithms fail. They did not, however, implement such a technique. The two best known and most frequently cited papers are those of Goldstein and Stanfield (1970) and Naphtali and Sandholm (1971), the latter providing more details of an application of Newton's method described by Naphtali at an AIChE meeting in May 1965.

To the best of our knowledge, a method to solve all the MESH equations for all stages at once using Newton's method was first implemented by Whitehouse (1964) (see, also, Stainthorp and Whitehouse, 1967). Among other things, Whitehouse's code allowed for specifications of purity, T, V, L or Q on any stage. Interlinked systems of columns and nonideal solutions could be dealt with even though no examples of the latter type were solved by Whitehouse. Since the pioneering work of Whitehouse, Naphtali and Sandholm and Goldstein and Stanfield, many others have employed Newton's method or one of its relatives to solve the MESH equations.

Simultaneous correction procedures have shown themselves to be generally fast and reliable. Extensions to the basic method to include complex column configurations, interlinked columns, nonstandard specifications and applications to column design result in only minor changes in the algorithm. In addition, simultaneous correction procedures can easily incorporate stage efficiencies within the calculations (something that is not always possible with other algorithms). Developments to about 1980 have been described in a number of textbooks (see, for example, Holland, 1963, 1975, 1981; King, 1980; Henley and Seader, 1981) and a recent review by Wang and Wang (1980). Seader (1985) has written an interesting history of equilibrium stage simulation.

Seader (1986) lists a number of things to be taken into consideration when designing a simultaneous correction method; a revised and extended list follows and is discussed in more detail below.

  1. What equations should be used?
  2. What variables should be used?
  3. How should the equations be ordered?
  4. How should the variables be ordered?
  5. How should the linearized equations be solved?
  6. Should the Jacobian be updated on each iteration or should it be held constant for a number of iterations or should it be approximated using quasi-Newton methods. Should derivatives of physical properties be retained in the calculation of the jacobian (J)?
  7. Should flexibility in specifications be provided and, if so, how?
  8. What criterion should we use to determine convergence?
  9. How should the initial guess be obtained?
  10. What techniques should we use to improve reliability?

4.2 Equations

Each equilibrium stage in the column has a vapour entering from the stage below and liquid from a stage above. They are brought into contact on the stage together with any fresh or recycle feeds. The vapour and liquid streams leaving the stage are assumed to be in equilibrium with each other. A complete separation process is modeled as a sequence of s of these equilibrium stages. Each stage can have optional sidedraws where part of the vapour or liquid stream leaving the stage is leaving the column.

The equations that model equilibrium stages are termed the MESH equations, MESH being an acronym referring to the different types of equations that form the mathematical model. The M equations are the Material balance equations, the E equations are the Equilibrium relations, the S equations are the Summation equations and the H equations are the entHalpy balances:

where we have vapour and liquid leaving flows from stage j, V_j and L_j, with mole fractions y_ij and x_ij, vapour and liquid sidedraws, W_j and U_j, Feeds F_j with mole fractions z_ij, equilibrium K-values K_ij, enthalpies H_j, and stage heat duty Q_j.

If we count the equations listed, we will find that there are 2c+4 equations per stage. However, only 2c+3 of these equations are independent. These independent equations are generally taken to be the c component mass balance equations, the c equilibrium relations, the enthalpy balance and two more equations. These two equations can be the two summation equations or the total mass balance and one of the summation equations (or an equivalent form). The 2c+3 unknown variables determined by the equations are the c vapour mole fractions, y; the c liquid mole fractions, x; the stage temperature, T and the vapour and liquid flowrates, V and L. For a column of s stages, we must solve s(2c+3) equations. The table below shows how we may easily end up having to solve hundreds or even thousands of equations.

c s s(2c+3)
2 10 70
3 20 180
5 50 650
10 30 690
40 100 8300

The first entry in this table corresponds to a simple binary problem that could easily be solved graphically. The second and third are fairly typical of the size of problem encountered in azeotropic and extractive distillation processes. The last two entries are typical of problems encountered in simulating hydrocarbon and petroleum mixture separation operations.

4.3 Condenser and Reboiler

The MESH equations can be applied as written to any of the interior stages of a column. In addition to these stages, the reboiler and condenser (if they are included) for the column must be considered. The MESH equations shown above may be used to model these stages exactly as you would any other stage in the column. For these special stages it is common to use some specification equation instead of the enthalpy balance. Typical specifications include:

  1. the flowrate of the distillate/bottoms product stream,
  2. the mole fraction of a given component in either the distillate or bottoms product stream,
  3. a component flow rate in either the distillate or bottoms product stream,
  4. a reflux/reboil ratio or rate,
  5. the temperature of the condenser or reboiler,
  6. a heat duty to the condenser or reboiler.
ChemSep includes all of these specifications as well as a few others that have not been listed.

In the case of a total condenser, the vapour phase compositions used in the calculation of the equilibrium relations and the summation equations are those that would be in equilibrium with the liquid stream that actually exists. That is for the total condenser, the vapour composition used in the equilibrium relations is the vapour composition determined during a bubble point calculation based on the actual pressure and liquid compositions found in the condenser. At the same time, these compositions are not used in the component mass balances since there is no vapour stream from a total condenser.

4.4 "Nonequilibrium" Stages

In actual operation the trays of a distillation column rarely, if ever, operate at equilibrium despite attempts to approach this condition by proper design and choice of operating conditions. The degree of separation is, in fact, determined as much by mass and energy transfer between the phases being contacted on a tray or within sections of a packed column as it is by thermodynamic equilibrium considerations. The usual way of dealing with departures from equilibrium in multistage towers is through the use of stage and/or overall efficiencies. The Murphree stage efficiency is most often used in separation process calculations because it is easily combined with the equilibrium relations:

         MV                        MV
  E   = E    K   x   - y   - (1 - E   ) y      = 0
   ij      j  ij  ij    ij           j   i,j+1

where E^MV_j is the Murphree vapour efficiency for stage j. If the efficiency is unity we obtain the original equilibrium relation from above. The Murphree efficiency must be specified for all stages except the condenser and reboiler which are assumed to operate at equilibrium.

4.5 Solution of the MESH Equations

Almost every one of the many numerical methods that have been devised for solving systems of nonlinear equations has been used to solve the MESH equations. However, as mentioned earlier, ChemSep uses mostly the Newton's method to solve the nonlinear algebraic MESH equations. Here we will discuss how ChemSep uses the Newton's method.

4.5.1 How to Order the Equations and Variables?

Separation problems in ChemSep result in stages with each a set of various types of equations. There are basically two ways to group the equations and variables: by type or by stage. Grouping the equations and variables by stage is preferred for problems with more stages than components (practically all distillation and many absorption and extraction problems) while grouping by type is preferred for systems with more components than stages (some gas absorption problems). ChemSep employs a by-stage grouping of the equations and variables. We define a vector of variables for the j-th stage as:

  (X )  = (V , y  , y   ... y  , x  , x   ... x  , L )
    j       j   1j   2j      cj   1j   2j      cj   j

and a vector of functions for the j-th stage (F_j):

      T     T                                            L-V
  (F )  = (M  , M  , M   ... M  , H , E  , E   ... E  , S    )
    j        j   1j   2j      cj   j   1j   2j      cj      j


   L-V       c
  S    = sum      (x   - y  ) = 0
              i=1   ij    ij

If the equations and variables are grouped by stage we have

     T        T      T         T
  (F)  = ((F ) , (F )  ... (F )  )
            1      2         s

     T        T      T         T
  (X)  = ((X ) , (X )  ... (X )  )
            1      2         s

4.5.2 The Jacobian

To evaluate the Jacobian matrix, one must obtain the partial derivative of each function with respect to every variable. Part of the appeal of the grouping by stage approach is that for single columns at least the Jacobian matrix is block tridiagonal in structure:

. . .
. . .
[Am] [Bm]

in which each entry [A], [B], [C] is a matrix in its own right. The [A] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for stage j-1. The [B] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for the j-th stage. Finally, the [C] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for stage j+1. The structure of the submatrices [A], [B] and [C] is shown below.

The elements of [B] include partial derivatives of K-values with respect to temperature and composition. Since it is rather a painful experience to differentiate, for example, the UNIQUAC equations with respect to temperature and composition, in some SC codes this differentiation is done numerically. This can be an extremely time consuming step. However, neglect of these derivatives is not recommended unless one is dealing with nearly ideal solutions, since, to do so, will almost certainly lead to an increase in the required number of iterations or even to failure.

Almost all of the partial derivatives needed in ChemSep are computed from analytical expressions. The exception is the temperature derivatives of the excess enthalpy which requires a second differentiation with respect to temperature of the activity and fugacity coefficient models. Coding just the first derivatives was bad enough.

If the column has pumparounds extra matrices will be present which are not on the diagonal and the use of block tridiagonal methods becomes less straight forward. Similar problems arise with non-standard specifications that are not on the variables of the condenser and reboiler. When we solve multiple interlinked columns (currently not supported by ChemSep) a special ordering is required to maintain the diagonal structure of the jacobian. Therefore, ChemSep now uses a sparse solver to solve the system of equations involved. In the case of non-standard specifications (which allow "wild" specification such as the temperature on a stage in the middle of a column, instead at one of its ends) the specification equation must be calculated using finite differencing.

4.5.3 How Should the Linearized Equations be Solved?

It is absolutely essential to take account of the sparsity of the Jacobian when solving the linearized equations; straightforward matrix inversion is totally impractical (and probably numerically impossible). Linear systems with a block tridiagonal structure may reasonably efficiently be solved using a generalized form of the Thomas algorithm. The steps of this algorithm are given by Henley and Seader (1981).

Still further improvements in the block elimination algorithm for solution of separation process problems can be effected if we take advantage of the special structure of the submatrices [A], [B] and [C]. In fact, [A] and [C] are nearly empty. ChemSep uses the sparse matrix solver NSPIV for solving the sparse linear system effectively.

As ChemSep calculates the Jacobian mostly analytically, the Jacobian is computed for each iteration of the Newton's method. This has shown to be the most effective method until now. Alternatively, the Jacobian could be held constant for a number of iterations or approximated using quasi-Newton methods. ChemSep currently does not implement these options. Another possibility is to selectively retain the derivatives of physical properties in the calculation of the Jacobian. For dynamic simulations this has shown to be counterproductive (Kooijman, 1995) and therefore ChemSep doesn"t support such a technique.

4.5.4 The Convergence Criterion

Convergence of the Newton's method in ChemSep is determined from the residuals of the equations which make up the model. These function-values are squared and summed to give one convergence parameter: the sum of squares (SS). In ChemSep we compare the square root of this sum of squares with a user specified convergence criterion. Typical values for this convergence criterion vary between 10^-6 for equilibrium column and flash calculations to 10^-3 for nonequilibrium column calculations.

As ChemSep does not scale the equations it is possible that this sum of squares remains too large for such small convergence criterions. As this is not always forseeable, we extended the convergence test in ChemSep with a check on the change in the variables. If the absolute relative value of the change of each variable is smaller than the convergence criterion, ChemSep also considers the problem solved.

Finally, in order to prevent calculations to continue endlessly, there is a maximum number of iterations after which the Newton's method exits. ChemSep will then ask the user whether to continue the calculations, and for how many additional iterations. Typically, we found that 30 Newton iterations suffices to solve most separation problems at hand. Often, when more iterations are needed, there are other problems which prevent the method to converge. For example, specifications can be impossible to attain, or the specifications define a problem with multiple solutions, where the Newton's method oscillates between two or more solutions (such problems have only been found and documented recently).

An exception to this observation of convergence form the simulations of packed columns using the nonequilibrium with plug flow model approximations, where convergence tends to be slow.

4.5.5 The Initial Guess

In order to obtain convergence, Newton's method requires that reasonable initial estimates be provided for all s(2c+3) independent variables. ChemSep uses an automatic initialization procedure where the user does not need to make any guesses. Flow rates are estimated assuming constant molar flows from stage to stage. If the bottoms flow rate and reflux ratio are NOT specified and cannot be estimated from the specifications that are supplied then the bottoms flow rate is arbitrarily assigned a value of half the total feed flow and the reflux ratio is given a default value of 2. This, of course, could cause serious convergence problems. In the future, optional guesses might be added to the specifications to circumvent this problem.

The next step is to estimate the compositions and temperatures. This is done iteratively. We start by estimating the K-values assuming ideal solution behavior at the column average pressure and an estimate of the boiling point of the combined feeds. (This eliminates the normal requirement of estimates of end stage temperatures).

Mole fractions of both phases are estimated by solving the material balance and equilibrium equations for one component at a time.

We define a column matrix of discrepancy functions

     T      c           c           c
  (F)  = (ME    , ... ME    , ... ME     )
             i,1         i,j         i,s

where ME^c_ij is the component material balance equation combined with the equilibrium equations to eliminate the vapour phase mole fractions. Each equation depends on only three mole fractions. Thus, if we define a column matrix of mole fractions (X) by

  (X)  = (x   , ... x     , x  , x     , x    )
           i,1       i,j-1   ij   i,j+1   i,s

we may write
  (F) = [ABC](X) - (R) = (0)
With the equations and variables ordered in this way, the coefficient matrix [ABC] has three adjacent diagonals with coefficients:
  A  = L   
   j    j-1

  B  = - (V  K j + L  )
   j       j  i     j

  C  = V  K j
   j    j  i

The right hand side matrix (R) has elements
  R  = - F  z  
   j      j  ij

This linear system of equations can be solved for the mole fractions very easily using Gaussian elimination. Temperatures and vapour compositions are computed from a bubble point calculation for each stage. The bubble point computation provides K-values for all components on all stages. So we solve the tridiagonal system of equations again using the old flow rates and the new K-values. Temperatures are recomputed as before. This procedure is repeated a third time before proceeding with the main simulation.

For columns without condenser and reboiler a different initialization is used where the compositions of the liquid are set equal to the top liquid feed compositions and the compositions of the vapour equal to the bottom vapour feed compositions. Temperatures in the whole column are set equal to the temperature of the first feed specified.

For columns with either a condenser or a reboiler the compositions are initialized to the overall compositions of all feeds combined, and the temperatures to the bubble point at the column pressure and the overall feed compositions.

Currently there is no special initialization routine that will handle pumparounds if they are present. The result is that columns with large pumparoundflows will require flow initialization by the user, or, to repetitively solve the problem using the old results as the initialization and increasing pumparoundflow (see below).

4.5.6 Reliability

SC methods are far more reliable and versatile than most other methods. The same method will solve distillation, gas absorption and liquid extraction problems. It must be admitted though, that although the probability that Newton's method will converge from the automatic initial estimates is quite high, there is no guarantee of convergence. The difficulty of supplying good initial estimates is particularly severe for problems involving strongly nonideal mixtures, interlinked systems of columns and nonstandard specifications.

Several methods have been used to improve the reliability of Newton's method; damping the Newton step, use of steepest descent (ascent) formulations for some of the iterations, and combination with relaxation procedures; none of which has proven to be completely satisfactory. The methods most recently proposed for assisting convergence of Newton's method are continuation methods.

In default mode, ChemSep does NOT use any of these techniques, other than a check to make sure that all quantities remain positive. Mole fractions, for example, are not permitted to take on negative values. The user does have the option of supplying damping factors.

4.5.7 Damping factors

The Newton's method in ChemSep has some extra features that will enhance the convergence to the solution. The Newton's method computes a new solution vector based on the current Jacobian and function vector. However, the new solution vector might be physically meaningless, for example if a composition becomes smaller than zero or larger than unity. Also, the new solution vector might represent a too large a change in stage temperatures or flows for the method to be stable. To eleviate these problems ChemSep uses damping of the newton's iteration changes. Temperature changes are limited to a maximum (default to 10 K) and flow changes up to a maximum fraction of the old flows (default set to 0.5). The compositions require a special type of damping. If a composition is becoming negative or larger than unity, the change is limited to half the distance to the extreme. Also, if a damping factor is specified, the maximum change in composition equals the factor (the default factor is 1 allowing a change over the whole mol fraction range). This type of damping has turned out to be very effective. The damping factors can be found under the Options - Solve options menu.

If for some reason your column simulation does not converge, changing the damping factors might help. If you know the iteration history (by either limiting the number of iterations or by printing out the intermediate answers. Both can be done in the Options - Solve options menu) you can adapt the factors so the column simulation might converge. Note that convergence is mostly slower when you start to apply extra damping by making the factors smaller, the Newton method looses its effectiveness when damped. Nor does damping guarantee convergence.

4.5.8 User Initialization

For difficult problems it might be necessary for the user to provide initial temperature, pressure, or flow profiles. In case the stage temperatures or pressures are not calculated user initialization is a way to define these profiles.

The user can specify either temperature or flow profiles, or both. The only requirement is that values for the first and last stages are provided. Missing temperatures on intermediate stages are computed by linear interpolation, missing flows are computed on a constant flow from stage to stage basis. Therefore, it is better to specify the flows of the first and last two stages in case a condenser and reboiler are present. Composition profiles are computed through the method described above, however, temperatures are not computed using the bubble point calculations. If both user specified temperature and flow profiles are incomplete ChemSep switches to the automatic initialization method.

4.5.9 Initialization with Old Results

In some cases it might prove advantageous to use the converged results of a previous run as the initial guess for a new problem (for example, when bottoms flowrate and reflux rate are not specified and cannot be estimated, and the automatic initialization always uses a reflux ratio of 2). This is a very straight forward way of specifying the initial guess as long as the number of components remains the same. Care must be taken when feed or product specifications or locations are changed. The results are interpolated if the number of stages is changed.


J.R. Friday, B.D. Smith, "An Analysis of the Equilibrium Stage Separations Problem - Formulation and Convergence", AIChEJ, Vol 10, 698 (1964).

R.P. Goldstein, R.B. Stanfield, "Flexible Method for the Solution of Distillation Design Problems using the Newton-Raphson Technique", Ind. Eng. Chem. Process Des. Dev., 9, 78 (1970).

E.J. Henley, J.D. Seader, Equilibrium-Stage Separation Operations in Chemical Engineering, Wiley (1981).

C.D. Holland, Multicomponent Distillation, Prentice Hall Inc., NJ (1963).

C.D. Holland, Fundamentals and Modelling of Separation Processes, Prentice Hall Inc., NJ (1975).

C.D. Holland, Fundamentals of Multicomponent Distillation, McGraw-Hill Inc.; New York (1981).

C.J. King, Separation Processes, Second Edition, McGraw Hill (1980).

L.M. Naphtali, "The Distillation Column as a Large System", presented at AIChE 56-th National Mtg., San Francisco, May 16, (1965).

L.M. Naphtali, D.P. Sandholm, "Multicomponent Separation Calculations by Linearization", AIChE J., Vol 17 (1), 148 (1971).

J.D. Seader, "The BC (Before Computers) and AD of Equilibrium Stage Operations", Chem. Eng. Ed., Spring, 88, (1985a).

J.D. Seader, Chapter on Distillation in Chemical Engineers Handbook, (Green D. Editor), 6th Edition, McGraw Hill, New York, (1986)

J.D. Seader, Computer Modelling of Chemical Processes, AIChE Monograph Series, No. 15, 81 (1986).

B.D. Smith, Design of Equilibrium Stage Processes, McGraw-Hill, New York, (1964).

F.P. Stainthorp, P.A. Whitehouse, "General Computer Programs for Multi Stage Counter Current Separation Problems - I: Formulation of the Problem and Method of Solution", I. Chem. E. Symp. Ser., Vol 23, 181 (1967).

J.C. Wang, Y.L. Wang; "A Review on the Modeling and Simulation of Multi-Stage Separation Processes" in Foundations of Computer-Aided Chemical Process Design, vol. II; R.S.H. Mah and W.D. Seider, eds.; Engineering Foundation; 121 (1981).

S.M Walas, Phase Equilibria in Chemical Engineering, Butterworth Publishers, (1985).

P.A. Whitehouse, A General Computer Program Solution of Multicomponent Distillation Problems, Ph.D. Thesis in Chem.Eng., University of Manchester, Institute of Science and Technology, Manchester, England (1964).

Chapter 5. Nonequilibrium Columns

The nonequilibrium model and the model equations are introduced. Models that describe the mass transfer, the flow type, pressure drop, entrainment, and weeping are discussed. The design method which enables the simultaneous design of the the column layout and column simulation is explained.

5.1 The Nonequilibrium Model

A second generation nonequilibrium model was developed by Taylor and coworkers and is described in detail by Taylor et al. (1994) and Taylor and Krishna (1993). It can be used to simulate trayed columns as well as packed columns. Packed columns are simulated with stages representing a discrete integration over the packed bed. The more stages are used the better the integration, and the more accurate the results will be (to check if the specified number of stages in a packed column simulation was sufficient, increase the number of stages and repeat the column simulation, the results should be similar). A schematic diagram of a nonequilibrium stage is shown in Figure 5.1. This stage may represent one (or more than one) tray in a trayed column or a section of packing in a packed column. The vertical wavy line in the middle of the diagram represents the interface between the two phases which may be vapor and liquid (distillation), gas and liquid (absorption) or two liquids (extraction).

Figure 5.1. Schematic diagram of a nonequilibrium stage (after Taylor and Krishna, 1993). Figure 5.1 also serves to introduce the notation used in writing down the equations that model the behavior of this nonequilibrium stage. The flow rates of vapor and liquid phases leaving the j-th stage are denoted by V_j and L_j respectively. The mole fractions in these streams are y_ij and x_ij. The N_ij are the molar fluxes of species i on stage j. When multiplied by the area available for interphase mass transfer we obtain the rates of interphase mass transfer. The temperatures of the vapor and liquid phases are not assumed to be equal and we must allow for heat transfer as well as mass transfer across the interface.

If Figure 5.1 represents a single tray then the term phi_j^L is the fractional liquid entrainment defined as the ratio of the moles of liquid entrained in the vapor phase in stage j to the moles of downflowing liquid from stage j. Similarly, phi_j^V is the ratio of vapor entrained in the liquid leaving stage j (carried down to the tray below under the downcomer) to the interstage vapor flow. For packed columns, this term represents axial dispersion. Weeping in tray columns may be accounted for with a similar term.

The component material balance equations for each phase may be written as follows:

     V            V       V                               V                V            n       V
  M    == ( 1 + r   + phi   ) V y   - V    y      - phi    V   y      - f    - sum        G       + N  
   ij            j       j     j ij    j+1  i,j+1      j-1  j-1 i,j-1    ij        nu =1   ij nu     ij

  = 0 i = 1,2,...,c

     L            L       L                               L                L            n       L
  M    == ( 1 + r   + phi   ) L x   - L    x      - phi    L   x      - f    - sum        G       - N  
   ij            j       j     j ij    j-1  i,j-1      j+1  j+1 i,j+1    ij        nu =1   ij nu     ij

  = 0 i = 1,2,...,c
where G_ij nu is the interlinked flow rate for component i from stage nu to stage j, and n is the number of total stages (trays or sections of packing). The last terms in Equations () and () are the mass transfer rates (in kmol/s), where mass transfer from the ``V" phase to the ``L" phase is defined as positive. At the V/L interface we have continuity of mass and, thus, the mass transfer rates in both phases must be equal.

The total material balances for the two phases are obtained by summing Equations () and () over the component index i.

        V            V       V                     V          V          c         n      V
  M       == ( 1 + r   + phi   ) V  - V    - phi     V    - F   - sum      sum       G      + N  
   t_{j}            j       j     j    j+1      j-1   j-1    j        i=1      nu=1   ijnu     tj

  = 0

        L          L       L                    L          L          c         n      L
  M       == ( 1+r   + phi   )L  - L    - phi     L    - F   - sum      sum       G      - N  
   t_{j}          j       j    j    j-1      j+1   j+1    j        i=1      nu=1   ijnu     tj

  = 0
F_j denotes the total feed flow rate for stage j, F_j=sum _i=1^c f_ij. Here total flow rates and mole fractions are used as independent variables and total as well as component material balances are included in the set of independent model equations. In the nonequilibrium model of Krishnamurthy and Taylor (1985a) component flow rates were treated as variables.

The nonequilibrium model uses two sets of rate equations for each stage:

     V           V
  R    == N   - N    = 0 i = 1,2 ...,c-1
   ij      ij     ij

     L           L
  R    == N   - N    = 0 i = 1,2 ...,c-1
   ij      ij     ij

where N_ij is the mass transfer rate of component i on stage j. The mass transfer rate in each phase is computed from a diffusive and a convective contribution with

     V    I   V
  N    = a   J    + y   N  
   ij      j   ij    ij  tj

     L    I   L
  N    = a   J    + x   N  
   ij      j   ij    ij  tj

where a^I_j is the total interfacial area for stage j and N_tj is the total rate on stage j (N_tj=sum ^c_i=1 N_ij). The diffusion fluxes J are given by (in matrix form):

    V     V   V    V    I
  (J ) = c  [k ]((y  - y ))

    L     L   L    I    L
  (J ) = c  [k ]((x  - x ))

where (y^V - y^I) and (x^I - x^L) are the average mole fraction difference between the bulk and the interface mole fractions (Note that the fluxes are multiplied by the interfacial area to obtain mass transfer rates). How the average mole fraction differences are calculated depends on the selected flow model. The matrices of mass transfer coefficients, [k], are calculated from

    P      P -1      P
  [k ] = [R ]  [Gamma ]
where [Gamma^P] is a matrix of thermodynamic factors for phase P. For systems where an activity coefficient model is used for the phase equilibrium properties the thermodynamic factor matrix Gamma (order c-1) is defined by
  Gamma   = delta   + x  ( (d ln gamma  / d x ) )                        
       ij        ij    i              i      j   T,P,x_k,k != j=1 ... c-1

If an equation of state is used gamma_i is replaced by phi_i. Expressions for the composition derivatives of ln gamma_i are given by Taylor and Kooijman (1991). The rate matrix [R] (orderc-1) is a matrix of mass transfer resistances calculated from the following formulae:

     P            P        c                    P
  R    = (z  / k   ) + sum             (z  / k   )
   ii      i    ic          k=1,k != i   k    ik

     P                 P            P
  R    = -z  ( (1 / k   ) - (1 / k   ) )
   ij      i         ij           ic

where k_ij^P are binary pair mass transfer coefficients for phase P. Mass transfer coefficients, k_ij, are computed from empirical models (Taylor and Krishna, 1993) and multicomponent diffusion coefficients evaluated from an interpolation formula (Kooijman and Taylor, 1991). Equations () and () are suggested by the Maxwell-Stefan equations that describe mass transfer in multicomponent systems (see Taylor and Krishna, 1993). The matrix of thermodynamic factors appears because the fundamental driving force for mass transfer is the chemical potential gradient and not the mole fraction or concentration gradient. This matrix is calculated from an appropriate thermodynamic model. The binary mass transfer coefficients are estimated from empirical correlations as functions of column internal type as well as design, operational parameters, and physical properties including the binary pair Maxwell-Stefan diffusion coefficients. Thus, the mass transfer coefficient models form the basis of the nonequilibrium model and it is possible to change the behavior of a column by selecting a different mass transfer coefficient correlation.

Note that there are c times c binary pair Maxwell-Stefan diffusion coefficients, but only c-1 times c-1 elements in the [R^P] and [k^P] matrices and, therefore, only c-1 equations per set of rate equations. This is the result of the fact that diffusion calculations only yield relative transfer rates. We will need an extra equation that will "bootstrap" the mass transfer rates: the energy balance for the interface. Note also that in this model the flux correction on the mass transfer coefficients has been neglected.

The energy balance equations on stage j are written for each phase as follows:

    V            V       V        V            V         V          V     V   VF            n       V     V
  E   == ( 1 + r   + phi   ) V  H   - V    H     - phi     V    H     - F   H    - sum        (G)     H    
   j            j       j     j  j     j+1  j+1       j-1   j-1  j-1     j   j         nu =1     jnu   jnu

      V     V
  + Q   + e   = 0
     j     j

    L            L       L        L            L         L          L     L   LF           n       L     L
  E   == ( 1 + r   + phi   ) L  H   - L    H     - phi     L    H     - F   H    - sum       (G)    H     
   j            j       j     j  j     j-1  j-1       j+1   j+1  j+1     j   j         nu=1     jnu  j nu

      L     L
  + Q   - e   = 0
     j     j

where G_jnu is the interlink flow rate from stage nu to stage j. The last term in the left-hand-side of Equations () and (), e_j, represents the energy transfer rates for the vapor and liquid phase which are defined by

    V    I   V   V    I        c       V       V
  e   = a   h  (T  - T ) + sum      N    ( H)   
   j      j                     i=1  ij      ij

    L    I   L   I    L        c       L       L
  e   = a   h  (T  - T ) + sum      N    ( H)   
   j      j                     i=1  ij      ij

where H_ij are the partial molar enthalpies of component i for stage j. We also have continuity of the energy fluxes across the V/L interface which gives the interface energy balance:

    I      V     L
  E   == e   - e   = 0
   j      j     j

where h^V and h^L are the vapor and liquid heat transfer coefficients respectively, and T^V, T^I, and T^L the vapour, interface, and liquid temperatures. For the calculation of the vapour heat transfer coefficients the Chilton-Colburn analogy between mass and heat transfer is used:
  Le = (lambda / D C  rho) = (Sc / Pr)

   V              2/3
  h  = k rho C  Le   

For the calculation of the liquid heat transfer coefficients a penetration model is used:

  h  = k rho C  sqrt(Le)

where k is the average mass transfer coefficient and D the average diffusion coefficient.

In the nonequilibrium model of Krishnamurthy and Taylor (1985a) the pressure was taken to be specified on all stages. However, column pressure drop is a function of tray (or packing) type and design and column operating conditions, information that is required for or available during the solution of the nonequilibrium model equations. It was, therefore, quite straightforward to add an hydraulic equation to the set of independent equations for each stage and to make the pressure of each stage (tray or packed section) an unknown variable. The stage is assumed to be at mechanical equilibrium so p_j^V = p_j^L = p_j.

In the second generation model, the pressure of the top tray (or top of the packing) is specified along with the pressure of any condenser. The pressure of trays (or packed sections) below the topmost are calculated from the pressure of the stage above and the pressure drop on that tray (or over that packed section). If the column has a condenser (which is numbered as stage 1 here) the hydraulic equations are expressed as follows:

  P  == p  - p  = 0
   1     c    1

  P  == p     - p  = 0
   2     spec    2

  P  == p  - p    - (Delta p   ) = 0 j = 3, 4, ..., n
   j     j    j-1           j-1

where p_c is the specified condenser pressure, p_spec is the specified pressure of the tray or section of packing at the top of the column, and Delta p_j-1 is the pressure drop per tray or section of packing from section/stage j-1 to section/stage j. If the top stage is not a condenser, the hydraulic equations are expressed as
  P  == p     - p  = 0
   1     spec    1

  P  == p  - p    - (Delta p   ) = 0 j = 2, 3, ..., n
   j     j    j-1           j-1

In general we may consider the pressure drop to be a function of the internal flows, the fluid densities, and equipment design parameters.

                                   V        L
  Delta p    = f(V   , L   , rho    , rho    , Design)
         j-1      j-1   j-1     j-1      j-1

The pressure drop term, Delta p_j-1, is calculated from liquid heights on the tray (from various correlations, see Lockett, 1986, and Kister, 1992) or specific pressure drop correlations for packings (see the section below on pressur drop models).

For bubble cap trays the procedures described by Bolles (1963) can be adapted for computer based calculation. Kister (1992) also covers methods available for estimating the pressure drop in dumped packed columns. The pressure drop in structured packed columns may estimated using the method of Bravo et al. (1986).

Phase equilibrium is assumed to exist only at the interface with the mole fractions in both phases related by:

     I          I      I
  Q    == K  x    - y    = 0 i = 1,2, ...,c
   ij      ij ij     ij

where K_ij is the equilibrium ratio for component i on stage j. The K_ij are evaluated at the (calculated) temperature, pressure, and mole fractions at the interface.

The mole fractions must sum to unity in each phase:

    V           c
  S   == sum      y   - 1 = 0
   j         i=1   ij

    L           c
  S   == sum      x   - 1 = 0
   j         i=1   ij

as well as at the interface:

    VI           c    I
  S    == sum      y    - 1 = 0
   j          i=1   ij

    LI           c    I
  S    == sum      x    - 1 = 0
   j          i=1   ij

Table 5.1 lists the type and number of equations for the nonequilibrium model. To solve the model we have 5c+6 equations and variables, where c is the number of components. They are solved simultaneously using Newton's method.

Table 5.1. Nonequilibrium model equations type and number

Equation Number
Material balances 2c+2
Energy balances 3
transfer Rate equations 2c-2
Summations equations 2
Hydraulic equation 1
interface e Quilibrium relations c
Total MERSHQ 5c+6
Nonequilibrium and equilibrium models require similar specifications. Feed flows and their thermal condition must be specified for both models, as must the column configuration (number of stages, feed and sidestream locations etc.). Additional specifications that are the same for both simulation models include the specification of, for example, reflux ratios or bottom product flow rates if the column is equipped with a condenser and/or a reboiler. The specification of the pressure on each stage is necessary if the pressure drop is not computed; if it is, only the top stage pressure needs be specified (the pressure of all other stages being determined from the pressure drop equations that are part of the model described above).

If we solve the nonequilibrium model with Newton's method, we also require initial guesses for all the variables. ChemSep uses the same automatic initial guess routine for the nonequilibrium model as the equilibrium model. The temperatures of the vapour, interface, and liquid are initialized all equal to the temperature from the automatic guess. Mass and energy transfer rates are initialized as zero and the interface mole fractions are set equal to the bulk mole fractions which are also provided by the initial guess. Pressure drops are initially assumed to be zero.

The nonequilibrium model, in comparison with the equilibrium model, requires the evaluation of many more physical properties and of the heat and mass transfer coefficients. In addition, a nonequilibrium simulation cannot proceed without some knowledge of the column type and the internals layout. Tray type and mechanical layout data, for example, is needed in order to calculate the mass transfer coefficients for each tray. For packed columns the packing type, size and material must be known. Libraries with standard tray and packing data are available on-line. Table 5.2 lists the currently supported types of column internals.

Table 5.2. Currently supported column internals for the nonequilibrium model

Bubble-cap trays
Sieve trays
Valve trays (including double weight valves)
Dumped packings
Structured packings
Equilibrium stage (with Murphree stage efficiency)
Rotating Disk Contactor (RDC) compartment (for extraction)
To avoid the problem that during the design of a column no column layout is available, the nonequilibrium column simulator has an optional design mode to automatically assign layout parameters. The user just needs to select one of the types of internals (for each section in the column). The design-mode is activated by not specifying the column diameter (leaving it as a "default" with "*") for a specific section. With the design-mode "on" each tray or stage is automatically adapted during each iteration while keeping the layout within each section the same.

For the evaluation of the heat and mass transfer coefficients, pressure drop, and the entrainment/weeping flows a nonequilibrium simulation needs the following:

Each of these are discussed in separate sections below.

5.2 Mass Transfer Coefficient Correlations

Mass transfer models are the basis of the nonequilibrium model. The models incorporated in ChemSep are all from the published literature. It is possible to change the behavior of a column by selecting a different mass transfer correlation. Therefore, we have tried to document the origin of the data of each method in order to guide you in selecting models. Table 5.3 gives a summary of the available correlations per type of internals. The various correlations are discussed below.

Table 5.3. Available mass transfer coefficient correlations per internals type

Bubble-Cap Sieve Valve Dumped Structured
tray tray tray packing packing
AIChE AIChE AIChE Onda 68 Bravo 85
Hughmark Chan-Fair Bravo 82 Bravo 92
Zuiderweg Billet 92 Billet 92
Harris Nawrocki 91
Binary mass transfer coefficients (MTC's) can be computed from Number of Transfer Units (NTU's = N) by:

   V      V       V
  k  = (N)  / t  a 

   L      L       L
  k  = (N)  / t  a 

where the vapor and liquid areas are calculated with

  a  = a  / epsilon h 
        d            f

  a  = a  / alpha h 
        d          f

the interfacial area density is computed according Zuiderweg (1982).

5.2.1 Trays

5.2.2 Random Packings

5.2.3 Structured packings

5.3 Flow Models

For the calculation of the diffusion fluxes the average mole fraction difference between the bulk and the interface mole fractions were required (see Equations and ). How these average mole fraction differences are computed depends on the selected flow model. Here three flow models are discussed: mixed flow, plug flow, and dispersion flow (which is only applied to the liquid phase).

5.3.1 Mixed flow

If we assume both phases are present in a completely mixed state, we can use

     V    I       V    I
  ((y  - y )) = (y  - y )

     I    L       I    L
  ((x  - x )) = (x  - x )
this keeps the rate equations (relatively) simple and only a function of the mole fractions leaving the current stage. However, on a tray where the vapour bubbles through a liquid which flows from one downcomer to the opposite downcomer this model is not accurate. Indeed, only for very small diameter columns will the mixed flow model give reasonable results. The mixed model is the most simple flow model and is the easiest to converge. For packed columns the convergence to the true column profiles by using increasing number of stages can be quite slow using the mixed flow model.

5.3.2 Plug flow

In the plug-flow model we assume that the vapour or liquid moves in plug flow (thus, without mixing) through the froth. This complicates the rate equations so much that no exact solution is possible. The mass transfer actually needs to be integrated over the froth. To approximate the total mass transfer an average mole fraction difference is computed. Kooijman and Taylor (1994) derived expressions for the average vapour and liquid compositions assuming constant mass transfer coefficients and that the interface compositions is constant (it isn't, but its "average" value is obtained):

     V    I               V    V    I
  ((y  - y )) = Omega[-(N) ] (y  - y )

     I    L               L    I    L
  ((x  - x )) = Omega[-(N) ] (x  - x )
where the mole fractions are of the leaving streams and the number of mass transfer units (N) for the vapour and liquid are defined as:

     V    V   V  V
  (N)  = c   k  a  h  A  / V
           t        f  b

     L    L   L  L
  (N)  = c   k  a  h  A  / L
           t        f  b

and Omega[M] is a matrix function defined as

                                -1          -1                        -1
  Omega[M] = [exp [M] - [I]] [M]   [exp [M]]   = [exp [-M] - [I]] [-M]  
Using this model predicted efficiencies for tray column experiments can be more accurately described. The plug flow model can also be used for packed columns, providing a much faster convergence to the true column profiles compared to the mixed flow model.

Currently no correction terms is applied to the plug flow model to correct for the change in mole fractions over the integration (as is discussed by Kooijman and taylor, 1994).

5.3.3 Dispersion flow

In the dispersion-flow model we assume the liquid to flow over the tray in plug flow with dispersion. Kooijman and Taylor (1994) also derived a formula to compute the average mole fraction difference for the liquid phase for this case. However, it is rather involved:

     I    L                              -1          -1                          -1          -1      -1
  ((x  - x )) = [ [p] [exp [m] - [I]] [m]   [exp [m]]   - [m] [exp [p] - [I]] [p]   [exp [p]]   ] [b]   ((X   ) / 2)

where we have defined

  a = Pe/2 = (L Z / D  W h   c )
                     e    cl

                    L             1/2
  [ b ] = a [ 2 [(N) ] / a + [I] ]   
  [ p ] = a [I] + [b]
  [ m ] = a [I] - [b]
Currently only a binary implementation is working for the dispersion model. Eddy dispersion coefficient are computed from Zuiderweg's (1982) correlation (this model is recommended by Korchinsky, 1994).

Results of the dispersion flow model are close to the plug flow model. How close depends on the eddy dispersion coefficient. Expect that the dispersion coefficient is larger for smaller diameter columns and trays with small weirs (or low liquid heights). We intend to extend the number of correlations predicting this coefficient. For now, it is advised to not use this flow model and it is not available.

5.4 Pressure Drop Models

Nowadays, there are many models and ways to compute tray pressure drops. For packings we see a shift from generalized pressure drop charts (GPDC) to more theoretically based correlations. We have chosen to employ the most recently published models. For packings we have in total 7 methods available (see table 5.4). For packing operating above the loading point (FF>0.7) we advise the use of models that take the correction for the liquid holdup into account, such as SBF-89, BS-92, and BRF-92. BRF-92 has the advantage of requiring very few fitted parameters, but is limited to structured packings.

Table 5.4. Pressure drop correlations per internals type

Bubble-Cap Sieve Valve Dumped Structured
tray tray tray packing packing
Fixed Fixed Fixed Fixed Fixed
Estimated Estimated Estimated Ludwig 79 Billet 92
Leva 92 Bravo 86
Billet 92 Stichlmair 89
Stichlmair 89 Bravo 92
Pressure drop can also be fixed to the pressure at the top of the section. However, this will can have an important effect on the designed column diameter, especially at very low pressures.

5.4.1 Tray pressure drop estimation

The liquid heights on the trays are evaluated from the tray pressure drop calculations. The wet tray pressure drop liquid height is calculated with:
  h   = h  + h 
   wt    d    l

where h_d is the dry tray pressure drop liquid height and h_l the liquid height:
  h  = h   + h  + (h   / 2)
   l    cl    r     lg

The clear liquid height, h_cl, is calculated with
  h   = alpha h  + h  
   cl          w    ow

where the liquid fraction alpha is computed with the Barker and Self (1962) correlation:
  alpha = (0.37 h  + 0.012 F  + 1.78 Q  / W  + 0.024 / 1.06 h  + 0.035 F  + 4.82 Q  / W  + 0.035)
                 w          s         L    l                 w          s         L    l

The choice of correlation for the liquid fraction turns out to be important as certain correlations are dynamicly unstable. The height of liquid over the weir, h_ow, is computed by various correlations for different types of weirs (see Perry) and a weir factor (F_w) correction (see Smith, pp. 487) is employed. For example for a segmental weir:

  h   = 0.664 F  ((Q  / W ))   
   ow          w    L    l

  w = (W  / D )
        l    c

    3    2                           2.5  2/3           2  2
  F   = w  / 1 - (F  w ((1.68 Q  / W    ))    + sqrt(1-w )) 
   w               w           L    l

where Q_L is the volumetric flow over the weir per weir length. The residual height, h_r, is only taken into account for sieve trays. Bennett's method (see Lockett, pp. 81) is:

                                      2/3                     1/3
  h  = ((6 / 1.27 rho )) ((sigma / g))    ((rho  - rho  / d ))   
   r                 L                         L      V    h

Dry tray pressure, h_d, is calculated with:

  h  = K (rho  / rho ) u  
   d         G      L   h

  K = (xi / 2 g)
where the orrifice coefficient xi for sieve trays is computed according to Stichlmair and Mersmann (1978). For valve trays we use the method of Klein (1982) as described in Kister (1992, pp. 309--312) where K is given for the cases with the valves closed or open. It is extended for double weight valve trays as discussed by Lockett (1986, pp. 82--86). The dry tray pressure drop is corrected for liquid fractional entrainment.

The froth density is computed with

  h  = (h   / alpha)
   f     cl

The liquid gradient, h_lg, is computed according to Fair (Lockett, 1986, pp. 72):
  R  = (W h  / W + 2 h )
   h       f          f

  U  = (Q  / W h  )
   f     L      cl

  Re  = (R  U  rho  / eta )
    f     h  f    L      L

          4       -1.06
  f = 7 10  h  Re      
             w   f

  h   = (Z f U   / g R )
   lg         f       h

where W is the average flow-path width for liquid flow, and Z the flow path length. The height of liquid at the tray inlet is:

                                                  2                                                         2
  h = sqrt( (2 \over g) \left((Q  \over W )\right)  \left((1 \over h  ) - (1 \over h )\right) + (2 \alpha h   \over 3) )
   i                            L        l                          cl              c                      f

where h_c is the height of the clearance under the downcomer. The pressure loss under downcomer (expressed as a liquid height) is

  h    = ((1 / 2 g)) ((Q  / C  W  h )) 
   udc                  L    d  l  c

where C_d=0.56 according to Koch design rules. The height of liquid in the downcomer can now be calculated with the summation:
  h   = h   + h  + h   
   db    wt    i    udc

Bubble-cap liquid heights are done according to Perry's (1984) and Smith (1963). Additionally the liquid fraction of the froth is computed according to Kastanek (1970).

5.4.2 Random packing pressure drop correlations

For packings the vapour and liquid mass flow per cross sectional area (kg/m^2 s) and velocities (m/s) are:

  L  = L M  / A 
   a           t

  V  = V M  / A 
   a           t

  u  = L  / rho 
   L    a

  u  = V  / rho 
   V    a

5.4.3 Structured packing pressure drop correlations

5.5 Entrainment and Weeping

Entrainment and weeping flows (for trays only) change the internal liquid flows and influence the performance of the column internals. ChemSep currently does not support the handling of these flows. This is due to the fact that few entrainment models behave properly. Neither is the effect of the entrainment and weepflows on the mass transfer properly taken into account.

Entrainment is computed from the fractional liquid entrainment which is computed from Hunt's correlation and from figure 5.11 of Lockett (1986) for sieve trays:

     L          -5                                                3.2
  phi  = 7.75 10   ( (0.073 / sigma) ) M  ( (U  / T  - 2.5 h   ) )   
                                        v     v    s        cl

The weeping factor is estimated from a figure from Smith (1963, plot on page 548), which was fitted with the following correlation
  WF = (0.135 phi ln (34 (H +H  ) + 1) / (H  + H ))
                           w  ow           d    r

where phi is the open area ratio.

5.6 The Design Mode

The initial layout is determined after the flows are known from the initial guess. Each stage in the column is designed separately and independently of adjacent stages. Then the sections in the column are rationalized so that trays or stages within a section have the same layout. During each iteration (that is, an update of the flows) each stage is re-designed only if the flowrates have changed more than by a certain fraction (which can be specified). Only sections with re-designed trays or stages are rationalized again. After convergence a complete design of any trayed or packed section in the column is obtained. In this manner trayed and packed sections can be freely mixed in a column simulation/design.

Different design methods can be employed:

The methods generate a column-design that might not be optimal from an engineers viewpoint. They must be seen as starting points for the actual design layouts. Also, the design does not include constructional calculations to determine tray support constructions or thicknesses of trays or the column. Design mode is automatically triggered if the column diameter is not specified. Other layout parameters can be specified but they may be changed by the design mode. Each of these methods behaves differently and they are discussed in more detail below. An additional and very important de-rating factor is the system factor (SF). It represents the uncertaincy in design correlations with regard to phenomena which are currently still not properly modeled, such as foaming.

Tray layout parameters that specify a complete design (for the calculation of mass transfer coefficients and pressure drops) are shown in Table 5.5. For packings only the column diameter and bed height are design parameters, other parameters are fixed with the selection of the type of packing (such as void fraction, nominal packing diameter, etc.). The packed bed height must be specified since it determines the desired separation and the capacity.

Table 5.5. Tray layout data

General (sieve) tray layout data:
Column diameter Active area
Number of flow passes Total hole area
Tray spacing Downcomer area
Liquid flow path length Weir length
Hole diameter Weir height
Hole pitch Deck thickness
Downcomer clearance
Additional data for bubble caps:
Cap diameter Slot area
Slot height Riser area
Skirt clearance Annual area
Additional data for valves:
Closed Loss K Open Loss K
Eddy Loss C Ratio Valve Legs
Valve Density Valve Thickness
Fraction Heavy Valves Heavy Valve Thickness

5.6.1 Tray Design: Fraction of flooding

The first task in this approach to tray design is to assign all layout parameters to consistent values corresponding to the required capacity defined by the fraction of flooding and current flowrates. These defaults function as starting points for subsequent designs.

The initial free area ratio is taken to be 15 % of the active area. The active area is determined with capacity factor calculation with internals specific methods (for sieve and bubble-cap trays the default is Fair's correlation by Ogboja and Kuye (19), and the Glitsch method for valve trays). The tray spacing is initially set to the default value (of 0.5 m) and the downcomer area is calculated according the Glitsch manual (limited by a minimum time residence check). From the combined areas the column diameter is computed. The number of liquid passes on a tray is initially set by the column diameter; under 5 ft one pass, under 8 ft two, 10 ft three, under 13 ft four, else five passes. With the number of passes and the column diameter the total weir length is computed. Once the weir length is determined the liquid weir load is checked, if too high the number of passes is incremented and a new weir length is evaluated until the weir load is below a specified maximum.

Initial weir height is taken as 2", but limited to a maximum of 15 % of the tray spacing. For notched or serrated weirs the notch depth is a third of the weir height. For serrated weirs the angle of serration is 45 degrees. Circular weirs have diameters 0.9 times the weir length. Hole diameter is set to 3/16" for sieve trays and tray thickness 0.43 times the hole diameter (or 1/10"). The hole pitch is computed from the free area ratio and hole diameter according to a triangular pitch. The default downcomer clearance is 1.5" but is limited by the maximum allowed downcomer velocity according to the Glitch method de-rated with the system factor. The clearance is set to be at least half an inch lower than the weir height to maintain a positive liquid seal but is limited to a minimum of half an inch.

For bubble-cap trays the cap diameter is 3" for column diameters below 4.5 ft and 4" for above. The hole diameter can vary between 60 % to 71 % of the capdiameter, and default taken as 70 %. Default skirt clearance is 1" with minimum of 0.5" and maximum of 1.5". slot height can vary in between 0.5" and 1.5", default 1" for cap diameters below 3.5" and 1.25" for larger cap diameters. The pitch can vary from 1.25" to half the flow path length (minimum number of rows is two), default set to 1.25".

Valve trays are initialized to be Venturi orifice uncaged, carbon steel valves of 3 mm thick with 3 legs (see Kister, 1992, p312). The hole diameter is 1" for column smaller than 4.5 ft, otherwise 2". No double weight valves are present.

The second task in the fraction of flooding method consists of finding the proper free area ratio (beta = A_h / A_b = hole area / active area) so that no weeping occurs. This ratio can vary between a minimum of 5% (for stable operation) and a maximum of 20%. To test whether weeping occurs, we use the correlation by Lockett and Banik (1984): Fr_hole>2/3. The method requires all liquid heights to be evaluated at weep rate conditions. This task is ignored for bubble-cap trays. The weep test is done at weeping conditions, with a weep factor at 60 % (this can be changed). Calculating liquid heights is done by adding various contributions with correlations from Lockett (1986) and Kister (1992), see Appendix A. If weeping occurs at the lower bound for the free area ratio, a flag is set for the final task to adapt the design. The final task consists of evaluating all liquid heights at normal conditions and to do a number of checks:

If a check fails the design is adapted to correct the problem, according to the adjustments shown in Table 5.6 after which new areas are calculated with capacity correlations. Part of this task is also to keep the layout parameters that are adjusted within certain lower and upper bounds to maintain a proper tray design. Finally the number of iterations for the design method is checked against a maximum (default 30) to prevent a continuous loop.

Table 5.6. Tray design checks and adjustments

Problem Test Adjustments
Bubble cap vapor distribution h_lg/h_d > 0.5 p+f_1,
Weeping Fr_h / (2/3) < 1-f_a
free < 0.05
A_b < A_bf: A_b=A_bf
else: A_b-f_1
t_v+f_2 (vt)
Hydrodynamic (downcomer) flooding) T_s < h_db/FF T_s+f_1
Excessive liquid entrainment A_b+f_1
Froth height limit h_f > 0.75 T_s A_b+f_1
Excessive pressure drop g rho h_wt > Delta p_max A_b+f_1
p+f_1 (bc)
h_skirt+f_2 (bc)
h_slot+f_3 (bc)
Excessive vapor entrainment A_d+f_1
The adjustment factors f_1, f_2, and f_3 are percentual in/decrements, normally set at 5, 2, and 1 %. These factors -- together with all the default, lower, and upper settings that are used in the design routine -- are stored in a ``design file'' (TDESIGN.DEF) that can be tailored to handle specific kinds of designs and columns. This allows the selection of different methods for capacity and hydrodynamic calculations as well. Also the fraction that the flows need to change before a re-design is issued can be changed in this manner together with other design criteria. The design file must be in the current directory for the nonequilibrium program to use it, otherwise the normal defaults will be used.

Here we discuss the most important parameters of the file. The file starts with a comment on the first line. The second line specifies the factors f_1, f_2, and f_3 for adapting the design layout parameters. The third line specifies the fraction of change allowed in the flows before a redesign occurs. It also specifies the fraction of deviation allowed in downcomer and bubbling area between current and design values. Line 11 specifies the volumetric weir load after which the number of passes is incremented. Line 14 specifies the maximum froth height as fraction of the tray spacing as is used in the froth height limit check. Line 15 specifies the criterium to which the free area ratio has to conmverge. Line 16 sets the maximum allowed pressure drop for the excessive pressure drop check. Line 20 specifies the maximum number of loops for the design method. Line 21-23 specify the methods to calculate capacity factors for bubble-cap, sieve, and valve trays. Line 24-25 set the downcomer are method and velocity check. Line 45 sets a flag to generate tray parameter output and line 46 sets a flag for intermediate design messages

5.6.2 Packing Design: Fraction of flooding

For packed columns only the column diameter is a design parameter to be evaluated. Default packing data are used for all packing parameters that are not specified; values of 1" inch metal Pall rings for random packed sections and of Koch Flexipack 2 (316ss) for structured sections.

To determine the packed column diameter, the diameter that gives rise to the flooding pressure drop (as specified) is computed using the selected pressure drop model. The resulting diameter is corrected for the fraction of flooding and the system factor:

  D  = (D        / sqrt(FF~SF))
   c     c,flood

This does make the resulting column diameter depend on the selected pressure drop model. If no pressure drop model is selected the Leva (1992) model is selected (which is only a function of the packing factor). If no pressure drop at flood is specified, it is estimated with Kister's correlation (1992) (which is only a function of the packing factor). Thus, as long as the packing factor is known, this method will not fail.

5.6.3 Pressure drop

Tray design on pressure drop works as discussed above but with a default fraction of flooding of 75 %. However, the specified pressure drop functions as a maximum allowed pressure drop per tray. No adjustment is done if the pressure drop is below this specified pressure drop.

Packing design automatically finds the diameter resulting in the specified pressure drop (with the selected pressure drop model). This is done by using a linear search technique as the different packing pressure drop correlations can behave quite irregularly. The maximum allowed pressure drop is the flooding pressure drop as specified or computed from Kister's correlation and the packing factor. If the pressure drop is specified to be very low the column diameter might converge to unrealistic diameters. A zero or larger than flooding pressure drop specification results in a 70 % fraction of flooding design.

5.6.4 Optimizing

This tray design only method tries to optimize the tray design for the following four aspects: However, this partical design mode is not yet available.

Symbol List

a_d      Interfacial area density (m^2/m^3)
a^I      Interfacial area (m^2)
A_h      Hole area (m^2)
A_b, A_bub  Bubbling area (m^2)
A_d      Downcomer area (m^2)
c        Number of components,
           Molar concentration (kmol/m^3)
d_h      Hole diameter (m)
D        Binary diffusivity coefficient (m^2/s)
D_c      Column diameter (m)
D_e      Eddy dispersion coefficient (m^2/s)
e        Energy transfer rate (J/s)
f_ij   Component i feed flow to stage j (kmol/s)
f_1, f_2, f_3  Design adjustment factors
F_j      Total feed flow rate to stage j (kmol/s)
F_p      Packing factor (1/m)
F_s      F factor F_s=U_v sqrt(\rho_V) (kg^0.5/m^0.5s)
FF       Fraction of flooding
FP       Flow parameter FP=M_L/M_V sqrt(\rho^V_t/\rho^L_t)
 Fr  Froude number
g        Gravitational constant, 9.81 (m/s^2)
G        Interlinked flow rate (kmol/s)
h        Heat transfer coefficient (J/m^2 K s)
h_c      Clearance height under downcomer (m)
h_cl   Clear liquid height (m)
h_d      Dry tray pressure drop height (m)
h_db   Downcomer backup liquid height (m)
h_f      Froth height (m)
h_i      Liquid height at tray inlet (m)
h_lg   Liquid gradient pressure drop height (m)
h_l, h_L  Liquid pressure drop height (m)
h_ow   Height of liquid over weir (m)
h_r      Residual pressure drop liquid height (m)
h_wt   Wet tray pressure drop liquid height (m)
h_w      Weir height (m)
h_udc  Liquid height pressure loss under downcomer (m)
H        Molar enthalpy (J/kmol)
 H_i  Partial molar enthalpy of component i (J/kmol)
J        Molar diffusion flux (kmol/m^2 s)
k        Binary mass transfer coefficient (m/s)
K_i      K-value or equilibrium ratio component i: K_i=y_i/x_i
L        Liquid flow rate (kmol/s)
 Le  Lewis number ( Le= Sc/ Pr)
M        Mass flow rate (kg/s)
N        Mass transfer rate (kmol/s)
n        Number of stages
p        Hole pitch (m),
           Pressure (Pa)
Delta p  Pressure drop (Pa)
Delta P_max  Maximum design pressure drop (Pa/tray or Pa/m)
 Pr  Prandtl number
Q        Heat input (J/s)
Q_L      Volumetric flow over the weir (m^3/s)
r        Ratio sidestream to internal flow
[R]      Matrix defined by () and ()
 Sc  Schmidt number
SF       System derating factor
t        Residence time (s)
t_v      Valve thickness (m)
T        Temperature (K)
T_s      Tray spacing (m)
V        Vapor flow rate (kmol/s)
We       Weber number
W_l      Weir length (m)
x        Liquid mole fraction
y        Vapor mole fraction
z        Mole fraction
alpha   Fraction liquid in froth
beta    Fractional free area beta=A_h/A_b,
phi     Fractional entrainment
rho     Density (kg/m^3)
sigma   Surface tension (N/m)
eta     Viscosity (Pa s)
[Gamma]  Thermodynamic matrix
lambda  Heat conductivity (W/m/K)
I        Interfacial
L        Liquid
P        Phase P
V        Vapor
flood    at flooding conditions
i        component i
j        stage j,
           component j
spec     specified
t        total
nu      from interlinking stage nu


P.E. Barker, M.F. Self, ``The evaluation of Liquid Mixing Effects on a Sieve Plate using Unsteady and Steady-State Tracer Techniques'', Chem. Eng. Sci., Vol. 17 (1962) pp. 541.

S.D. Barnicki, J.F. Davis, ``Designing Sieve-Tray Columns, Part 1: Tray Design'', Chem. Engng., Vol. 96, No. 10, (1989) pp. 140--146.

S.D. Barnicki, J.F. Davis, ``Designing Sieve-Tray Columns, Part 2: Column Design and Verification'', Chem. Engng., November, pp. 202--212 (1989).

D.L. Bennett, R. Agrawa, P.J. Cook, ``New Pressure Drop Correlation fo Sieve Tray Distillation Columns'', AIChE J., Vol. 29, No. 3 (1983) pp. 434.

R. Billet, M. Schultes, "Advantage in correlating packed column perfomance", IChemE. Symp. Ser. No. 128, B129 (1992).

R. Billet, Distillation Engineering?, Heyden (1979?).

W.L. Bolles, in B.D. Smith, Design of Equilibrium Stage Processes, Chap. 14, McGraw-Hill (1963).

J.L Bravo, J.R. Fair, Ind. Eng. Chem. Proc. Dev., 21, 163 (1982).

J.L. Bravo, J.A. Rocha, J.R. Fair, "Mass transfer in gauze packings", Hydrocarbon Processing, January, 91 (1985).

J.L. Bravo, J.A. Rocha, J.R. Fair, "Pressure drop in structured packings", Hydrocarbon Processing, March (1986).

J.L. Bravo, J.A. Rocha, J.R. Fair, "A comprehensive model for the performance of columns containing structured packings", IChemE. Symp. Ser. No. 128, A439 (1992).

Chan, J.R. Fair, ``Prediction of point efficiencies on sieve trays'', Ind. Eng. Proc. Des. Dev., Vol. bf 23, 814 (1984)

Chen and Chuang, Ind. Eng. Chem. Res, Vol. 32, 701--708 (1993).

Gerster et al., AIChE J., (1958).

I.J. Harris, `` Optimum Design of Sieve Trays'', Brit. Chem. Engng, Vol. 10, No. 6 (1965) pp. 377.

G.A. Hughmark, ``Mdels for Vapour Phase and Liquid Phase Mass Transfer on Distillation Trays'', AIChE J., Vol. 17, No. 6 (1971) pp. 1295.

F. Kastanek, ``Efficiencies of Different Types of Distillation Plate'', Coll. Czech. Chem. Comm., Vol. 35 (1970) pp. 1170.

H.Z. Kister, Distillation Design, McGraw-Hill, New York (1992).

G.F. Klein, Chem. Engng, May 3, 81 (1982).

H.A. Kooijman, R. Taylor, ``On the Estimation of Diffusion Coefficients in Multicomponent Liquid Systems'', Ind. Eng. Chem. Res., Vol 30, No. 6, (1991) pp. 1217--1222.

W.J. Korchinsky, ``Liquid Mixing in Distillation Trays: Simultaneous Measurement of the Diffusion Coefficient and Point Efficiency'', Trans. I. Chem. E., Vol. 72, Part A, (1994) 472-478.

R. Krishnamurthy, R. Taylor, ``A Nonequilibrium Stage Model of Multicomponent Separation Processes. Part I: Model Description and Method of Solution'', AIChE J., Vol. 31, No. 3 (1985), pp. 449--455.

Leva, (1953?).

M. Leva, ``Reconsider Packed-Tower Pressure-Drop correlations", Chem. Eng. Prog. January, 65 (1992).

M.J. Lockett, Distillation Tray Fundamentals, Cambridge University Press (1986).

M.J. Lockett, S. Banik, ``Weeping from Sieve Trays'', AIChE Meeting, San Francisco, Nov. (1984).

E.E. Ludwig, Applied Process Design for Chemical and Petrochemical Plants, Vol. 2, 2nd Ed., Gulf Pub. Co., Houston, TX, (1979).

O. Ogboja, A. Kuye, ``A procedure for the design and optimisation of sieve trays'', Trans. I. Chem. E., Vol. 68, Part A, Sept. (1990) pp. 445-452.

K. Onda, H. Takeuchi, Y. Okumoto, "Mass transfer coefficients between gas and liquid phases in packed columns", J. Chem. Eng. Jap., Vol. 1, No.1, 56 (1968).

R.H. Perry and D. Green, Perry's Chemical Engineering Handbook, 6th edition, section 18, Liquid-Gas System, 18-8 -- 18-12 (1984).

M. Prado, The bubble-to-Spray Transition on Sieve Trays: Mechanisms of the Phase Inversion, Ph.D. thesis, University of Texas, Austin (1986).

B.D. Smith, Design of Equilibrium Staged Processes, McGraw-Hill, New York (1963)

J. Stichlmair, A. Mersmann, ``Dimensioning Plate Columns for Absorption and Rectification'', Chem. Ing. Tech., Vol. 45, No. 5 (1978) pp. 242.

J. Stichlmair, J.L. Bravo, J.R. Fair, Gas. Sep. Purif., Vol. 3, 19 (1989).

R. Taylor, H.A. Kooijman, ``Composition derivatives of Activity Models (for the estimation of Thermodynamic Factors in Diffusion)'', Chem. Eng. Comm., Vol. 102 (1991) pp. 87--106.

R. Taylor, H.A. Kooijman, J-S. Hung, ``A second generation nonequilibrium model for computer simulation of multicomponent separation processes'', Comput. Chem. Engng., Vol. 18, No. 3, pp. 205--217 (1994).

R. Taylor, R. Krishna, Multicomponent Mass Transfer, Wiley, New York (1993).

P.C. Wankat, Separations in Chemical Engineering - Equilibrium Staged Separations, Elsevier, 420--428 (1988).

F.J. Zuiderweg, ``Sieve Trays - A View of the State of the Art'', Chem. Eng. Sci., 37, 1441--1461 (1982).

Chapter 6. Nonequilibrium Extraction

This chapter deals especially with the application of the nonequilibrium model to the modelling of extraction columns. In such operations the two phases present are both liquids instead of a liquid and a vapor as in the case of distillation, stripping, or absorption. This requires fundamentally different mass transfer coefficients and flow models, as well as a completely new design method, that an entire chapter is devoted to the subject.

6.1 Introduction

Nonequilibrium extraction uses the same model as described in the nonequilibrium section, with the exception that there is no vapor. Instead we have a light and a heavy liquid phase, where the light liquid behaves as the vapor with, of course, liquid-like properties. If the heavy phase (L) is lighter than the light phase (V) the program stops. However, either phase (that is, L or V) can be the disperse phase. The user must specify which is the disperse phase, since this changes the design and the calculation of MTC's. Currently sieve trays, structured and random packed columns, rotating disc contacters, and spray columns are supported as internals (as well as equilibrium stages with a specified stage efficiency). The K-values must be the Liquid-Liquid model, which uses activity coefficients. The energy balance can be ignored (Enthalpy=None) or included. In case it is ignored the column temperature is dictated by that of the feeds, and linear interpolation is used to provide a column temperature profile. A specific temperature profile can be imposed if the energy balance is ignored and user temperature initialization is supplied. Default values for the total interfacial area and mass transfer coefficients are: A_i=100 m^2, k_d=10^-5 and k_c=10^-4 m/s. Mass transfer in coelescencing layers and jet zones are neglected (they could be modeled by a special stage for packed/RDC columns). Thus, only the drop rise zone is taken into account for mass transfer.

Current limitations of the nonequilibrium extraction model are:

6.2 Sieve trays

ChemSep will attemp to design the extraction column if no design is specified, this design method is adapted from the notes by R. Krishna.

6.2.1 Design

The default free area ratio is 5 %, tray spacing is 0.4 meter, and the clearance under the downcomer a quarter of the tray spacing. There is no weir. The hole diameter is set by default to:
  x = sqrt(\sigma \over \Delta \rho g)
  d  = 1.8 x

but d_h is limited (if supplied) by:
  0.5 x < d  < pi x

and the practical limits (overriding):
  3 mm < d  < 8 mm

The hole velocity is computed with:

  (Eo) = (Delta rho g d   / sigma)

  (We) = 4.33 (Eo)     
  U  = sqrt( We \sigma \over \rho  d )
   h                             d  h

If the hole velocity is less than 0.15 (m/s) then its design value is kept at 0.15 (m/s). The Froude number is computed from

  (Fr) = (U   / g d )
            h      h

For Eo is less than 0.4 the Sauter mean droplet diameter is computed by:

           -0.4                              0.67
  d  = (Eo)     ( 2.13 ( (Delta rho / rho ) )     + exp (-0.13 (Fr)) ) d 
   p                                     d                              h


           -0.42                    0.42
  d  = (Eo)      ( 1.24 + exp (-(Fr)    ) ) d 
   p                                         h

The hole area is
  A  = (Q  / U )
   h     d    h

The ratio of the hole area over the active area (free area ratio, f) is limited between 1 and 20 %.
  A  = A  / f
   a    h

The hole pitch can be computed if the hole diameter and free area ratio are known. The downcomer velocity can be computed if a minimum droplet diameter, d_min, is assumed which will not be entrained. The downcomer velocity equals the velocity of the continuous phase, U_c:

                                  2                0.33
  U  = 0.249 d    ( ((g Delta rho)  / rho  eta ) )     
   c          min                        c    c

This droplet diameter is taken to be 0.5 mm. With U_c known we can compute the downcomer area:
  A  = Q  / U 
   d    c    c

The total area is equal to two downcomer areas plus the active area and 0.5 % area for support etc.:
  A  = (A  + 2 A ) / 0.995
   t     a      d

With the total tray area known the column diameter can be computed. The net area for the disperse phase, A_n, and the disperse velocity, U_d, are:
  A  = A  + A 
   n    A    d

  U  = (Q  / A )
   d     d    n

Next the dispersed phase velocity holdup and slip velocity are computed. The slip velocity (V_s) is guessed at one tenth of the disperse phase velocity, making the disperse phase holdup equal to a tenth since it is defined as
  phi  = (U  / V )
     d     d    s

The slip velocity (which is a function of the dispersed phase holdup and needs to be obtained iteratively) can be calculated from:

                                                                                                      0.33           1.834
  V  = sqrt( 2.725 g d  \left( (\Delta \rho \over \rho ) \right) \left( ( (1 - \phi ) \over (1 + \phi     ) ) \right)      )
   s                  p                               c                            d                 d

After the dispersed phase holdup is computed (it depends on V_s) it is checked to be within 1 and 20 % for standard operation conditions. If it is too small the free area ratio is increased, if it is too large the free area ratio is decreased (each by 5 %) till it is within the desired range.

The Weber number

  (We) = rho  U   d  / sigma
            d   h  p

must be larger than 2 to ensure all holes produce drops (i.e. to avoid inactive holes, see Seibert and Fair, 1988).

The height of the coalesced layer is (according to Treybal, 1980):

          2     2               2
  h  = ((U   - U  ) rho  / 2 g C   Delta rho) + (4.5 U  rho  / 2 g Delta rho) + (6 sigma / d  g Delta rho)
   c       h     d     d         d                    c    c                                p

(with C_d=0.67). The first term is height to overcome flow through the orrifices, the second for friction losses due to contraction/expension on entry/exit (0.5 + 1.0) and change of direction (2 times 1.5 velocity heads), and the third term for the interfacial tension effects at the holes. The height needs to be larger than 4 cm (to ensure safe operation). If not, then the hole diameter is decreased by 5 % and we repeat the procedure from the hole velocity calculation (). This design is for a one pass sieve tray, and flow path length, L_f is computed from geometric relations. The weir length is (segmental downcomer):
  W  = A  / L 
   l    a    f

The tray thickness is defaulted to a tenth of an inch. To prevent entrainment of droplets, the flow under the downcomer is only allowed to be 50 % higher than the downcomer velocity. If higher, then the downcomer clearance is enlarged until this requirement is met. The tray spacing is adjusted so that the coeleseced layer and coalescence zone divided over the length of the downcomer equals the fraction of flooding (multiplied with the system factor).

6.2.2 Report

The reported fraction of flooding equals to the ratio of the height of the coelesced layer over the height of the downcomer (according to Seibert and Fair the flooding calculation is within 20 %). The lower operating limit is the ratio of two over the Weber number (to guarantee proper droplet formation).

6.2.3 Mass Transfer Coefficients

The "Handlos-Baron-Treybal" method is used. The hole velocity U_h, Eo, Fr, net area A_n, Sauter mean drop diameter d_p, disperse velocity U_d, slip velocity V_s, disperse phase holdup phi_d, h_c, and h_z are computed as above (but with fixed design parameters). The interfacial area per unit of volume is
  A  = (6 phi  / d )
   i         d    p

and the drop rising zone (where mass transfer is assumed to take place):
  h     = t  - h 
   drop    s    c

where t_s is th etray spacing. The volume for mass transfer on a tray is
  V  = A  h    
   i    n  drop

The mass transfer coefficients for transport from the disperse phase are (Handlos and Baron, 1957):
  k  = (0.00375 V  / ( 1 + eta  / eta  ) )
   d             s            d      c

and for transport from the continuous phase are (Treybal, 1963):

                     -0.43                    -0.58
  k     = 0.725 (Re)       (1 - phi ) V  (Nu)      
   c,ij             c              d   s     c

  (Re)  = d  V  rho  / eta 
      c    p  s    c      c

  (Nu)  = eta  / rho  D   
      c      c      c   ij

Note that k_d is not a function of the diffusion coefficient and, thus, is the same for all components.

6.3 Packed columns

Column design and calculation of mass transfer coefficients is done the same way for structured packed column and random packed columns, following the methods and correlations as outlined by Seibert and Fair (1988).

6.3.1 Design

For mass transfer from the continuous phase to the disperse phase we have x=1 for the calculation of the Sauter mean drop diameter:
  d  = 1.15 x sqrt( \sigma \over \Delta \rho g)

The slip velocity of a single droplet at zero disperse phase holdup is given by

  V   = sqrt( 4 \Delta \rho g d  \over 3 \rho  C )
    s                          p             c  d

where Cd=0.38 (for high values of Reynolds). Static disperse holdup is:
  phi   = (0.076 a  d  / xi)
     ds           p  p

where a_p is the packing area and xi the packing void fraction. The static holdup area and total area are:
  a  = 6 0.076 a 
   s            p

  a = a  + a 
       p    s

The tortuosity is defined as
  zeta = (a d  / 2)

The superficial velocity of the continuous phase at the flood point is
  e = cos ((pi zeta / 4))

                     o                         2
  U   = (0.192 xi * V   / (1.08 + (Q  / Q ) / e ))
   cf                 s             d    c

This needs to be corrected for the fraction of flooding (and system factor):
  U  = SF FF U  
   c          cf

to give the net area
  A  = (Q  / U )
   n     c    c

from which the packed column diameter (D_c) can be calculated.

6.3.2 Report

The reported fraction of flooding is the quotient of computed U_c to U_cf as discussed above. The dispersion coefficients are given by (Vermeulen et al., 1966):
  log ( E  / V  d  ) = 0.046 (V  / V ) + 0.301
         d    d  p             c    d

  log ( E  / V  d  ) = 0.161 (V  / V ) + 0.347
         c    c  p             c    d

where d_p is the packing diameter.

6.3.3 Mass Transfer Coefficients

The method of Seibert and Fair (1988) is used. The phase velocities are computed by
  U  = (Q  / A )
   c     c    n

  U  = (Q  / A )
   d     d    n

The drop diameter d_p, slip velocity V^o_s, area a, static holdup area a_s, and tortuosity zeta are calculated as above. Then the disperse phase holdup, phi_d, is determined iteratively (starting at 0.1.) from:
  f(phi ) = exp ( (-6 phi  / pi))
       d                 d

                    o                 2
  phi  = (U  / xi (V   f(phi ) - U ) e )
     d     d         s      d     c

Then the slip velocity is

  V  = V   f(phi ) e + (1 - e) U 
   s     s      d               c

since U_d=V^o_s f(phi_d). The mass transfer coefficient for the disperse phase is computed by:
  phi = (sqrt(Sc ) / ( 1 + eta  / eta  ) )
                d             d      c

  phi > 6 : k     = (0.023 V  / sqrt((Sc) ))
             d,ij           s            d

  phi < 6 : k  = (0.00375 V  / ( 1 + eta  / eta  ) )
             d             s            d      c

  (Sc)  = ( eta  / rho  D     )
      d        d      d  d,ij

If phi is larger than 6 the Laddha and Degaleesan correlation is used otherwise the Handlos-Baron method. For the mass transfer coefficient in the continuous phase:

                    0.5      0.4
  (Sh)  = 0.698 (Re)     (Sc)     (1-phi )
      c                c        c       d

  k     = ( (Sh)  D     / d  )
   c,ij         c  c,ij    p

  (Re)  = ( rho  V  d  / eta  )
      c        c  s  p      c

  (Sc)  = ( eta  / rho  D     )
      c        c      c  c,ij

The interfacial area per unit volume is:
  a  = ( 6 xi phi  / d  )
   i             d    p

The total interfacial area in a stage is the stage height times the net area times the interfacial area per unit volume:
  a      = a  A  h     
   i,tot    i  n  stage

6.4 Rotating Disk Contactors

This design method is based on the Handbook of Solvent Extraction (chapter 13.1) and notes by R. Krishna.

6.4.1 Design

The phase ratio alpha is
  alpha = ( Q  / Q  )
             d    c

The maximum stable drop diameter is

                            5/21      6/21       10/21     1/21
  u  = 0.9 ( ( g Delta rho )     sigma     / rho       rho      )
   0                                            c         d

  d      = (sigma / rho  u  )
   p,max               c  0

A stable drop diameter is selected as half of the maximum diameter
  d  = 0.5 d     
   p        p,max

and the require power input (P_i=N^3 R^5 / H D^2) is computed

                               0.6        2.5
  e = ( (0.25 ((sigma / rho ) )    / d ) )   
                           c          p

  P  = (pi e / 4 C )
   i              p

(C_p=0.03 for Re>10^5). If no column diameter is known, an estimate is made from assuming a cross-sectional area for a combined velocity of 0.05 m/s with:
  A  = (Q  + Q ) / 0.05
   c     c    d

The required rotational speed (using these standard ratios) is then

                          5        2    0.33
  N = ( ( P  ((0.1 / (0.6) ) ) / D   ) )    
           i                      c

Now the slip-velocity can be calculated using a correlation from Kung and Beckman (1961):

                                             0.9          2.3          0.9          2.6          2
  (V  * eta  / sigma) = (( Delta rho / rho ))    ((S / R))    ((H / R))    ((R / D))    ((g / R N ))
    s      c                              c

The disperse holdup at flood is determined from

  phi  = (sqrt(\alpha +8 \alpha) - 3 alpha / 4 ( 1 - alpha))

from which the continuous phase velocity at flood can be determined with

  U    = V  (1 - phi )  (1 - 2 phi )
   c,f    s         d             d

Correction for fraction of flooding (and system factor) gives
  U  = SF FF U  
   c          cf

from which the column area and diamater can be calculated
  A  = ( Q  / U  )
   c      c    c

The rotor diameter R, stator diameter S, and the height of the compartment have standard ratios with respect to the column diameter (D_c)
  R = 0.6 D 

  S = 0.7 D 

  H = 0.1 D 

so the size of the column is determined. Below a Renolds number of 10^5 C_p becomes a function of the Renolds number. Normally RDC's are operated in the regime above 10^5 so the Renolds number is computed by

  Re  = (rho  N R  / eta )
    d       d           d

and a smaller diameter is selected (and the calculations repeated) if necessary. On re-design the layout of the stage with the largest diameter is used for the entire section.

6.4.2 Report

The reported fraction of flooding is the quotient of computed U_c over U_cf as discussed above. The operating velocity is proportional to the slip velocity and so inverse proportional to the square of the rotation speed. One of the design rules was to keep the disperse Reynolds number larger than 10^5 so the lower operating limit is defined as:

       5         2
  (( 10  / Re  )) 

Stemerding et al. (1963) gave a correlation for the axial dispersion coefficient for the continuous phase

  (E  / V  H) = 0.5 + 0.012 N R (S/D)  / V 
    c    c                                c

The disperse dispersion coefficient is set to twice this number.

6.4.3 Mass Transfer Coefficients

The method of "Kronig-Brink-Rowe" is used. Phase ratio alpha, energy input P_i (from N, R, H, and D_c) are computed as above. The drop diameter is computed from
  C  = 0.03

  e = (4 C  P  / pi)
          p  i

                         0.6       0.4
  d  = 0.25 ((sigma / e))    / rho    
   p                              c

The dispersed holdup phi_d is calculated iteratively as above and the slip velocity is determined as described above (with ). The mass transfer coefficients are:
  (Sh)  = 10.0

  k     = ( (Sh)  D     / d  )
   d,ij         d  d,ij    p

                       0.62      0.36
  (Sh)  = 2 + 0.42 (Re)      (Sc)     
      c                    c         c

  k     = ( (Sh)  D     / d  )
   c,ij         c  c,ij    p


                   1.33  0.33
  (Re)  = ( rho  d      e     / eta  )
      c        c  p                c

  (Sc)  = ( eta  / rho  D     )
      c        c      c  c,ij

The interfacial area per unit volume is
  a  = ( 6 phi  / d  )
   i          d    p

Alternatively the "Rose-Kintner-Garner-Tayeban" method can be used:
  (Sh)  = 0.6 sqrt((Re) ) sqrt((Sc) )
      c                c           c

  b = d       / 1.242

  w = (8 sigma b / dp) (n (n+1) (n-1) (n+2) / (n+1) rho  + n rho )
                                                       d        c

  k  = 0.45 sqrt(D     \omega)
   d              d,ij

where n=2, and d_p is in cm for the calculation of b and w.

6.5 Spray columns

This design method is adapted from Jordan (1968) and Lo et al. (1983).

6.5.1 Design

The height of a stage in a spray column is set to the default value of 0.4 m and the hole diameter in the distributor to 0.005 m. The hole velocity (U_o) in the distributor is set to 0.1 m/s from which the total hole area is then:
  A  = Q  / U 
   o    d    o

The droplet diameter can be calculated from (Vedaiyan et al, 1972):

                  2            -0.067
  d  = 1.592 (( U   / 2 g d  ))       sqrt( \sigma \over g \Delta \rho)
   p             o         o

The flood velocity of the continuous phase is (Treybal, 1963):

                         0.28                0.075                        0.056                     2
  U   = (0.3894 Delta rho     / [ 0.2165 eta       sqrt(\rho ) + 0.2670 d       sqrt(\rho  \alpha) ] )
   cf                                       c               c            p               d

where alpha=Q_d/Q_c. The disperse holdup at flood is

  phi   = (sqrt(\alpha  + 8 \alpha) - 3 alpha / 4 ( 1 - alpha ))

The velocity of the continuous phase is then
  U  = FF SF U  
   c          cf

and the column area
  A  = Q  / U 
   c    c    c

from which the column diameter can be calculated (The column area must also be larger then the total hole area, if not, the column area is set to four times the hole area).

6.5.2 Report

The fraction of flooding reported is calculated as
  FF = (U  / SF U  )
         c       cf

where U_cf is computed as in the spray column design and U_c=A_c/Q_c. No lower operating limit is calculated. The dispersion coefficient for the continuous phase is (Vermeulen et al., 1966):
  (E  / V  H) = 7.2 sqrt(U  D )
    c    c                d  c

Since the dispersion coefficient for the disperse phase is unknown it is set equal to that for the continuous phase.

6.5.3 Mass Transfer Coefficients

The transition drop size below which droplets become stagnent is calculated from

           2      4         4
  P = (rho   sigma  / g eta   Delta rho)
          c                c

  d    = 7.25 sqrt( \sigma \over g \Delta \rho P     )

The drop terminal velocity is (Satish et al., 1974):

                  2            -0.082                            2   1/4
  V  = 1.088 (( U   / 2 g d  ))       (( sigma g Delta rho / rho   ))   
   t             o         o                                    c

With the continous operating and flood velocities the fraction of flooding is calculated and then the disperse phase holdup is
  phi  = FF phi  
     d         df

and the slip velocity
  V  = (1 - phi ) V 
   s           d   t

If the drops are stagnent (d_p < d_p,t) the disperse MTC is computed from
  k     = 18.9 D     / d 
   d,ij         d,ij    p

else the Handlos-Baron correlation (1957) is used:
  k  = (0.00375 V  / (1 + eta  / eta ))
   d             s           d      c

For the continuous phase MTC we use (Ruby and Elgin, 1955)

                -0.43    -0.58
  k  = 0.725 Re       Sc       (1 - phi ) V 
   c           c        c              d   s

  Re  = (d  V  rho  / eta )
    c     p  s    c      c

  Sc  = (eta  / rho  D    )
    c       c      c  c,ij

The interfacial area for mass transfer per unit of volume is
  A  = (6 phi  / d )
   i         d    p

6.6 Modeling Backflow

The backflows in the column are computed from the dispersion coefficients with:
  alpha  = (E  / V  H) - 0.5
       d     d    d

  alpha  = (E  / V  H) - 0.5
       c     c    c

where alpha is the fractional backflow ("entrainment") in the stage, and H is the stage height.

For spray columns (Perry, 198x):

  E  = 7.2 sqrt( V  D  )
   c              d  c

For packed columns:
  log (E  / V  d ) = 0.046 (V  / V ) + 0.301
        d    d  F            c    d

  log (E  / V  d ) = 0.161 (V  / V ) + 0.347
        c    c  F            c    d

For a RDC:
  E  = 0.5 H V  + 0.012 R N H ( ( S / D  ) )
   c          c                        c

  E  = F E 
   d      c

where F is calculated by

             5     2                3.3
  F = (4.2 10  / D  ) ( ( V  / h ) )   
                  c        d

and must be larger or equal than one. Krishna uses:
  E  = (0.5 H U  / (1-Phi )) + 0.012 R N H ( ( S / D  ) )
   c           c         d                          c

  E  = (0.5 H U  / Phi ) + 0.024 R N H ( ( S / D  ) )
   d           d      d                         c

Symbol List

a_p      Packing area per unit volume (m^2/m^3)
a_s      Static holdup area per unit volume (m^2/m^3)
A_a      Total tray active area (m^2)
A_d      Downcomer area (m^2)
A_i      Interfacial area per unit volume (m^2/m^3)
A_h      Total tray hole area (m^2)
A_n      Netto tray area (m^2), A_n=A_a+A_d
A_t      Total tray area (m^2)
C_d      Drag coefficient
D        Binary diffusion coefficient (m^2/s)
D_c      Column diameter (m)
d_e      Effective drop diameter (m) ?
d_h      Hole diameter (m)
d_min  Minimum droplet diameter (m)
d_p      Sauter mean drop diameter (m)
Eo       Eotvos number ( Delta rho g d_h / sigma)
f        Free area ratio (A_h/A_a)
F        Molar flow (kmol/s)
Fr       Froude number (U^2_h / g d_h)
FF       Fraction of flooding
g, g_c  Gravitational constant, 9.81 (m/s^2)
H        RDC compartment height (m)
h_c      Height of coalesced layer (m)
h, h_drop  Height of drop rising zone (m)
h_stage  Stage height for packed column (m)
k        Binary mass transfer coefficient (m/s)
M_w      Molecular weight (kg/kmol)
N        Rotation speed (rad/s)
Nu       Nusselt number
Pe       Peclet number
P_i      Power input (?)
Q        Volumetric flow (m^3/s)
R        Rotor diameter (m)
Re       Reynolds number
S        Inner stator diameter (m)
Sc       Schmidt number
Sh       Sherwood number
SF       System derating factor
t        Contact time (s)
t_s      Tray spacing (m)
U_c,U_d  Continuous, disperse velocity (m/s)
U_cf     Continuous phase superficial velocity at flood (m/s)
U_h      Hole diameter (m/s)
V_i      Tray volume for interfacial mass transport (m^3)
V_s      Slip velocity (m/s)
V^o_s    Slip velocity at zero disperse phase holdup (m/s)
We       Weber number ( rho_d U^2_h d_p / sigma)
W_l      Weir length (m)
alpha   Phase ratio (Q_d / Q_c)
rho     Mass density (kg/m^3)
phi_d   Disperse phase holdup fraction
phi_ds  Static disperse phase holdup fraction
sigma   Interfacial tension (N/m)
eta     Liquid viscosity (Pa.s)
mu      Kinematic viscosity (eta / rho)
zeta    Tortuosity
c     Continuous phase
d     Disperse phase,
i     Interface,
        Component i
j     Component j


F.H. Garner, M. Tayeban, Anal. Real Soc. Espan. Fis. Quim. (Madrid), Vol. B56 (1960) pp. 479.

R.M. Griffith, Chem. Eng. Sci., 12, 198 (1960).

A.E. Handlos, T. Baron, ``Mass and Heat Transfer from Drops in Liquid-Liquid Extraction'', AIChE J., 3 (1957) pp. 127--136.

A.E. Handlos, T. Baron, AIChE J., 6, 145 (1957).

Hughmark, Ind. eng. Chem. Fundam., 6, 408 (1967).

D.G. Jordan, Chemical Process Development, Part 2, John Wiley, New York (1968).

W.J. Korchinsky, "Liquid-Liquid Extraction Column Modelling: Is the Forward Mixing Influence Necessary?", Trans. I. Chem. E., Vol. 70, Part A, 333--345 (1992).

R. Krishna, S.M. Nanoti, A.N. Goswami, "Mass-Transfer Efficiency of Sieve Tray Extraction Columns", Ind. Eng. Chem. Res., Vol. 28 (1989) 642-644.

R. Krishna, Design of Liquid-Liquid Extraction Columns, University of Amsterdam (NL), (1993).

R. Kronig, J.C. Brink, Appl. Sci Res., A2, 142 (1950). A. Kumar, S. Hartland, "Prediction of Axial Mixing Coefficients in Rotating Disc and Asymmetric Rotating Disc Extraction Columns", Can. J. Chem. Eng, Vol. 70, 77--87 (1992).

A. Kumar, S. Hartland, "Prediction of drop size, dispersed-phase holdup, slip velocity, and limiting throughputs in packed extraction columns", Trans. IChemE., 72, Part A, 89--104 (1994).

M. Lao et al., "A Nonequilibrium Stage Model of Multicomponent Separation Processes VI: Simulation of Liquid-Liquid Extraction", Chem. Eng. Comm., 86, p73--89 (1989).

G.S. Laddha, T.E. Degaleesan, Transport Phenomena in Liquid Extraction, McGraw-Hill (1978).

T.C. Lo, M.H.I. Baird, C. Hanson, Handbook of Solvent Extraction, John Wiley, NY (1983).

J.A. Rocha, J.L. Humphrey, J.R. Fair, "Mass transfer Efficiency of Sieve Tray Extractors", Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4 (1986) pp. 862--871.

P.N. Rowe, K.T. Claxton, J.B. Lewis, Transact. Inst. Chem. Eng., 43, T14 (1965).

P.M. Rose, R.C. Kinter, ``Mass Transfer from Large Oscillating Drops'', AIChE J., Vol. 12, No. 3 (1966) pp. 530.

E.Y. Kung, R.B. Beckman, ``Dispersed-Phase Holdup in a Rotating Disk Extraction Column'', AIChE J., Vol. 7, No. 2 (1961) pp. 319-324.

A.M. Rozen, A.I. Bezzubova, Theor. found. Chem. Eng., 2, 715 (1968); translated from Teor. Osnovy Khim. Tekh, 2, 850 (1968).

Ruby, Elgin, Chem. Eng. Prog., 51, Sump. Ser. 16, 17 (1955).

L. Satish, T.E. Dagaleesan, G.S. Laddha, Indian Chem. Eng., Vol. 16 (1974) pp. 36.

A.F. Seibert, J.R. Fair, "Hydrodynamics and Mass Transfer in Spray and Packed Liquid-Liquid Extraction Columns", Ind. Eng. Chem. Res., 27, No.3, 470--481 (1988).

A.H.P. Skelland, R.M. Wellek, AIChE J., 10, 491 (1964).

A.H.P. Skelland, W.L. Conger, Ind. Eng. Chem. Process De. Devel., 12, 448 (1973).

A.H.P. Skelland, D.W. Tedder, Handbook of Separation Processes, Ed. R.W.Rousseau, Wiley (1987).

S. Stemerding, E.C. Lumb, J. Lips, ``Axiale Vermischung in einer Drehscheiben-Extraktions Kolonne'', Chem. ing. Tech, Vol. 35 (1963) pp. 844--850.

G. Thorsen, S.G. Terjesen, Chem .Eng. Sci., 17, 137 (1962).

R.E. Treybal, Mass Transfer Operations, 3rd ed., McGraw-Hill, New York (1980)

R.E. Treybal, Liquid Extraction, 2nd ed., McGraw-Hill, New York (1963).

S. Vedaiyan, T.E. Degaleesan, G.S. Laddha, HE. Hoelscher, AIChE J., Vol. 18 (1972) pp. 161.

Vermeulen et al., Chem. Eng. Prog., Vol. 62, No. 9 (1966) pp. 95.

M.E. Weber, Ind. Eng. Chem. Fund., 14, 165 (1975).

Mass Transfer Coefficient correlations

Mass Transfer Coefficients correlations for the continuous phase (chapter 3.4, Handbook of Solvent Extraction): Mass Transfer Coefficients correlations for the disperse phase (chapter 3.4, Handbook of Solvent Extraction): where Reynolds, Schmidt, and Sherwood numbers are defined as:
  (Re) = ( d  rho V  / eta )
            p      s

  (Sc) = ( eta / rho D   )

  (Sh)   = ( k   d  / D   )
      ij      ij  p    ij

  mu = ( eta / rho )
Table 1 of chapter 10 in the Handbook of Solvent Extraction supplies us with three more models for the drop rise zone. One for stagnant drops (Skelland and Conger, 1973):

                                                         0.5    0.5
  k  = - (( d  / 6 t )) (( rho  / M  ))   ln ( 1 - ( pi D      t    / 0.5 d  ) )
   d         e                d    d   av                   vd             e

                                                                   0.5                     0.333
  k  = 0.74 (( D   / d  )) (( rho  / M  ))   (( d  V  rho  / mu  ))    (( mu  / rho  D   ))     
   c            vc    e          c    c   av     e  s    c     c            c      c  vc

for circulating drops (Treybal, 1963):

                                                           2    -0.34                     -0.125        2                 0.37
  k  = 31.4 (( D   / d  )) (( rho  / M  ))   (( 4 D   t / d   ))      (( mu  / rho  D   ))       (( d  V   rho  / sigma ))    
   d            vd    e          d    d   av       vd       e              d      d  vd              e   s    c

                                                     -0.34                     -0.58
  k  = 0.725 (( rho  / M  ))   (( d  V  rho  / mu  ))      (( mu  / rho  D   ))      V  ( 1 - phi  )
   c               c    c   av     e  s    c     c              c      c  vc          s          d

and for oscillating drops (Skelland and Conger, 1973):

                                                           2    -0.14                       0.68         3    2      4               0.10
  k  = 0.32 (( D   / d  )) (( rho  / M  ))   (( 4 D   t / d   ))      (( d  V  rho  / mu  ))     (( sigma  rho   / mu  g Delta rho ))    
   d            vd    e          d    d   av       vd       e             e  s    c     c                      c

                                                              -0.34                                      1.0                     0.7
  k  = (( D   / d  )) (( rho  / M  ))   (( d  V  rho  / mu  ))      [( 50 + 0.0085 (( d  V  rho  / mu  ))    (( mu  / rho  D   ))    )]
   c       vc    e          c    c   av     e  s    c     c                            e  s    c     c            c      c  vc

where t=h/V_s (with h as the height of the drop rise zone) and V_s=V_t(1-phi_d). Perry's also supplies us with some more correlations. There we find that () is from Ruby and Elgin (1955) and is to be applied for circulating drops. Another correlation for the continuous mass transfer coefficients for circulating drops is by Hughmark (1967):

                                 0       0.339         1/3    1/3    0.072
  (k  d  / D ) = [ 2 + 0.463 (Re) .484 Sc       (( d  g    / D     ))      ] F
    c  p    c                                 c     p            c

                                      2             3
  F = 0.281 + 1.615 kappa + 3.73 kappa  -1.874 kappa 

            1/8                1/4                        1/6
  kappa = Re    (( mu  / mu  ))    (( mu  V  / sigma g  ))   
                     c     d            c  s          c

where Re is the droplet Reynolds number. A correlation for the disperse mass transfer coefficient for oscillating droplets by Rose and Kinter (1966) is:

  k  = sqrt( (4 D  \omega \over \pi) (1 + \delta + (3 \over 8) \delta ) )
   d             d

  w = (1 / 2 pi) sqrt( 192 \sigma g  b \over d   (3\rho +2\rho ) )
                                   c           p       d      c

  b = 1.052 d      

where delta can be taken as 0.2 if unknown.

RDC's: Korchinsky

Korchinsky (1992) summarizes correlations for RDC's from literature and adivises on to use the Kumar and Hartland correlations (1986). They use the following dimensionless groups:

  N  = (( N D   rho  / mu  ))
   1          r    c     c

  N  = (( N  D  / g ))
   2          r

  N  = (( mu  / sqrt( \sigma \rho  D ) ))
   3        c                    c  r

  N  = (( rho  / rho  ))
   4         d      c

  N  = (( D   rho  g / sigma ))
   5        r    c

  N  = (( H / D  ))
   6           r

           2   2           2
  N  = (( D   H  rho  g / D   sigma ))
   7        s       c       c

  N  = (( Delta rho / rho  ))
   8                     c

               0.25      0.25       0.75
  N  = (( mu  g     / rho      sigma     ))
   9        c                c

  N   = (( V   rho  / g sigma ))
   10        d    c

  N   = (( g sigma / rho  ))
   11                   c

  N   = (( D  / D  ))
   12       r    c

            4      0.25    0.25      0.25
  N   = (( V   rho      / g     sigma     ))
   13        d    c

  N   = (( N D  / V  ))
   14         r    c

  N   = (( N D  / V  ))
   15         r    d

  N   = (( D   g Delta rho / sigma ))
   16        r

  N   = (( V  D  rho  / mu  ))
   17       c  r    c     c

  N   = (( D  / H ))
   18       c

  N   = (( D  / D  ))
   19       s    c

The Sauter droplet size is computed by the high Reynolds formula from Kumar and Hartland (1986):

                  0.55                    -1.3   0.75   -0.3   0.28
  (d   / D ) = k N      exp ((-0.23 N )) N      N      N      N     
    32    r           1              2        3      4      5      6

with 10^3 k=7.01 for no mass transfer. The Kumar and Hartland disperse holdup:

                    n_1     n_2   n_3   n_4   0.22                      0.35
  phi  = [ k  + k  N     ] N     N     N     N       (( 1 + (V  / V ) ))    
     d      1    2     2       7     8     9      10          c    d

where (using all data) k_1=65.73, k_2=74.20, n_1=1.24, n_2=-0.34, n_3=-.049, and n_4=0.53. The slip velocity is computed by

                                     0.52   0.25    -0.45   0.08   1.03   0.51    0.28
  V  = [ k  + k  exp ((-1.28 N )) ] N      N       N       N      N      N       N      
   s      6    7              2          8      11       9      5      6      12      13

with (for all data) 10^2 k_6=-5.11 and k_7=0.20. the continuous phase dispersion coefficient is given by

                                                                                 -0.08    -0.16    0.1    2
  (E  / V  H) = 0.42 + 0.29 (V  / V ) + ( 0.0126 N   + ( 13.38 / 3.18 + N   ) ) N      7 N        N      N   
    c    c                    d    c              14                     14           1        12     18   19

and the disperse phase coefficients

                                                   -0.64    -0.7      -0.9
  (E  / V  H) = 0.3 (( V  + V  / V  )) + 9.37 N   N        N       phi     
    d    d              c    d    d            15       16      12        d

Packed columns: Kumar and Hartland

Kumar and Hartland (1994) developed new correlations for the drop diameter, dispersed phase holdup, slip velocity, and flooding velocities for packed extraction columns using a large database. The Sauter mean droplet size is

                  1/4                 1/4      3/4        0.19
  d  = C  [( mu  g    rho  / Delta rho    sigma    rho  )]     sqrt( \sigma \over g \Delta \rho )
   p    1      w         w                            d

where C_1 is 2.54, 2.24, or 3.13 for no mass transfer, transfer from c to d, and from d to c. The dispersed phase holdup is

             -1.11                       -0.50                   2        2    1/3   -0.72                0.10                        1/3  1.03                               1/3
  phi  = C  e      (( Delta rho / rho  ))      [( (1 / a ) (( rho   g / mu   ))    )]      (( mu  / mu  ))     [ Vd (( rho  / g mu  ))    ]     exp [ 0.95 V  ( rho  / g mu  )    ]
     d    2                          c                  p         c        c                    d     c                   c       c                         c      c       c

where C_2 is 5.34, 6.16, or 3.76. The slip velocity is

              -0.11                       0.40                   2        2    1/3   0.61                -0.10                   -1/3
  V     = C  e      (( Delta rho / rho  ))     [( (1 / a ) (( rho   g / mu   ))    )]     (( mu  / mu  ))      (( rho  / g mu  ))    
   slip    3                          c                 p         c        c                   d     c               c       c

where C_3 is 0.24, 0.21, or 0.31. Or, as function of the dispersed phase holdup:

              -0.17                       0.41                   2        2    1/3   0.59                -0.10                   -1/3
  V     = C  e      (( Delta rho / rho  ))     [( (1 / a ) (( rho   g / mu   ))    )]     (( mu  / mu  ))      (( rho  / g mu  ))     (1-phi )
   slip    4                          c                 p         c        c                   d     c               c       c              d

where C_4 is 0.30, 0.27, or 0.38. The flooding velocity is:

                  2                          1.54                       0.41                   2        2    1/3   0.30                                          0.15
  V    (1+sqrt(R))  sqrt(a  / g) = alpha C  e     (( Delta rho / rho  ))     [( (1 / a ) (( rho   g / mu   ))    )]     (( mu  / sqrt(\Delta \rho \sigma / a ) ))    
   c,f                    p               1                         d                 p         c        c                   c                              p


                                1.54                       0.41                   2        2    1/3   0.30                                          0.15
  V    sqrt(a  / g) = alpha C  e     (( Delta rho / rho  ))     [( (1 / a ) (( rho   g / mu   ))    )]     (( mu  / sqrt(\Delta \rho \sigma / a ) ))    
   d,f       p               1                         d                 p         c        c                   c                              p

where alpha is 1.0 for continuous phase packing wetting and 1.29 for dispersed phase packing wetting.

Chapter 7. Interface and Technical Issues

In this chapter we will briefly discuss ChemSep's internals. We discuss which programs make up ChemSep and explain how they cooporate. All the supporting libraries and files are identified and explained. More information on printing from within ChemSep and other technical issues can be found in this chapter.

7.1 ChemSep Commandline Parameters

The following commandline parameters are optional when you start ChemSep by typing cs at the commandline: Any other parameter will be handled as a Sep-file, which ChemSep will try to load after the introduction screen (the default "sep" extension does not need to be added). The interface of ChemSep is too large to fit in conventional memory. Therefore, it uses an "overlay" technique to switch parts in and out of conventional memory from XMS, EMS, or the hard disk (in this order). Use the /v options to manipulate what type of memory is used for the overlays. You can not disable the hard disk, as there would be possibly no place for the overlay manager to store unused parts. Normally you should not need to use these options. To see which type of overlay is used type Ctrl-Y in the interface (which shows also the DOS version, coprocessor type, and memory status).

7.2 ChemSep Environment Variables

In order to solve column problems ChemSep requires more than the standard DOS memory of 640 kilobytes. Either EMS ( expanded) or XMS ( extended) memory can be used.

7.2.1 CauseWay DOS extender

CauseWay (Devore software) is currently our default DOS extender. Executables linked with this extender start with the "CW" characters. This DOS extender supports memory up to 4 GB (although you won't need that much to run ChemSep!). If physical memory is limited it will use the disk as virtual memory. All options can be set with one environment variable. The format is
   SET CAUSEWAY=[setting_1;] [setting_2;] [setting_n;]
Seven options are available:

7.2.2 Rational

The DOS Rational extender was our previous choice of DOS extender. Executables linked with this extender start with the "4G" characters. The DOS extender selects EMS and then XMS. You can force the DOS extender to use XMS above EMS by specifying it to use a block of memory larger than the available EMS. For example, the environment variable dos16m=:4m requests a block of memory of 4 MB (use the DOS set command to specify environment variables). If there is more than 4 MB of EMS available the DOS extender will use EMS, otherwise it will try to use 4 MB of XMS. If the physical memory in your machine does not allow you to run a large problem you can use virtual memory which is swapped to your hard disk by using the following environment variable setting: dos4gvm=deleteswap which allows up to 16 MB of virtual memory and deletes the swap file after the run is completed. The virtual memory is swapped to the DOS4GVM.SWP file which is placed in the root directory of the current drive.

7.2.3 SVGA drivers

In case the automatic selection of the Super VGA drivers XVGA16 and XVGA256 select a wrong video chipset, you can set the CHIPSET environment variable to change the detection test. For example SET CHIPSET=VESA,CRRS will set the CHIPSET environment variable so that the Super VGA driver will first test for VESA, if this fails for Cirrus, and if this also fails it will use a generic VGA mode. The chipset codes are (in the standard detection order): EVRX (Everex), CMPQ (Compaq), V7 (Video 7), C&T (Chips & Tech), CRRS (Cirrus), ATI, TSNG (Tseng), OAK, (Oak Technologies), GNOA (Genoa), TRID (Trident), PRDS (Paradise), NCR, AHED (Ahead Systems), S2, VESA.

7.2.4 Printer drivers

The printer drivers also can make use of extended/expanded memory for generating temporary raster files for the printouts. At the moment, no options can be set through environment variables for the memory use of these drivers.

7.3 ChemSep's Programs

ChemSep is not one program but is split into several executables and associated data files. Before we can explain in detail how ChemSep uses the programs lets first see what program names are present in the executables directory and what they do:

The drivers acts as the glue between all the interfaces and simulators. It also gives us flexibility in the way we run our programs and hides details from the average user. ChemSep and ChemProp are normally started by running their drivers: CS.EXE and CP.EXE, respectively. However, both ChemProp and ChemLib can be run from within the ChemSep interface through menu options under the input menu.

The "CW-" in front of the calculation programs denotes the DOS extender type that is linked to the executables. CauseWay executables include the DOS extender, the Rational DOS-Extender uses a separate file (DOS4GW.EXE). The executables often need at least two MegaByte of Extended memory in order to run. They also require a 386-based system (minimum) to run. The different DOS extender executables can be selected under Options/DOS-Extender. Check this setting when you get error messages stating that executables with the "CW-" or "4G-" are not found. CauseWay is currently our default.

As we support the simulators on other platforms as well, they must also run without the ChemSep interface. That is why we store all information (in- and output) in the problem files, with the .SEP extension. As the simulator must know what problem file to run it checks for the existance of the CHEMSEP.FIL file. If it exists, it will read this name from this file. If the file doesn't exist or no CHEMSEP.FIL file was found in the current directory, the simulator will prompt you for a SEP file. DOS and Windows executables of the simulators also allow specification of the SEP file on the commandline, making the CHEMSEP.FIL unnecessary.

To build in maximum flexibility we use the CS and CSW drivers. Both these drivers can run the simulators upon cooperation with the CS2.EXE interface. The default configuration makes the interface run the CSW driver which in turn runs the necessary simulator. The CSW driver insures that all the output of the simulator gets written into a window, giving the illusion that the simulator is an integral part of the interface. While solving the problem the interface is swapped out of DOS memory, and swapped back in upon termination of the simulator. In this case the CS driver is not doing anything and the interface can be invoked directly as well.

Alternatively, the interface can write a CHEMSEP.FIL file, exit, and the CS driver reads it. It determines which simulator is requested to run and start it. The simulator will read the CHEMSEP.FIL again and read the name of the SEP file to solve. Upon termination the CS driver restarts the ChemSep interface, telling it which SEP file it was solving. The interface reads the SEP file and deletes the CHEMSEP.FIL file. Again, the CS driver also ensures the output to be written into a window on the screen. If the interface is loaded directly, without the driver, this is not possible. However, the interface will still run the simulator by swapping itself out of memory, and back into memory after the simulator returns control. It will call the simulator directly without writing a CHEMSEP.FIL file. Without any driver the interface can't redirect the ouput to a window and thus the screen is cleared before and restored after the simulator runs.

Switching between these modes of solving SEP files is done by specifying the user program under the solve options. The default setting is to call the CSW driver (the interface will locate it, don't specify its path!). If the user program is left empty, the interface will write a CHEMSEP.FIL and quit if the CS driver was loaded, else it will call the simulator directly. Only in this case the screen will be cleared (it also requires the least amount of DOS memory). If the interface is swapped out of memory, it is written to extended memory or disk. If, for whatever reason, the swapfile on disk is deleted the interface will not be able to recover and abort to DOS.

If you are running the simulator programs by your self (that is not abnormal) you might like the small utility MAKEFIL to generate your CHEMSEP.FIL file. Since the simulator programs do not delete the CHEMSEP.FIL file you have only to create it once for each new problem. The MAKEFIL program takes as first commandline argument the SEP-file and as (optional) second argument the scrap file name:

    MAKEFIL  [Temporary file]
The CHEMSEP.FIL file will be written in the current directory. If you run under DOS or Windows, you can also specify the SEP-file pathname as commandline parameter to the simulator, and avoid the use of the CHEMSEP.FIL file all together.

7.4 Running ChemSep - Advanced Use

As noted in the previous section, ChemSep has separate programs for the calculations and interfacing with the user. In order for the Driver to know what is going on, the interface saves information in the CHEMSEP.FIL such as the name of the SEP-file, the program to run, the temporary scrap file name, the user program, the run window coordinates and color as well as the starttime.

In order to make ChemSep as versatile as possible, we implemented the User Program entry under the "User program" option in the "Solve Options" menu. If you want to run your own program you can enter its full path and name (with the extension!) there and that program will be run no matter what kind of operation is selected. In order for a user program to have access to the SEP-file, the interface provides the user program with its name as the second commandline parameter. In case it also wants to use the temporary file, that is supplied as the third parameter. Running a user program will clear the screen before it starts executing the user program.

It is logical to suppose that the user just might want to run several programs or his own program(s) before/after calling ChemSep's original calulation programs. To allow this you can use a batch file as the user program (use complete path, name and extension !). Look up in your DOS manual how to make batch files. In order to run the original calculations program we provide it to the batch file as the first commandline parameter and the SEP-file as the second. Even running the user program, the Interface generates a CHEMSEP.FIL file to be read by the calculation program. Here's an example of such a batch file that will type the problem SEP-file first before running (note how the parameters are accessed with %1 and %2):

  @Echo off
  Rem -----------------------------------------------------------
  Rem Echo is set off to avoid to show this batch file is running
  Rem use the "@" to suppress echoing of commands to the screen
  Rem -----------------------------------------------------------
  Rem Type the SEP-file:
  Rem ------------------
  Type %2
  Rem Pause for the user to strike a key:
  Rem -----------------------------------
  Rem Run the appropriate ChemSep calculation:
  Rem ---------------------------------------
  Rem Done!
Of course you can leave out all the Rem(arks) if you type this batch file but it is a good habit to document your code ! You might create a complete set of batch files for different problems. Programs that might be run after the calculation might be cost estimation or design (in case of equilibrium simulation) programs. Using command line parameters (or the CHEMSEP.FIL file) your program might append its results to the SEP-file, so all results will be collected in there. Here is an example of a batch file that will automatically run ChemProp to generate physical property information in the sep-file:
   @Echo off
   echo Running ChemSep with physical property information generation
   rem Run executable
   rem Run ChemProp to generate physical property information
   c:\chemsep\bin\cp /c %2
   rem Done!
The ultimate freedom is allowed by typing "DODOS" as User Program in the Interface. The driver will automatically locate DOS and run it. All Dos commands will be available to you. The only way to access the SEP-file and other information is to read the CHEMSEP.FIL file. When you are ready to go back to the Interface you type "EXIT" and press Enter.

By allowing you to shell to DOS or run batch files from within ChemSep we have created the maximum flexibility. Although ChemSep takes a lot of effort to prevent your system or ChemSep from crashing, it is possible to do so, using batch files or the DOS shell. Avoid deleting crucial files (such as executable files), or changing system parameters while the ChemSep-Driver is loaded. Do not load any TSR (Terminate and Stay Resident) programs or device drivers, since these programs will be removed from memory after the Driver takes over again. However, the interrupt to trigger these programs usually remains active. Changing directories is allowed, but remember the current CHEMSEP.FIL is written in the current directory ! The driver will change the directory back to what it started running from, when returning to the Interface. If you want to change the current directory use the Directory option in the File menu ! Do not load ChemSep again by invoking the driver, as multiple copies will be loaded into memory.

7.5 ChemSep Libraries and Other Files

A number of libraries has been added to the ChemSep package. We can divide the data libraries into the following groups:

The Pure Component Data (PCD) files contain the information of pure components like molecular weight, critical temperature, acentric factor, vapour pressure correlation constants, UNIFAC group ID's and number of groups etc. To access this data we have developed ChemLib which is a completely menu driven data manager - very similar to ChemSep, in fact - that will allow you to search PCD-file(s) and select component data records to edit. It can also move components from one PCD-file to another, or to text files. ChemSep v3.5 and higher can also use components from text files instead of PCD files. We prefer the faster (to search) binary PCD format for distribution, however, component data information in text format can have additional information as long as you append this information after the regular component's data items in the text file ( ChemSep and ChemLib will stop reading after the fixed set of items and look for the start of the next component). The components ID numbers (Library Index) are based on the system developed at Penn State University and adopted by DIPPR.

The Interaction Parameter Data (IPD) files are ASCII text files with the interaction parameters for activity coefficient models and equations of state. Currently we have the following IPD files:

All these files are in plain text format file so the user can extend the data files (it is probably better to backup the original files or to use a new name for the extended files). Edit these files with any ASCII editor, or the built-in editor in ChemSep. The data in these files comes after the line with the [IPD] keyword. The next line contains the file-comment after the "=", which will be used as a header in the interface. Next comes a line (only for activity coefficient models) with the interaction parameters units. Next comes the data, one line per binary pair. The first two numbers are the component library indices, then the interaction parameters. Appended text is optional but will also be displayed. As an example we include the first relevant lines of the NRTL.IPD file:

Comment=DECHEMA NRTL data @ 1atm.
  1101   1921  -189.0469   792.8020 0.2999 Methanol/Water p61 1/1a
Lines starting with a "#" are comment lines which may appear anywhere. Since the interface will only start reading the file from the [IPD] keyword on, you can start the file with some text describing where the data was obtained and remarks on who/when/how changed the file. Most of the IPD files contain information from the DECHEMA series, a very extensive collection of interaction parameters. The Hayden O'Connell virial parameters are from Prausnitz et al. (1980).

Polynomial K-value and enthalpy correlation coefficients as well as extended Antoine coefficients are stored as component LIBraries (LIB files), which are ASCII files as well. The default LIB files are

Here the interface starts reading the file after the [LIB] keyword. Again a comment is read from the next line (after the "=") and then the data starts with a line for each component (first the library index followed by the coefficients). For example the extended Antoine file (with data from Prausnitz et al., 1980) starts like:

Comment=Extended Antoine Prausnitz et al.
#  ID  A            B          C   D   E          F   G
  902  3.15799e+01 -3.2848e+2  0.  0. -2.5980e+0  0.  2.0  Hydrogen
Group Component Data (GCD) files are binary files containing data for group contribution method such as UNIFAC or ASOG. Currently only UNIFAC files are present (UNIFACRQ.GCD, UNIFACVL.GCD and UNIFACLL.GCD) for Vapour-Liquid and Liquid-Liquid systems. These files need not be changed, unless for example the UNIFAC group tables have changed. Several GCD-files are used for the estimation of pure component data in ChemLib and UserPcd.

Internal layout data (ILD) files are text files which store tray or packing layouts for use by the nonequilibrium model. This way a specific design can be saved and reloaded upon demand. Two formats are used, which resemble eachother very much. The difference lies in the use of an equal sign in the internals type. Without the equal sign we are dealing with one of the fixed internal layouts. If an equal sign is present then the two following numbers define the type of operation (VL=1 or LL=2) and the number of parameters that describe the internalstype. For example the variable internals type for a baffle tray column section is:

10 Discrete Baffle = 1 3
1.2 Column diameter (m)
0.5 Baffle area (%)
0.7 Baffle spacing (m)
Note that the parameters are -always- written in base units, not in the units that are used for display/input. See also the discussion of the model and internals definitions in ChemSep.

Packing Data (CHEMSEP.*PD) files contain many physical and model parameters for various random and structured packings. They are text files that might be modified by the user with an ASCII editor, though there is one restriction, namely that the first line should not be changed! A shortened version of the structured packing data file CHEMSEP.SPD is shown below.

# CHEMSEP SPD Structured Packing Data
#   Type (Name):        Specific Equiv. Channel Packing Void    ...
#                       packing  diam.  flow    factor: fract:  ...
#                       surface:        angle:
Koch Flexipac 1       M   558/m  0.00897m   45    98/m   0.91  ...
Koch Flexipac 2      SS   223/m  0.01796m   45    43/m   0.95  ...
Koch Flexipac 3       M   135/m  0.03592m   45    26/m   0.96  ...
Koch Flexipac 4       M    69/m  0.07183m   45    20/m   0.98  ...

Glitsch Gempak 1A     M   131/m  0.03592m   45    30/m    *    ...
Glitsch Gempak 2A    SS   223/m  0.01796m   45    52/m   0.95  ...
Glitsch Gempak 2AT   SS   223/m  0.01796m   45     *     0.96  ...
Glitsch Gempak 3A     M   394/m  0.01346m   45    69/m    *    ...
Glitsch Gempak 4A     M   525/m  0.00897m   45   105/m    *    ...
Lines starting with "#" are comment lines and are ignored (except for the first line). A line starting with "!" is used to set the length of the packing type identifiers which is set equal to the length of that line. Lines starting with an "@" will insert a separator in the list with packings and blank lines will be ignored. As you can see units may be added as long as there is no space between the number and the unitstring (otherwise errors will occur in reading this file!).

Miscellaneous ChemSep files include:

The CHEMSEP.UDF file contains the definitions for the units and the unit conversions in ChemSep. The first line must have the number of following lines with on each a unit definition. Such a unit definition consists of 15 characters (from column 1) with the unit abreviation (take care, these are case sensitive !) followed by 15 characters (from column 16) with the full unit name (not case sensitive). Then, from column 31 the offset-factor (fo) and multiplication factor (fm) come and finally the reference unit. The conversion is done according the following formula:

  Number (in Reference units) = fm * ( Number (in Units) - fo )
For example 22 C = 1.0 ( 22 - -273.15) = 295.15 K. You can inspect ASCII text files with ChemSep's file viewer ( F7) or make simple modifications with edit-file. However, we strongly suggest you do not change the original data files that come with ChemSep. We carefully selected and typed the data into these files and other users might use your changed data and obtain erroneous results. We urge you first to copy the file to another name before you change or add anything in these files. Errors in the unit conversion can be very frustrating so it is good to check some results if you have changed or added a unit definition. To encourage you to do so we made CHEMSEP.UDF a read-only file !

The CHEMSEP.SYN is a file containing synonyms for over 1000 compounds. While searching for a special component name you can use synonyms if you have selected to do so in the options interface spreadsheet. You must select this file as your synonyms file. The synonym search does not work while typing in a searchlist for a synonyms name. You will have to issue a search under the synonyms name ! The synonyms file is an ASCII file you can modify to your needs.

The CHEMSEP.HLP file contains the information to provide you with help when you press ( F1) for help. It is a binary file that can not be changed. It is generated with the MAKEHELP utility from help source (HSR) files. This utility and the source files are not part of the ChemSep distribution. If you find errors or shortcomings in the help please notify the authors.

The CHEMSEP.TXT file contains the ASCII text of this book. Normally ChemSep will be configured with function key ( F9) assigned to automatically load this file. There are also other formats available on our ftp-site (see Author information). The ASCII version is somewhat limited in available characters, sub- or superscripts, and equations.

The CHEMSEP.CNF file is used to store the default configuration, such as the macro definitions, directory structure, selected video and printer devices, and solve options. CHEMSEP.CNF is the default configuration file which will be loaded upon startup of the interface. If none is found in the current directory, the file supplied in the original distribution is used. This allows one to have multiple CHEMSEP.CNF files in different directories, which automatically configure the interface to (a) specific problem(s).

Finally, the CHEMSEP.SCR file contains the additional introduction screen(s) that are shown on startup of the program. These are used to stipulate conditions of the use of the program, but could can be adapted to suit the users needs (for example if ChemSep is installed on a network, the operator can place important notes here). If the file contains more lines than can be shown on the screen, the user will have to press multiple times to go through the various screens sequentially.

7.6 The SEP-file format

The SEP-files are written in a format which is readible by a human as well as by the calculation program. However, there are some strict rules ! The SEP-file is constructed using delimiters in square brackets: []. The order of the delimeters is not directly of importance (although ordering delimiters enhances the speed of reading a SEP-file). The sections under the delimiters are ordered and have a fixed format. Usually they consist of lines with a selected value and a comment. ChemSep uses two major sections: the INPUT and RESULTS sections. The INPUT starts with the delimiter [ChemSep] and ends with [End of Input]. The Results starts with [Results] and ends with [End of Results]. Within these sections there are sub sections. Note that a "*" denotes that the value is not yet set by the user. In some cases the Interface might create cryptic "* *" lines where the first star denotes the (not yet known) value of a selector and the second the description of the (not yet) selected item.

Remember that normally the Interface will not read any comments you have added yourself to the SEP-file and thus, not save them again when saving from within the Interface. To overcome this problem we have added the [User-Data] and [End User-Data] delimiters. When loading a SEP-file it is the last section that is looked for. If found, it is read into a buffer and written back to the file if saved again. Note that the delimiters each have to be on a separate line. You can use this data block to save information about the problem or to store parameters for your own programs that process SEP-files. You can edit User Data within the interface of ChemSep under the Solve Options ( F6). The following keywords are used to switch several hidden features (with explanation in "()"):

Diffusivity model            (0=Maxwell-Stefan,  1=Effective)
Liquid MS-diffusivity model  (0=Kooijman-Taylor, 1=Wesselingh-Krishna)

[No user interaction]        (if present the user is not asked for
                              more iterations if maximum number has
                              been reached, but the program exits)

[Sensitivity]                                   (sensitivity factors)
Vapour/Light-liquid Mass Transfer Coefficient
Liquid/Heavy-liquid Mass Transfer Coefficient
Interfacial Area
For a column problem the INPUT subsections are:

Version number and SEP file name

Current directories

Current set of units

Number of components, for each component library offset (in the PCD file),
Index, Name, and PCD-library filename

Operation type and kind, condenser and reboiler types, number of stages,
feeds, sidestreams, and pumparounds.

Property selections:

  K model, Activity coefficient, Wilson model, UNIQUAC model, Equation of
  State, Cubic EOS, Virial EOS, Vapour pressure models

  Enthalpy model

  [Physical Properties]
  Physical Properties model selections

  [Property Data]
  Pure component data or interaction parameters if required


  Number and stage with duty, if specified.

  Number and for each section: section number, begin and end stages, model
  selections, and tray/packing layout data.

  default efficiency, number of exceptions and stage with value, if

  Type of pressure specification, condenser, top, bottom pressures, pressure

  Number and for each feed: feed state, stage, temperature, pressure, vapour
  fraction, number of componentflows, molar component flows.

  Number and for each sidestream: stage, phase, specification type and value.

  Specification type and value(s) if present

  Specification type and value(s) if present

[Solve options]
Initialization type, solving method, damping factor (if present), accuracy,
maximum number of iterations, and print options.

Temporary file, and user program

Here the user data text is written.
[End User-Data]

[End of Input]
For an equilibrium column the [Sections] subsection will be absent, for a nonequilibrium model the [Efficiencies] subsection is not needed. For a flash a different set of specifications is present consisting of the [Feeds] subsection and a [Flash] specifcation subsection where flash type and specifactions are made. The RESULTS section for a nonequilibrium problem looks like:

[Vapour phase compositions]
[Liquid phase compositions]
[Interface vapour mole fractions]
[Interface liquid mole fractions]
[Murphree efficiencies]
[Mass transfer rates]
[Condenser Heat Duty]
[Reboiler Heat Duty]
[Feed streams]
[Top product]
[Bottom product]
[Designed Sections]
[Operating Limits]
[End of Results]
For equilibrium problems the Interface mole fractions, Murphree efficiencies, mass transfer rates, designed sections, and operating limits will be missing from the above list.

7.7 Printing graphs in ChemSep

ChemSep supports dot matrix printers, laser printers, inktjet printers, or HP plotters to print/plot its graphs. Besides printing directly to a printer, ChemSep can also write the output to a file. ChemSep supports nine DeskTop Publishing (DTP) file formats as well. In order to print to a device or write to an output file ChemSep needs to know what type of device you have. Select device, mode, port (file) and work path (for temporary files written while generating output) in the printer setup under the graphs or the output setup in the options. Of course each printer has usually several different modes to print. By default none of the graphs are printed in color, unless a color device is selected. ChemSep lets you choose between three page formats:

However, these page formats vary slightly from printer to printer. Plotters usually plot only at FULL size. Besides these three page formats there are (usually) additional modes. Select "Other" and type the mode number you want (see table below). The lowest mode number is zero and will always work.

ChemSep supports the following printers and plotters and desktop publishing file formats:

Dot matrix printers: Max.Mode Desktop Publishing: Max.Mode
Epson 9-pin dot matrix 8 Zsoft PCX 1
Color Epson 9-pin dot matrix 5 Windows 3 BMP 1
Epson 24-pin dot matrix 8 Gem IMG 2
Color Epson 24-pin dot matrix 5 TIFF compressed 2
IBM Proprinter X24 8 TIFF uncompressed 2
IBM Quietwriter 8 ANSI CGM 1
Toshiba 24-pin dot matrix 2 AutoCad DXF 0
OkiData ML-92 dot matrix 2 Video Show 0
Word Perfect Graphic 1
Laser/Inktjet printers: Max.Mode Hewlett-Packard plotters: Max.Mode
LaserJet II 8 HP 7090 3
LaserJet III 8 HP 7470 1
DeskJet 8 HP 7475 7
Color DeskJet 8 HP 7550 7
PaintJet 14 HP 7585 9
Postscript 11 HP 7595 9

7.8 Model Definition and Selection

ChemSep reads a definitions file (CHEMSEP.DEF) at startup, where models for the mass transfer coefficients, pressure drop, flow models, entrainment, and holdup are defined. This alleviates us from adapting the ChemSep interface upon any addition or modification of a model. The defintion file also includes the definitions of any non-standard type of internal. In case no definitions file is found, the nonequilibrium part of ChemSep is disabled.

The definitions file must start with "[ChemSep Definitions]" followed by a line which defines the version (the current version is 3). Only after these two lines content may appear. Everything after a "mbox char35" is regarded as comment. Comments may appear on any content line, and may also be appended after any content. There are five different definitions types that are read from the definition file: [InternalType], [Operation], [ModelType], [Internal], and [Model]. Other entries are ignored. The five different types of definitions may be mixed throughout the definitions file. A definition may have no blank lines in it as entries are ended by a blank line!

The [InternalType] defines the nature of an internal, that is whether it is continuous or discrete. The [Operation] defines the nature of the operation, that is whether we are contacting a vapor and a liquid (VL) or two liquids (LL). The [ModelType] defines the type of models, that is whether the model describes mass transfer coefficients, pressure drops, etc. Each of these has the following fields: ID, Name, and Short, for example:

The [Internal] definitions link these types together: its "Type" must be one of the defined internal types, and its "Operation" one of the defined operation types. It also defines which "Models" need to be specified for the internal as well as which layout "Parameters" the internal has. These parameters and models are separated by semicolons (;) or comma's (,). The parameters might include a units string in parenthesis ("(",")") and a default value (with different units) after a colon (:). Note that multiple "Parameters" entries are allowed but that that its total lentgth is limited to 255 characters (these rules also apply for the "Models" and "Internals" entries).

A parameter entries may also contain a picklist in brackets ("[","]"), separated by vertical lines. Then instead of a numeric entry the interface will ask the user to select the parameter from a picklist of items. The first item in the list will result in the value 1, the second in 2, etc. A default item may also be specified using the colon and the numeric value of the parameter. For the model entries, either ID numbers or short notation may be used, as long as they are defined ModelTypes. For example the definition for a baffle tray may be:

Name=Baffle tray
# multiple parameter lines:
Parameters=Column diameter (m):120cm;Baffle area (%):50;Baffle spacing (m)
Parameters=Baffle type[Disk & donut|Single flow|Multiflow]:2 #default=Single
The [Model] definitions has instead of a "Models" entry the "Internals" entry, where the types of internals for which the model can be used are defined. It also has the "Parameters" where model parameters can be defined. For the definition of a baffle tray above the parameters are a column diameter (with meters as units and a default value of 1.2m), the baffle area (as percentage, default 50%), the baffle spacing (in meters, without any default), and the baffle type (with a pick list between 3 different types where the default is the single flow).

The models that are required for the baffle tray are a mass transfer model (MTC), a pressure drop model (PD), models for the vapor and liquid flow (VF and LF), as well as a design model (DSGN). Note that the design method is also considered as a model. A sample design method definition is:

Name=Fraction of flood
Internals=Bubble cap;Sieve;Valve
Parameters=Fraction of flooding:75%;Fraction of weeping:70%
Model parameters can be loaded from libraries (*.PAR) where for each internal different parameter values are stored. Note that a model defintion may also be repeated (with the same id, name, type, and operation) so to add another internalstype to the list to which the model applies. Short fields are optional, and have a maximum length of ten characters, used for displaying selected models etc. ID fields associate a unique number to the definition. Only for the internal and model definitions non-unique numbers are allowed. When the interface reads the definitions file it uses the Short descriptions to assign the ID's for the Operations, InternalTypes and ModelTypes, see Table 7.1 (the Short descriptions used by the interface are in parenthesis). Thus, you will have to use these Short descriptions but are free to change the ID numbers or names.

Table 7.1. Interface Assigned Types

Operation InternalType ModelType
Vapor-Liquid (VL) Discrete (DIT) Mass transfer coefficient (MTC)
Liquid-Liquid (LL) Continuous (CIT) Pressure drop (PD)
Vapor flow (VF)
Liquid flow (LF)
Entrainment (ENTR)
Holdup (HOLD)
Light liquid flow (LLF)
Heavy liquid flow (HLF)
Backmixing (BACK)
Design method (DSGN)
There are standard definitions for internal types that are built into the interface. These do not require the parameters entry to be specified. If they are specified, they -override- the internal buildin types. For example, the definition of the bubble cap tray internal is:
Name=Bubble cap tray
Short=Bubble cap
The builtin types are listed in Table 7.2. This table also lists defined vapor and liquid flow models and models for entrainment.

Table 7.2. Assigned Internals and Vapor flow, Liquid flow, and Entrainment Models

Internal Vapor Flow / Liquid Flow / Entrainment /
Light liquid flow Heavy liquid flow Backmixing
Bubble cap tray (1) Mixed (1) Mixed (1) None (1)
Sieve tray (2) Plug flow (2) Plug flow (2) Estimated (2)
Valve tray (3)
Dumped packing (4)
Structured packing (5)
Equilibrium stage (6)
RDC compartment (7)
Spray column stage (8)
An example of a model definitions is:

Internals=Bubble cap,Sieve tray,Valve tray
Assigned models for Mass Transfer Coefficient and Pressure Drop models are listed in Table 7.3. Model parameters may be read from parameter libraries (*.PAR) that have a format like the packing data files. The first line of such a ASCII text file must start with "# CHEMSEP xxxx" where xxxx must be replaced by the full name of the model (as defined in CHEMSEP.DEF). Upon choosing the library option in the interface the library will be automatically pre-selected if the name of the library file is the same as the short name of its model.

Table 7.3. Assigned MTC and PD Models

Mass Transfer Coefficient Pressure Drop
AIChE (1) Fixed (1)
Chan Fair (2) Estimated (2)
Zuiderweg (3) Ludwig 1979 (3)
Hughmark (4) Leva GPDC (4)
Harris (5) Billet-Schultes 1992 (5)
Onda et al. 1968 (6) Bravo-Rocha-Fair 1986 (6)
Bravo-Fair 1982 (7) Stichlmair-Bravo-Fair 1989 (7)
Bravo-Rocha-Fair 1985 (8) Bravo-Rocha-Fair 1992 (8)
Bubble-Jet (9)
Bravo-Rocha-Fair 1992 (10)
Billet-Schultes 1992 (11)
Nawrocki et al. 1991 (12)
Chen-Chuang (13)
Handlos-Baron-Treybal (14)
Seibert-Fair (15)
Kronig-Brink-Rowe (16)
Rose-Kintner-Garner-Tayeban (17)
Sherwood (20)

7.9 Author and program information

The authors can be reached through
  Regular mail:     Ross Taylor / Harry Kooijman
                    Department of Chemical Engineering
                    Clarkson  University
                    Potsdam, NY 13699, USA

  Telephone:        Ross Taylor    (315) 268 6652
                    Harry Kooijman (908) 771 6544

  Electronic mail:

  WWW page:
We typed most of the code with the Multi-Edit text editor. It allows us to switch between the source code of the drivers, interfaces, and column simulation executables, as well as the text files for the help and documentation. Each different type of file has its own commands associated with it (compile source code, run latex on the documentation, etc.). As the project is over several hundred thousands lines of source and text, the editor has proven to be very valuable to us.

The drivers and interfaces are written in Turbo Pascal (version 7.0), the column and flash programs are written in standard Fortran 77, which we compile with WATCOM Fortran and link with the CauseWay DOS extender. The source code for the simulation executables has also been successfully compiled and executed on a range of operating systems and platforms .


WordPerfect is a product of the WordPerfect Corporation. Microsoft Word is a product of Microsoft Incorporated. DOS4GW is a product of Rational. CauseWay is a product of Devore software Incorporated. Multi-Edit is a product of American Cybernetics. WATCOM F77 is a product of WATCOM Incorporated. Turbo Pascal is a product of Borland International Incorporated.


J.M. Prausnitz, T.F. Anderson, E.A. Grens, C.A. Eckert, R. Hsieh, J.P. O'Connell, Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria, Prentice-Hall, Englewood Cliffs, NJ (1980).

Chapter 8. FlowSheeting

By combining different ChemSep models you can actually simulate (small) flowsheets. Several other common unit operations (reactor, make-up stream, and stream splitter) are available to make this possible. Simulating flowsheets with the utility program fs is illustrated with several examples in this chapter.

8.1 Flowsheet Input File

The flowsheet utility uses a text file as input. You will have to make this file yourself (with, for example, the edit option in ChemSep's file menu). The ChemSep distribution contains various examples (in the fs directory). To explain the format of this file we will use an example where we simulate the production of diethylether from an ethanol-water (85%-15%) feed, as shown in Figure 8.1. Figure 8.1. Diethylether Production The conversion of the ethanol in the reactor is only 50% and the unreacted ethanol has to be recycled. Pure ether (99.5%) and water (with 1% ethanol) are the products. The reaction is:
  2 C H OH -> (C H ) O + H O
     2 5        2 5 2     2

In the flowsheet input file (ep.fs) you define the units and streams in your flowsheet. The input file consists of four parts. The first part consists of one section where all the components, units, streams, feeds, and stream estimates are declared as well as the output file, executable directory, and the method and accuracy used in solving the flowsheet.
Comment=Ethylether Production
! this is a comment
The section starts with the [Flowsheet] identifier. The next lines are part of this section until an empty line is found. You can add (non empty) comment lines in the input file by, for example, using any punctuation character to start the line. The flowsheet program looks for the specific keywords and if none is found it regards the line as a comment. However, comment lines do not get copied to the output file. In each section the keywords can be different. The keywords are the identifiers left of the equal signs in the above flowsheet section. They are not case-sensitive and you may enter them in any order (some restrictions do apply, though).

Comma's are used to separate the various components, streams, units, feeds, or stream estimates (therefore, they can not be part of a name. This is important as IUPAC component names sometimes use comma's. In that case omit the comma's in the component name. Currently the component names are not checked and it is assumed that the same components in the same order are used in all the units! This is very important. The case of the names you specify is not important but the one used in the flowsheet section is used in the rest of the report.

In our example we have the following streams we have given our streams numbers from 1 to 8 (instead of numbers names may also be used). We have only one feed, 1, and we will make a stream estimate for the recycle ( 7). Stream 2 is a stream that is not used (we need it since the flash program that simulates the mixer produces a vapor product stream which we have to assign). Our unit operations are a mixer, reactor, and two columns ( Sep-1 and Sep-2) where we separate the products. Our output file is ep.fs (the same name as our input file which will be overwritten) and the flowsheet method is Direct substitution with a convergence criterium of 0.001. The maximum number of iterations is specified here as 20.

The next part describes the units and, therefore, can consist of multiple sections. We have four units in our example. Each unit section starts with the [Unit] identifier and has five keywords: Name, File, Type, Inlets, and Outlets.




The specified unit name must match the one given in the flowsheet section. The unit file is the associated data file, for most units it is the Sep-file which was generated by ChemSep. The unit type is the program that must be run to simulate this specific unit. For a ChemSep flash that is 4g-flash.exe, for an equilibrium column 4g-col2.exe, and for a nonequilibrium column 4g-neq2.exe. This allows you to make your own unit simulation program that reads a datafile with feed and product section as in a Sep-file. In the case of the reactor this is done by the reactors.exe program (which is described below).

The units are connected by in- and outlets streams, which were declared in the flowsheet section above. Be sure to use the same names (case is not important) as otherwise the flowsheet will be incorrect or incomplete. If you just want to analyze a flowsheet or not to simulate the unit it, ommit the unit file and type (they will not be executed).

The next part consists of the [Feed] sections where the feed streams are defined. The specified feed name must match the one given in the flowsheet section (again, case is not important). Furthermore, the feed stream pressure, temperature, flowrate, and composition must be specified. Append units when specifying temperatures, pressures, and flows. The default units (which can be omitted) are temperatures in degrees Celcius, pressure in bar (absolute), and flows in mol per second. The flowrate and compositions can also be set by specifying the component flows, for example for the definition of the feed above we could use as well
where F is the list of component flows (default units mol/s). If all the component flows are specified, the total flowrate and compositions are computed. A partial list (like Z=*,0.85) can be specified as well, which is useful for supplying stream estimates.

All streams are reset at the start, and when the inlets are written to the unit-file only the specifications supplied in the input file are written to the unit file. If stream values are reset the values already present in the unit-file will be used.

The last part consists of [Estimate] sections where estimates of stream variables (flow, temperature, pressure, or composition) can defined. This part is optional and only required if you specified streams under the estimates keyword in the flowsheet section. Stream estimates have the same input format as feeds have, except for the different identifier, of course. In our case we estimate the recycle stream 7 to start with a better value of the flowrate to the reactor:

The files that make up our example all start with EP (ep.fs, ep-mix.sep, ep-reac.sep, ep-sep1.sep, ep-sep2.sep) and are included in the ChemSep distribution (in the fs directory).

8.2 Flowsheet execution

If the flowsheet input is correct the flowsheet program starts the execution of each units program with the specified file written to the ChemSep.Fil file in the default directory ( ChemSep programs read this file to obtain the Sep-file to simulate). The order of execution is the order as was specified in the flowsheet section (so you can manipulate it). If a unit file or type is missing execution of that unit is skipped.

A unit evaluation takes the unit inlet streams and puts them (in the same order!) into the file under the [Feeds] section. Then it runs the associated program and reads the [Top product] section as the first outlet stream, the [Bottom product] section as the second outlet stream, and each following [Sidestream section as the third, fourth, etc. outlet streams.

In the case that the flowsheet analysis finds one or more cycles it will inform you and set the recycle flag to start an iterative run. The criterium is the maximum relative (or absolute) differences in stream variables when any stream gets updated. For each iteration this is initially set to zero and computed over the simulation of all the units. Convergence is obtained if the criterium is below the specified accuracy (default is 10^-2) or if the maximum number of iterations (default is 20) is attained (In the case that an absolute criterium is used the temperature difference is divided by 10 and for the pressure by 10^4 to scale the differences).

Sometimes it may be handy to abort the simulation or to stop it and inspect certain streams during the simulation. This can be accomplished by holding the Shift keys or toggleing Caps Lock on. The program will beep and display the next unit to be simulated, the current attained convergence, and the current iteration number. A simple menu allows you to change the maximum number of iterations or the accuracy, display a stream, continue, or to quit the simulation. It will also allow you to swap to DOS to do other work, like loading Sep-files into the ChemSep interface to make changes or to evaluate the intermeadiate results. Typing "exit" will then return you to the flowsheet program (don't forget this).

Once the flowsheet is convergenced the output file is written, consisting of four parts: the input file (generated from the information read in), a flowsheet analysis, the mass balances, and a report of all the streams.

8.3 Flowsheet Analysis

The flowsheet analysis of our depropanizer example is rather simple. Flowsheet reports the incidence, adjacency, distance, and cycle matrices for the specified flowsheet.
Incidence Matrix
Unit\Flow:  1  2  3  4  5  6  7  8
Mixer       +  -  -           +
Reactor           +  -
Sep-1                +  -  -
Sep-2                      +  -  -

Adjacency Matrix
Unit:      Mixer Reactor   Sep-1   Sep-2
Mixer                  3
Reactor                        4
Sep-1                                  6
Sep-2          7

Distance Matrix
Unit:      Mixer Reactor   Sep-1   Sep-2
Mixer          4       1       2       3
Reactor        3       4       1       2
Sep-1          2       3       4       1
Sep-2          1       2       3       4

Cycle Matrix
Cycle:  Rank:  Streams:
C1      4      3,4,6,7

Stream:  Frequency:
3        1
4        1
6        1
7        1

Cycle:  Rank:  Units:
C1      4      Mixer,Reactor,Sep-1,Sep-2

Evaluation order from analysis = Reactor Mixer,Sep-1,Sep-2
These matrices can be useful in assessing the structure of the flowsheet and the flowsheet evaluation order of the units. The flowsheet analysis will also give an ordering of the units that might be better than the specified order.

8.4 Convergence

This short section displays whether the flowsheet was converged and, if iteration was required, the attained convergence criterium and number of iterations.
Flowsheet is solved.
Attained convergence = 0.000826326
Number of iterations = 10
Used relative differences.

8.5 Mass Balances

In the balances section total mass balance is given for each unit, and for the overall flowsheet. If a flowsheet is converged, each balance should equal zero or some small number (relative to the in- and outlet flows of the specific unit). The balances are first written as the names of the units and associated streams, then in the total molar flowrates:
Balances in mol/s

Total Molar Balances:
Mixer   = 1+7-2-3 = 20+19.7097-0-39.694 = 0.0157021
Reactor = 3-4 = 39.694-39.694 = 0
Sep-1   = 4-5-6 = 39.694-8.45862-31.2354 = -0.0000204891
Sep-2   = 6-7-8 = 31.2354-19.7097-11.5257 = 0
Overall = 1-2-5-8 = 20-0-8.45862-11.5257 = 0.0156797

Water Balances:
Mixer   = 3+2.799822-0-5.797864 = 0.00195764
Reactor = 5.797864-14.21037 = -8.412508
Sep-1   = 14.21037-0.0000697587-14.21033 = -0.0000232831
Sep-2   = 14.21033-2.799822-11.41044 = 0.0000614673

Ethanol Balances:
Mixer   = 17+16.66376-0-33.65003 = 0.0137314
Reactor = 33.65003-16.82501 = 16.82502
Sep-1   = 16.82501-0.0460432-16.77897 = 0
Sep-2   = 16.77897-16.66376-0.115257 = -0.0000505315

Diethylether Balances:
Mixer   = 0+0.246113-0-0.246111 = 1.920853E-06
Reactor = 0.246111-8.658619 = -8.412507
Sep-1   = 8.658619-8.412504-0.246113 = 1.949957E-06
Sep-2   = 0.246113-0.246113-3.203753E-08 = -3.203753E-08
Care must be taken that the individual component balances are also satisfied, they are written after the total molar balances. Also, a reactor total molar balance will not be zero if the reaction changes the number of moles in the mixture. A reactor component balance will not be zero if that component is involved in one of the reactions that has a nonzero conversion. In the case of our example we see that the reaction rate is 8.4 mol/s dietheylether.

8.6 Stream Report

The stream report lists the streams after the execution has stopped.
Stream                        1
Temperature (C)               40
Pressure (bar)                1.01325
Flowrate (mol/s)              20
Zwater                        0.15
Zethanol                      0.85
Zdiethylether                 0

Stream                        2 is zero

Stream                        3
Temperature (C)               58.68402
Pressure (bar)                1.01325
Flowrate (mol/s)              39.694
Zwater                        0.146064
Zethanol                      0.847736
Zdiethylether                 0.00620021

Stream                        4
Temperature (C)               49.85001
Pressure (bar)                1.01325
Flowrate (mol/s)              39.694
Zwater                        0.357998
Zethanol                      0.423868
Zdiethylether                 0.218134

Stream                        5
Temperature (C)               34.76102
Pressure (bar)                1.01325
Flowrate (mol/s)              8.45862
Zwater                        8.24705E-06
Zethanol                      0.00544335
Zdiethylether                 0.994548

Stream                        6
Temperature (C)               77.918
Pressure (bar)                1.01325
Flowrate (mol/s)              31.2354
Zwater                        0.454943
Zethanol                      0.537178
Zdiethylether                 0.0078793

Stream                        7
Temperature (C)               76.11902
Pressure (bar)                1.01325
Flowrate (mol/s)              19.7097
Zwater                        0.142053
Zethanol                      0.84546
Zdiethylether                 0.0124869

Stream                        8
Temperature (C)               96.62302
Pressure (bar)                1.01325
Flowrate (mol/s)              11.5257
Zwater                        0.99
Zethanol                      0.01
Zdiethylether                 2.77966E-09

8.7 Commandline Options

The flowsheet program has several options which can be specified on the commandline when you start it. These are described below in some examples followed by explanation (after the equal).
\x        = Skip execution
\r        = Force iteration
\utF      = Set default temperature units to degrees Fahrenheit
\utpsia   = Set default pressure units to psia
\utkmol/s = Set default flow units to kmol/s
\da       = Use absolute difference criterium
\dr       = Use relative difference criterium (default)
Although the streams in the input file can be defined by any units, the output file will be written with the default units (including the part with the input file!). Unit definitions are read from the FS.UDF file, or if that file is not found, from the CHEMSEP.UDF file.

8.8 Other Unit Operations

If you want to simulate a complex flowsheet you need several other unit operations models besides flash and column operations. The flash unit can also be used to model heaters, coolers, or pumps. However, a nonsense stream has to be added as most streams are either a vapor or liquid. Here we discuss several other unit operations which are commonly used in flowsheets. They have there own little programs for which we supply the pascal code (to be compiled with Borland Pascal, v7 or later). Most of them are limited to only a couple of hundred lines which mostly cover the in- and output. If you require other unit operations you could code them yourself (we welcome your pascal code for unit operations to enhance the flowsheet capabilities).

8.8.1 Simple Reactor

A unit operation used in almost any flowsheet simulating a chemical plant is the reactor. We have implemented a simple reactor model which can handle multiple reactions with specified conversion(s). The conversion is adapted in case one of the component flows is constraining. The input file has a similar style as that of the sep-files and flowsheet input file. The reactor specifications (conversion, outlet temperature, pressure drop, and the reaction stoichiometry) are given under the [Reactor] section which is followed by a [Feeds] section following the sep-file format (only one feed is allowed). The first line in the [Reactor] section contains the reactor's name, followed by the outlet temperature, pressure drop, number of components, and number of reactions. Then, for each reaction, lines for the base component, the conversion based on the base component feedflow, and the reaction's stoichiometry coefficients (coefficients for each component; negative for reactants, positive for products).
Production of benzene by hydrogenation of toluene
   C5H5CH3 + H2 ->   CH4 + C6H6
 2 C5H5CH3 + H2 -> 2 CH4 + (C5H5)2
 Toluene, Hydrogen, Methane, Benzene, Diphenyl

400 K outlet temperature
10000 Pa pressure drop
5 components
2 reactions
1 base component r1
0.9 conversion r1
-1 -1 1 1 0 stoichiometry coef. r1
1 base component r2
0.2 conversion r2
-2 -1 2 0 1 stoichiometry coef. r2

1 number
400 K temperature
101325 Pa pressure
5 components
0.1 kmol/s toluene
0.1 kmol/s hydrogen
0 kmol/s methane
0 kmol/s benzene
0 kmol/s diphenyl
After running the reactors program the output is written to the same file and contains a small [Results] section which reports whether the specified convergence(s) was attained followed by a [Top product] section with the reactors outlet stream.
Conversion(s) on base components was limited by component 1
to 90% of specified conversion(s)
If there are multiple reactions and one (or more) of the component feedflows is constraining the reactions, the conversions are all equally affected. For the ether reactor of our example the reactors input file looks like:
323 K outlet temperature
0 Pa pressure drop
3 components
1 reactions
2 base component r1
0.5 conversion r1
1 -2 1 stoichiometry coef. r1

1 Number
1 Feed state T & P
331.834 Temperature
101325 Pressure
* Vapour fraction
3 componentflows
0.00579786 Component 1 flow
0.03365 Component 2 flow
0.000246111 Component 3 flow

8.8.2 Make-Up Feeds

Often, when a flowsheet contains a recycled solvent or base material, a (small) make-up feed is required. Since the exact flowrate of the make-up feed is unknown, we made a unit operation that will add a make-up feed (the first feed) to another stream (the second feed) to obtain a specified total flowrate.
Make-up Unit
.06 Total flowrate (always in kmol/s)

2 Number
1 Feed state T & P
1 stage
373.15 Temperature
100000 Pressure
* Vapour fraction
3 componentflows
0 Component 1 flow
0 Component 2 flow
0.001 Component 3 flow
1 Feed state T & P
1 stage
373.15 Temperature
101325 Pressure
* Vapour fraction
3 componentflows
5.041514E-09 Component 1 flow
0.0000207905 Component 2 flow
0.0598143 Component 3 flow
When the make-up program runs it appends a small [Results] section displays the total flowrate of the make-up feed followed by a [Top product] section containing the resulting stream.
0.0001649 = Make-up flow
Since the make-up feed has to be specified at the beginning this could cause the mass balances to be incorrect. Therefore, if the flowsheet encounters a make-up unit, and the first inlet stream is declared as a feed, its flowrate is adapted on writing the output.

8.8.3 Stream Splitter

Similarly to the Make-up unit, flowsheets containing recycles sometimes need to purge some part of the recycle to prevent build-up of various components in the cycle. This means that the recycle stream must be split into two parts. The splits program implements the stream split operation. The input file for the splitter consists of a [Splitter] section, containing the name of the unit and a splitfactor (which can range from zero to unity):
Recycle purger
This is followed by the regular [Feeds] section (see for example, the reactor description), where only one feed may be specified. After running the splits program top and bottom product streams will be appended to the input file, where the top product stream has a flowrate equal to the feedflowrate multiplied with the splitfactor, and the bottom product contains the rest of the feedflow.

8.9 Examples

Here we discuss several examples which have been run with the flowsheet program. We will not supply all the details as they can be found in the various files that are distributed with ChemSep. Note that all these examples can be solved using the nonequilibrium column models.

8.9.1 Extractive Distillation (PH)

We need to separate a equimolar mixture of methylcyclohexane (MCH) and toluene, and do this by extractive distillation with phenol as solvent. The flowsheet is shown in Figure 8.2. Valve trays are used for both the columns using the design mode nonequilibrium simulator. The Phenol recycle is cooled to 100 C. For a high purity of the products the solvent feed to MCH/Toluene feed ratio as well as the reflux ratio need to be sufficiently high (for the extractive column). Figure 8.2. Extractive Distillation

8.9.2 Distillation with a Heterogeneous Azeotrope (BW)

Water and butanol form an azeotrope, so that they cannot be seaparated by conventional distillation. However, at not too high temperature, they form two liquid phases; an aqueous phase with little butanol and a butanol phase with a large mole fraction of water. We can use a decanter to separate the two liquid phases. A feed which is predominatly water but contaminated with butanol (1 mole%) can be separated into two pure products using such a decanter. The flowsheet is shown in Figure 8.3. The first column produces a water bottom product, and the vapor is fed to a condenser/cooler and then to the decanter. There we obtain a water- and a butanol-rich stream which get recycled to the columns. In a second column we can then produce butanol as the butanol-rich recycle contains much more butanol than the butanol-water azeotrope.

Figure 8.3. Distillation with Decanter

8.9.3 Distillation of a Pressure Sensitive Azeotrope (MA)

Methanol and acetone form also an azeotrope. However, the composition of the azeotrope is sensitive to the pressure. We can make use of this to separate the two components into pure products by operating two columns at different pressures to change the azeotrope composition. The separation of an equimolar feed of methanol and acetone is shown in Figure 8.4. This type of azeotropic distillation is rare as the azeotrope composition needs to be quite sensitive to the pressure in order to obtain a recycle stream which is not unreasonably large, and that the pressures are such that no special columns or equipment is required.

Figure 8.4. Azeotropic Distillation at Two Pressures

8.9.4 Petyluk Columns (PETYLUK)

To lower the energy consumption of separation trains, two columns separating 3 components can be replaced by one column with a condenser and reboiler plus one column without condenser and reboiler. The feed is fed to this (first) column, which receives a vapor to the bottom and a liquid at the top from the (second) column with the condenser and reboiler. In turn, the products of the first column are fed into the second column at the sidestream stages. Figure 8.5 shows such a configuration for the separations of three alcohols. To get this flowsheet to run requires that you can solve both the columns separately, which is not easy. The second column needed some specific user initialization information to run. Afterwards convergence can be obtained more easily by using the old results as initialization.

Figure 8.5. Petyluk Columns

8.9.5 Extraction with Solvent Recovery (BP)

Extraction can be used to separate aromatics from parrafins. This is a common type of separation in the crude oil refining. a simplified example is shown here where we separate benzene from n-pentane, using an extraction with sulfolane as solvent. The extractor is a rotating disk contactor (RDC) operating at 50 C. The solvent is recovered from the extract by "ordinary" distillation. UNIQUAC parameters for the components are given in Table 8.1. The extraction is dependent on the temperature as can be seen from the interaction parameters, which can be used to manipulate the separation. See Figure 8.6 for the flowsheet.

Table 8.1. UNIQUAC LLE interaction parameters (K)

Components i-j A_ij A_ji
$25^o C$:
pentane-benzene 295.38 -171.63
pentane-sulfolane 532.04 187.84
benzene-sulfolane 151.02 -36.08
$50^o C$:
pentane-benzene 179.86 -95.70
pentane-sulfolane 375.93 247.77
benzene-sulfolane 131.51 -6.28
Figure 8.6. Extraction with Solvent Recovery