Some aspects of rate-based modelling and simulation of three-phase distillation columns

Computers and Chemical Engineering 25 (2001) 603– 612 www.elsevier.com/locate/compchemeng

Some aspects of rate-based modelling and simulation of three-phase distillation columns Egon Eckert *, Toma´sˇ Vane˘k Department of Chemical Engineering, Institute of Chemical Technology, Technicka´ 5, 166 28, Prague 6, Czech Republic Received 10 May 2000; accepted 5 January 2001

Abstract The application of the bate-based approach to the dynamic and stationary modelling of three-phase distillation columns is presented. The main problem can be seen in finding adequate models for the mass and heat transfer over phase interfaces to be utilised in the Maxwell-Stefan equations. While for the vapour– (continuous) liquid interface, there are a number of methods and data available, the interface between continuous and dispersed liquid phases has not been studied for the case of distillation processes. A reasonable alternative approach to modelling of a three-phase distillation column can be found where the liquid–liquid interface is treated as equilibrium. The adequacy of this modification is supported by phenomena observed on distillation trays, e.g. strong agitation of liquid phases. The resulting combined non-equilibrium and equilibrium model is compared to the classical equilibrium model and also to experimental data for ethanol– water– cyclohexane separation in a number of examples. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Rate-based approach; Three-phase distillation; Non-equilibrium tray model; Simulation

1. Introduction The rate-based approach (Krishnamurthy & Taylor, 1985a,b,c) has become a standard in the steady-state and dynamic modelling and simulation of separation equipment in last 15 years, partially due to the availability of powerful computer systems and effective model solving simulation programs. While the modelling principles for the rate-based approach had been established and verified quite a long time ago, the problem of a suitable description of phenomena taking place at interfaces between various phases is in most cases still open. The non-equilibrium, i.e. rate-based, concept is undoubtedly superior to the concept of an equilibrium stage, which has to use various types of corrections to fit the real performance of a separation equipment (Seader & Henley, 1998). At the same time, it is well known that the equilibrium concept produces more or less unreliable results and, moreover, it is not fully suited for the investigation of the dynamic behaviour of separation columns. * Corresponding author. E-mail address: [email protected] (E. Eckert).

The rate-based approach to the modelling and simulation of separation columns is suitable for any type of separation processes (Taylor & Krishna, 1993; Wesselingh, 1997; Lao & Taylor, 1994). As the most common types we can recognise, for example, the V–L separation processes (i.e. two-phase distillation), G–L (absorption), L1 –L2 (extraction) or V–L1 –L2 separation (three-phase distillation), the last case being in the focus of this paper. The problem is in finding suitable models for mass and heat transfer coefficient calculations for a sufficient number of phase interfaces to be accordingly inserted into the model of the mass and heat transfer between appropriate two interacting phases. Widely used for this purpose are the Maxwell-Stefan equations (Taylor & Krishna, 1993). To describe the mass transfer between phases the two-film theory is applied here, i.e. thermodynamic equilibrium is postulated at the interface and it is assumed that mass transfer resistance is limited to boundary layers on both sides of the interface. The case of a three-phase distillation requires to decide what kind of mass and heat transfer models should be used and to find out if some sources of data are available. In this contribution these aspects are

0098-1354/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 6 4 0 - 8

604

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

discussed and on examples based on data from a real two- or three-phase distillation column it is shown, how reasonable simplifications and modifications of the pure rate-based model can be introduced. Some of the numerical problems with standard solvers, e.g. SPEEDUP, will be also discussed.

2. Modelling concepts For ordinary distillation, there is only one interface but for the V–L1 – L2 separation we recognise three interfaces. First of all, we must differ between the two liquid phases. The prevailing phase is called continuous (LC) while the other is expected to be present in the form of drops dispersed in the continuous phase and accordingly is called the dispersed liquid phase (LD). For the V–LC interface the problem can be usually sufficiently analysed, but there are almost always difficulties in other cases. The overall rate-based model must be then simplified or abandoned. Generally, the possibilities for the modelling of a phase interface are as follows: “ full heat and mass transfer model — no simplification;

Fig. 1. Simplified scheme of a three-phase (V–LC–LD) stage.

“

