entropy

ISSN 1099-4300

© 2007 by MDPI www.mdpi.org/entropy

Full Research Paper

Distillation of a Complex Mixture. Part I: High Pressure Distillation Column Analysis: Modeling and Simulation

Douani Mustapha*, Ouadjenia Fatima and Terkhi Sabria

Department of Chemical Engineering, University of Mostaganem, Mostaganem, Algeria E-Mails: douani_mustapha@yahoo.com; fatouadj@yahoo.fr.

* Author to whom correspondence should be addressed. E-Mail: douani_mustapha@yahoo.com Received: 7 December 2006 /Accepted: 14 April 2007 /Published: 25 April 2007

Abstract: In this analysis, based on the bubble point method, a physical model was established clarifying the interactions (mass and heat) between the species present in the streams in circulation in the column. In order to identify the externally controlled operating parameters, the degree of freedom of the column was determined by using Gibbs phase rule. The mathematical model converted to Fortran code and based on the principles of: 1) Global and local mass conservation balance, 2) Enthalpy balance, and 3) Vapour-liquid equilibrium at each tray, was used to simulate the behavior of the column, concentration distributions, temperature and streams for each phase along the column at high pressure in each tray. The energy consumption at the condenser and the boiler was also evaluated using the Starling equation of state.

Keywords: High pressure distillation, complex mixture, modeling, simulation, bubble point method.

1. Introduction

Numerous methods are cited in the literature regarding the resolution of the complex mixtures distillation problem. However; the engineer often resorts to the shortcut method which enables the computation of the minimum number of trays (Nmin) [1, 2-3] and to the formalism of Underwood to determine the minimum reflux rate Rmin [4, 5-6]. As for the optimal reflux rate, the Malkanov's equation which requires the preliminary definition of the key components (light and heavy) is called

upon [1, 4, 5-6]. For the complex mixtures, other light components are combined with the light key component according to the vapour-liquid equilibrium constants K. In crude distillation, it is clear that the definition of the key component varies from one tray to another according to the values of K. Obviously, this remark cannot be applied for strongly non ideal complex mixtures without introducing a certain number of errors where distillation with a minimum reflux exhibits a pinching point in the two sections of the column [7]. The quantity of minimal energy necessary for a given separation is specified purely by thermodynamic analysis of the process and depends only on the quality of the feed, the composition of the desired products and the pressure in the column. Although rigorous, the shortcut method remains a universal procedure of design presenting sufficient preliminary information to initiate complex mixtures distillation simulation [8]. In exchange, the simulation programs exploited in the industry require a significant number of iterations carried out manually which is reflected by a long computation period of time. Bausa et al [9] proposed a geometrical method of non-ideal separation while calculating the rate of minimum reflux. In the petrochemical field, the related industries require products of specific quality meeting some standards. The complex mixtures distillation does require a series of columns whose arrangement depends on several factors. The synthesis of processes of complex mixtures distillation is very significant for an optimal design. Indeed, due to the complex interactions between components in the mixture, the research has been focused on ternary mixtures [10]. Reversible distillation of ternary mixtures implies a consequent minimization of the potential difference of transfer (mass and heat) at the tray level. This minimization results graphically in the reduction of the difference between the operating line and the equilibrium curve on an x-y diagram.

2. Column Modeling

Schematically, a distillation column is composed of a cascade of trays between which liquid and vapour phases flow in counter-current directions according to hydrodynamic diagrams depending on tray model. These interactions lead to a mass transfer so that the less volatile components are recoverable at the lower trays, whereas the lightest are recovered mainly in the distillate. To take account of physical interactions, we used the rigorous method based on the bubble point calculation without definition of the light and heavier keys in the column. [11] In the absence of a chemical reaction, the most complex tray with ramified interactions can be depicted in Figure 1. Mj and Uj are the liquid and vapour side streams and Fj the feed flow rate for tray J. For a complex mixture made up of C chemical species, the mathematical model related to a tray is written: Mass balance for component i: (2 < i < c):

Lj-1 xXi,j-1 + V/+1X Y,j+1 + Fj xZ,j -(Lj +Mj) x Xt,j -(Vj + Uj) x Y,j =0 (1)

