A discrete multicomponent fuel evaporation model with liquid turbulence effects

International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A discrete multicomponent fuel evaporation model with liquid turbulence effects O. Samimi Abianeh, C. P. Chen ⇑ Department of Chemical and Materials Engineering, University of Alabama in Huntsville, Huntsville, AL 35899, United States

a r t i c l e

i n f o

Article history: Received 27 January 2012 Received in revised form 27 June 2012 Accepted 2 July 2012 Available online 2 August 2012 Keywords: Spray droplet evaporation Turbulent atomization Multi component fuel

a b s t r a c t A new approach to simultaneously account for finite thermal conductivity, finite mass diffusivity and turbulence effects within atomizing multicomponent liquid fuel sprays has been developed in this study. The main contribution of this paper is to incorporate the liquid turbulence effect in modeling the boundary layer heat and mass resistance during multi-component droplet evaporation. The finite conductivity model is based on an existing two-layer film theory, where the turbulence characteristics of the droplet are used to estimate the effective thermal conductivity. The present paper extends the two-layer film theory formulation to include multi-component mass diffusivities within the droplet liquid phase. In this model four regions are considered: the interior region of the droplet, droplet-side interface, gas-side interface, and the surrounding gas phase. Approximate solutions to the quasi-steady energy and mass transfer equations were used to derive an explicit expression for the heat and mass flux from the surrounding gas to the droplet–gas interface, and within the multi-component droplet. Extension of the model to high pressures using the Peng–Robinson equation of state is also considered. The validation study was carried out for a bi-component decane/hexadecane fuel, followed by application studies of complex gasoline–ethanol blended fuels evaporating in hot gas environments. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction In liquid fueled combustion systems such as gasoline direct injection and diesel engines, combustion performance depends strongly on the fuel type, liquid atomization, spray vaporization, and the mixing process. The droplet size history during evaporation in the combustion chamber influences the dynamic behavior of the droplets, whereas the variation of the composition determines the distribution of the fuel compounds within the combustion chamber. The fundamental understanding of these processes is essential for the modeling of evaporating fuel sprays. A comprehensive overview of droplet evaporation can be reference in [1]. Commercial fuels are typically composed of hundreds of complex compounds with different physical properties. Especially for fuels derived from bio-mass, the range of relative volatility is very large. The highly volatile (light) components, such as water and simple acids, evaporate relatively fast at the liquid–vapor interface compared to the internal diffusion, leading to non-uniform concentration distributions and a build-up of lighter components in the center of the fuel droplet [2]. Thus, the evaporation of multi component droplets is a complex process. To efficiently predict the real fuel evaporating sprays, surrogate fuels involving a small number of components matching distillation curves are usually used to

⇑ Corresponding author. Tel.: +1 (256) 824 6194; fax: +1 (256) 824 6839. E-mail address: [email protected] (C. P. Chen). 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.07.003

represent the real fuels. In [3], a seven-component surrogate was used as gasoline simulant, and a four-component surrogate was proposed for simulating kerosene fuel [4]. Numerous theoretical and modeling studies have been carried out for describing droplet evaporation [1–11], and these studies can be classified into general categories such as: discrete or continues thermodynamic modeling [2–8], and how the finite thermal/ mass transfers within droplets was modeled [2–5,9–11]. The Discrete Multi-Component (DMC) approach has high computational overhead because additional transport equations are to be solved for each species in order to track the fuel composition and the vaporization behavior. On the other hand, the Continuous Multicomponent Model (CMC) represents the fuel composition as a continuous distribution function of molecular weight or other appropriate parameters; however, modeling the consumption of each species is complicated, when using this model in combustion simulation [3]. Rapid Mixing Model (RMM), Thin Skin Model (TSM) and Diffusion Limit Model (DLM) are different approaches for heat/mass rate description within fuel droplet [1]. The heat/mass diffusion resistance within the droplet is zero for the RMM approach. This model has reasonable results for slow evaporation processes, when droplet internal heat conduction and diffusion do not have major effect on the internal temperature and concentration profiles [12]. In the TSM approach, only a thin layer on the droplet surface is heated and evaporated. Hence the transient heat transport within the droplet is completely neglected. It is assumed that the droplet

6898

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

Nomenclature a A ca cD cp D h J k Kc Ke L Lin Lt Lw _ m M Nu P Pr Q r rp rt rw R Re Sc s Sh t T v

radius of blob or parent drop droplet surface area constant (=0.0163) discharge coefficient of injection nozzle specific heat at constant pressure binary diffusion coefficient heat transfer coefficient mass transfer coefficient thermal conductivity loss coefficient due to nozzle inlet geometry turbulence model constant (=0.27) latent heat of evaporation length of injection nozzle characteristic length scale associated with turbulence characteristic length scale associated with surface wave motion mass transfer rate molecular weight Nusselt number pressure Prandtl number heat transfer rate instantaneous drop radius radius of parent drop used in the secondary break up formulation radial length scale associated with surface wave motion instantaneous drop radius universal gas constant Reynolds number Schmidt number nozzle contraction area ratio Sherwood number time temperature velocity

surface adapts instantaneously to the changes of ambient conditions [13]. The validity of the model is limited to cases in which the droplet reaches its boiling temperature immediately after the start of the evaporation process [12]. Using DLM approach, the temperature and concentrations vary temporally as well as spatially within the liquid droplet. The temperature and concentration distributions inside the droplet can be determined either by theoretical/modeling approaches or using a numerical discretization scheme [14] to solve the equations of heat conduction and mass diffusion. The major disadvantage of discretized DLM model is the enormous computational effort involved [11,14]. In this study, a new finite thermal-conductivity/mass-diffusivity model for heat and mass transfer within a multicomponent liquid fuel droplet was developed. The finite heat and mass effects are phenomenologically modeled through a thermal boundary layer and mass diffusion boundary layer within the droplet. The current model was based on the hybrid primary/secondary atomization model, the T-blob/T-TAB model, of [15,20]. In this model, the liquid droplets were assumed to carry turbulence from issuing liquid injector flow field. The importance of liquid turbulence within atomizing droplets was described in [16,17]. This liquid turbulence was used to model a thermal resistance layer for single component fuel as described in [5]. In [5], both one-way and two-way coupled single component evaporating spray studies were carried out to validate the two-layer film theory model. The current study extents this two-temperature model for multicomponent fuels.