infinite product of mass and heat transfer coefficients with interfacial area on one or both sides of the interface — corresponding to the equilibrium concept if on both sides; “ no mass nor heat transfer — corresponding to the case, when the phases are not in a direct contact or the contact time is negligible. When designing a mathematical model of a distillation stage with a potential occurrence of two liquid phases (V–LC –LD) using the non-equilibrium concept the V–LC interface can be treated in the same way as for two-phase distillation, except the problem of physical properties of the liquid phase as a whole. Here we need to consider a unique liquid phase for which a set of average properties must be delivered depending on the tray model chosen. The main uncertainty is in the LC –LD interface. In this case the flow of liquid phases across and between column plates is co-current and the estimation methods for mass transfer coefficients and interfacial area employed in case of liquid –liquid extraction are useless for three-phase distillation. We expect the second liquid phase to be present in the form of drops dispersed in the prevailing continuous liquid phase (Fig. 1). If we would want to develop a more reliable rigorous model of the process, it would be necessary to know the appropriate mass and heat transfer coefficients and splitting and/or coalescing rate constants (Eckert, 1995). Unfortunately, no estimation methods are known for this purpose and we could only choose some compromise until the knowledge of such phenomena increases. Obviously, certain possibilities would appear if we take into account some model simplifying assumptions, i.e. (a) The volume occupied by any phase is always ideally mixed, i.e. the phase temperature and composition are uniform in the entire volume and at the same time identical to outlet values. (b) In the case when both liquid phases are present, their composition is in equilibrium and their temperatures are identical. (c) No heat nor mass transfer between dispersed liquid and vapour phases takes place as the contact between vapour bubbles and dispersed liquid drops is incidental only and the contact time is likely to be very short. (d) The work of expansion is neglected. Assumption (a) enables to simplify the model of the stage but we can imagine the usage of more sophisticated approaches (e.g. plug flow model). Generally in liquid –liquid extraction processes the mass transfer resistance on both sides of the interface is considerably higher than in vapour –liquid distillation. However, Herron, Kruelskie, and Fair (1988) claimed that the case of three-phase distillation is diametrically different. First, Priestly and Ellis (1978) showed that the gas agitation in liquid –liquid extraction can decrease the mass transfer resistance significantly. In distillation the

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

vapour phase flow rates are even 100 times or more higher than used by Priestly and Ellis (1978). Moreover, in the real experimental column modelled in examples presented further it was observed that the second liquid phase if existing is in form of very small drops. This is the base for the acceptance of assumption b). These assumptions served for the development of a combined non-equilibrium and equilibrium model of a separation stage depicted in Fig. 1 and in detail described further. 3. Dynamic non-equilibrium model of a three-phase distillation stage

605

Energy balance: du Lj C C C D D D D =L C j − 1H j − 1 − L j H j + L j − 1H j − 1 − L j H j dt I

