Simultaneous Momentum, Mass, and Energy Transfer Analysis of a

Ind. Eng. Chem. Res. 2010, 49, 6599–6611

6599

Simultaneous Momentum, Mass, and Energy Transfer Analysis of a Distillation Sieve Tray Using CFD Techniques: Prediction of Efficiencies Dirceu Noriler,*,†,‡ Antonio A. C. Barros,† Maria R. Wolf Maciel,‡ and Henry F. Meier*,† Laboratory of Computational Fluid Dynamics, Chemical Engineering Department, Regional UniVersity of Blumenau, 89030-000, Blumenau, Brazil, and Laboratory of Separation Processes DeVelopment, Chemical Engineering School, State UniVersity of Campinas, 13083-970, Campinas, Brazil

Conventional models for distillation columns are based on equilibrium and nonequilibrium stage concepts. Although equilibrium and nonequilibrium stage models provide useful results, they neglect the fluid dynamics phenomena by assuming a perfect mixture on the plates in each phase. However, the flow pattern on a distillation tray is of great importance in terms of the mass and energy transfer mechanisms, and this influence can only be analyzed by carrying out a fluid dynamics study. The main objective of this study is to obtain experimental data of clear liquid height in a distillation sieve tray apparatus operating with air and water, and to apply a computational fluid dynamics (CFD) model in a Eulerian-Eulerian framework for gas-liquid flows that is able to predict the momentum, mass, and thermal phenomena of multiphase flows. A two-phase, three-dimensional and transient model with chemical species, energy, and momentum conservation balances was applied to predict the volume fraction, velocity, pressure, temperature, and concentration fields of twophase flows on distillation sieve trays. The mathematical model was applied in the CFD commercial code for numerical studies, with the construction of a particular numerical grid and with its own subroutines in Fortran language for the closure equations of the model. The hydrodynamic model was compared with experimental data and presented good agreement. The results of the model show the volume fraction, velocity, temperature, and concentration profiles as a function of time and the position on the distillation sieve tray. The model implemented in this study allows direct application to predict efficiencies in distillation plates, more specifically point, plate, component, and global efficiencies. Introduction A better understanding of the mechanisms that occur in largescale industrial processes is important in order to improve equipment, design, and process development. Distillation columns are one of the most important separation techniques, and equilibrium and nonequilibrium stage models are available in the open literature to study and predict the behavior of distillation columns.1-8 However, despite the useful results obtained with these models, they assume a perfect mixture of the phases on the plates and it has been recognized that the flow pattern on a distillation tray greatly affects the mass and energy transfer mechanisms,9,10 and this influence can only be analyzed through a fluid dynamics study. Recent advances and interest in the use of computational fluid dynamics (CFD) techniques have allowed the study of fluid dynamics in processes and equipment.11-19 Contributions have been made toward the modeling and simulation of gas-solid flows12-16 and gas-liquid flows,17-39 for example, in ozonation towers.17 Bubble columns were initially considered for modeling bubbling flows,19-30 where gas bubbles cross the liquid column, introducing, at the same time, advances for modeling gas-liquid flow on a distillation tray. Some studies on gas-liquid flow on a distillation tray have used CFD techniques.31-38 Mehta et al.31 analyzed the liquid phase flow patterns by solving the time-averaging equations of mass and momentum continuity only for the liquid phase. Liu et al.32 attempted to model two-phase flow considering that the liquid phase is the most important, and modeled the gas action * To whom correspondence should be addressed. E-mail: noriler@ furb.br (D.N.); [email protected] (H.F.M.). † Regional University of Blumenau. ‡ State University of Campinas.

with an empirical equation. Krishna et al.33 and van Baten and Krishna34 suggested a three-dimensional and multiphase model to represent the hydrodynamics on a distillation sieve tray. The authors use the Eulerian-Eulerian framework and the Reynolds averaged Navier-Stokes (RANS) equations to model turbulent gas-liquid flows assuming that momentum exchange considers only the bubble-liquid interaction; i.e., the bubble-bubble interaction was not considered. Soares et al.35 modeled the twophase flow behavior considering that the drag force between phases is very high for distillation cases. This implies that the gas and liquid velocities are the same. The authors showed a comparison between heterogeneous34 and homogeneous models for free surface flow and bubble columns. Gesit et al.36 applied the heterogeneous model to a commercial-scale sieve tray and compared the results with experimental data obtained by Solari and Bell.10 However, energy and mass conservation are neglected in most studies. Component mass conservation is applied in bubbling flow by some authors.17,20,30,37,38 Wang et al.37 and Sun et al.38 developed a new model for turbulent mass diffusivity to represent the turbulence effects on the mass equation. The model is an improvement of the model proposed by Liu et al.32 that considers the effects of the gas phase under liquid phase using empirical equations. The concentration profiles were corroborated by the authors when compared with experimental results. The energy equation was made available by Noriler et al.,39 who used a two-phase model with empirical equations for momentum and energy interphase transfer. The momentum transfer was correlated using the equation proposed by Krishna et al.,33 and the energy interphase transfer was correlated using the Nusselt number for a sphere particle. The results showed that the temperature fields have an unsteadystate behavior, as also occurs with the volume fraction and velocity fields. The authors demonstrated that the liquid tem-

10.1021/ie9013925  2010 American Chemical Society Published on Web 06/21/2010

6600

Ind. Eng. Chem. Res., Vol. 49, No. 14, 2010

perature decreases from the inlet up to the weir and the gas temperature increases from the holes up to the top of the tray. Rahimi et al.40 used the model that considers energy and chemical species conservation to predict the temperature and concentration distributions on a distillation sieve tray. In this study, a comparison was made between experimental data obtained in a cold system and the CFD model in an Eulerian-Eulerian framework for gas-liquid flows, able to predict the main phenomenological aspects of multiphase flows. In association, two-phase, three-dimensional and transient models, with mass continuity, energy, momentum, and chemical species conservation, were applied to predict the volume fraction, velocity, pressure, temperature, and composition fields of the two-phase flow on distillation sieve trays. The model was implemented in CFD commercial code, which uses the finite-volume method with variables located in a generalized coordinate system, for operating conditions taken from the nonequilibrium stage model simulation for an ethanol-water mixture with the DESTIL code.5