We x y

Weber number mole fraction mass fraction

Greek letters density surface tension turbulence kinetic energy in j–e model e dissipation rate of kinetic energy in j–e model s time associated with surface wave motion sw characteristic time scale associated with surface wave motion st caracteristic time scale associated with liquid turbulence d thermal layer or mass transfer layer thickness l viscosity

q r j

Superscripts g gas phase l liquid phase o initial condition p at the constant pressure s at the droplet surface t parameter associated with turbulence 1 free stream condition Subscripts b associated with boiling condition c associated with critical condition d associated with droplet E associated with properties of ethanol G associated with properties of gasoline i the ith species m associated with mixture  first derivative with respect to time

The inherent turbulence in the injected fuel spray affects the heat and mass transfer rates of the vaporization process. The effects of these changes in the rates have to be accounted for in the numerical models for spray evaporation. Detailed model description and validations can be found in [5,15,18,20], and it is suffice to state that within each numerical droplet, turbulence characteristics such as fluctuating velocity level, length and time scales are supplied by the model. Effective mass and thermal liquid diffusivities are implemented in this model to account for internal circulation within the liquid caused by relative motion between the liquid and the gas and also inherent turbulence in the injected fuel drop. Contrary to [11], detailed temperature and concentration profiles (described by partial differential equations) within the droplet are not solved; thus the current modeling approach is much more efficient. The effective turbulent mass diffusion coefficient proposed in this research encompasses the method of laminar binary diffusion coefficient of Sazhin et al. [19] for calculating mass transfer inside the drop. The current model provides new contribution in modeling finite-conductivity, finite-mass-transfer and liquid turbulence simultaneously for multicomponent fuel evaporation. The model development will be described in the following sections.

2. Theoretical formulations For the multi-component fuel droplet, four regions as shown in Fig. 1 are considered: droplet interior region (or droplet core),

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

"

o

k ¼

o

e

Fig. 1. Schematic of the modeled four regions of an evaporating droplet: droplet interior region, droplet-side interface, gas-side interface, and surrounding gas.

U 1  K c  ð1  s2 Þ ; 8L=Dnozzle C 2D

" # U3 1 2 ¼ Ke  K c  ð1  s Þ : 2Linj C 2D

6899

#

ð1Þ

ð2Þ

U is the liquid velocity at the injection nozzle, which has length Linj, and diameter Dnozzle. The discharge coefficient, the loss coefficient due to the nozzle entrance sharpness and the down stream-to-upstream contraction area ratio of the injection nozzle are presented by CD, Kc, and s, respectively. After the first ‘‘blob’’ is issued at the injector exit, the k–e model was used to predict turbulence level within each droplet. This was achieved by calculating the k and e within the droplet in Lagrangian tracking method. These k and e values supply the required turbulent characteristic length and time scales within the droplet. The mean droplet size after the first primary breakup is [20]:

liquid-side interface of the droplet, gas-side interface, and the surrounding gas phase (infinity). An approximate solution to the quasi-steady energy equation was used to derive an explicit expression for the heat flux from the surrounding gas to the droplet–gas interface, with inter-diffusion of fuel vapor and the surrounding gas taken into account. Normal evaporation is considered in this study. In the normal evaporation regime, the toP tal mass fraction of fuel vapor at the surface ( Y g;s i ) is less than unity and the Spalding’s mass transfer number is used [1]. In this study, the fuel droplets were assumed to have originated from an atomization/process. A short overview of the specific atomization model will be provided here both for convenience and to highlight the different treatment of the equations.

rw, associated with the wave motion generated droplet radius, can be determined from the Reitz model [21]. The value of rt is estimated from turbulence related breakup droplet radius [20]. ct is a weighting coefficient and is determined by the kinetic energy ratio of the turbulent motion and wave perturbation. The detailed T-blob breakup model can be found in Trinh and Chen [20]. Eq. (3) is used for predicting the initial droplet radius at the beginning of the evaporation calculation or after the breakup.

2.1. Turbulent model within droplet

2.2. Heat transfer formulation

The atomization model used in this study is the T-blob/T-TAB model [20] that extended two widely used models, the Kelvin– Helmholtz (KH) instability of Reitz (the ‘‘blob’’ model) [21] and the Taylor-Analogy Breakup (TAB) secondary droplet breakup of O’Rourke and Amsden [22] to include liquid turbulence effects. Fig. 2 schematically depicts the droplet breakup process. Based on integral analysis of a straight injector, initial liquid turbulent kinetic energy and its dissipation rate at the injector nozzle exit are [18,20]:

A spherical liquid droplet with multi-component vaporizing without chemical reaction in a gaseous environment is considered. Conservation of energy at the droplet surface requires:

rd ¼

rw rt ; ð1  ct Þr t þ ct rw

_ Q g þ Q r  Q l ¼ mL; g

ð3Þ

ð4Þ

where Q is the heat transfer rate from the external gas to the droplet, Ql is the rate of heat transfer within the liquid droplet. The radiative heat transfer Qr may be obtained from a relevant radiation model for a given combustor. In the following analysis, it’s assumed

Fig. 2. Schematic of injection, breakup and evaporation modeling of the fuel droplet.

6900

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

to be negligible. The heat transfer to the droplet is expressed in terms of heat transfer coefficient, as follows: g

g

Q g ¼ h AðT 1  T s Þ ¼ Nug k AðT 1  T s Þ=2r d :

ð5Þ

Different correlations have been proposed for estimation of gas side Nusselt number [1,11,23–25]. In this study, two different correlations (model I and II) were used. The more recent Nusselt number correlation that is widely used (e.g. [1]) for spray calculation is: 0:33 Nug ¼ ð2 þ 0:454Re0:5 Þ lnð1 þ BH Þ=BH : d Pr

ð6Þ

The second correlation (model II) is taken from Haywood et al. [27] as given below: 0:33 Nug ¼ ð2 þ 0:57Re0:5 Þ  ð1 þ BH Þ0:7 ; d Pr !   l p;g Q C ðT 1  T s Þ : BH ¼ 1  g  L Q