C + % N*i, jH C *− T Lj )+ Q C i, j + h j a*(T j j j

(8)

i=1

Internal energy hold-up: C D D L 0= u Lj − n C j H j − n j H j + s j pj

(9)

where the definitions of molar hold-ups and of liquid volume are: I

nC = % nC j i, j

(10)

i=1

We shall present model equations separately for bulk vapour phase, equilibrium mixture of liquid phases and the interface between vapour and continuous liquid phases (Eckert and Vane˘k, 1999a,b).

(11)

i=1 C D D s Lj = n C j /z j + n j /z j

1.5 D C D L 0 = s Lj (L C j + L j )− (n j + n j )[(s j /Aj − Wj )/Cj ]

The impact of mass and heat transfer through the phase interface is given by component molar fluxes and overall heat flux. Component balances: dn V i, j =Vj + 1yi, j + 1 − Vj yi, j −N*i, j, i= 1, …, I dt

(1)

Energy balance: I

du V V = Vj + 1H V j + 1 −Vj H j − % N* i, jH i, j dt i=1 V j

− h a*(T − T*) j +Q

V j

(13)

Relation for the dispersed liquid flow: 0=





nD LD j j − LC LD s Lj j j + D zC zj j

(14)

Alternative equation in case of single liquid phase: 0=LD j

V j

(12)

Tray hydraulics – Francis weir formula:

3.1. Bulk 6apour phase

V j j

I

D nD j = % n i, j

(14a)

Stage volume balance (2)

0= sj −

nV nC nD j j j − − zV zC zD j j j

(15)

Internal energy hold-up: V V 0= u V j −nj Hj +

nV j pj zV j

3.3. Vapour/liquid interface (3)

where the definitions of molar hold-up and fractions are as follows: I

V nV j = % n i, j

(4)

i=1

yi j =

nV i, j nV j

(5)

Pressure drop equation (simplified for the case of a low pressure drop): 0= Vj F /FB − z (pj −pj − 1 −Dp0)/nj V j

V j

(6)

For generalised Maxwell-Stefan equations the binary mass transfer coefficients are calculated according to the appropriate stage type. Using I-1 by I-1 dimensional multicomponent mass transfer coefficient matrices [k] (Taylor and Krishna, 1993; Eckert and Vane˘k, 2000b) the model of the interface between the vapour phase and the continuous liquid phase can be written in the following matrix form: Mass transfer rates: I

[0]= [N*]− % N*i, j [yj ]− h[k V, j j ]a*[y j j − y*] j

(16)

i=1 I

, C % N*i, j [x C ]a*[x *− xC [0]= [N*]− j j ]− [k j j j j ]

3.2. Continuous and dispersed liquid phases

(17)

i=1

Component balances: L i, j

dn C D D C C =L C j − 1x i, j − 1 +L j − 1x i, j − 1 +N* i, j − L j x i, j dt D −LD j x i, j, i= 1, …, I

(7)

Eq. (16) is modified by introducing a scalar multiplicative factor h which is used for the correction of the predicted mass transfer coefficients (see Section 5). Value h= 1 corresponds to uncorrected values. Energy balance:

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

606 I

V V C [0] =h V *j −T C j a*(T j j − T*) j + % N* i, jH i, j −h j a*(T j j ) i=1 I

− % N*i, jH C i, j

(18)

i=1

Thermodynamic equilibrium: 0 = K*i, jx*i, j − y*i, j,

i= 1, …, I

(19)

Interface mole fraction summations: I

0= % y*i, j −1

(20)

i=1 I

0= % x*i, j −1

(21)

i=1

3.4. O6er6iew of model equations and unknowns Eqs. (1) –(3), (6) –(9), (14) – (21) constitute a system of 5I + 9 simultaneous equations with 5I + 9 unknowns V V L L L C D nV x*, N*, j , u j , T j , Vj, nj , u j , T j , L j , L j , pj, y*, j j j T*. Additional 2I + 2 equations are delivered by the j flash procedure (Heinichen, 1994) to calculate remainD C D ing unknowns x C j , x j , n j , n j . This special procedure was designed to overcome troubles that might appeared when the overall liquid composition approached the critical point at the liquid – liquid binodal envelope. In case of a single liquid phase detected this procedure yields only the following trivial set of equations: D 0=x C i, j −x i, j,

i= 1, …, I

n Li, j

0 =x C i, j −

I

, i= 1, …, I

(22) (23)

% n Lk, j k=1

0=n

D j

(24) I

L 0=nC j − % n i, j

(25)

i=1

4. Solution methods and tools The non-equilibrium concept has been incorporated into the commercial simulation system ASPEN PLUS in the form of a standard steady state two-phase separation module RATEFRAC. Up to now we are not informed about other attempts to design such modules in other commercial simulation programs. When developing more complex, i.e. three-phase non-equilibrium steady-state and/or dynamic, models it is necessary to have powerful and reliable numerical tools for their solution. Chemical engineers can profit from the existence of a number of equation-oriented simulation programs as SPEEDUP, Aspen Dynamics, gPROMS, etc. which usually contain the following important subsystems:

“

A thermodynamic database and a library of thermodynamic functions for the calculation of properties of pure components and their mixtures. “ A library of robust and stable numerical methods capable to solve large systems of algebraic or differential and algebraic equations (DAE) even in case when the model has discontinuities or some other numerically difficult features. For the construction of the model we can usually use a modelling language combined with standard modules and parts coded in higher programming language, e.g. FORTRAN 77. The modelling language typically serves for the design of an overall tray and column model while the external FORTRAN modules are well suited for detailed description of the mass and heat transfer processes. Nevertheless, our experience with the SPEEDUP system has shown some numerical and technical problems, which might be general when solving models possessing various local discontinuities or singularities. For example: “ SPEEDUP usually requires an extremely good starting point which is in most cases unavailable. In case of complex processes it is necessary to build the model gradually and to utilise previous solutions as starting points. “ If we use the homotopy method to move from one steady-state solution to another it is impossible to continue the solution through the binodal envelope. On the other hand, this problem is possible to overcome if we formulate the dynamic model and use integration methods for obtaining the new stationary solution. “ The bounds for variable values can have a negative influence on the convergence when for some variable the upper or lower bound is reached. “ Usage of standard modules with an internal iteration procedure usually leads to a failure in finding the overall solution. The decision about the number of phases serves as an example. It must be incorporated into the model along with appropriate equations. Build-in modules for three phase flash calculations (FLASH3 family) could not be used because the overall iteration process diverges.

5. Examples The suite of problems we have used for testing of the proposed approach concerns the dehydration of ethanol using cyclohexane as the entrainer in various continuous or batch columns. This system exhibits two - as well as three - phase regions when changing the composition and thermodynamic conditions. There has been an extensive experimental work done using a real laboratory column (Mu¨ller et al., 1997; Klein, 1996; Du¨x,

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

Fig. 2. Scheme of the model of a batch distillation column operating at the total reflux conditions. Saturated vapour corresponds to the output from the reboiler (stage 9).

607

1996). The scheme of the batch laboratory column with eight real stages equally used as the structure of the column model is depicted in Fig. 2. In case of rate based modelling the number of stages in the model and the column are of course the same. Above stage 1 is a total condenser from which the oneor two-phase liquid is entirely returned back to the column, i.e. during the experiment the column was performing at total reflux conditions. For modelling purposes the condenser is treated as a stage at equilibrium conditions without any hold-up at which the incoming vapour condenses and instantly forms one or two liquid phases at the bubble point. The specifications for all examples are presented in Table 1. The vapour flow rate entering the column was estimated from the measured heat duty of the total condenser. Examples 1 and 2: These examples correspond to experiments No. 281195 and 201295B documented by Klein (1996) and also have been studied by Eckert and Vane˘k (2000a). They differ in the feed composition and accordingly in the presence of second liquid phase on stages. While in Example 1 the results prove the existence of two liquid phases on all stages in Example 2 only some of upper stages appear to be three-phase. Example 2 is closer to conditions used for industrial dehydration of ethanol. It could be observed from Figs. 3 and 4 that the results of neither approach gave a satisfactory match with experimental data.

Table 1 Specifications for Examples 1–3 Parameter type

Value or source Example 1

Example 2

Example 3

Components

Ethanol (1)

Water (2)

Cyclohexane (3)

Flow rate, mol s−1 Temperature, K Pressure, Pa Composition–mole fractions (1) (2) (3) Tray Type Exit weir height, m Active tray bubbling area, m2 Flow path length, m Constant C in Eq. (13), s 2/3 m−1 Dry tray pressure drop, Pa Pressure drop parameter n, Pa s−1 m−3 Density FB, kg m−3 Thermod. model Activity coefficients

8.09×10−3

Vapour feed

Mass transfer Heat transfer

Binary mass-transfer coefficients Heat-transfer coefficients, W m−2 K−1

6.3×10−3 10.98×10−3 Dew point (calculated) 1.013×105 0.007 0.926 0.982 0.239 0.049 0.005 0.754 0.025 0.013 Cap tray 0.011 0.00155 0.0865 38 140 54000 1.20325 Modified NRTL method, data published by Connemann, Gaube, Karrer, Pfennig, and Reuter (1990) Estimated by modified AIChE method for cap trays (AIChE, 1958) Estimated by Chilton–Colburn analogy (Taylor and Krishna, 1993)

608

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

Fig. 3. Example 1. Mole fraction profiles in vapour phase-experimental data (Exp), simulation results using the rate-based approach (NE) and the equilibrium approach (EQ).

Fig. 4. Example 2. Mole fraction profiles in vapour phase-experimental data (Exp), simulation results using the rate-based approach (NE) and the equilibrium approach (EQ).

The experimental data are situated approximately in the middle between equilibrium and non-equilibrium profiles. Generally, the results of the equilibrium modelling usually correspond to more ‘optimistic’ behaviour of individual stages while the rate-based approach is more conservative in the prediction of the separation on stages. There are a number of potential reasons, e.g. “ The presence of the second liquid phase on some stages. “ Non-adequacy of the assumption of ideal mixing of phases on the stage. “ Non-adequacy of AIChE estimation methods for the laboratory column in consideration.

Example 3: In order to exclude the possibility that the observed discrepancies are caused by the presence of second liquid phase we tried to model the experiment No. 171095 with equal components and performed in the same column as in previous cases but entirely in the two-phase (V–L) region (Du¨x, 1996). The comparison between rate-based approach and the equilibrium approach with experimental data is given in Fig. 5. It is obvious that there is a similar discrepancy between computed and experimental profiles as in previous results. As a potential reason we could see that the mass transfer coefficients provided by the AIChE method are for the small laboratory column in consid-

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

eration underestimated. We decided simply to modify the AICHE method using a scalar multiplicative empirical correction factor h for vapour phase, see Eq. (16). The determination of the value for this factor needed some numerical experiments. The results for value h= 1.8 proved that the hypothesis presented above might be realistic, see Fig. 6 for the comparison of the profiles of mole fractions in the vapour phase and Fig. 7 for temperature profiles. With the same value of the correction factor we were able to receive a very nice match with experimental data also for Examples 1 and 2, see Figs. 8 and 9. The

609

temperature profiles were very similar to those of Example 3 and are not presented here.

5.1. Dynamic beha6iour In order to illustrate the complexity of the numerical solution of dynamic models where qualitative changes (e.g. discontinuities) occur we add also Fig. 10, which shows the response of mole fraction of ethanol in both liquid phases on stage 2 for Example 2 when the composition of the vapour feed changes from Ž0.926; 0.049; 0.025 to Ž0.95; 0.049; 0.001 during a short

Fig. 5. Example 3. Mole fraction profiles in vapour phase-experimental data (Exp), simulation results using the rate-based approach (NE) and the equilibrium approach (EQ).

Fig. 6. Example 3. Mole fraction profiles in vapour phase-experimental data (Exp) and simulation results using the rate-based approach (NE) with the correction factor h=1.8.

610

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

Fig. 7. Example 3. Temperature profiles-experimental data (Exp) and simulation results using the rate-based approach with the correction factor a= 1.8. I = V – LC interface, L = Continuous liquid phase, V = Vapour phase.

Fig. 8. Example 1. Mole fraction profiles in vapour phase-experimental data (Exp) and simulation results using the rate-based approach (NE) with the correction factor h=1.8.

time interval (200 s). After a while this caused the disappearance of the dispersed liquid phase. Setting back the composition of the feed at time 3000 s resulted in the reappearance of the second liquid phase.

6. Conclusion The more detailed rate-based steady state as well as dynamic modelling of separation processes can be useful for the analysis and design of real industrial separation columns. The resulting solution of the model, if

reached, can influence the column design and performance with direct economical effects. Moreover, for this type of column model the use of stage or column efficiencies is not needed, thus overcoming any ambiguous determination of global column efficiency, usually taken as 70% (Mu¨ller et al., 1997). There are certain difficulties resulting from our insufficient knowledge how to calculate certain physical and transport properties of complex mixtures and how to estimate the performance of column trays with various constructions. However, this contribution shows that even in case of more complicated type of separa-

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

611

Fig. 9. Example 2. Mole fraction profiles in vapour phase-experimental data (Exp) and simulation results using the rate-based approach (NE) with the correction factor h= 1.8.

Fig. 10. Example 2. Time dependency of mole fractions of ethanol in liquid phases on stage 2 — disappearance and reappearance of the second (dispersed) liquid phase.

tion process the rate-based modelling principles can be preserved despite the fact that some modifications are necessary.

Acknowledgements The authors appreciate the support of the Volkswagen-Stiftung (Germany) (Project no. I/71413Ak), of the Grant Agency of Czech Republic (Project no. 104/97/ 0916) and fund MSM 223400007.

Appendix A. Nomenclature a A C h H I k K L

interfacial area effective cross sectional area height over weir coefficient heat transfer coefficient molar or partial molar enthalpy number of components mass transfer coefficient equilibrium constant flow rate of the liquid phase

E. Eckert, T. Vane˘k / Computers and Chemical Engineering 25 (2001) 603–612

612

n N p Dp0 Q s t T u V W x y h F z n

molar holdup mass transfer rate pressure dry tray pressure drop heat flux volume time temperature energy hold-up flow rate of the vapour phase weir height mole fraction in liquid phase mole fraction in vapour phase empirical correction factor, see Eq. (16) density molar density pressure drop parameter

Subscripts B air i component i j stage j Superscripts C continuous liquid phase D dispersed liquid phase L mixture of liquid phases V vapour phase * vapour-continuous liquid interface

References AIChE, 1958. Bubble Tray Design Manual. New York: AIChE. Connemann, M., Gaube, J., Karrer, L., Pfennig, A., & Reuter, U. (1990). Measurement and representation of ternary vapour– liquid– liquid equilibria. Fluid Phase Equilibrium, 60, 99–118. Du¨x, A. (1996). Verification of an equilibrium model in dynamic simulation of three– phase distillation. Diploma project. RWTH Aachen (in German). Eckert, E. (1995). A nonequilibrium model of distillation stage with possible occurrence of two liquid phases. Collec. Czech. Chem. Commun., 60, 2085– 2096.

.

Eckert, E., and Vane˘k, T. (1999a). Simplified rate-based modelling and simulation of three-phase distillation columns. Proceedings of the 18th IASTED Int. Conference MIC’99, Anaheim: ACTA Press, pp. 223– 226. Eckert, E., and Vane˘k, T. (1999b). Rate-based modelling and simulation of three-phase distillation. Proceedings of the 45th Conference CHISA’99 (CD ROM form). Srnı´, October 18–21, Czech Republic. Eckert, E., and Vane˘k, T. (2000a). Some aspects of rate-based modelling and simulation of three-phase distillation columns. Computer-Aided Chemical Engineering, vol. 8, Elsevier Science, pp. 313– 318. Eckert, E., and Vane˘k, T. (2000b). Advances in rate-based modelling and simulation of three-phase distillation columns. Proceedings of the 14th International Congress CHISA’2000 (CD ROM form). August 27 – 31, Prague, Czech Republic. Heinichen, H. (1994). Dynamic negative liquid– liquid flash. Diploma project. ENSIGC Toulouse/RWTH Aachen. Herron, C. C., Kruelskie, B. K., & Fair, J. R. (1988). Hydrodynamics and mass transfer on three-phase distillation trays. American Institute of Chemical Engineering Journal, 34, 1267– 1274. Klein, W. (1996). Experimental investigation of the stationary and dynamic behaviour of heteroazeotropic distillation. Diploma project. RWTH Aachen (in German). Krishnamurthy, R., & Taylor, R. (1985a). A nonequilibrium stage model of multicomponent separation processes I. American Institute of Chemical Engineering Journal, 31, 449– 456. Krishnamurthy, R., & Taylor, R. (1985b). A nonequilibrium stage model of multicomponent separation processes II. American Institute of Chemical Engineering Journal, 31, 456– 465. Krishnamurthy, R., & Taylor, R. (1985c). A nonequilibrium stage model of multicomponent separation processes III. American Institute of Chemical Engineering Journal, 31, 1973– 1985. Lao, M., & Taylor, R. (1994). Modeling mass transfer in threephase distillation. Ind. Engineering Chem. Res., 33, 2637–2650. Mu¨ller, D., Marquardt, W., Hauschild, T., Ronge, G., & Steude, H. (1997). Experimental validation of an equilibrium stage model for three-phase distillation. IChemE Symposium Series, 142, 149– 159. Priestly, R., & Ellis, S. R. M. (1978). Gas agitated liquid extraction columns. Chem. Ind., 19, 757– 760. Seader, J. D., and Henley, E. J. (1998). Separation process principles. New York: Wiley. Taylor, R., and Krishna, R. (1993). Multicomponent Mass Transfer. New York: Wiley. Wesselingh, A. (1997). Non-equilibrium modelling of distillation. Trans. IChemE, 75, 529– 538.