+ Fkj is the momentum transfer between phases, M+ kj vj - Mkj vk represents the momentun transfer between the phases due to the mass transfer between the phases, and g is the gravity vector. Energy Equation.

∂ (f F h ) + ∇ · (fkFkvkhk + fkqeff k ) ) ∂t k k k

∑

∂ (f F v ) + ∇ · (fkFkvk) ) ∂t k k k

nf

∑M

(1)

kj

k*j

k ) 1, ..., np In eq 1 fk, Fk, and vk represent the volume fraction, macroscopic density, and velocity vector for the kth phase, respectively, Mkj represents the changing mass rate between the phases, and “np” is the number of phases. The gas and liquid volume fractions are related through the summation constraint: np

∑f

k

)1

k)1

where fk )

Qk Qk +

∑ Sl Q j

j ) 1, ..., np and j * k

;

j

j*k

Slj is the velocity relation, i.e. Slj )

|vj | |vk |

Momentum Equations. ∂ (f F v ) + ∇ · (fkFkvkvk + fkTeff k ) ) -fk∇pk + ∂t k k k

np

∑F

kj

+

k*j

nf

Fkfkg +

∑ (M

+ kj vj

- M+ kj vk) (2)

k*j

k ) 1, ..., np pk is the pressure field, having the same value for the gas and liquid phases, i.e., pk ) pj. Teff k represents the effective tensor,

Qkj +

k*j

∑ (M

+ kj hjs

k*j

-

M+ kj hks) (3)

k ) 1, ..., np hk represents the enthalpy of the kth phase, and qeff k represents an effective heat flux for the kth phase composed of the diffusion flux plus an additional flux due to the velocity and enthalpy fluctuations. Qkj represents the energy transfer between phases + and M+ kj hjs - Mkj hks represent the energy transfer between the phases due the mass transfer between the phases. The subscript “s” represents the value of interface. Chemical Species Equation.

Mathematical Modeling The model considers the gas and liquid flows in a EulerianEulerian framework, where each phase is treated as an interpenetrating continuum having separate transport equations. Therefore, the Reynolds averaged Navier-Stokes equations (RANS) for the kth phase can be written as shown in eq 1. Mass Continuity Equations.

nf

np

∂ (f F y ) + ∇ · (fkFkvkyi,k + fkJeff ik ) ) ∂t k k i,k

np

∑Y

i,kj

+

k*j nf

∑ (M

+ kj yi,j

- M+ kj yi,k) (4)

k*j

k ) 1, ..., np;

i ) 1, ..., nc

yi,k and Jieff k represent the mass fraction and effective mass flux for component i in the kth phase, respectively. Yi,kj represents + the chemical species transfer between phases, M+ kj hjs - Mkj hks represents the chemical species transfer between the phases due to the mass transfer between the phases, and “nc” is the number of components. Model Assumptions. The following assumptions are imposed on the model developed for a distillation sieve tray: • A two-phase system is considered. • The gas phase is the dispersed phase. • The binary mixture ethanol-water was used. • The net rate of the mass transfer between the phases is considered to be so small that it can be neglected. In this case, the right-hand side of eq 1 vanishes and the terms associated with it in eqs 2, 3, and 4 are nulls. • No turbulence model was used for dispersed phases. Closure Equations. In order to solve eqs 1-4, additional equations are required relating the effective fluxes of momentum, eff eff Teff k , energy, qk , and chemical species equations, Jik , and the np interphase momentum transfer, ∑k*j Fkj, the interphase energy np transfer, ∑k*j Qkj, and the interphase chemical species transfer, np ∑k*j Yi,kj. Effective Fluxes. The effective fluxes are related to the diffusion and the dispersion of the properties due to the turbulence effects. Thus, the effective fluxes for momentum, energy, and chemical species equations can be written as a combination between diffusion flux and turbulent flux [For the sake of convenience, the diffusion flux is represented without the superscript letter and the turbulent flux with the “t” superscript.]: Teff ) Tk + Ttk k qeff ) qk + qtk k t Jeff ) Jik + Jik ik

(5)


Considering a Newtonian fluid and the eddy viscosity hypothesis, the effective fluxes can be written as T Teff ) -µeff k k [∇vk + (∇vk) ]

qeff ) -λeff k k ∇Tk

(6)

eff Jeff ) -Dijk ∇yik ik

encountered in the prediction of gas-liquid bubbling flow. This strategy was incorporated in the model by implementation of the routines written in the Fortran language as subroutines for the CFD commercial code used in this study. Interphase Energy Transfer. The rate at which energy crosses the interface between phases for gas-liquid two-phase flow can be described as

where µkeff, λkeff, and Dijeffk are the effective viscosity, thermal conductivity, and diffusivity for species i f j in the kth phase, respectively, expressed by

) λk +

λtk

t ) Dijk + Dijk

np

kj

) Fgl ) -Flg )

kj

) Qgl ) -Qlg ) hglAg(Tg - Tl)

βglApg(vg

- vl)

(8)

k*j

where subscripts “g” and “l” represent the properties in the gas and liquid phases, respectively. βgl is the momentum transfer coefficient, and Apg is the projected bubble area in the flux direction per unit volume. For spherical bubbles

Ag )

[( )(

(Ug√Fg)20hcl√Fg /Fl 43 σgl Ahb0.3 hf

hgl )

(14)

0 e Reg e 200 and 0 e Prg e 250 with particle Reynolds and Prandtl numbers, Reg and Prg, i.e. Fl |vg - vl |dg µl

Prg )

Cp,gµg λg

where Cp,g is the heat capacity of the gas phase. Interphase Chemical Species Transfer. Analogous to the momentum and energy interphase transfers, the rate at which the chemical species mass crosses the interface between phases for gas-liquid two-phase flow can be described as np

(9) 2