2.3. Mass transfer formulation For multi-component droplets, vaporization is governed not only by the volatility of each component but also by the rate of species diffusion, droplet surface regression, and the nature of the fluid motion within the droplet. The general continuity equation for species i with Fickian diffusion is written in the form:

  @Y i þ v  rY i ¼ r  ðqDi rY i Þ þ si ; @t

q

where v is convective reference velocity. si is a source term or rate of production or rate of dissipation of mass of i per unit volume. For a control volume centered around the interface of the droplet with no species accumulation (si = 0):

ZZZ ð7Þ

ð14Þ

q

@Y i dV d ¼ 0: @t

ð15Þ

By integrating the equation of continuity over the control volume: The second model (Eq. (7)) was mainly used in this study for comparison with the detailed numerical method of [28]. The mixture latent heat of evaporation is:

Lmixture ¼

X

yl;s i Li :

ð8Þ

Latent heat of evaporation, as a function of temperature, for each component (Li) was employed and was updated in each time step. Heating of the liquid droplet is treated by considering effective convection effect inside the droplet [5]: l

l

s

keff

d

Q ¼ h AðT  T Þ ¼

dthermal

s

d

AðT  T Þ:

ð9Þ

The averaged temperature of the interior region of the droplet is determined from the conservation of energy applied to the droplet, which is:

dT ¼ Q l: dt

ð10Þ

The thermal boundary layer thickness is calculated based on the effective thermal diffusivity (see [5]):

dthermal ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

v is velocity. Because of the no accumulation condition:

_ ¼ ql ðv l;s  r_ Þ ¼ qg ðv g;s  r_ Þ; m

ð17Þ

and for spherical symmetric droplet [29]: g;s g;1 g;1 g g qg Dgi rY g;s  Y g;s  Y g;s i ¼ J i ðY i i Þ ¼ Shi q Di ðY i i Þ=2r d :

ð18Þ

The species mass transfer within the liquid droplet, by analogy, is: l;s l l;s ql Dli rY l;s i ¼ J i ðY i  Y i Þ:

ð19Þ

Mass transfer rate between the surface and inside of the droplet is modeled following the analogous heat transfer rate formulation [5]: l J l;s i ¼ q

Deff : ddiffusion

ð20Þ

ð11Þ

This is a resistance layer characterization for accounting the finite conductivity effect [5]. Internal recirculation effect was not considered in the present model, but can be added as additional resistance. In the current model, both turbulent and molecular thermal diffusivity are taken into account for evaluating the effective thermal diffusivity: l

k

ql C

Deff ¼ Dl þ Dt ;

ð21Þ

l

in which, D is the molecular diffusivity coefficient [30]; and the turbulent diffusivity is modeled by: 2

paeff st :

aeff ¼ al þ at ¼

where

ð16Þ

The effective mass transfer coefficient is calculated as:

d

mC p;l

g;s l;s l l g g;s g g ql ðv l;s  r_ ÞY l;s  r_ ÞY g;s i  q Di rY i ¼ q ðv i  q Di rY i ;

p;l

þ Cl

ðkÞ2 ; ePrt

ð12Þ

Dt ¼ C l

k ; eSct

ð22Þ

here Sct is the turbulent Schmidt number. It should be noted the second term on the right hand side of Eq. (21) would vanish when the turbulence is low at the end of the droplet lifetime. The effective diffusivity gradually transitions to laminar diffusivity at the end of the droplet lifetime. The individual species mass transfer rate thus becomes (see also [11,29]): g;s g;s l;s l l;s 1 _ l;s _ g;s _ i ¼ mY m i þ J i ðY i  Y i Þ ¼ mY i þ J i ðY i  Y i Þ:

ð23Þ

Due to the high turbulent diffusivity with respect to the laminar diffusivity, the mass transfer coefficient can be estimated by simplifying Eq. (20):

where Cl = 0.09. Prt is the turbulent Prandtl number and is set to be 0.9 according to a previous study [5]. It should be noted that the second term on the right hand side equation (12) would vanish when the turbulence is low at the end of the droplet lifetime. The effective diffusivity gradually transition to molecular (laminar) diffusivity at the end of the droplet lifetime. The effective thermal conductivity is calculated by the effective thermal diffusivity:

Following the film resistance approach as described in [5], the mass transfer boundary layer thickness is calculated by:

keff ¼ aeff ql C p;l :

ddiffusion ¼

ð13Þ

The liquid droplet turbulence quantities k and e are obtained from the T-blob/T-TAB atomization spray model of Ref. [20] as discussed in Section 2.1.

l Deff  Dt ! J l;s i ¼ q

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pDeff st :

Dt ddiffusion

2

¼ ql C l

k : eSc ddiffusion t

ð24Þ

ð25Þ

st in Eqs. (11) and (25) is the turbulence characteristic time scale and can be estimated based on the analytical transient solution of the k–e turbulence model [20]. It should be noted that both k and

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

e are time dependent, thus the mass transfer boundary layer thickness is also time dependent. Assuming air is insoluble in the liquid phase and summing Eq. (23) over all fuel species results in:

PN

g;s g;s g;1 Þ i¼1 J i ðY i  Y i : PN g;s 1  i¼1 Y i

_ ¼ m

ð26Þ

The instantaneous droplet radius during the evaporation process is calculated from:

r2 ¼

 1=3  1=3 _ m2 ql1 m1  mdt ql  l ¼  1l r1 ; m1 m1 q2 q2

ð27Þ

where the subscripts 1 and 2 denote the beginning and ending states of the droplet, respectively. The average species mass fraction in the interior region of the droplet is calculated from:

Y li;2

_ i dt mi;2 mi;1  m ¼ ¼ : _ m1  mdt m2

ð28Þ

Mass transfer from the droplet surface into the gas environment occurs by the bulk motion and diffusion, and is modeled by: g

g;s g J g;s i ¼ Shi q Di =2r d ;

ð29Þ

For model I, Sherwood number is taken from [26]: g

0:33 Shi ¼ ð2 þ 0:39Re0:5 Þ lnð1 þ Bm Þ=Bm ; d Sc i

ð30Þ

where

l1 Sci ¼ 1 g : q Di

