Computers and Chemical Engineering 24 (2000) 1843 – 1849 www.elsevier.com/locate/compchemeng
A new approach to simulation of distillation in packed columns Lida Yang, K.T. Chuang * Department of Chemical & Materials Engineering, Uni6ersity of Alberta, Edmonton, Alta, Canada T6G 2G6 Received 11 July 1997; received in revised form 10 April 2000; accepted 10 April 2000
Abstract A modified equilibrium stage model for simulating and designing packed distillation columns has been developed. The new model incorporates mass transfer efficiency and convective heat transfer along the packed bed. This results in a set of simultaneous (MESHEC) equations to replace the classic MESH equations. By solving the MESHEC equations, detailed temperature, flowrate and composition profiles in the column can be determined. Furthermore, the new model is based on the well established methods of calculation for the equilibrium stage and mass transfer efficiency. Thus, it is more practical and efficient than the nonequilibrium model, which requires the precise input of mass transfer coefficients. The model presented in this paper can also be extended to simulate the plate columns and even absorption and extraction operations. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Distillation design and simulation; Equilibrium stage model; Nonequilibrium model; Efficiency
1. Introduction Packed columns are widely used in the chemical and petroleum industry for distillation operations. Shortcut and rigorous procedures are the two calculation methods used for the design and simulation of such operations. There are also two fundamentally different kinds of rigorous models to describe such operations (i.e. the equilibrium stage model and the nonequilibrium model) (Kister, 1992). Conventional equilibrium stage models are based on the concept which assumes the streams leaving any particular stage are in equilibrium with each other. Component Material balances, phase Equilibrium, Summation equations and Heat balances for each stage (the so-called MESH equations) are solved simultaneously to give product distributions, temperature profiles, and other operating parameters (Holland, 1975). A shortcoming of the equilibrium stage model is that, during operation, stages rarely operate at equilibrium, despite attempts to approach this condition by proper * Corresponding author. Tel.: +1-780-4923321; fax: + 1-7804922881. E-mail address:
[email protected] (K.T. Chuang).
design and choice of operating conditions (Krishnamurthy & Taylor, 1985). The equilibrium stage model assumes that the contacting time of the streams with each other is infinite in each stage. The product distributions and temperature profiles calculated using the equilibrium stage model do not correspond to those of any real stages and, therefore, cannot describe the actual operating conditions in the packed column. Since the late 1970s, several nonequilibrium models (Krishnamurthy & Taylor, 1985; Sivasubramanian, Krishnamurthy & Taylor, 1987) have been proposed to overcome this shortcoming of the equilibrium stage model. All of these models abandoned the assumption that each stage operates at equilibrium. Instead, they introduced mass and heat transfer coefficients to describe the stage operation. These nonequilibrium models offer a reasonable approach to distillation calculations. However, techniques to reliably predict the mass transfer coefficients, interfacial areas and diffusion coefficients are unavailable, and thus their application is restricted to academic exercises (Kister, 1992). This is the reason that the equilibrium stage model is still widely used in the design and simulation of distillation processes, despite the criticisms. In our opinion, the equilibrium stage model is so well established in the commercial computer software packages that it is more
0098-1354/00/$ - see front matter © 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 0 ) 0 0 5 6 2 - 7
1844
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
pragmatic to modify the model by incorporating the concept of mass transfer efficiency, rather than to establish a totally new design procedure. The advantage of this approach is that the required input data for the stage efficiency or height equivalent to a theoretical plate (HETP) are readily available for a large number of separation systems and equipment. The objective of our research was to develop a reliable design and simulation model for distillation calculations. Our approach was to modify the traditional MESH equations and incorporate the mass transfer efficiency into the equilibrium stage model. Unlike the conventional simulation methods, which only generate operating conditions at discrete theoretical stages, this model will give operating conditions at any bed height of a packed tower. Since the new model can describe the detailed picture of a real process, it provides an efficient method for reliable process and equipment design.
2. Model development The equilibrium stage model assumes that the streams leaving any particular stage are in equilibrium with each other. Although this is not the case in real operations, it is reasonable to assume that one phase leaving a certain stage is in equilibrium with another phase leaving another stage, the location of which is determined by the mass transfer efficiency. In the case of packed columns, one phase leaving a certain crosssection of the column could be assumed to be in equilibrium with another phase leaving another crosssection of the column, the separation between which is an HETP. This is the basic principle for our model development. The model described herein was developed from an investigation into the causes of coking, which frequently occurs in the packed section (wash oil section) at the bottom of a vacuum tower for crude oil distillation (Fig. 1). Since there are thousands of towers of this type, which were designed using the conventional equilibrium stage model, and which have the same problem, it is believed that the cause is to be found in the process design methodology. Further, in this model for a packed distillation column, the impact of effects such as entrainment and conduction heat of the column internal, which are of more importance in dealing with packed columns than plate columns, will be considered. The model described can be easily extended to simulate the plate columns and even absorption and extraction operations. The calculating procedures for our model are summarized as follows: 1. The separation process is simulated using the conventional equilibrium stage model, and product distri-
butions, temperature profiles, etc., corresponding to the equilibrium stages are found. The results calculated using our model are therefore, in principle, consistent with those of the well established equilibrium stage model. The entrainment is usually not included in the conventional calculations; however, for packed columns the entrainment should be considered in most cases. 2. A theoretical stage is equivalent to a portion of the packed tower with a height of a HETP called a sub-bed. Several such successive sub-beds join together to form the whole pack tower. The vapor leaving the top of this sub-bed is considered to be in equilibrium with the liquid leaving the bottom of the sub-bed. A differential element of our model is a cross-sectional element with a height of dH (Fig. 2). Each phase is assumed to be uniformly distributed within the differential element, and the composition of the liquid entrainment is assumed to be identical to the composition of the liquid phase. For the ith differential element of each sub-bed, the following equations can be described: a) Component Material equations: Li Xi, j + Gi Yi, j + Ei Xi, j = Li + 1, j Xi + 1, j + Gi − 1Yi − 1, j + Ei − 1Xi − 1, j
(1)
b) It is assumed that the liquid phase is at its bubble point and the vapor phase is at its dew point. However, it should be noted that the liquid and vapor leaving the differential element are not at the same temperature, although the vapor leaving the top of the sub-bed is at the same temperature as the liquid leaving the bottom of the sub-bed (Fig. 2). The temperature profiles of the vapor and liquid phase are determined by the mass
Fig. 1. Wash oil section at the bottom of a vacuumtower for crude oil distillation.
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
1845
Gi − 1H(TYi − 1)+ Li + 1h(TXi + 1)+ Ei − 1h(TXi − 1) − lAd
TPi − TPi − 1 dH
(6)
e) For a vacuum distillation operation, the F-factor for the vapor flow (defined as uy ry ) is usually small. In such cases, the total mass transfer rate is controlled by the mass transfer in the vapor phase (van Winkle, 1967). Therefore, the efficiency equation for the vapor phase is written as: Yi, j − YN -1, j i · dH = YN, j − YN − 1, j HETP
(7)
If the mass transfer resistance lies in the liquid phase, the Efficiency equation is written as: Xi, j − XN, j i · dH = XN + 1, j − XN, j HETP
Fig. 2. Diagram of a differential element for the MESHEC equations.
transfer rate (i.e. the efficiency-see the following Efficiency equation). This modification is one of the major improvements of our model over other suggestions in the literature, which also incorporate the efficiency into the equilibrium stage model. In those models, although the efficiency is applied to the equilibrium constant Ki,j, and appears as the product of oi,jKi,j (Kister, 1992), the temperatures of the liquid and vapor phases leaving each ‘real’ stage or each bed height of a ‘real’ packed bed are assumed to be the same. This assumption is, however, not generally valid. Therefore, to describe our modification more clearly, the modified phase Equilibria are expressed in our model as: Xi, j =fi (TXi )
(2)
Yi, j =Fi (TYi )
(3)
c) The Summation equations remain as: (4)
j=1 n
% Yi, j =1
(5)
j=1
d) The Heat balance which includes axial conduction heat through the packing under adiabatic conditions, is expressed as: Gi H(TYi )+Li h(TXi ) +Ei h(TXi ) −lAd
However, Eqs. (7) or (8) are valid only when there is no concentration extrema (maxima or minima) in the subbed. The existence of an extrema can be determined by plotting Yj versus TY in the temperature range TN B TY B TN − 1. The Yj − TY relationship can be obtained by solving equilibrium equations (Eqs. (2) and (3)). If there is a concentration extrema which occurs at Ym,j, the sub-bed must be divided into two parts (i.e. one part of the sub-bed, say the first k differential elements, for the concentration to reach Ym,j from YN − 1,j, and the remainder of the sub-bed for the concentration to reach YN,j from Ym,j ). Therefore, Eq. (7) can be revised for the first part of the sub-bed as: i · dH Yi, j − YN − 1, j = Ym, j − YN − 1, j k/s HETP
TPi + 1 −TPi = dH
(9)
and for the remainder of the sub-bed as: (i− k) · dH Yi, j − Ym, j = YN, j − Ym, j (1−k/s)HETP
(10)
(f) The packing in the tower is assumed to be uniformly ‘coated’ with liquid, and only the liquid phase has direct heat exchange with the packing. Therefore, the heat convection equation for the liquid phase is written as: −lAd
n
% Xi, j =1
(8)
TPi − TPi − 1 T − TPi + lAd Pi + 1 dH dH
= bAV(TPi + 1 − TXi )
(11)
The two terms on the left hand side of the above equation represent conduction heat transfer through solid packing. 3. The above set of equations are termed MESHEC equations. By solving these equations simultaneously for each differential element, the actual component distributions, temperature profiles for each phase and other operating parameters, can be determined. It can be seen that both newly proposed efficiency and con-
1846
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
vection equations are linear. The rest of the MESHEC equations are in form nearly identical to those in the classic MESH equations. Therefore, solving of the MESHEC equations is not significantly more difficult than that for the original MESH equations. The existing computational methods, such as the bubble-point methods and the global Newton methods, can well be applied to solve the MESHEC equations. According to our experience, by solving at first the liquid temperature TXi through the dew-point equation (Eq. (7)), or the vapor temperature TYi through the bubble-point equation (Eq. (8)), the total MESHEC equations can be solved without difficulties.
3. Model demonstration and discussion To demonstrate the validity of the model, the ‘wash oil section’ of the vacuum tower at a refinery, shown in Fig. 1, was simulated. This section normally contains grid or structured packing. The liquid from the HVD section, which is directly above the wash oil section, is collected and split into three streams. One stream is drawn as side product heavy vacuum distillate (HVD), one is fed back to the liquid distributor of the HVD section, and one is fed to the liquid distributor of the wash oil section. The inlet liquid flow rate at the top of the wash oil section is adjustable. In our simulation, the conditions of the feed, which is a liquid – vapor mixture, as well as the conditions of the inlet stream at the top of the wash oil section, are specified. The entrainment flow along the packed bed was assumed to have the function of E =E0e − 1.0H (E0 was assumed to account for 10% of the bottom product of the vacuum tower (vacuum residue)). The grid packing used in our simulated tower has an HETP of 1.2 m, and the packed section has a height approximately equal to two theoretical stages. The thermodynamic properties of the system were evaluated using the commercial simulation package PROII (supplied by SimSci, CA), which was based on the Braun K10 model. The Braun K10 model is generally recommended for heavy hydrocarbon systems at low pressure (Watkins, 1979).
be encountered. Since coking (from dry-out) has been observed on many occasions, a design of the packed bed based on the conventional simulation results may be in error.
3.2. Results from the MESHEC model 3.2.1. Without heat conduction (l= 0) In the case where the packing material is a poor heat conductor (e.g. ceramic), the l value would approach zero. Under such conditions, the heat convection equation (Eq. (11)) becomes TPi + 1 = TXi, and the heat balance equation (Eq. (6)) is reduced to a simple enthalpy balance equation. The results from the MESHEC model (Fig. 3) show that the liquid flow rate was as low as 2.5 × 103 kg h − 1 at the bed height of 1.9 m. This portion of the bed may not operate properly with such a low liquid rate. Consequently, ‘coking’ of the heavy oil on the packing seems to be inevitable. The concentration distributions of the distillate fraction of heavy oil with NBP= 440–468°C in the packed bed is shown in Fig. 4. By examining the concentration distribution shown in the figure, a concentration maximum in the bed is observed. Because of the existence of the maximum, Eqs. (8) and (9), instead of Eq. (7), were used to describe the mass transfer efficiency of the equilibrium stage III (Fig. 4). It can be seen that the minimum liquid flow rate occurs at the maximum concentration. As a matter of fact, it is quite common that in non-ideal multicomponent distillation systems the composition distributions along the distillation column show extrema at certain column height(s) (King, 1980; Holland, 1981). However, calculations using the conventional equilibrium stage model can only reveal the operating conditions at the equilibrium stages and, quite naturally, could mislead to such an expectation
3.1. Con6entional method To compare the results obtained by our model with those obtained by the conventional method, the operation of the wash oil section was also simulated using the commercial simulation package PROII. The results are shown in triangle in Fig. 3. The results form the PROII package suggest that the liquid flow rate along the packed bed is always greater than 17× 103 kg h − 1. This value is above the minimum wetting rate for the packing. Therefore, it is expected that the packed bed should receive sufficient liquid and dry spots should not
Fig. 3. Liquid flow rate distributions along the packed bed.
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
1847
using the proposed method, are different from that calculated using the conventional theoretical stage method. In the latter method, the temperatures of the liquid and vapor phase of each calculated element are assumed to be equal. For a two-stage packed bed, there would be only four calculated temperatures (one liquid and one vapor temperature for each theoretical stage), which result in over-simplification of the liquid flow for a packed bed with temperature difference of 130°C. In our model, both the liquid and vapor temperature profiles are continuous, and they are consequently different from the temperatures calculated using the conventional theoretical stage method at the stage boundaries. Therefore, our method showed that different loads at the stage boundaries to the conventional theoretical stage method and our results should be closer to reality. The detailed temperature and flow rate profiles of the packed bed are shown in Table 1. To investigate the sensibility of the differencing parameter (dH) on the solution of the MESHEC equations, we calculated liquid flow rates under the same conditions, but with dH values of 0.05 and 0.02 m. The results are reported in Fig. 5, along with the results of dH=0.01 m. No significant difference was observed for using different dH values.
Fig. 4. Concentration distributions of the distillate fraction of heavy oil with NBP= 440 468°C.
that the operating parameters along the packed bed always change in manners which involve no extrema between equilibrium stages as shown by the dashed line in Figs. 3 and 4. Consequently, this method could lead to an improper design of the equipment. The proposed MESHEC equations can overcome these shortcomings. We can also see from Fig. 3 that the liquid loads at the stage boundaries (H =0.0 and 1.2 m), calculated Table 1 Detailed temperature and flow rate profiles of the packed bed Bed height H m
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4
Liquid phase T °C
Vapor phase T °C
Liquid flow kg h−1
Vapor flow kg h−1
New method
Pro II
New method
Pro II
New method
Pro II
New method
Pro II
356.4 350.7 346.0 340.8 335.6 330.5 325.1 319.1 313.2 307.6 302.0 296.6 291.2 286.8 282.3 278.9 274.4 272.2 269.8 267.8 265.8 262.8 258.9 254.9 250.5
356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 250.5
379.0 377.0 375.0 373.0 371.1 369.1 367.2 365.3 363.4 361.5 359.7 357.8 356.4 350.4 344.9 339.5 334.2 328.9 323.6 315.0 303.7 299.6 296.3 293.4 291.2
379.0 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 356.4 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2 291.2
32.1 32.1 30.7 29.6 28.5 27.6 27.2 27.2 27.3 27.4 27.5 27.6 27.8 24.5 21.0 17.2 13.6 9.7 6.0 2.5 7.3 13.0 20.0 28.4 38.6
18.0
105.1 105.1 105.1 104.6 104.1 103.7 104.1 104.6 104.5 105.0 105.4 105.8 106.2 103.1 99.9 96.3 92.9 89.1 85.5 82.1 87.0 92.8 99.9 108.4 118.6
105.1
33.1
38.6
109.0
119.0
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
1848
Fig. 6. Influences of conduction heat on the liquid flow rate and the outlet liquid composition along the packed bed. Fig. 5. Effect of dH on liquid flow rate.
3.2.2. With heat conduction To determine the degree of the contribution due to conduction heat through the packing on the NESHEC simulation results, parameters relevant to the heat transmission in the packed tower were examined. Thermal conductivity of metallic materials is generally smaller than 100 kcal h − 1·m·°C (Perry & Green, 1984). The equivalent area for heat conduction of a tower with packing of moderate to large nominal sizes should be equal to only a few percent of the total tower cross-section At. This should be true because the packing elements have only ‘point contact’ with each other. For the grid type packing considered in our simulation, it is reasonable to assume that Ad 50.01At. Considering the flow rates in the tower and the possible temperature gradient of the packed bed, it can be seen that the terms in Eq. (6) which account for the contribution of the heat conduction through packing can be negligible. The other heat transfer equation, Eq. (11), can be rearranged as: TPi + 1 −TXi =
C (T +TPi − 1 −2TPi ) dH Pi + 1
(12)
where C is defined as a heat transmission number: C=
lAd bAV
Since the purpose of the investigation was to determine the influence of conduction heat on the liquid flow rate along the packed bed, the flow rate of the inlet liquid (at the top of the packed bed) was specified to be the same as the case where l=0 (instead of the composition of the outlet liquid at the bottom of the bed). In addition, it was assumed that the conduction heat had no influence on the temperature profile (i.e. the composition profile) of the vapor phase. This assumption was justified since the vapor phase had no direct heat exchange with the packing and, more importantly, the conduction heat through the packing was much smaller than the enthalpy of the vapor phase. The results for the case C= 1.0, together with the results for C= 0 (l=0) are shown in Fig. 6. It can be seen that the influence of the conduction heat on the distillation process is insignificant. The evaporation or condensation on the packing caused by the conduction heat is shown to be practically negligible. Therefore, although the wetting and liquid holdup on the packing could be improved by replacing metallic packing with ceramic packing, which may alleviate the ‘coking’ problem, the extremely low liquid flow rate zone will continue to be a concern to the operator. Since the coking problem is mainly caused the low liquid flow rate in the bed, a change in the inlet liquid flow rate at the top of the wash oil section should be considered (Fig. 1).
(13)
The convective heat transfer coefficient for falling film is generally greater than 10 kcal h − 1·m2·K (Chapman, 1967). The equivalent area for heat convection of the differential element dH of the packed tower can be estimated as: AV = At · dH · a where a is the specific surface area of the packing. C values are, in most cases, smaller than 0.1, a conservative C value of 1.0 was used in our simulation.
4. Conclusions The model presented in this paper for simulating packed distillation columns has advantages over conventional methods. By modifying the classic MESH equations of the equilibrium stage model, the MESHEC model can simulate actual operating conditions. The use of the MESHEC equations avoids the possible mistakes caused by obtaining only the operating conditions at each ‘imagined’ theoretical stage, us-
L. Yang, K.T. Chuang / Computers and Chemical Engineering 24 (2000) 1843–1849
ing the conventional equilibrium stage methods. Compared with other proposed modifications which incorporate efficiency into the equilibrium stage model, the MESHEC model describes a distillation process with greater precision. Since the model is based on the well established understanding of equilibrium stage and mass transfer efficiency, the required parameters for calculating a given process are the basic equilibrium and efficiency data. These data are usually available or can be estimated with fair reliability, and any effort to improve the reliability of these basic data will lead to a more precise application of our model. The demonstration in the present study shows clearly that our model offers a new approach to reliable simulation and design of separation processes.
Appendix A. Nomenclature A C
m2 m−1
dH
m
E
kg h−1
G
kg h−1
H
kcal kg−1
H HETP h
m m kcal kg−1
K L
kg h−1
NBP
°C
S
T u X
°C m s−1 kg kg−1
area heat transmission number defined as in Eq. (13) height of the differential element flow rate of entrainment flow rate of the vapor phase enthalpy of the vapor phase packed bed height HETP enthalpy of the liquid phase equilibrium constant flow rate of the liquid phase normal boiling point a theoretical stage is divided into S differential elements temperature velocity mass concentration of the liquid phase
.
1849
Y
kg kg−1
mass concentration of the vapor phase
Greek letters a
m−1
b
kcal h−1·m2·°C
o l r
kcal h−1·m·°C kg m−3
specific surface area of the packing convective heat transfer coefficient efficiency heat conductivity density
Subscripts d i j m N P t v X Y
conduction ith differential element component j extrema (minima or maxima) theoretical stage number packing total convection liquid phase vapor phase
References Chapman, A. J. (1967). Heat transfer. New York: Macmillan. Holland, C. D. (1975). Fundamentals and modeling of separation processes. New Jersey: Prentice-Hall. Holland, C. D. (1981). Fundamentals of multicomponent distillation (p. 358). New York: McGraw-Hill. King, C. J. (1980). Separation processes (pp. 326 – 327). New York: McGraw-Hill. Kister, H. Z. (1992). Distillation design. New York: McGraw-Hill. Krishnamurthy, R., & Taylor, R. (1985). A nonequilibrium stage model of multicomponent separation processes, part I. American Institute of Chemical Engineering Journal, 31, 449 – 456. Perry, R. H., & Green, D. (1984). Perry’s chemical engineers’ handbook (pp. 3 – 261). New York: McGraw-Hill and pp. 3 –262. Sivasubramanian, M. S., Krishnamurthy, R., & Taylor, R. (1987). A nonequilibrium stage model of multicomponent separation processes, part VI. American Institute of Chemical Engineering Journal, 33, 325 – 327. van Winkle, M. (1967). Distillation (p. 545). New York: McGrawHill. Watkins, R. N. (1979). Petroleum Refinery Distillation (second ed.). Houston, TX: Gulf.