where Ug is the superficial velocity which can be calculated from Ug ) Qg/Ab. For the gas average holdup, jfg, the Bennett et al.7 correlation was used:

[ ( )] Fg Fl - Fg

∑Y

i,kj

) Yi,gl ) -Yi,lg ) KcglAg(KcFi,g - Fi,l)

(15)

k*j

fg

fg ) 1 - exp -12.55 Ug

(13)

for

3fg Apg ) 2dg

| |

λgNugl dg

Nugl ) 2.0 + 0.6Reg0.5Prg0.3

and

1 4 Fl - Fg gdg Ug 3 Fl

(12)

where λg is the gas thermal conductivity and Nugl is the Nusselt number. The Nusselt number for spherical bubbles in a turbulent regime can be calculated using the Ranz-Marhan41 equation:

Reg )

where dg is the bubble diameter and CD is the drag coefficient. The drag coefficient, CD, was estimated using the drag correlation of Krishna et al.,33 proposed for the rising of a swarm of large bubbles in the churn turbulent regime:

)]

where Ahb is the perforated fractional area, hcl is the clear liquid height, σgl is the surface tension, and hf is the froth height. The clear liquid height and froth height were determined by reconstitution of the volume fraction fields in each time step. The global heat transfer coefficient is defined as

1 βgl ) FlCD |vg - vl | 2

CD )

(11)

k*j

(7)

where µk, λk, and Dijk are the molecular viscosity, thermal conductivity, and diffusion for species i f j in the kth phase, respectively. The µtk, λtk, and Dijt k are the additional viscosity, thermal conductivity, and diffusion for species i f j in the kth phase, respectively, due to the turbulence effects. The former are physical properties of the kth phase, and the latter are the dispersion coefficients that need to be predicted using a turbulence closure or turbulence model. Interphase Momentum Transfer. For interphase momentum transfer, it was considered that momentum transfer is only due to the drag force. The drag force per unit volume for gas-liquid two-phase flow can be written as

∑F

np

∑Q

where hgl is the global heat transfer coefficient and Ag is the interfacial area per unit volume, which can be calculated through an empirical correlation for bubbling on a sieve tray, i.e.

µeff ) µk + µtk k λeff k Dijeffk

6601

0.91

(10)

Note that when eq 10 is substituted in eq 8 the bubble diameter is canceled, eliminating, thus, the main problem

where Kcgl is the global mass transfer coefficient, Ag is the interfacial area per unit volume, Kc is the thermodynamic constant, and Fi,gand Fi,l are the bulk mass concentrations in the gas and liquid phases, respectively. The global mass transfer can be calculated using the empirical equations presented by Zuiderweg2 and specially developed for a sieve tray: Kcgl )

2.6 × 10-5 µ0.25

(16)

6602


The thermodynamic constant can be calculated employing the γ-φ procedure, i.e. Kc )

γipvap i Φip

(17)

where γi is the activity coefficient in the liquid phase, pivap is the component i vapor pressure, Φi is the fugacity coefficient of the gas phase, and p is the total pressure. The activity coefficients are dependent on the composition of the liquid phase and they are calculated, in this study, employing the UNIQUAC activity coefficient model,42 i.e. ln γi ) ln γCi + ln γRi

(18)

where the first term on the right side represents the combinatorial part, which describes the entropic contribution, and the second part is the residual term, due to the chemical species interaction. The combinatorial part is determined through the composition, size, and form of the species and can be written as ln γCi ) ln

()

Φi Φi θi jz + qi ln + li Yi 2 Φi Yi

nc

∑Yl

j j

(19)

j)1

The thermodynamic constant can be converted into units appropriate for eq 15. The equations for the global mass transfer coefficient, interfacial area, vapor pressure, and UNIQUAC model were incorporated into the model by routines written in the Fortran language. Turbulence Equations. The standard k-ε turbulence model can be used to represent the turbulence on a distillation sieve tray. For the standard k-ε model the turbulent viscosity is calculated as follows: µtk ) CµFk

)

φi )

jz (r - qj) - (rj - 1) 2 j riYi

(24)

where kk is the turbulent kinetic energy in the kth phase and εk is the dissipation rate of turbulent kinetic energy in the kth phase. The conservation equation for turbulent kinetic energy and its dissipation rate can be written as

{[ {[

( (

) ]} ) ]}

µtk ∂ (fkFkkk) + ∇ · fk Fkvkkl - µk + ∇kk ∂t σk µtk ∂ (fkFkεk) + ∇ · fk Fkvkεk - µk + ∇εl ∂t σε

with lj

kk2 εk

) fk(Pk - Fkεk) εk ) fk (C1Pk kk C2Fkεk) (25)

k ) 1, ..., np with

nc

∑rY

Pk ) µtk∇vk:[∇vk + (∇vk)T]

qiYi

where Cµ, C1, C2, σk, and σε are the model constants. The turbulent viscosity, µtk, predicted through the k-ε turbulence model and with the concept of the turbulent Prandtl number, σh, and the turbulent Schmidt number, σm, it is possible to predict λtk and Dijt k with the following equations:

j j

j)1

θi )

nc

∑qY

j j

j)1

(20)

where Yi is the mole fraction of the ith component in the liquid phase; Φi and θi represent the area fraction and the segment fraction, which is similar to the volume fraction, respectively; jz is the coordination number defined as 10; ri and qi are obtained from the van der Waals group volume and surface areas of the ith component, respectively. The residual part can be written as

[

nc

ln γRi ) qi 1 - ln(

∑

nc

θjτji) -

j)1

∑ j)1

θjτij nc

∑θ τ

k kj

k)1

]

(21)

with

[

τij ) exp -

uij - uii RTl

]

(22)

where uij and uii are the model parameters and must be determined using experimental data. The Antoine model, eq 23, was used for the vapor pressure calculation. ln pvap ) Ai i

Bi Tl + Ci

(23)

where Ai, Bi, and Ci are the model parameters and Tl is the temperature of the liquid phase. The gas phase was considered ideal; i.e., Φig ) 1.

λtk )

Dijt k

µtk σh

µtk ) σm