ð31Þ

Again, for model II, Sherwood number is taken from Haywood et al. [27]: g

0:33 Shi ¼ ð2 þ 0:87Re0:5 Þ  ð1 þ Bm Þ0:7 : d Sc i

ð32Þ

Drag coefficient is also taken from Haywood et al. [27]:

 Cd ¼

 24 1 0:63 Þ  ð1 þ BH Þ0:2 : 1 ð1 þ 0:2ðRed Þ Red

ð33Þ

To characterize turbulence within the droplets, the T-blob/TTAB atomization model as described in [20] will be used in this study. The liquid turbulent kinetic energy and its dissipation rate at the injector nozzle exit were specified as the initial conditions for Lagrangian tracking of the droplet. A set of ODEs were derived to track the evaluation of k and e within the droplet according to the T-blob/T-TAB model. The values obtained from the evolution of k and e are used in the estimate of thermal and mass boundary layer thickness to facilitate fine-rate heat and mass transfer within the liquid droplet for calculating the droplet evaporation. In this study, one-way coupled methodology was used for spray calculation with a prescribed gas environment.

6901

formulations can be found in [36,37]. For a multicomponent mixture, the latent heat of vaporization of each species is defined as the difference between the partial molar enthalpy of that species in the vapor and liquid phases. The following thermodynamic relation then gives the partial molar enthalpy of i species [31]: 0

hi  hi ¼ RT 2

@ ðln /i Þ; @T

ð35Þ

where the superscript 0 denotes the quantity in an ideal state. The fugacity concept will be used for the binary component mixture of the decane and hexadecane at high pressure in this study (Section 3.1). If the vapor phase is assumed to be an ideal gas mixture, the gas side surface mole fraction of the fuel vapor can be determined using the modified Raoult’s law [31]: v ap ci X l;s ¼ X g;s i Pi i P system ;

ð36Þ

where ci is the activity coefficient. The modified Raoult’s law is valid for low to moderate pressure [31]. For a given set of free stream conditions, the global conservation equations can be solved for mass transfer and heat transfer. However, all these equations are nonlinear and strongly coupled. It is therefore necessary to adopt an iterative procedure at each time step to find the desired solution. The solution algorithm used in this study is given below: 1. Calculate the droplet diameter after the break up at time step from Eq. (3). 2. Calculate the thermal and mass transfer boundary layer thickness from Eqs. (11) and (25). 3. Guess surface temperature Ts. 4. Guess surface mass fraction Y l;s i : 5. Evaluate fuel vapor mole fraction at gas-side interface X g;s i from Eqs. (34) or (36) and convert it to the mass basis ‘‘Y g;s i ’’. _ from Eq. (26). 6. Evaluate the total mass transfer m _ i from right 7. Evaluate mass transfer for each component m hand side of Eq. (23). 8. Solve left hand side of Eq. (23) to find new Y l;s i : 9. If change of Y l;s relative to the previous iteration is within i specified tolerance proceed to the step 10. Otherwise, return to the step 4. 10. Evaluate the droplet radius from Eq. (27). 11. Evaluate droplet mass fraction at interior region Y li from Eq. (28). 12. Evaluate ‘‘Qg’’ from Eq. (5). 13. Evaluate ‘‘Ql’’ from Eq. (4). 14. Evaluate ‘‘Ts’’ from Eq. (9). 15. Calculate ‘‘Td’’ from Eq. (10). Proceed to next time step. 16. If changes of ‘‘Ts’’ and ‘‘Td’’ relative to the previous iteration are within specified tolerance proceed to the next time step. Otherwise, return to the step 3.

2.4. Models at the droplet–gas interface At the liquid droplet/vapor interface, multi-component phase equilibrium is established to account for the jump conditions [1]. The chemical potential of the liquid phase and the vapor phase at the interface of the droplet are equal for each species at the equilibrium condition, and can be written as [31]:

X gi /gi ¼ X Li /Li ;

ð34Þ

where /i represents fugacity coefficients. For a high-pressure environment, the thermodynamics model needs to consider real gas and liquid behaviors, and the effect of pressure on thermo-transport properties. The real gas and liquid behavior is determined by the Peng–Robinson equation of state [31] in this study. Further detailed

Drag coefficient and thermal properties are updated in each time step. 3. Results and discussion 3.1. Model validation The present phenomenological model introduces a model coefficient, turbulent Schmidt number Sct, which requires tuning. Due to the lack of detailed benchmark experimental data, the multidimensional direct numerical simulation results of [28] were used for determining Sct and for model validation. In their work [28], multi-dimensional transport equations for temperature and con-

6902

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

centrations in gas-phase and within the liquid droplet were discretized and numerically solved. A bi-component decane/hexadecane fuel with homogeneous species concentration of 50% decane and 50% hexadecane by mass was used. The ambient air temperature and pressure are specified to be 1000 K and 10 bar. The initial droplet radius and temperature are 31.8 lm, and 300 K. To be consistent, the same heat and mass transfer coefficient correlations as used in [28,29] (the model II described in Sections 2.2 and 2.3) were used for the gas-side heat and mass convective transport in the validation study. All the thermal properties at the gas side were evaluated based on the free stream ambient temperature and pressure. Within the liquid side, all of the properties are updated at the start of each time step using the database from [30]. Latent heat of evaporation and fuel vapor pressure were calculated based on the droplet surface temperature and the other droplet properties such as thermal conductivity and density were calculated based on droplet interior temperature. Eq. (34) and the Peng–Robinson equation of state were used to account for high pressure (10 bar) effects. Initially the liquid droplet is at 300 K and the initial Reynolds number, based on the free stream condition, droplet diameter and initial relative velocity between droplet and surrounding gas, is 100. The final turbulent Schmidt number for present model is tuned to 0.1, which produced the most favorable comparisons. It should be noted that in the simulation of [28] the Reynolds number range may not coincide with fully turbulent regimes. For the purpose of tuning, the current turbulent Schmidt number may not be optimal. In general, increasing turbulent Schmidt number would decrease the rate of mass transfer between surface and interior region of the droplet due to reduction of the mass transfer coefficient (Eq. (22)). Further effects of different turbulent Schmidt numbers on droplet mass transfer rates can be found in [36]. The turbulent Prandtl number was set at 0.9 as previously determined [5]. For this bi-component fuel, the volatilities of both decane and hexadecane are low; the evaporation rate at the beginning of the evaporation is low. Thus, more energy is transferred into the droplet interior to heat up the droplet, which leads to a decrease in the fuel droplet density, resulting in noticeable swelling of the droplet. This is seen from the profile of the normalized droplet radius that exceeds one as shown in Fig. 3. Also shown is the comparison of the present model predictions with the detailed multi-dimensional numerical results of [28]. Droplet size varies in the manner of the well-known d2-law after the initial swelling period. The averaged evaporation constant for this case is 0.0375117 mm2/s. Maximum differences between the present model and the benchmark numerical results in determining the droplet radius and remaining mass