For each component, the vapour-liquid equilibrium is expressed by:

Yi, j - Ki, j X Xi, j = 0 (2)

The condition related to the composition, expressed in molar fraction, is:

X Xi-j -1 = 0

Fj—H

Tray J

-►Mi

Figure 1. Theoretical tray with heat exchange.

The enthalpy balance for a tray J is:

L. ,x HI., + V. x Hv., + F x Hf -(l. +M.) x HI - Vj +Ui) x m - Q. =0 (4)

j-1 j-1 j+1 j+1 J Jj V j J! J w J! J ^J (4)

Hfj is the feed enthalpy at the tray J. Proceeding by elimination of Lj and Yij, the mathematical model for a column of N trays can be reduced to:

Aj x Xi, J-1+Bi, J x Xi, J+Ci, J x Xi, J+1=Di, J

where:

B i, J = -

A J = VJ + X (Fm - Um - Mm )- Vl r^ < J < N

m = 1 J

V j + 1 + X ( Fm - U m - Mm )- V1 + M j + (V j + U j )x K i, j

Ci, j = Vj+1 x Ki, j+1 1 < j < n - 1

Di, j = - Fj x Zi, j 1 < j < N

with m refers to number of the trays. It is noticed that the mass balance based on Equation (5), for the various components, takes the N x C form of linear algebraic equations where the unknowns are the compositions Xi,j. For each component, the Thomas numerical method can be used to solve the matrix form (Equation (10)) in order to determine the composition profiles as function of the tray position. To avoid the redundancy of some operating variables, analysis of the degree of freedom of the column is carried out.

i,1 i1 a B,2

3. Analysis of Column Variance

r ^ D i,1

D i ,2

V i' V

The method of determination of the degree of freedom of the installation is based on the Gibbs phase rule, and on the one hand, and the laws of mass and energy conservation, on the other [12]. If Nv is the number of unknown variables and Ne that of independent equations, the number of degrees of freedom ND (variance) is:

Nd = Nv - Ne

For the theoretical tray depicted in Figure 1, its interaction with seven mass streams and one thermal stream leads to:

Nv = 7 X (C + 3) + 1

Ne = 4 X (C + 3)

which results in the following number of degrees of freedom:

Nd = 3 x C + 10

As a consequence, it is important to specify the list of these design variables, in particular, the following:

Specification Number

- Total liquid flow rate (L) 1

- Molar fraction of component (Xij) C-1

- Total vapour flow rate (V) 1

- Molar fraction of component (Yi,j) C-1

- Total feed flow rate (F) 1

- Molar fraction of component of F (Zij) C-1

- Temperature and pressure of L, V and F 2 x 3 = 6

- Tray temperature (of each stream) (Tj ) 1

- Tray pressure (of each stream) (Pj ) 1

- Liquid side stream (M) 1

- Vapour side stream (U) _1_

3 X C +10

The presence of additional elements essential to column operation (heat exchangers, reflux system, etc.) on the one hand and the redundancy of the liquid and vapour streams interacting between two successive trays should also be taken into account. The study of complex units has led to the following relations [6-13]:

(Nv )unit = XNV )e - Nr(C + 3)+ Na (15)

all elements