In this study, no turbulence models were used for the vapor phase. This approach was used by Krishna et al.,33 Rahimi et al.,40Gesit et al.,36 Noriler et al.,39 and others. It can be shown that turbulent kinetic energy associated with the vapor phase represents 0.1% of the turbulent kinetic energy. Therefore, there are no conservation equations for k and ε in this phase. Boundary and Initial Conditions. Due to the elliptical characteristics of the partial differential equations of the model, boundary conditions on all frontiers of the physical domain are necessary: at the inlet, a uniform profile of velocities and turbulent properties are imposed; there are no slip conditions on the wall for either phase; pressure conditions at the outlet were also applied for the two phases. Figure 1 shows the physical domain used in the simulations. The liquid phase enters the domain through the downcomer clearance area, identified by “Inlet Liquid” in Figure 1, and vapor enters the domain by square holes located at the bottom of the tray, identified by “Holes” in Figure 1. The vapor and


6603

Figure 1. Flow geometry and boundary conditions.

liquid phases leave the domain by the downcomer clearance area and by the top of the tray, identified by “Outlet” in Figure 1. At the beginning of the simulation, the procedure consists of filling with liquid up to the weir height, and with gas up to the weir height at a homogeneous temperature equal to T0. The velocity fields and the turbulent properties were also considered as initial conditions to close the model. Details regarding the initial conditions for each simulation will be shown below. Methodology of Solution. Due to the model complexity in three-dimensional space-time, an efficient strategy is required to guarantee the stability of the numerical solution. Therefore, a methodology was formulated to obtain the solution composed of four steps. The first step consists of solving the model without the energy and chemical species equations. This solution was carried out for 15 s of real time. This result was then used as an initial condition for the second step of the solution. In the second step, the energy equation, together with the energy interphase transfer, was incorporated into the model and the solution was carried out for 35 s of real time. From this solution, the next step began, which considered the chemical species equations in the model, but without interphase chemical species transfer. In this step, the solution was carried out for 50 s of real time. In the last step, the interphase chemical species transfer was incorporated into the model and the solution was carried out until

Figure 2. View of the numerical grid.

the system reached the quasi-steady state condition, i.e., about 80 s of real time for all cases analyzed in this study. All of the variables were reconstituted from the statistical analysis in the last period, from 50 to 80 s, according to eq 26, seeking to obtain average variables over the time period:

∫ ∫

t+∆t

φ¯ )

t

t

φ(t) dt

t+∆t

dt

)

1 ∆t

∫

t+∆t

t

φ(t) dt

(26)

This filter is necessary because the system has no steady state. Rather, the behavior is oscillatory, characteristic of a quasi-steady state, where the variables oscillate around a mean value. Noriler43 presents the solution analysis for the partial differential equations of the model. Based on this, the finite volume method was applied on a generalized coordinated and collocated grid to solve the model. A highorder interpolation scheme was used for the hydrodynamics equations and upwind for the turbulence equations. The pressure-velocity coupling was obtained using the SIMPLEC algorithm. The improved Rhie-Chow algorithm was used to calculate the velocity at the cell faces to avoid numerical problems like checkboarding and zigzag. The relaxation factors were not used. The grid was generated with a maximum edge for the finite volumes of 5 mm to ensure that the solution is independent of the grid size. The temporal

6604


Figure 3. Schema and photograph of the experimental apparatus.

discretization was available by first-order approximation. The total cell number, in this study, is around 70 000 cells, which gives a solution independent of the grid. Figure 2 shows a view of the numerical grid used in this work.

where dh is the perforation diameter and σ is the surface tension. The pressure drop due to the passage of the vapor through the perforation can be assumed to be equal to the vapor flow resistance of the dry tray. Thus, this contribution can be correlated through

Experimental Section Materials. The materials used in the experiments were air under ambient conditions, which simulate the vapor phase, and tap water obtained fron the public supply. Sieve Tray Apparatus. The sieve tray apparatus used in this work consisted of a distillation column with three sieve trays made of acrylic for better visualization. Air entered the column from the first tray driven by a fan, and the water recirculated due to action of a pump. Water entered at the top of column, on the third tray, and was collected at the base of column. The central plate was used for the measurements. The sieve tray was 0.350 m in diameter, D, with 0.234 m weir length, W. The liquid entered the plate through a rectangular entrance of 0.234 m × 0.040 m. The weir height, hw, was 0.060 m, and the fractional hole area to bubbling area ratio, Ah/Ab, was 0.0227. The distance between the liquid inlet and the weir was 0.260, with 66 holes in the tray. Further details of the apparatus and experimental procedure can be found in Noriler.43 Figure 3 shows the schema views and an image of the apparatus. Experimental Procedure and Measurements. The total sieve tray pressure drop is attributed to three factors:2,6,7 the liquid inventory on the tray, the passage of vapor through the perforations, and the formation of vapor bubbles, i.e. ht ) hcl + hv + hσ

(27)

where ht is the total sieve tray pressure drop, hcl is the pressure drop due to the liquid height, hv is the pressure drop due to the passage of the vapor through the perforation, and hσ is the pressure drop due to the formation of vapor bubbles. The last term in eq 27, hσ, can be correlated using the Bennett equation:6 hσ )

6σ2/3

(

dh 1.27g Fl Fl - Fg 2/3

)

1/3

(28)

hv )

2 a FgVh cv2 Flg

(29)