Fig. 3. Droplet radius and remaining mass history for a decane–hexadecane bi-component droplet; ambient temperature 1000 K and pressure 10 bar.

Fig. 4. Temporal history of the Reynolds number; ambient temperature 1000 K and pressure 10 bar.

are within 2.5% and 4% respectively. The discrepancy is attributed to the difference used in evaluating thermo-physical properties. The temporal history of Reynolds number, which is a combined effect of droplet acceleration and shrinking, is shown in Fig. 4. The predicted results shown in Figs. 3 and 4 are in close agreement with detailed numerical data.

3.2. Multi-component fuel evaporating spray Next we apply the current multi-component model within evaporating spray. The ethanol–gasoline blend evaporating spray within a quiescent nitrogen environment at temperature of 600 K and pressure of 1 bar was considered. Development of the alternative energy resources to replace fossil fuels has great importance at present. A possible alternate fuel is ethanol since it can be made from biomass. In the past few decades, 10% ethanol by volume has been added to automotive gasoline due to its ability to boost the fuel’s octane rating. However, with the increased focus on reducing dependence on fossil fuels, auto manufacturers have been manufacturing cars that are capable of running on gasoline containing 85% ethanol; this fuel known as E85. In this study, different gasoline–ethanol blended fuels, E5 (95% gasoline and 5% of ethanol by volume), E15, E40, E60, and E85, were considered. Typical commercial gasoline composed of hundreds of hydrocarbon species. To be practical in numerical modeling, surrogate fuels with evaporation target is required to mimic real gasoline fuel. Surrogate mixtures are designed with the purpose of numerical simulation of complex mixtures using a small number of components. In this study, surrogate mixtures are chosen to reproduce distillation curves of the gasoline. Here we used the seven-component surrogate proposed by Ra and Reitz [3]. Other choices involving less component number are possible (see [41]). The averaged molecular weight of this gasoline surrogate fuel is 108.5 g/mol, and Ra and Reitz [3] had verified that these seven components can predict the distillation curve of gasoline fuel. Ethanol as a new component is added to the gasoline surrogate. The mole fractions of E5, E15, E40, E60, E85, and gasoline are plotted in Fig. 5. Although the components used in the model may not characterize all of the actual chemical components found in practical pump gasoline–ethanol blended fuels, the selected paraffin fuels were considered adequate to demonstrate the performance of the present multi-component evaporation model.

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

6903

Fig. 5. Composition of surrogate ethanol–gasoline blend fuel.

Ethanol has a highly polar nature and when around 20% ethanol is mixed with gasoline, it forms a positive azeotrope having a higher vapor pressure than either pure gasoline or ethanol [32]. To account for this behavior, the activity coefficients are calculated using the Wilson equation formulation [30]:

lnðcE Þ ¼

 lnðX l;s E

þK

l;s EG X G Þ

þ

X l;s G

! Keg Kge  l;s ; l;s X l;s X G þ KGE X l;s E þ KEG X G E ð37Þ

l;s l;s lnðcG Þ ¼  lnðX l;s G þ KGE X E Þ  X E

! KEG KGE :  l;s l;s X l;s X l;s E þ KEG X G G þ KGE X E ð38Þ

The parameters KEG and KGE have values of 0.1665 and 0.3527 respectively for mixture of ethanol and gasoline [32,33]. The gasoline fuel was treated as a single pseudo-component with respect to ethanol; therefore, the activity coefficients for individual components within gasoline group (seven components) are represented by one. The activity coefficients in Eqs. (37) and (38) were used between the group of seven components in gasoline, or one pseudo component, and ethanol. The modified Raoult’s law (Eq. (36)) was used for multicomponent mixture of ethanol and gasoline due to the lower pressure environment. During the course of the study, it was found that after secondary breakup regime, the droplet diameter is very small and it heats up rapidly. Thus, the current model was only implemented within the T-blob (i.e. the primary breakup) model. Initial droplet diameter is equal to the injector nozzle and after the first break up the radius is equal to 16.81 lm and the Reynolds number, based on the free stream condition, droplet diameter and initial relative velocity between droplet and surrounding gas, is 78. Latent heat of evaporation and fuel vapor pressure were calculated based on droplet surface temperature and the other droplet properties like thermal conductivity and density were calculated based on droplet interior region temperature. Initially droplet temperature is equal to 300 K. Thermal and transport properties of the liquid phase are calculated as described by Reid et al. [30] and updated at each time step. Again, turbulent Schmidt number is set to 0.1 following the validation study described in Section 3.1, and the turbulent Prandtl number is set to 0.9 following [5]. Fig. 6 shows the behavior under normal evaporation conditions of a multi-component gasoline–ethanol, 15% by volume of ethanol or E15, droplet in terms of its variation of droplet mass and droplet surface area. One of the important results of the atmospheric drop-

Fig. 6. Normalized droplet radius and remaining mass versus time; ambient temperature 600 K and pressure 1 bar.