(NE Lit = X (NE )e- NR

all elements (16)

where Nr and NA indicate the constraints of redundant molar fractions and the additional variables respectively. Application to a distillation column (Figure r) gives:

- the total variables number of the unit as:

(n v)unit = 16 x n + 5 x n x c + 2 x c + 7 (17)

- the number of independent relations as:

(Ne )umt = 4 x N x C + 10 x N + 2 (18)

Using the equation (11), the number of the degrees of freedom of the column can be expressed by:

(nd )unit =(N + 2)x(C + 2)+ 4 x N + 1 (19)

The external variables can be easily specified:

Specification Number

- Feed streams rates (F) N x (C + 2)

- Vapour streams (V) (C+2)

- Liquid streams (L) (C+2)

- Pressure at the various trays (P) N

- Temperature at the various trays (T) N

- Liquid side stream (M) N

- Vapour side stream (U) N

- Total number of trays _1

(N + 2) x (C + 2) + 4 x N +1

Figure 2. Diagram of the distillation column with liquid and vapour streams.

4. Simulation of the column

After identifying the operating parameters, simulation of the column operation consists in solving the equations of mass and heat balances by using an iterative calculation in order to determine the:

1. Composition distribution of vapor and liquid phases in the whole column.

2. Temperature of each tray.

3. Flow rates of the phases in circulation between trays.

The calculation steps are illustrated by the following algorithm depicted in Figure 3: • Initiate values for Tj(0) and Vj (0) so that T0(1) corresponds to the dew point temperature and T0(N) to the bubble point temperature. For intermediate trays, the temperatures of initialization are supposed to vary linearly as:

7 »(■) -7 "(1) + (j-1)x ^'(N>-f'(1)

Evaluate the equilibrium coefficients Kij (ideal mixture)

□ Calculate all the compositions by solving the system of equations in the form of C tridiagonal

matrix (Equation (10)) by using the Thomas's method. [14] • If all the equations of balances converge with Xij =1 for any tray, then it is imperative to standardize the values of Xij by taking:

X i, j ■

Figure 3. Design procedure of a column distillation.

Calculate the bubble point Tj(k) at each tray by using the Newton-Raphson's method [14]. • Calculate the vapour concentrations:

yi, j = ki, j x x i, j

Calculate the vapour flow rate values Vj (from the heat balance):

V( j) = (CH ( J -1) - AH (J -1) X V ( J -1)) / BH ( j - 1) (23)

with: AH (J-1)=H/ (J-2)-Hv (J-1) (23.a)

BH (j-1)=Hv (J) - HZ (/-1)

(23.b)

CH( J-1) =

YJF (m)- U(m)-MU(m)-V1) m=1

X HZ (J-1)

+ HZ (j-2)+F (J-1)x(HZ (J-1)-Hf (j-1))

+ U ( J - 1) x (Hv ( J - 1) - HZ ( J - 1))+ g ( j - 1) (23.c)

where Hlj and Hvj are respectively the enthalpies of liquid and vapour phases for tray j [15].

• Check if the calculated values Tj(k) and Vj(k) are of the same order of magnitude as the values Tj(k-1) and Vj(k-1) used for the preceding iteration. For this purpose, the following criterion of convergence is adopted:

n ft.k tk - n f^/k .rk - 1 \ 2