where Vh is the vapor velocity in the perforations and a/cv2 is the constant of the tray. The constant a/cv2 must be determined for each sieve tray and can be obtained through experimental runs without liquid. For this, the downcomer regions were closed to air rising only through perforations. In this study, the constant was determined applying a nonlinear regression using the leastsquares method. For this, data on the pressure drop were obtained with a dry tray for a superficial velocity of air in the range 0.200-0.615 m/s. The constant a/cv2 had a value of 0.507 36. For the purpose of measuring the clear liquid height on the tray, the total pressure drop was measured and eq 27 was applied. The pressure drop in both cases, dry tray and normal operation, was measured using a digital differential manometer, Model GTPD-A (Samrello Inc.). The vapor and liquid flow rates were measured with a turbine flow meter model SVTG 2′′ for air and SVTL 1/2′′ for liquid (Contech Inc.). The sensors were connected to a programmable logic controller connected to a computer for data acquisition and control. Results and Discussion Validation of the Model. The validation of the model was carried out simulating the hydrodynamics of the sieve tray presented in the Experimental Section. The operation conditions were the same as those used in the experiments; i.e., the rate of the liquid per weir length, Ql/W, was 2.14 × 10-3 m3/ms and the superficial gas velocity was varied between 0.170 and 0.567 m/s. The simulation results in terms of the volume fraction are shown in Figure 4 together with a photograph of the dispersion on the tray to allow a comparison between the real sieve tray and that predicted by the model. It can be noted in Figure 4 that, at the top of the plate, the liquid volume fraction is close


6605

Figure 4. Snapshots of the volume fraction in a central plane and photographs of the experiments.

to zero and increases to a constant value close to the weir. In the region close to the bottom tray, the average liquid volume fraction then begins to increase. This profile is due to the fact that, close to the bottom of the plate, the bubbles rise in jet form. In this region, the liquid is in the continuous phase. Above this, the liquid and the vapor phases are mixed, forming the froth region. This region extends to above the weir where the froth region disappears and the vapor behaves as a continuous phase. This behavior has been previously described by Zuiderweg2 and Bennett et al.,7 and confirms that the model predicts well the flow on a distillation sieve tray. A qualitative comparison can be carried out using the height of the continuous liquid above the bottom of the tray. It is possible to perceive that in both the simulation and the experimental situation the continuous liquid height decreases with an increase in the superficial velocity of the same order of magnitude. This is better seen by the dashed line in Figure 4. The agreement between these data allows the model to be qualitatively validated. For the quantitative validation, the pressure drop due to liquid height was compared between the experimental and CFD models. In the experimental setup, the total pressure drop and dry tray pressure drop were measured and bubble formation was calculated by eq 28. Then, the experimental liquid height pressure drop was determined by eq 27. The predicted liquid height pressure drop was reconstituted by the volume fraction fields. Note that it is not possible to compare the total pressure drop or the dry tray pressure drop because the model does not consider the passage of the vapor through the perforation and its effects. The superficial gas velocity was varied from 0.170 to 0.454 m/s, and the rate of the liquid per weir length, Ql/W,

Figure 5. Pressure drop as a function of superficial gas velocity at Ql/W ) 2.14 × 10-3 m3/ms.

was kept at 2.14 × 10-3 m3/ms. Figure 5 shows the pressure drop for the experiments at Ql/W ) 2.14 × 10-3 m3/ms and the simulation as a function of the superficial gas velocity. There was a good agreement between the simulation and experimental data, as can be seen in Figure 5. In both cases, the clear liquid height decreased with an increase in the superficial gas velocity. Through the figure it is possible to calculate the error between the clear liquid height predicted by the model and that of experimental data. The average error [%

6606


Table 1. Operational Conditions Used in Simulations wall

volume fraction normal velocity [m/s] temperature [°C] ethanol mass fraction turbulent kinetic energy [m2/s2] dissipation rate of turbulent kinetic energy [m2/s3] a

vapor inlet

liquid inlet

vapor

liquid

1.000 10.650 83.750 0.812 -

1.000 0.080 81.650 0.785 1.5(ivR)2 k01.5/(0.3Ls)

impermeable 0 insulate impermeable -

impermeable 0 insulate impermeable 0 dεR/dξ ) 0a

ξ, orthogonal wall direction.

Table 2. Physical Properties of Liquid and Vapor Phases

liquid vapor

viscosity [kg/ms]

density [kg/m3]

thermal conductivity [W/m K]

heat capacity [J/kg K]

diffusivity [m2/s]

molecular weight [kg/kg mol]

1.03 × 10-3 1.01 × 10-5

853.800 1.292

0.203 0.019

3253.200 1843.400

7.14 × 10-9 1.68 × 10-5

33.238 36.857

error ) (ht,simulation - ht,simulation)/ht,experimental; average error ) ∑(% error)/number of data points; mean absolute error ) ∑(|% error)|/number of data points] is -0.6% and the mean absolute error is 8.9%. These values show that the model prediction represents the behavior of the liquid-vapor flow on a distillation sieve tray satisfactorily. Prediction of Tray Efficiency. Based on the physical domain presented in Figure 1, a distillation sieve tray of 0.30 m diameter, D, and 0.180 m weir length, W, was used. The column has a rectangular entrance for the liquid of 0.180 m × 0.015 m. The weir height, hw, is 0.080 m, and the fractional hole area to bubbling area ratio, Ah/Ab, is 0.065. The length between the liquid inlet and the weir is 0.240, with 180 holes in the tray. The operating conditions were obtained through a nonequilibrium stage model simulation for an ethanol-water mixture with the DESTIL code.5 The column simulated using DESTIL had 21 trays with the feed located in tray 11. The feed stream was at 40 °C with an ethanol molar fraction of 0.35. The specifications of this column were a reflux rate of 2.0 and an ethanol molar fraction of 0.88 at the top of the column. From the simulation, tray 12 was defined, and the parameters for the CFD simulation were then obtained. Table 1 gives the operational conditions used in the CFD simulation where liquid (with volume fraction of liquid equal to 1) enters in the liquid inlet section and vapor (vapor volume fraction equals 1) enters through the holes. Table 2 presents the physical properties of the liquid and vapor phases used in the numerical simulation. In agreement with the methodology previously described, the simulation was carried out in four steps. In the last step, three bulk variables were monitored in order to verify the stability of the solution and the steady state regime. These bulk properties, clear liquid height, bulk temperature of the liquid phase, and bulk mass fraction of ethanol in the liquid phase, were calculated applying a volumetric average, that is φ¯ k(t) )

∫ ∫ ∫ φ (x, y, z, t) dx dy dz ∫ ∫ ∫ dx dy dz z

y

x

z

k

y