let evaporation theory is that K (derivative of droplet surface area with respect to time) is constant and independent of droplet size, being function mainly of the fuel properties [1,34]. This result is based on the quasi-steady state theory yielding the well-known d2 law; the droplet size varies in the manner of this law. The averaged evaporation constant for this case is 0.1191 mm2/s. The K value for a gasoline surrogate with seven components is 0.1315. Therefore the evaporation rate of E15 ethanol–gasoline blend decreases when compared with gasoline. Due to the finite-conductivity of the droplet heating [3,5], the temperature difference between the surface and interior region of the droplet increases with time initially and finally start to decrease as seen in Fig. 7 for E15 fuel. For the single component fuel, the temperature deference between the droplet surface and interior monotonously decreased after the initial periods. For multicomponent fuel, this temperature difference persists until 20% of the mass remains (see Fig. 6). At this point, the mass transfer rate _ and droplet side mass transfer coefficient (J l;s (m), i ) in Eq. (23) are _ i ¼ mY _ l;s the same, and left hand side of this equation (m i þ l;s l l;s J i ðY i  Y i Þ) could not be solved for finding an implicit equation for component mass fraction at droplet-side interface (Y l;s i ). At this time the droplet size and turbulence level inside of droplet is small l and can be considered as perfectly mixed (Y l;s i ¼ Y i ). The infinite conductivity/diffusivity model [1] can be utilized, or the secondary break up could happen at this time based on the modeling methodologies of [20,21]. In these models, some smaller droplets will be generated and infinite conductivity/diffusivity model could be uti-

6904

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

transfer causes the droplet to evaporate lighter components instead of heating up. The droplet interior may have the same temperature as the surface temperature depending on the rate of the heat transfer, evaporated mass and diameter of the droplet. For the present case, the surface and interior temperature are the same and droplet encounters mostly evaporation during the initial time period because the boiling temperature of pentane is near the drop initial temperature (time = 0 in Fig. 7). Following this time period, part of the heat transfer causes evaporation of more components on the surface of the droplet, and part of it causes heating up of the droplet interior as shown in Fig. 7. With passing time, the density and specific heat of droplet will increase due to the decrease in more volatile components within the droplet. The temperature difference between droplet bulk and the surface would increase as a result. After 0.5 ms, most ethanol evaporates from surface of the

Fig. 7. History of mixture boiling temperature, temperature of surface and interior region of the droplet, and temperature difference between droplet surface and interior region; ambient temperature 600 K and pressure 1 bar.

lized as a result of very small droplets size and consequently very quick droplet evaporation. As shown in Figs. 9 and 10, the mass fraction at the droplet-side interface and droplet interior region are not the same at this point (20% remaining mass) for some components (e.g. C9H20 and C10H22). Therefore considering droplet breakup after this point (20% remaining mass) and utilizing the infinite conductivity/diffusivity model for these smaller droplets would be acceptable in this case. Typical droplet thermal dynamics encounters both heat up and evaporation process. During the first phase of the droplet life time, the droplet may encounter heating up without significant evaporation. It depends on the difference between boiling temperatures of components and initial temperature of the droplet. At initial evaporation phase the more volatile components with boiling temperature near the droplet initial condition temperature (such as pentane) evaporate faster than other components and the heat

   l _ Fig. 8. Normalized Mass transfer rate mm and relative heat transfer rate QQg vs. _o time; ambient temperature 600 K and pressure 1 bar.

Fig. 9. E15 fuel component mass fraction at droplet-side interface; ambient temperature600 K and pressure 1 bar.

Fig. 10. E15 fuel component mass fraction at interior region of the droplet; ambient temperature 600 K and pressure 1 bar.

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

6905

presence of ethanol at the surface of the droplet. As shown in Fig. 9, the mass fraction of ethanol at the surface of the droplet is almost zero at about 0.50 ms. Therefore, less energy is needed for evaporation and more energy will be transferred to the droplet  l interior after 0.50 ms. As a result, the rate of heat transfer QQg in-

Fig. 11. E15 fuel vapor mass fraction at gas-side interface; ambient temperature 600 K and pressure 1 bar.

droplet (see Fig. 9, this point will be explained later). Therefore, after this time, the temperature of the droplet (surface and interior region of the droplet) and mixture boiling temperature increase with higher rate as shown in Fig. 7. This phenomenon was also observed experimentally and numerically by Hallett and BeauchampKiss [40]. The surface temperature of the droplet is always less than the boiling temperature of the mixture during droplet life time. Therefore, the droplet never experienced boiling condition.   _ The normalized mass transfer rate mm and comparative heat _o  l transfer rate QQg are shown in Fig. 8. The ratio of heat transfer  l Q increases until 0.1 ms. In this period, most of the energy is Qg transferred to the surface for heating up the droplet. The ratio of  l heat transfer QQg decreases between 0.1 and 0.5 ms due to the

creases at this time as shown in Fig. 8. Mass transfer rate also increases after evaporation of ethanol from the surface of the droplet as shown in Fig. 8. This increase of evaporation rate is due to the strong effects of droplet density variation. Droplet density is a function of temperature and composition and those two factors have opposite effects as time goes under the conditions considered (see Figs. 6–8). Preferential evaporation of lighter component tends to increase the droplet density while droplet heat up tends to decrease the liquid density. The overall droplet evaporation rate increase with time, as seen in Fig. 8, indicates the droplet density increases much rapidly than the decreasing rate of drop size, since the K (derivative of droplet surface area with respect to time) value is almost constant (from Fig. 6). This phenomenon is in contrast to single-component (no composition change) evaporation [5], in which the droplet density decreases with time, and thus evaporation rate is likely to decrease with time. As shown in Figs. 9 and 10, due to the preferential vaporization of the more volatile components, the composition of the droplet changes continuously, thus E15 gasoline droplets do not reach an equilibrium state. Surface and interior mass fractions of the more volatile components are high at the early period of the droplet lifetime and are much lower near the end of the droplet lifetime as expected. The resulting trends compared favorably with existing models in the literature (see, e.g. [38]). The total mass fraction of the vapor fuel components at the gas-side interface is less than unity during the droplet lifetime as shown in Fig. 11, which validates that vaporization occurs in the normal evaporation regime. At specific temperature and pressure, fuel vapor pressure of each component in the mixture is a strong function of component latent heat of evaporation and also activity coefficient of the component with respect to the other components. Ethanol in discrete mixture of gasoline–ethanol blend has higher latent heat of evaporation (energy/mol) than the first four more volatile components: iC5H12, iC6H14, iC7H16, and iC8H18. However, ethanol life time in

Fig. 12. Droplet radius squared, r 2d , history for different fuels; ambient temperature 600 K and pressure 1 bar.

