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

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 (N


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 (N min ) [1,[2][3] and to the formalism of Underwood to determine the minimum reflux rate R min [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.

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.M j and U j are the liquid and vapour side streams and F j 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: For each component, the vapour-liquid equilibrium is expressed by: 0 The condition related to the composition, expressed in molar fraction, is: Figure 1.Theoretical tray with heat exchange.
The enthalpy balance for a tray J is: H fj is the feed enthalpy at the tray J. Proceeding by elimination of L j and Y ij , the mathematical model for a column of N trays can be reduced to: (5) where: 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 C N × 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.

Analysis of Column Variance
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 N v is the number of unknown variables and N e that of independent equations, the number of degrees of freedom N D (variance) is: For the theoretical tray depicted in Figure 1, its interaction with seven mass streams and one thermal stream leads to: (12) ( ) C N E (13) which results in the following number of degrees of freedom: (14) As a consequence, it is important to specify the list of these design variables, in particular, the following:  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][7][8][9][10][11][12][13]:

Specification
where N R and N A indicate the constraints of redundant molar fractions and the additional variables respectively.Application to a distillation column (Figure 2) gives: -the total variables number of the unit as: ( ) -the number of independent relations as: ( ) Using the equation ( 11), the number of the degrees of freedom of the column can be expressed by: ( ) ( ) ( ) The external variables can be easily specified:

Specification
Number -Feed streams rates (F) -Vapour streams (V) -Liquid streams (L) -Pressure at the various trays (P) -Temperature at the various trays (T) -Liquid side stream (M) -Vapour side stream (U) -Total number of trays

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 T j (0) and V j (0) so that T 0 (1) corresponds to the dew point temperature and T 0 (N) to the bubble point temperature.For intermediate trays, the temperatures of initialization are supposed to vary linearly as: • Evaluate the equilibrium coefficients K i,j (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 X i,j =1 for any tray, then it is imperative to standardize the values of X i,j by taking: No Yes

Initially
Calculate the bubble point Tj(k) at each tray by using the Newton-Raphson's method [14].
• Calculate the vapour concentrations: • Calculate the vapour flow rate values Vj (from the heat balance): ) ) where H lj and H vj 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 T j (k-1) and V j (k-1) used for the preceding iteration.For this purpose, the following criterion of convergence is adopted: • 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 ρ L or ρ V are given by the Starling's equation [16].It is interesting to notice that only the lowest value and the greatest value of ρ have a physical significance, which correspond to ρ V and ρ L respectively.
iii.Equilibrium coefficients K at the tray temperature and pressure.

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.5mRegarding energy consumption, the maximum power provided to the condenser and to the reboiler is equal to 25.23 10 6 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 14 th tray.The unit is fed by a load of LPG whose composition is given by Table 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 C 1 / C 2 / C 3 /iso-C 4 /n-C 4 /iso-C 5 /n-C 5 into a distillate rich in C 3 a residue fairly rich in n-C 4 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 3 /h (2548.36kmol/hr ≡100%), the results compared with the actual values are presented in Table 2, a and 2, b.

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 (Q r -Q c ) 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 (C 7 , C 8 , C 9 , C 10 , C 15 ), 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 T c =140.4 o C and T r =207.4 o C. [18] Power consumption in the condenser and the boiler was respectively equal to 37.72 and 82.52 MM Btu/h.

Temperature profile in the column:
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.

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.

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.

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 (C 3 ) whose boiling temperature is 54.6 o C, whereas the residue is rather rich in i-C 4 and n-C at a temperature of 110.7 o C.Moreover, for an adiabatic column (except at the ends), the energy exchanged at the condenser and boiler is so Number -Total liquid flow rate (L) -Molar fraction of component (Xi,j) -Total vapour flow rate (V) -Molar fraction of component (Yi,j) -Total feed flow rate (F) -Molar fraction of component of F (Zi,j) -Temperature and pressure of L, V and F -Tray temperature (of each stream) (Tj ) -Tray pressure (of each stream) (Pj ) -Liquid side stream (M) -Vapour side stream

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

Figure 3 .
Figure 3. Design procedure of a column distillation.

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

Figure 5 .Figure 6 .
Figure 5. Profile of the temperature along the column.

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

Table 1 .
Average composition of LPG.