values are around a defined mean value. It is possible to verify that energy and chemical species transfers have little influence on the clear liquid height, which reaches the quasi-steady state around 56 s of real time. The energy and chemical species balance needs more time to stabilize in a quasi-steady state,

Figure 6. Clear liquid height as a function of real time.

(30)

x

Figures 6, 7, and 8 show the clear liquid height, bulk temperature, and bulk mass fraction of ethanol, respectively, as a function of real time as a progress variable of the quasisteady state regime. It can be seen in Figures 6-8 that the flow does not present a well-defined steady state, but rather a quasi-steady state, as has been verified in previous studies,33-36 where the property

Figure 7. Liquid bulk temperature as a function of real time.


φ¯ k(x, t) )

∫ ∫ φ (x, y, z, t) dy dz ∫ ∫ dy dz z

y

k

z

Figure 8. Bulk mass fraction of ethanol in the liquid phase as a function of real time.

Figure 9. Average volume fraction of the liquid phase as a function of dispersion height.

that is, 68 s of real time to obtain a bulk temperature of the liquid phase and 78 s for the bulk mass fraction of ethanol in the liquid phase, presenting a regime similar to that of the clear liquid height, characteristic of the quasi-steady state. Therefore, the variables analyzed in this paper were calculated from eqs 26 and 30 in the quasi-steady state regime, observed after 78 s of real time. Figure 9 shows the average profile for the volume fraction of the liquid phase calculated in planes along the dispersion of the liquid-vapor flow in the sieve tray. It can be observed that, at the top of the plate, the average liquid volume fraction is close to zero, increasing to values of around 0.5 near the weir and remaining constant around this value until 0.020 m above the tray, when the average liquid volume fraction begins to increase. This behavior is in agreement with that discussed in the previous section. The average temperature of the liquid phase in the y-z plane calculated along the x coordinate, i.e., in the liquid flow direction, can be seen in Figure 10. The instantaneous average temperature is calculated from eq 31, and then eq 26 is applied in order to determine the average temperature of the liquid phase as a function of the x coordinate.

6607

(31)

y

The profile shows the average temperature of the liquid phase from the inlet of the liquid up to the weir. The dependence of the temperature on the position can be observed, mainly in the liquid flow direction. For the other direction, this dependence does not appear. The profile shows a lower temperature at the inlet of the liquid which increases quickly near the liquid inlet. From this point onward, the rate is constant up to near the weir, where there is no increase in the temperature. This profile is due to the fact that the liquid enters the plate with a temperature lower than that of the vapor phase, and crosses the plate exchanging energy with the vapor that enters the plate with a higher temperature. Figure 11 shows a snapshots of the central plane from the liquid inlet up to the weir of the sieve tray for the liquid volume fraction (Figure 11a), liquid temperature (Figure 11b), ethanol mass fraction in the vapor (Figure 11c), and ethanol mass fraction in the liquid (Figure 11d). In Figure 11a, a behavior similar to that discussed in relation to Figure 9 can be seen, with a higher liquid volume fraction region near the bottom, a froth region, and a region of continuous vapor. Figure 11b shows that the liquid temperature field is not influenced by the volume fraction field, with higher temperature gradients in the liquid flow direction. This is due to the high global heat transfer coefficient, calculated using eq 13, that results in a high interphase energy transfer rate and, consequently, in a small influence of the volume fraction. In Figure 11c and 11d, the ethanol mass fraction in the vapor and liquid phases, respectively, can be observed. The white regions in these figures represent an absence of the respective phase. This can be noted through a comparison with Figure 11a. The vapor mass fraction of ethanol increases in the vapor flow direction, from the base of the tray up to the top of the tray, while the liquid ethanol mass fraction decreases in the liquid flow direction, from the liquid inlet up to the weir. Moreover, there is a small variation in the vapor ethanol mass fraction at the top of the tray in the flow direction. This indicates gradients to the tray, in both phases, which are important in the performance of the tray. These results can be used to develop equations for the distillation

Figure 10. Profile of average temperature of the liquid phase from the inlet of the liquid up to the weir.

6608


Figure 11. Snapshots of the central plane from the liquid inlet up to the weir of the sieve for the liquid volume fraction, liquid temperature, and ethanol mass fraction in the vapor and liquid phases.

sieve tray efficiencies based, for example, on the fundamentals proposed by Murphree44 and West et al.45 The global tray efficiency can be calculated based on the Murphree44 efficiency: ηnv

)

xAn g - xAn-1 g xAn g* - xAn-1 g

(32)

where xnAg* is the vapor ethanol mole fraction that would be in equilibrium with the liquid that leaves the tray. Equation 32 is the Murphree44 efficiency based on the vapor phase and relates the vapor stream composition that enters and leaves the tray with that which would leave the tray if it were in equilibrium with the liquid that leaves the tray. In this regard, Figure 12 gives the ethanol mass fractions in the vapor phase (Figure 12a) and in the liquid phase (Figure 12b) that leave the tray as a function of time. It can be seen that the ethanol mass fractions in the two phases have a chaotic behavior similar to the other conservation properties. The

Figure 12. Ethanol mass fraction in the streams that leave the tray as a function of time: (a) vapor phase; (b) liquid phase.

Figure 13. Vapor and liquid ethanol mass fraction profiles and vapor ethanol mass fraction that would be in equilibrium with the liquid phase as a function of the position.


Figure 14. Point efficiency distribution.

average of these properties can be calculated applying the time average, eq 26. In the vapor phase, the mean value for the ethanol mass fraction is 0.8465, and in the liquid phase it is 0.7605. From the results in Figure 12, it is possible to determine the global tray efficiency. For the operational and geometric

6609