6906

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

Fig. 13. Total fuel vapor mass fraction at gas-side interface for different blended fuels; ambient temperature 600 K and pressure 1 bar.

Fig. 14. Evaporation of ethanol for different fuels at droplet-side interface; ambient temperature 600 K and pressure 1 bar.

the mixture is shorter than iC7H16 and iC8H18 as shown in Figs. 9 and 10. This is due to the non-ideal thermodynamic behavior of the mixture and the individual component does not evaporate just based on the volatility. This phenomenon is modeled by using the activity coefficient (see Eqs. (37), (38), and (36)). Due to this nonideality correction, ethanol has higher vapor mass fraction at the gas-side interface and reaches its peak of vapor pressure faster than iC7H16, and iC8H18 as shown in Fig. 11. Surface mass fraction of the light component (iC5H12) is high during the early period of the evaporation and then it rapidly decreases as the vaporization continues. In the later stage of evaporation mostly the heavy-end component (dodecane, C12H26) remains in the droplet, as shown in Fig. 11. Time history change of droplet radius squared (r2d ) for different blends are shown in Fig. 12. The rate of droplet radius change for pure gasoline is larger than the other blends. Moreover, by adding more ethanol to the mixture this rate decreases further. This phenomenon is in agreement with experimental observation of Saha [35] and Basu et al. [39]. Total fuel vapor mass fractions at gas-side interface for different blended fuels are shown in Fig. 13. As can be seen in this figure, increasing the amount of ethanol in blended fuel from 5% to 85% by volume reduces total fuel vapor mass fraction at gas-side interface. Fig. 14 shows the evaporation history of ethanol at the droplet-side interface for different blended fuels. Comparison between

Figs. 13 and 14 shows that adding more ethanol to the mixture is the source of increasing the latent heat of evaporation and decreasing the total fuel vapor mass fraction at gas-side interface with respect to gasoline. From the transient droplet evaporation history, it is also observed that after evaporation of ethanol from the surface of the droplet (Fig. 13), the total fuel vapor mass fraction at gasside interface starts to increase as shown in Fig. 12. In gasoline– ethanol blends, pure gasoline is more volatile than pure ethanol as shown in Figs. 11 and 12.

4. Conclusions A new model for unsteady vaporization of multicomponent sprays using the discrete multicomponent approach was developed. Explicit equations for thermal film thickness and mass-transfer film thickness within the droplet based on liquid turbulence characteristics were derived in this paper. The model was validated against direct numerical simulation results of a decane/hexadecane fuel within high pressure environment with very good comparisons. The only tunable model coefficient, the turbulent Schmidt number liquid phase species, was also determined based on the validation study. Due to the low volatilities of the bi-component decane/hexadecane fuel, noticeable swelling

O. Samimi Abianeh, C. P. Chen / International Journal of Heat and Mass Transfer 55 (2012) 6897–6907

of the droplet during the initial stage of droplet heat-up was correctly predicted. For application, this model was applied to simulate the evaporation process of complex multi-component blended fuels at atmospheric pressure. Four different ethanol/gasoline blended fuels (E5, E15, E40, E60, and E85) were modeled and their evaporation characteristics were compared with gasoline fuel. Pure gasoline is more volatile than gasoline ethanol blend fuel. In this study, the transient behaviors are present during the entire droplet lifetime of multicomponent fuels. Liquid and gas-phase properties change considerably during the droplet lifetime. Consequently, an accurate calculation of the thermo-transport properties should be an essential part of spray computations. The present model is shown to be capable of taking into account many important physical effects including variable thermophysical properties, turbulence and transient heat and mass transfer inside the droplet. Compared to the discretization method [11] used within droplet to account for finite heat/mass transfer effects, only a small fraction (about 10%) of CPU time was required for the present modeling approach. Hence, this model is suitable for two-way coupled spray codes in which a large amount of droplets with changing ambient conditions need to be tracked trough the spray field. Acknowledgements The financial support of the BP is gratefully acknowledged. This research was partially funded by Grant DISL T2-008-UAH from BP/ Gulf of Mexico Initiative. References [1] W.A. Sirignano, Fluid Dynamics and Transport of Droplets and Sprays, Cambridge University Press, Cambridge, 1999. [2] J. Brett, A. Ooi, J. Soria, The effect of internal diffusion on an evaporating bio-oil droplet- the Chemistry Free Case, J. Biomass Bioenergy 34 (8) (2010) 1134– 1140. [3] Y. Ra, R.D. Reitz, A vaporization model for discrete multi-component fuel sprays, Int. J. Multiphase Flow 35 (2) (2009) 101–117. [4] O. Samimi Abianeh, C. P. Chen, Multi-component turbulent droplet evaporation of kerosene fuel in hot gas environment, in: 50th AIAA Aerospace Sciences Meeting, Nashville, Tennessee, 2012, AIAA paper 20120345. [5] M.S. Balasubramanyam, C.P. Chen, H.P. Trinh, A new finite conductivity droplet evaporation model including liquid turbulence effect, J. Heat Transfer 129 (8) (2007) 1082–1087. [6] J. Tamim, W.L.H. Hallett, Continuous thermodynamics model for multicomponent vaporization, Chem. Eng. Sci. 50 (18) (1995) 2933–2942. [7] A.M. Lippert, R.D. Reitz, Modeling of multi-component fuels using continuous distributions with application to droplet evaporation and sprays, in: International Fuels & Lubricants Meeting & Exposition, SAE, Tulsa, OK, 1997, 972882. [8] G.S. Zhu, R.D. Reitz, A model for high pressure vaporization of droplets of complex liquid mixtures using continuous thermodynamics, Int. J. Heat Mass Transfer 45 (3) (2002) 495–507. [9] C.K. Law, W.A. Sirignano, Unsteady droplet combustion with droplet heating. ii: Conduction limit, Combust. Flame 28 (1977) 175–186. [10] A.Y. Tong, W.A. Sirignano, Analysis of vaporizing droplet with slip, internal circulation, and unsteady liquid phase and quasi-steady gas-phase heat transfer, in: ASME-JSME Thermal Engineering Joint Conference, Honolulu, HI, 1983, pp. 481–487. [11] D.J. Torres, P.J. O’Rourke, A.A. Amsden, A discrete multicomponent fuel model, Atomization Sprays 13 (2–3) (2003) 218–259. [12] K. Prommersberger, G. Maier, S. Wittig, Validation and application of a droplet evaporation model for real aviation fuel, in: 92nd Symposium on Gas Turbine Engine Combustion, Emissions and Alternative Fuels, Lisbon, Portugal, 1998.