5 (+ ( s10 " 7 x N

J = 1 v tj J j = 3 v vj J (24)

• Otherwise, calculations are redone starting from step (2) based on values obtained in iteration K.

Subroutines are needed for computing the following values:

i. Vapour and liquid enthalpies according to the temperature and the pressure whose expression is given in appendix.

ii. Density pL or p are given by the Starling's equation [16]. It is interesting to notice that only the lowest value and the greatest value of p have a physical significance, which

correspond to p and pj respectively.

iii. Equilibrium coefficients K at the tray temperature and pressure.

5. Operating mode

In order to test the model, we validated it by comparing the results of simulation with those obtained experimentally. For this, we worked on an installation made up primarily of a perforated tray column. The column diameter is a function of the cross section:

1. Enriching section: Diam=4.1m

2. Stripping section: Diam=5.5m

Regarding energy consumption, the maximum power provided to the condenser and to the reboiler is equal to 25.23 106 kcal/hr. The column, whose efficiency is 60%, is equipped with 55 real trays. Consequently, the number of corresponding theoretical trays is 33 with a feed at the 14th tray. The unit is fed by a load of LPG whose composition is given by Table 1.

Table 1. Average composition of LPG.

Component Ci C2 C3 Iso-C4 n-C4 Iso-C5 n-C5 C6+

Composition ( % mol.) 0.26 1.97 60.06 14.22 23.39 0.08 0.02 Traces

6. Results and Discussion

6.1. Validation of the column simulation model

In order to analyze the impact of the various operating parameters on the column performance, the model was merely validated as a first step. The column is designed to separate the load made up of C1 / C2 / C3/iso-C4/n-C4/iso-C5/n-C5 into a distillate rich in C3 and a residue fairly rich in n-C4 and iso-C4. Compared with the parameters provided by the ChemShare simulator controlling actually the column, the results obtained show an extremely low discrepancy (<2%). Consequently, the model presented describes perfectly the distillation column at high pressure in order to determine the operating variables with a high degree of accuracy. Thus, for a feed F = 240 m /h (2548.36 kmol/hr =100%), the results compared with the actual values are presented in Table 2, a and 2, b.

6.2. Variation of power exchanged with the column

This influence is graphically depicted in Figure 4. Two sections can be clearly distinguished with increasing slopes. This variation is linear since for a mixture of a given composition and rate of reflux, the energy required for vaporization (condensation) is proportional to the feed flow rate. The diameter of the column is also affected in an almost linear way by the feed flow rate. Takamatsu et al [17] analyzed the possibility of extending the method of Mc Cabe - Thiele to the column with an integrated source of energy while showing that this technology is possible without needing the presence of neither the condenser nor the boiler. Bandyopadhyay [18] showed that it is possible to reduce the energy consumption (Qr-Qc) and to locate the feed tray irrespective of the calculation method and by decreasing the reflux rate while exploiting the pinch temperature in the T-H diagram. For a mixture of five components (C7, C8, C9, Ci0, Ci5), the simulation of distillation in a column with 18 trays

(condenser and boiler included) was made with PROII simulator (1994-1995) using the RKS model (Redlich, Kwong and Soave) whereas the ninth tray was the feed tray. The Bandyopadhyay's simulation results gave Tc=140.4 oC and Tr=207.4 oC. [18] Power consumption in the condenser and the boiler was respectively equal to 37.72 and 82.52 MM Btu/h.

140 150

Operating rate

Figure 4. Power exchanged at boiler and condenser versus operating rate.

Table 2a. Column simulation results (F= 2548.36 kmol/hr).

Component Distillate composition ( molar fractions) Residue composition ( molar fractions)

Design Simulation values Deviation Design Simulation values Deviation

Methane 0.0007 0.0041 0.0034 4.34 10-17 0.417 10-19 0.0000

Ethane 0.0147 0.0310 0.0163 0.1698 10-10 0.127 10-10 0.0000

Propane 0.9659 0.9530 0.0129 0.0100 0.0037 0.0063

Iso-butane 0.0160 0.0099 0.0061 0.3647 0.3665 0.0018

N-butane 0.0025 0.0019 0.0006 0.6226 0.6270 0.0044

Iso-pentane 5.95 10-9 8.3 10-9 0.0000 0.0021 0.0021 0.0000

N-pentane 2.36 10-10 2.88 10-10 0.0000 0.0005 0.0005 0.0000

Temperature (°C) 58.07 58.42 0.6% 111.02 110.71 0.28%

Qc (kcal/hr) 15.85 106 16.17 106 2%

Qr (kcal/hr) 16.51 106 16.86 106 2%

Table 2b. Column simulation results ( F=3567.70 kmol/hr).

Qc=2,277 10+7 kcal/hr .

F = 3567.70 kmol/hr

Composition (%):

Ci = 0.26 C2 = 1.97 C3 = 60.06 i-C4 = 14.22 n-C4 = 23.29

i-C5 = 0.08 n-C5 = 0.02

iaaaAaa wvwyw

Distillate:

D =2231.88 kmol/hr Composition (%):

Ci = 0.413

C2= 3.134

C3 = 95.3306 i-C4 = 0.940

n-C4 = 0.183

i-C5 = 0.774 10-6

n-C5 = 0.270 10-7

v\/w\A/\ vwvw\

Qr=2,372 10+7 kcal/hr

N=33 T=110,7 oC

6.3. Temperature prof ile in the column:

f Residue:

R =1335.82 kmol/hr

Composition (%):

C1 = 0.412 10-17 C2 = 0.128 10-7 C3 = 0.386 i-C4 = 36.788 n-C4 = 62.556 i-Cs = 0.216 n-Cs = 0.054

Figure 5 represents the temperature profile versus the tray position for various feed flow rates. The profile exhibits a significant temperature increase near the bottom of the column. Indeed, this result is corroborated by thermal analysis of the column operation as, according to the model, the fluid on the tray is at its bubble temperature and since the low boiling point components remain at the higher trays (rectifying section), the temperature profile naturally takes that kind of shape. As for the operating rate, it can be noted that the three curves are practically identical. The low discrepancies stem from the disturbances and some unstable local regimes.

6.4. Distribution of the liquid flow rate in the column

The variation of the liquid flow rate versus the tray position for various rates of reflux is presented in Figure 6. A harmony in the liquid flow rate profiles is noticed for the various operating conditions studied, i.e., 120, 130 and 140%. Three distinct sections can be easily identified:

Zone I: corresponds to the liquid flow rate which varies slightly in the enriching section between tray 1 and tray 13.

Zone II: refers especially to the feed tray where a sudden increase of the liquid flow rate occurs because the feed is introduced in the form of boiling liquid.

Zone III: The liquid flow rate is fairly large in the stripping section, with a slight variation; but in the last trays, a decrease of the flow rate is noticed due to the high temperature which reigns at the bottom of columns. It should be noted that the shape of the second section would not be so sharp if the feed was a boiling liquid, a vapour or a biphasic mixture.

g 370 e

S 360 t

arpe 350

m £ 340

..........

♦ Operating rate: ■ Operating rate: A Operating rate:

0 10 20 30 40

Tray order in the column

Figure 5. Profile of the temperature along the column.

s e t ta r

w o lf

10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0

♦ Operating rate: 120% ■ Operating rate: 130% A Operating rate: 140%

20 30 40

Tray order in the column

Figure 6. Profile of the liquid flow rates along the column.

I /in«/

6.5. Distribution of the vapour flow rate in the column

This distribution is graphically represented in Figure 7. The shapes are almost flat with a slight convexity upwards for the first trays (< 14). This can be explained by the fact that at the first tray, all vapour is condensed. However, the vapour flow rate, composed essentially of the more volatile components, increases in the trays of the enriching section, then, it decreases at the feed tray. At lower trays, the vapour contains increasing amounts of less volatile components leading to a slight bending of the curve which accounts for the finding obtained.

9000 8000 7000

& 6000 I 5000

£ o q=!

S 3000

2000 1000 0

♦♦♦

♦♦♦♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦♦

♦ Operating rate: 120% ■ Operating rate: 130% 4 Operating rate: 140%

20 25 30 35 Tray order in the column

Figure 7. Profile of the vapour flow rates along the column.

Conclusion

The simulation of the behavior of a complex mixture distillation shows that for a flooding level close to 0.8, increasing the feed rate from 100 to 140% does neither affect significantly the liquid and vapour flow rates distribution nor the temperature profile along the column, but results in an increase in the column size. Indeed, the temperature follows a profile identical to a breakthrough curve without being significantly affected by the processing capacity of the column. The continuous stirred tank reactor model, based on the method of bubble point rather accurately represents the behavior of a theoretical tray. This method of analysis is a powerful tool for the study of a column operation providing the specific characteristics of the feed (composition, viscosity, surface tension and physical state). The distillate and residue analysis confirm the results predicted by the model. For a feed rate of 3567.70 kmol/hr, the distillate recovered is particularly rich in propane (C3) whose boiling temperature is 54.6 oC, whereas the residue is rather rich in i-C4 and n-C4 at a temperature of 110.7 oC. Moreover, for an adiabatic column (except at the ends), the energy exchanged at the condenser and boiler is so

large that it would be necessary to locally analyze the loss of energy. It can be concluded that the interactions between the different parameters are so strong so that a model sensitivity analysis is essential. In order to locate and to minimize the entropy creation, this study was completed by an exergetic analysis of the column, which is the subject of part II.

Nomenclature

C: Number of components of the mixture

D: Distillate flow rate (kmol/hr)

F: Feed flow rate (kmol/hr)

H: Enthalpy of a mass stream (kcal/kmol)

K: Vapour-liquid equilibrium coefficient (--)

L: Liquid stream flow rate (kmol/hr).

M: Liquid side stream flow rate (kmol/hr).

N: Number of trays in the column (--)

P: Operating pressure of the column (atm).

Qc: Calorific power exchanged with the condenser (kcal/hr)

Qr: Calorific power exchanged with the boiler (kcal/hr)

Tref: Reference temperature (273.15 K)

To: Temperature of the ambient conditions (298.15 K)

V: Vapour stream flow rate (kmol/hr)

U: Vapour side stream flow rate (kmol/hr).

X: Molar composition of the liquid phase

Y: Molar composition of the phase vapour

Z: Molar composition of the feed

Indices

f: Relating to the feed I: Relating to the component in the mixture J: Relating to the sequence number of the tray C, c: Condenser (J=1) R, r: Boiler (J=33).

References

1. Cicile, J.C. Distillation-Absorption, Techniques de l'Ingénieur; Ed. Technip: Paris, 1994; Sections J-2610, J-2611, J-2621, J-2622, J-2623.

2. Claudel, B.; Andrieu, J. ; Otterein, M. Bases de Génie Chimique; Techniques et Documentation: Paris, 1977 ; pp 164-178.

3. Walas, S.M. Textbook of Chemical Process Equipment-Selection and Design, 2nd Ed. Butterworths Series in Chemical Engineering: USA, 1988; pp 236-271.

4. Wuithier, P. Raffinage et Génie chimique. 2nd Ed.; Ed. Technip: Paris, 1972 ; pp 132-165.

5. Chittur, K. ChE 448 Chemical Engineering Design, MultiComponent Distillation. K. Chittur: Huntsville, 1998; http://www.che.uah.edu/courseware/che448-Spring1998/multcomp/ (online courseware).

6. Kwauk, M. A system for counting variables in separation processes. A. I. Ch. E. Journal 1956, 2, 240-248.

7. Koehler, J. Poellmann, P.; Blass, E. A review on minimum energy calculations for ideal and nonideal distillations. Ind. Eng. Chem. Res. 1995, 34, 1003-1020.

8. Rong, B.G.; Kraslawski, A.; Nystrom, L. Design and synthesis of multicomponent thermally coupled distillation flow sheets. Comp. Chem. Eng. 2001, 25, 807-820.

9. Bausa, J.; Vonwatzdorf, R.; Marquardt, W. Minimum energy demand for nonideal multicomponent distillations in complex columns. Comput. Chem. Engng. 1996, 20, 55-60.

10. Agrawal, R.; Fidkowski, Z.T. Are thermally coupled distillation columns always thermodynamically more efficient for ternary distillations. Ind. Eng. Chem. Res. 1998, 37, 34443454.

11. Perry, R.H.; Green, D.W.; Maloney, J.O. Perry's Chemical Engineer's Handbook, 7th Ed.; McGraw-Hill: New York, 1997; pp 1-54.

12. Smith, F. Analysis of the Equilibrium Stage Separation Problem Formulation and Convergence. Am. Inst. Chem. Eng. J. 1964, 10, 698-711.

13. Gilliland, E.R.; Reed, C.E. Degrees of freedom in multicomponent absorption and rectification columns. Ind. Eng. Chem. 1942, 34, 551-557.

14. Gourdin, A.; Boumahrat, M. Méthodes Numériques Appliquées; O.P.U: Algiers, 1993; pp 147202.

15. Holland, C.; Pendon, G. Solve more Distillation Problems. Hydrocarbon Process. 1974, 53, 148156.

16. Starling, K.E.; Han, M.S. Thermo data refined for LPG. Pt. 15. Industrial applications. Hydrocarbon Process. 1972, 51, 107-115.

17. Takamatsu, T.; Nakaiwa, M.; Huang, K.; Akiya, T.; Noda, H.; Nakanishi, T.; Aso, K. Simulation Oriented Development of a New Heat Integrated Distillation Column and its Characteristics for Energy Saving. Computers and Chemical Engineering 1997, 21, 243-247.

18. Bandyopadhyay, S.; Malik, R.K.; Shenoy, U.V. Temperature-enthalpy curve for energy targeting of distillation columns. Computers and Chemical Engineering 1998, 22, 1733-1744.

© 2007 by MDPI (http://www.mdpi.org). Reproduction is permitted for noncommercial purposes.