conditions used in this case, the efficiency for this plate and system is close to 64.16%. This value lies within the efficiency zone obtained by Barros et al.,5 who used the nonequilibrium stage model. When the point efficiency concept suggested by West45 is applied, where the tray is divided into regions, and then the Murphree44 efficiency is applied, the efficiency as a function of the position can be calculated. This information can be used to understand the tray behavior and then optimize the tray parameters. Figure 13 shows the vapor and liquid ethanol mass fraction profiles as a function of the position and ethanol mass fraction of the vapor that would be in equilibrium with the liquid phase. It can be observed that the liquid ethanol mass fraction decreases with greater intensity near the liquid inlet. However, it is clear that there is a gradient concentration at the tray outlet, showing that the thermodynamic equilibrium was not yet reached and, consequently, the efficiency is not 100%. For the vapor ethanol mass fraction, the greater decrease observed in the first part of the tray is due to the large concentration gradient in this region. In other regions, this rate is lower and constant due to the mixture in the vapor phase. The point efficiency is observed in Figure 14. It can be noted that the tray has higher efficiencies near the inlet and the weir. That happens due to the circulating zones that appear together at the inlet and the weir. On the other hand, these circulating zones may be considered as a sort of back mixing and have a

Figure 15. Vector plot and stream lines: (a) vapor phase in a central plane from the inlet to the weir; (b) liquid phase in a central plane from the inlet to the weir; (c) liquid phase in a cross-sectional plane at 10 mm above the tray bottom.

6610


negative effect on the tray efficiency. Near the weir, the efficiency reaches the maximum value because the circulating zones increase the liquid residence time and the ethanol in the vapor phase is at the minimum; consequently, the liquid gets more saturated. In the middle of the tray the efficiency reaches the minimum value due to the liquid phase having preferential flow paths in this region. In general, the net contribution to the efficiency of the tray due to the circulating zones is usually negative, because throughput of these zones is small and its occupation of the tray space is high. This behavior can be seen in Figure 15, which shows the vector plot and two-dimensional stream lines for the vapor and liquid phases in a central plane and a cross-sectional plane at 10 mm above the tray bottom. It can be noted in Figure 15 that the vapor phase rises through the froth without recirculating. However, the liquid phase presents two main recirculating regions: near the liquid inlet and near the weir. This elucidates the high efficiency noted in these regions. Near the base tray the stream lines suggest a flow from the wall to the center as can be seen in Figure 15c. From these results, it is possible to understand better the flow behavior in the sieve tray and its influences on energy and mass transfer. Moreover, it was possible to evaluate the point efficiency and the global efficiency of the sieve tray. Conclusions The methodology proposed in this was shown to be appropriate to describe the vapor-liquid flow on a distillation sieve distillation tray. The CFD model demonstrated good agreement with the experimental data for clear liquid height and the behavior of the sieve tray flow, as verified through qualitative analysis with an experimental sieve tray. The model was validated by the hydrodynamics point of view, and it is necessary to validate in terms of heat and mass. The main results of the simulations show temperature, concentration, and holdup profiles to be a function of time and position on the sieve tray, confirming the need for the analysis proposed herein. It was possible to calculate the tray efficiency based on models that consider spatial variations in all properties, not only pressure and velocity. This is important in terms of the design and optimization of distillation trays. The next step will be to use this model with other systems and geometries. Furthermore, it is intended to acquire the experimental data for temperature and concentration fields and to compare the results with predicted values in order to validate the model in terms of heat and mass transfer. Finally, the CFD tools presented and discussed allow a better understanding of the turbulent vapor-liquid flow on the sieve plates of distillation columns and they can be used to optimize the design and operating conditions of such processes. Acknowledgment The authors acknowledge CNPq and CAPES for financial support. Nomenclature a, cv ) constants A, B, C ) Antoine constants Ab ) active bubbling area [m2] Ag ) interfacial bubble area per volume [m2/m3] Apg ) projected bubble area per volume [m2/m3] Ah ) hole area [m2] Ahb ) perforated fractional area Cµ, C1, C2 ) turbulence model constants CD ) drag coefficient

Cp ) specific heat of fluid [J/kg K] D ) column diameter [m], diffusivity [m/s] dh ) hole diameter [m] dg ) bubble diameter [m] f ) volume fraction F ) interphase momentum exchange term [N/m3] g ) gravity vector [m/s2] h ) enthalpy [J/kg] hcl ) clear liquid height [m] hgl ) global heat transfer coefficient [J/s m2 K] hw ) weir height [m] ht ) sieve tray pressure drop [m] hv ) dry sieve tray pressure drop [m] hσ ) pressure drop due to bubble formation [m] i ) turbulence intensity J ) diffusion chemical species flux [kg/s m2] k ) turbulent kinetic energy [m2/s2] Kc ) thermodynamic constant c Kgl ) global mass transfer coefficient [m/s] l, r, q, z, uij ) UNIQUAC parameters Ls ) dissipation length scale [m] nc ) number of chemical species np ) number of phases Nu ) Nusselt number p ) pressure [N/m2] Pr ) Prandtl number q ) diffusion heat flux [J/s m2] Q ) interphase energy exchange term [J/s m3] Ql ) rate of liquid flow [m3/s] Qg ) rate of vapor flow [m3/s] Re ) Reynolds number Sl ) velocity relation t ) time [s] T ) temperature [K] Ug ) superficial gas velocity [m/s] v ) velocity vector [m/s] Vh ) hole velocity [m/s] W ) weir length [m] x ) coordinate [m] Y ) interphase chemical species exchange term [kg/s m3], mole fraction y ) coordinate [m] z ) coordinate [m] Greek Symbols β ) momentum transfer coefficient [kg/s m3] γ ) activity coefficient φ ) conservation properties Φ ) fugacity coefficient ε ) dissipation rate of turbulent kinetic energy [m2/s3] λ ) thermal conductivity [J/s m K] µ ) viscosity [kg/ms] F ) density [kg/m3] σ ) surface tension [N/m] σε, σk ) turbulence model constants ψ ) orthogonal direction to the wall γ ) stress tensor [N/m2] Subscripts g ) vapor phase i ) ith chemical species

Ind. Eng. Chem. Res., Vol. 49, No. 14, 2010 k ) kth phase l ) liquid phase Superscripts C ) combinatorial part eff ) effective property t ) turbulence property T ) transpose of vector R ) residual part