6907

[13] W.A. Sirignano, Fuel droplet vaporization and spray combustion, Prog. Energy Combust. Sci. 9 (1983) 291–322. [14] D.J. Torres, P.J. O’Rourke, A.A. Amsden, Efficient multicomponent fuel algorithm, Combust. Theory Model. 7 (1) (2003) 67–86. [15] H.P. Trinh, C.P. Chen, M.S. Balasubramanyam, Numerical simulation of liquid jet atomization including turbulence effects, J. Eng. Gas Turbines Power 129 (4) (2007) 920–929. [16] P.K. Wu, R.F. Miranda, G.M. Faeth, Effects of initial flow conditions on primary breakup of non-turbulent and turbulent round jets, Atomization Sprays 5 (2) (1995) 175–196. [17] K.A. Sallam, Z. Dai, G.M. Faeth, Liquid breakup at the surface of turbulent round liquid jets in still gases, Int. J. Multiphase Flow 28 (3) (2002) 427–449. [18] K.Y. Huh, E. Lee, J.Y. Koo, Diesel spray atomization model considering nozzle exit turbulence conditions, Atomization Sprays 8 (4) (1998) 453–469. [19] S.S. Sazhin, A. Elwardany, P.A. Krutitskii, G. Castanent, F. Lemoine, E.M. Sazhina, M.R. Heikkal, A simplified model for bi-component droplet heating and evaporation, Int. J. Heat Mass Transfer 53 (21–22) (2010) 4495–4505. [20] H.P. Trinh, C.P. Chen, Development of liquid jet atomization and breakup models including turbulence effects, Atomization Sprays 16 (8) (2006) 907– 932. [21] R.D. Reitz, Modeling atomization processes in high-pressure vaporizing sprays, Atomization Spray Technol 3 (4) (1987) 309–337. [22] P.J. O’Rourke, A.A. Amsden, The TAB method for numerical calculation of spray droplet breakup, in: SAE International Fuels and Lubricants Meeting and Exposition, Toronto, Ontario, 1987, SAE Technical Paper No. 872089. [23] M. Renksizbulut, M.C. Yuen, Experimental study of droplet evaporation in a high-temperature air stream, J. Heat Transfer 105 (2) (1983) 384–389. [24] B. Abramzon, W.A. Sirignano, Droplet vaporization model for spray combustion calculations, Int. J. Heat Mass Transfer 32 (9) (1989) 1605–1618. [25] G.F. Yao, S.I. Abdel-Khalik, S.M. Ghiaasiaan, An investigation of simple evaporation models used in spray simulations, J. Heat Transfer 125 (1) (2003) 179–183. [26] A.A. Amsden, P.J. O’Rourke, T.D. Butler, Kiva-P: a program for chemically reactive flows with sprays, Los Alamos National Laboratory Report No. LA11560-MS UC-96. [27] R.J. Haywood, R. Nafziger, M. Renksizbulut, A detailed examination of gas and liquid phase transient processes in convective droplet evaporation, J. Heat Transfer 111 (2) (1989) 495–503. [28] M. Renksizbulut, M. Bussmann, Multicomponent droplet evaporation at intermediate Reynolds numbers, Int. J. Heat Mass Transfer 36 (11) (1993) 2827–2835. [29] M. Renksizbulut, M. Bussmann, X. Li, A droplet vaporization model for spray calculation, Part. Part. Syst. Charact. 9 (1–4) (1992) 59–65. [30] R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, fourth ed., McGraw-Hill, New York, 1987. [31] J.M. Smith, H.C. Van Ness, M.M. Abbott, Introduction to Chemical Engineering Thermodynamic, seventh ed., McGraw-Hill, New York, 2005. [32] K. Kar, T. Last, C. Haywood, R.R. Raine, Measurement of vapor pressures and enthalpies of vaporization of gasoline and ethanol blends and their effects on mixture preparation in an SI engine, SAE Int. J. Fuels Lubr. 1 (1) (2009) 132– 144. [33] J.A. Pumphrey, J.I. Brand, W.A. Scheller, Vapor pressure measurements and predictions for alcohol–gasoline blends, Fuel 79 (11) (2000) 1405–1411. [34] J. Bellan, K. Harstad, An All-Pressure Fluid-Drop Model: Validation and Species System Dependency, California Institute of Technology, Pasadena, CA 91 1098099, 2002. [35] A. Saha, Vaporization Characteristic of Pure and Blended Bio-fuel Droplet Injected into Hot Stream of Air, MS Thesis, University of Central Florida, 2010. [36] O. Samimi Abianeh, Multi-component Droplet Evaporation Model, MS Thesis, University of Alabama in Huntsville, 2011. [37] O. Samimi Abianeh, C.P. Chen, Advanced model of bi-component fuel droplet heating and evaporating with liquid turbulence effects at high pressure, in: ASME Internal Combustion Engine, Division Fall, Morgantown, West Virginia, 2011, ICEF2011-60040. [38] S. Yang, Y. Ra, R. D. Reitz, A vaporization model for realistic multi-component fuels, in: ILASS Americas, 22nd Annual Conference on Liquid Atomization and Spray Systems, Cincinnati, OH, 2010, ILASS 2010-123. [39] S. Basu, B. Dubas, A. Saha and R. Kumar, Vaporization characteristics of pure and blended bio-fuel droplets injected into hot stream of air, in: 6th US National Meeting of the Combustion Institute in Ann Arbor, MI, May 2009. [40] W.L.H. Hallett, S. Beauchamp-Kiss, Evaporation of single droplets of ethanol– fuel oil mixture, Fuel 89 (2010) 2496–2504. [41] O. Samimi Abianeh, C.P. Chen, R.L. Cerro, Batch distillation: the forward and backward problems, I&EC Research, under revision.