Literature Cited (1) Henley, E. J.; Seader, J. D. Equilibrium-Stage Separation Operations in Chemical Engineering; Wiley: New York, 1981. (2) Zuiderweg, F. J. Chem. Eng. Sci. 1982, 37, 1441–1464. (3) Kister, H. Z. Distillation Design; McGraw-Hill: New York, 1992. (4) Krishnamurthy, R.; Taylor, R. AIChE J. 1985, 31, 456–465. (5) Barros, A. A. C.; Pescarini, M. H.; Maciel, M. R. W. Comput.-Aided Process Eng. 1996, 279–284. (6) Colwell, C. J. Ind. Eng. Chem. Process Des. DeV. 1979, 20, 298– 307. (7) Bennett, D. L.; Agrawal, R.; Cook, P. J. AIChE J. 1983, 26, 434– 442. (8) English, G. E.; van Winkle, M. Chem. Eng. Sci. 1983, 70, 241–246. (9) Bell, R. L.; Solari, R. B. AIChE J. 1983, 19, 688. (10) Solari, R. B.; Bell, R. L. AIChE J. 1986, 32, 640–649. (11) Bertodado, M. L. Nucl. Eng. Des. 1998, 179, 65–74. (12) Ropelato, K.; Meier, H. F.; Cremasco, M. A. Powder Technol. 2005, 154, 179–184. (13) Meier, H. F.; Ropelato, K.; Forster, H.; Iess, J. J.; Mori, M. ZKG Int. 2002, 22, 65–75. (14) Meier, H.; Mori, M. Comput. Chem. Eng. 1998, 22, 5641–5644. (15) Lun, C. K. K. Int. J. Multiphase Flow 2000, 26, 1707–1736. (16) Cheng, Y.; Wei, F.; Guo, Y.; Jin, Y.; Lin, W. Powder Technol. 1998, 98, 151–156. (17) Cockx, A.; Do-Quang, Z.; Ling, A.; Roustan, M. Chem. Eng. Sci. 1999, 54, 5085–5090. (18) Cheng, Y.; Wei, F.; Guo, Y.; Jin, Y. Chem. Eng. Sci. 2001, 56, 1687. (19) Boisson, N.; Malin, M. R. Int. J. Numer. Methods Fluids 1996, 23, 1289–1310. (20) Krishna, R.; van Baten, J. M. Catal. Today 2003, 79-80, 67–75. (21) Krishna, R.; van Baten, J. M.; Urseanu, M. I. Chem. Eng. Sci. 2000, 55, 3275–3286. (22) Krishna, R.; van Baten, J. M.; Urseanu, M. I.; Ellenberger, J. Chem. Eng, Sci. 2001, 56, 537–545.

6611

(23) Delnoij, E.; Kuipers, J. A.; Swaaij, W. P. M. Chem. Eng. Sci. 1997, 52, 3759–3772. (24) Delnoij, E.; Kuipers, J. A.; Swaaij, W. P. M. Chem. Eng. Sci. 1999, 54, 2217–2226. (25) Michele, V.; Hempel, D. C. Chem. Eng. Sci. 2002, 57, 1899–1908. (26) Pan, Y.; Dudukovic, M. P.; Chang, M. Chem. Eng. Sci. 1999, 54, 2481–2489. (27) Pfleger, D.; Becker, S. Chem. Eng. Sci. 2001, 56, 1737–1747. (28) Sokolichin, A.; Eigenberger, G. Chem. Eng. Sci. 1999, 54, 2273– 2284. (29) Spick, P.; Dias, M. M.; Lopes, J. C. M. Chem. Eng. Sci. 2001, 56, 6367–6383. (30) Wiemann, D.; Lehr, F.; Mewes, D. International Conference on Distillation & Absorption [CD-ROM]; 2002. (31) Mehta, B.; Chuang, K.; NandaKumar, K. Trans. Inst. Chem. Eng. 1998, 76, 843–848. (32) Liu, C. J.; Yuan, X.; Yu, K. T.; Zhu, X. J. Chem. Eng. Sci. 2000, 55, 2287–2294. (33) Krishna, R.; van Baten, J. M.; Ellenberger, J.; Higler, A. P.; Taylor, R. Chem. Eng. Res. Des. 1999, 77, 639–646. (34) van Baten, J. M.; Krishna, R. Chem. Eng. J. 2000, 77, 143–151. (35) Soares, C.; Noriler, D.; Barros, A. A. C.; Meier, H. F.; Maciel, M. R. W. International Conference on Distillation & Absorption [CDROM]; 2002. (36) Gesit, G.; Nandakumar, K.; Chuang., K. T. AIChE J. 2003, 49, 910–924. (37) Wang, X. L.; Liu, C. J.; Yuan, X. G.; Yu, K. T. Ind. Eng. Chem. Res. 2004, 43, 2556–2567. (38) Sun, Z. M.; Liu, B. T.; Yuan, X. G.; Liu, C. J.; Yu, K. T. Ind. Eng. Chem. Res. 2005, 44, 4427–4434. (39) Noriler, D.; Barros, A. A. C.; Meier, H. F.; Maciel, M. R. W. Chem. Eng. J. 2008, 136, 133–143. (40) Rahimi, R.; Rahimi, M.; Shahraki, F.; Zivdar, M. Inst. Chem. Eng. 2006, 220–229. (41) Ranz, W.; Marshall, W. R. Chem. Eng. Prog. 1952, 48, 141. (42) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. Properties of gases and liquids; Wiley: New York, 1987. (43) Noriler, D. Mathematical Modelling and Numerical Simulation of Vapor-Liquid Flow on a Distillation Sieve Tray. Ph.D. Thesis, State University of Campinas, 2007. (44) Murphree, E. V. Ind. Eng. Chem. 1925, 17, 747. (45) West, F. B. Ind. Eng. Chem. 1952, 10, 2470–2478.

ReceiVed for reView September 5, 2009 ReVised manuscript receiVed March 30, 2010 Accepted May 26, 2010 IE9013925

Simultaneous Momentum, Mass, and Energy Transfer Analysis of a

Recommend Documents