Superheated fuel injection modeling: An engineering approach

International Journal of Thermal Sciences 50 (2011) 1460e1471

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Superheated fuel injection modeling: An engineering approach S. Negro*, G.M. Bianchi University of Bologna, DIEM, Dipartimento delle Costruzioni Meccaniche, Nucleari, Aeronautiche e di Metallurgia, Viale Risorgimento, 2, 40138 Bologna, Italy

a r t i c l e i n f o

a b s t r a c t

Article history: Received 18 April 2010 Received in revised form 17 March 2011 Accepted 31 March 2011 Available online 6 May 2011

A liquid fuel undergoes a sudden decrease in pressure during an injection process. If the pressure outside the injector is below the saturation pressure of the fuel, superheated conditions are reached and the fuel undergoes a thermodynamic instability. Flash evaporation might take place inside or outside the injector nozzle, thus the liquid jet behaves following different mechanisms. Flashing conditions greatly influences atomization and vaporization processes as well as the mixture formation and combustion. This work presents a homogeneous one-dimensional model for the prediction of flash evaporation in superheated liquid fuel injections, able to deal with both internal and external flashing, useful for initializing 3D spray calculations in CFD codes. In this model a Volume Translated PengeRobinson Equation of State (VTPREOS) is used to calculate thermodynamic properties of hydrocarbon fuels. A Homogeneous Relaxation Model (HRM) is used to predict the evaporation rate during internal flashing which may lead to effervescent atomization outside the nozzle. Standard thermodynamic analysis of the CJ-point is used instead for the evaluation of the evaporation rate outside the nozzle, when an unbroken meta-stable liquid jet is predicted under superheated conditions. Two different sets of numerical simulations have been carried out and the results have been compared versus experimental data in literature. Ó 2011 Elsevier Masson SAS. All rights reserved.

Keywords: Flash evaporation Flash boiling Superheated liquid jet Fuel injection 1D homogeneous two-phase model

1. Introduction In recent years, injection processes involving superheated fuels are being intensively studied in aerospace and automotive sectors. Extensively investigated in nuclear engineering since the sixties, liquid flashing jets are not new in technical literature, but the complex physics behind them is still partially unclear. Superheated liquid jets are now the targets of new interest in modern aircraft and in internal combustion reciprocating engines, and other applications involving liquid fuel injection processes. In advanced aircraft engines a conventional hydrocarbon fuel is being considered for cooling the engine [1e3]. The increase in temperature of the fuel before it is injected and burned in the combustion chamber might lead to superheated conditions of the liquid jet. Flash vaporization might take place affecting the combustion performances and creating the danger of vapor lock within the fuel line [3]. Superheating conditions are likely to be found under partial load conditions in modern internal combustion reciprocating engines as well. In these engines the fuel goes through an increase in temperature due to the high pressure pump (i.e. GDI engines) and heat exchange within the fuel line or within the injector itself. In

* Corresponding author. E-mail address: [email protected] (S. Negro). 1290-0729/$ e see front matter Ó 2011 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.ijthermalsci.2011.03.028

fact, in direct injection engines the fuel injector faces the combustion chamber undergoing high thermal stresses. Also in these cases the fuel works as a coolant, increasing its own temperature when injected. If the pressure inside the combustion chamber is lower than the saturation pressure of the fluid at the injection temperature, the enthalpy of the liquid exceeds the saturated value and might promote flash evaporation, which in turn greatly influences atomization and vaporization processes, as consistently reported from technical literature. Compared to mechanical break-up processes, evidences of enhanced atomization, increased cone angle of the spray, and decrease of intact jet length with increasing superheating degree are some of the correlations that have been found and investigated by Reitz [4], Park and Lee [5] and others. When thermal effects supersede hydrodynamic influences, the jet’s behavior will result thermally driven, thus thermal effects need to be considered. This is due to the different time scale at which heat transfer related phenomena are taking place, compared to inertia related phenomena. Investigations in several working scenarios have shown that superheated fuels can alter mainly depending on nozzle geometry, pressure regimes, superheating degree and heat transfer. Oza and Sinnamon [6], identified two flashing modes in fuel injectors, described as internal and external flashing, depending on whether nucleation takes place within the nozzle or not. Nucleation was found to be promoted by the superheating degree and a high resident time of the fuel within the nozzle. Park and Lee [5]

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

analyzed the insurgence of two-phase flow in glass nozzles and proposed a critical value of length-to-diameter ratio (L/D) of 7 for the transition from one flashing mechanism to the other. Yldiz et al. [7] assumed that the two types of flashing were respectively promoted by a slow rate and a high rate conduction of the superheat from the liquid to the interface. The knowledge of the external-flashing regime has been drastically improved by the work of Vieira and Simoes-Moreira [8]. They studied a flashing mechanism that was taking place at very lowdischarge pressures, using iso-octane through a short, conical convergent nozzle. In their experiment the in-flow nucleation was absent and the jet was exiting the orifice with a liquid unbroken core. A fast evaporation was taking place from the jet surface, which was surrounded by a two-phase flow region and, in some cases, by a shock wave. In presence of shock wave the flow was also choked, indicating that a limit in the mass flow rate was reached. The surface evaporation can be seen as a discontinuity between the phases, or an ‘evaporation wave’. This concept was proposed by Simoes-Moreira and Shepherd [9]. Performing rapid decompression of saturated dodecane in a tube, they observed the propagation of a front, called wave of evaporation, starting on the liquid free surface and propagating into the undisturbed meta-stable liquid, with a steady mean velocity. The presence of the evaporation wave was correlated to high degree of superheat, the absence of activated nucleation sites, and the depressurization time scale. The evaporation wave was modeled using jump conditions as in case of deflagrations, in analogy with combustion fronts such as flames. According to these experimental results, two different flashing mechanisms can be identified: internal flashing with effervescent atomization outside the nozzle, and external flashing with surface evaporation outside the nozzle. The first regime is characterized by a second dispersed phase present within the nozzle. The bubble nucleation and growth along the nozzle are a consequence of the superheated degree of the liquid. Vigorous evaporation usually takes place and the jet behaves randomly outside the nozzle, exiting the orifice already broken up. The second regime is typical of low pressure injections but high injection to discharge pressure ratios, high degree of superheat in short convergent nozzles, with smooth entry. In external-flashing regime the in-flow nucleation is absent and the jet exits the orifice with a liquid unbroken core. The jet is usually surrounded by a two-phase region and the evaporation takes place on the liquid jet surface. Due to the different physical structure and research fields of study, these two different flashing mechanisms have been historically approached with different methods. Since nowadays they are both likely to be found in modern fuel injection systems, it is evident the necessity to develop a versatile framework able to deal with both, internal and external flashing.

Nucleation Model Is function of the superheating degree

2. Contribution of the work A 1D homogenous two-phase framework has been developed in order to detect flash evaporation of a superheated liquid inside and outside the injector nozzle. The 1D modeling provides quick responses and low computational efforts, as required by industrial applications. In internal flashing, a non-equilibrium model (HRM) has been coupled with a Volume Translated PengeRobinson Equation of State (VTPR-EOS). The PengeRobinson equation of state [10] is the most popular in petroleum calculations, and here it is used to calculate thermodynamic properties of the vaporizing fluid, granting the code the capability of dealing with hydrocarbon fuels. The saturation pressure of a substance is calculated using fugacity, and the calculation of enthalpy is carried out using departure functions, which account for real fluid behavior. Fig. 1 shows a schematic representation of the internal flashing model. An external-flashing model is presented and used instead when an unbroken meta-stable liquid jet is predicted outside the orifice. In this case, a kinetic closure for the model is adopted in order to calculate the superficial evaporation rate of liquid jet. Standard thermodynamic analysis of the CJ-point is used here for evaluate the evaporation mass flow rate. Fig. 2 shows the surface evaporation process assuming an axi-symmetric two-phase expansion. External and internal flashing models run simultaneously, using the same set of conservation equations within the nozzle. Validation is proposed versus two experimental test cases from literature [8,11], in which one flashing mechanism is predominant. Since flash vaporization is thermally driven, in both cases a 1D framework can be successful. The model is coded in FORTRAN and named 4-flash. 3. Two-phase nozzle-flow model 3.1. Governing equations The governing equations deal with a one-dimensional, compressible and adiabatic flow. The two phases are homogeneous and there is no slip between them. The conservation laws for mass, momentum and energy plus the transport equation for the vapor mass fraction, for a nozzle with area variation F and a flow with evaporation rate G, are respectively

1 vp vðruÞ ru dF þ þ ¼ 0 vx F vx a2 vt

(1)


 ru2 dF vðruÞ v ru2 þ p þ ¼ 0 þ F vx vt vx

(2)

The nozzle flow

Homogenous Relaxation Model (HRM) A relaxation time is required to model the transition from non equilibrium to equilibrium

1461

Nozzle-Spray Coupling

Spray Evaporation

Atomization Model Relates the bubble diameter to the void fraction and the superheated degree at nozzle exit

Fig. 1. Schematic representation of the flashing nozzle-flow coupled with the effervescent atomization model.

1462

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

Two phase region The nozzle flow

Intact superheated liquid core

Absence of in-flow nucleation

Surface evaporation front

Surface Evaporation Model Relates the nozzle-flow conditions to a superheated fast evaporating external jet Fig. 2. In absence of in-flow nucleation an intact liquid core is generated outside the nozzle and, in superheating conditions, a fast evaporation takes place on the surface of the liquid jet. The expansion can be considered axi-symmetric.

vðreo Þ vðruðe0 þ p=rÞÞ ruðe0 þ p=rÞ dF þ þ ¼ EV vt vx F vx

(3)

vðrYÞ vðruYÞ ruY dF þ þ ¼ G vt vx F vx

(4)

where r is the mixture density, rl is the density of the liquid, p is the mixture pressure and u is the mixture velocity. In the continuity equation (1) the mixture density time derivative has been replaced by the pressure time derivative, using the speed of sound squared. The source term in the energy equation (3) is contributions of the evaporation process. The term EV is computed from the knowledge of latent heat of vaporization and the vapor generation rate G.

where Hls and Hvs are respectively the saturated liquid specific enthalpy and the saturated vapor specific enthalpy at local vapor pressure, and H the local liquid specific enthalpy. The source term due to evaporation in equation (3) is calculated as

EV ¼ G$ðHvs  Hls Þ

(10)

In this work steady boundary conditions have been considered. It is important to point out that vapor predictions are subject to the relaxation time (equations (6) and (7)). Values of vapor fraction greater than zero require a relaxation time lower than the resident time of the liquid within the nozzle. Below this time limit, the contribution of equation (6) will generate void fraction within the nozzle.

3.2. Vaporization sub-model 3.3. Nucleation and atomization sub-models Based on the assumption of homogeneous two-phase mixture, the mixture density is given as a function of the vapor mass fraction

r¼

rv rl

Y$rl þ ð1  YÞ$rv

(5)

The vaporization sub-model is based on the Homogeneous Relaxation Model (HRM), a non-equilibrium model for two-phase flow, based on the concept that a relaxation time is required for the vapor mass fraction to adjust to its equilibrium value. The closure of the model is obtained computing the vapor mass fraction source term G

G ¼ r

Y  Ye

Qx

(6)

where Qx is a semi-empirical dryness fraction relaxation time, computed according to Downar-Zapolski and Bilicki [12]. Qx is a monotonically decreasing function of the void fraction a and the non-dimensional pressure difference f

Qx ¼ Q0 a0:257 f2:24

(7)

where

f¼

½Psat ðTin Þ  P Psat ðTin Þ

(8)

Q0 was set equal to 6.51 104 s, according to Ref. [12]. The vapor mass fraction at equilibrium Ye is function of local enthalpy and is given by

Ye ¼

H  Hls Hvs  Hls

(9)

Effervescent atomization is a particular break-up mechanism, in which thermal effects are stronger compared to aerodynamic and inertia effects. For a spray of this kind the atomization sub-model should consider in-flow nucleation, bubble growth and liquid jet break-up mechanisms. A detailed discussion on the thermodynamic break-up of the liquid jet has been presented by authors in Refs. [13] and [14], based on the work proposed by Senda et al. [15] and Kawano et al. [16]. A brief description will only be given in this section. Nucleation processes are divided into two categories, which are homogeneous and heterogeneous nucleation. Homogeneous nucleation takes place within the liquid when the depressurization through the orifice generates a tensile strength that overcomes the molecular forces of the liquid itself. This process seems to happen with more difficulty in the working conditions considered. In heterogeneous nucleation, activated clusters generate bubble nuclei within the liquid close to the boundaries, in presence of a second dispersed phase or surface crevices. In the present study, according to Ref. [15], heterogeneous nucleation due to solid surface is neglected because of the low L/D ratio of the injector nozzle geometry, thus only heterogeneous nucleation due to dissolved gas, such as air, is considered. The nucleation of vapor bubbles is correlated to the superheating degree Dq by this relation [15]

  5; 279 Dw

Nb ¼ CN exp

(11)

where CN is a constant taken equal to 5,757e12 according to Ref. [15]. Then, based on Eq. (11), and the knowledge of the volume occupied by the vapor on the outlet section, bubbles radius Rb can be determined

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

Rb ¼

sffiffiffiffiffiffiffiffiffiffiffiffi 3 3$Vv 4pNb

(12)

Finally, under the assumption that the number of liquid droplets formed because of flashing atomization is twice the vapor bubble number [15], a mean droplet radius Rd is determined. The atomization sub-model is triggered only if void fraction is predicted within the nozzle, because of the degree of superheat, and does not account for aerodynamic break-up that can be included using a secondary break-up model. 3.4. Boundary conditions A 1D-model has limited ways for describing the evolution of the local pressure with particular reference to the pressure drop because of three-dimensional effects within the nozzle. This is the main reason why the present model focuses on the thermal effects only, being aware that cavitation, induced by a pressure drop, might play a role in nucleation. The best possible care has been taken in order to match the boundary conditions in the two sets of simulations. In the entrance section it is assumed that a vena contracta is formed due to the turn, and the minimum pressure is calculated using the Bernoulli equation and loss coefficients, which are different for the two different nozzles (0.62e0.92). In this regard it is worthy to mention the work of von Kuensberg Sarre et al. [17], who proposed a methodology for investigating the effect of the fillet-radius-to-nozzle-diameter ratio (r/D) in the entrance section, on the spray characteristics in diesel injector orifices. However, with respect to the nozzle cross section, any attempt to model the pressure using a 1D nozzle-flow model will result in an average value of the pressure, and it would not account for local variations. Notwithstanding the existence of mentioned limits, the differences in hydrodynamic conditions, fluid temperatures and saturation properties in the two test cases proposed in our work were large enough to keep the average pressure above or below the relative saturation pressures, thus the nucleation triggering was consistent with the experiments without introducing any modification in the model. Initial void fraction and temperature were also defined in the first cell of the domain. The first one has been numerically low bounded to 1e-14. 3.5. Nozzle-flow numerical algorithm The system of four equations (1)e(4) can be written in symbolic vector form

vW vFðWÞ þ þ CðWÞ ¼ 0 vt vx

(13)

Variation Diminishing) scheme as proposed by Davis is used as flux limiter [18,19].

4. Thermodynamic state calculation 4.1. Equation of state In order to properly describe the behavior of an evaporating high temperature two-phase system, it is necessary to have an adequate description of the fluid state. When a pure real liquid is isothermally decompressed (compressed) while below its critical conditions, the actual behavior of the substance shows the evaporation (condensation) transition at constant pressure. So far the best solution in terms of accuracy and computational efficiency is represented by a semi-empirical generalized cubical equation of state, suitable for real fluids and in particular for hydrocarbons fuels [20,21]. In the present model the PengeRobinson equation of state has been chosen [10]

P ¼

RT aðTÞ  v  b vðv þ bÞ þ bðv  bÞ

    Z 3  ð1  BÞZ 2 þ A  3B2  2B Z  AB  B2  B3 ¼ 0

A ¼

aP R2 T 2

(16)

B ¼

bP RT

(17)

Z ¼

Pv RT

(18)

Equation (15) generates either one or three roots depending on the number of phases in the system. Dealing with physical quantities only real roots are of interest. In the case of a two-phase solution the largest root is the compressibility factor for the vapor phase while the smallest positive root is the compressibility factor for the liquid phase. The intermediate root is not computed as it is physically meaningless. The coefficient for equations (15)e(18) are obtained applying equation (14) at the critical point,

aðTÞ ¼ aðTC ÞaðTr ; uÞ

aðTr ; uÞ ¼ 1 þ k 1 

ru ru2



3

0 7 dðlnFÞ 6 0 7 6 7 6 6 þ4 7 CðWÞ ¼ 4 ruðeo þ p=rÞ 5 dx EV 5 ruY G In order to solve the discretized system of equations, a LaxWendroff two-step method has been adopted. This method is second order accurate in both space and time. A TVD (Total

(19)

R2 T 2 PC

3 2 1 3 2 ru p 2þp 2 r u 7 6a 7 6  6 ru 7 p 7 W ¼ 6 5 4 5 FðWÞ ¼ 4 ru eo þ reo r rY ruY 2

(15)

where

aðTC Þ ¼ 0:45724

3

(14)

The PR-EOS can describe the thermodynamic fluid state in subcritical, critical and supercritical conditions. Equation (14) can be rewritten in terms of compressibility factor Z

where

2

1463



(20) pffiffiffiffiffi2 Tr

(21)

kðuÞ ¼ 0:37464 þ 1:54226u  0:26992u2

(22)

RT b ¼ 0:07780 C PC

(23)

where a(Tr,u) is a dimensionless function of reduced temperature and acentric factor and is equal to unity at the critical temperature. The term k is a constant characteristic of each substance, which has been correlated with the acentric factor u [10]. Fig. 3 shows

1464

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

4.3. Enthalpy departure function

4 106 510 K 540 K

6

3.5 10

3 106

Dealing with real fluids it is necessary to consider the deviation from ideal behavior when calculating any thermodynamic extensive property. The enthalpy departure function gives the difference between the real fluid enthalpy and the ideal gas enthalpy. Using a superscript (*) to indicate an ideal gas state, the enthalpy change during a process as proposed by Peng and Robinson reads

supercritical critical

6

P [Pa]

569.4 K 580 K

2.5 10

6

vapour spinodal

2 10

Psat

6

1.5 10

subcritical

6

1 10

liquid spinodal

5

5 10

0 0

100

200

300 Rho [kg/m3]

400

500

Fig. 3. Pressure-density isotherm for n-octane at T ¼ 510 K, T ¼ 540 K, T ¼ 569.4 K and T ¼ 580 K calculated with the PR-EOS.

subcritical, critical and supercritical behavior of predicted isotherms. Table 1 shows the thermodynamic properties of the fluids at the moment implemented in 4-Flash code. 4.2. Fugacity and vapor pressure Any thermodynamic process requires a driving force in order to be carried out. A very convenient approach to equilibrium calculation is the one which uses fugacity. Fugacity, like the chemical potential, describes the tendency of a substance to prefer one phase over another. The phase with the lowest fugacity will be the most favorable. Fugacity also measures how far a real fluid is from its ideal behavior. The ratio between the fugacity and its ideal limit, the pressure, it is called the fugacity coefficient

F ¼

f P

(24)

and becomes equal to unity for an ideal gas. The following expression for fugacity of a pure substance was proposed by Peng and Robinson [10]

  f A Z þ 2:414B ln ¼ Z  1  lnðZ  BÞ  pffiffiffi ln P Z  0:414B 2 2B

(25)

At equilibrium the fugacity of the liquid phase must be equal to the fugacity of the vapor phase, for every component of the system

f

L

¼ f

V

(26)

This relation is satisfied using a convergence criterion, as done in Ref. [10], iterating from the critical point to the triple point. The pressure corresponding to this condition is the saturation pressure.

da T  a Z þ 2:414B pffiffiffi ln h  h* ¼ RTðZ  1Þ þ dT Z  0:414B 2 2b

(27)

The enthalpy per mass unit can be obtained multiplying the molar enthalpy h by the number of moles per mass unit, which is the inverse of the molar weight Mw,

da " # T  a Z þ 2:414B 1 * dT p ffiffiffi h ðTÞ þ RTðZ  1Þ þ H ¼ ln Mw Z  0:414B 2 2b

(28)

The ideal molar enthalpy is calculated as a function of the sole temperature using data from Irvine and Liley [22], Vargaftik [23], NIST [24], or McBride et al. [25]. The latter proposed a polynomial function as follow:

h* ðTÞ 1 1 1 1 b ¼ a1 þ a2 T þ a3 T 2 þ a4 T 3 þ a5 T 4 þ 1 T RT 2 3 4 5

(29)

The ideal enthalpy calculated with (29) does not account for the enthalpy of formation of the species which is already considered in data from Vargaftik. In the present work the enthalpy of liquid at triple point of the working fluid has been set to zero. In order to consider enthalpy of formation of the substance the following formula has been applied when using McBride * HPR

¼

h* ðTÞ MW

!

Z298  H298 þ



 CP dT þ LH TTp

(30)

Ttp

where H298 is the enthalpy of formation at 298 K, CP is the specific heat at constant pressure calculated between the triple point and 298 K and LH(TTp) is the latent heat of vaporization at triple point needed to bring the enthalpy of liquid at triple point to zero. Numerically this has been done simply pre-calculating the enthalpy of the substance of interest at triple point using PR-EOS with the McBride ideal enthalpy, and adding this value to translate the curves successively. Using data from Vargaftik the following expression is used instead


 * HVFTK ¼ CP ðTÞðT  298Þ þ CP 298 298  TTp

(31)

where CP298 is the average value of CP between the two temperatures TTp and 298 K. Figs. 4 and 5 show respectively enthalpy calculations for n-octane and for water at P ¼ 1 bar, using the PengeRobinson

Table 1 Thermodynamic fluid’s properties for the PengeRobinson equation of state.

Formula Boiling T @ P atm [K] Triple point T [K] Critical pressure Pc [MPa] Critical temperature Tc [K] Critical compressibility factor Zc [e] Molar weight [kg/mol] Acentric factor w [e]

n-dodecane

n-octane

n-heptane

Methanol

Ethanol

Water

C12H26 488.61 263.60 1.830 659.50 0.238 170.34e-3 0.575

C8H18 398.77 216.40 2.490 569.10 0.259 114.23e-3 0.398

C7H16 371.56 182.60 2.740 540.10 0.262 100.20e-3 0.349

CH4O 337.77 175.50 8.010 513.00 0.224 32.04e-3 0.556

C2H6O 351.45 159.10 6.360 515.80 0.240 46.07e-3 0.644

H2O 373.15 273.15 22.005 647.30 0.229 18.01e-3 0.344

Enthalpy [J/kg]

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

1.2 10

6

1 10

6

8 10

5

6 10

5

4 10

5

2 10

5

1465

Table 2 Parameters from the Gaussian-like volume correction function database, Monnery et al. [26].

Ideal enthalpy PR-EOS enthalpy Methanol Ethanol n-pentane n-heptane Water

C1

C2

C3

m

n

0.39584 0.39445 0.40451 0.40257 0.19963

1.02553 1.02040 0.97611 0.99236 2.03032

0.68135 0.66970 0.69291 0.67445 0.38268

1.23300 1.26260 0.74697 0.87071 0.91132

4.053  101 3.159  102 1.054  102 6.279  102 2.285  101

VCORR ¼ VPR þ C where

0 200

250

300

350

400

450

500

550

600

T[K] Fig. 4. Enthalpy calculation using the PengeRobinson EOS for n-octane at P ¼ 1 bar.

equation of state and reference ideal enthalpy from equation (29). The transition from liquid to vapor at the considered pressure occurs at 398.4 K for n-octane and at 374.2 K for water. Reference values are reported in Table 1. 4.4. Volume correction The PengeRobinson equation of state produces very accurate results dealing with hydrocarbons components but fails to generate satisfactory density values for the liquid phase for some substances, such as polar substances. In this work a Gaussian-like volume correction as developed by Monnery et al. [26] has been considered. Using the volume translation approach, equation (33) is used to calculate the temperature-dependent parameter a(T) of the PREOS instead of equation (21)



aðTr ; uÞ ¼

  pffiffiffiffiffi pffiffiffiffiffi2 2 1 þ m 1  Tr  n 1  Tr

(32)

The liquid volume is then corrected using a combination of a constant term plus a temperature-dependent Gaussian-like enhancement, which is only significant in the critical region

3.5 10

6

3 10

6

2.5 10

6

2 10

6

1.5 10

6

(33)

"  #  C4 Tr  C3 2 C ¼ C1 þ pffiffiffiffiffiffi exp  0:5 C2 2pC2

(34)

This model requires 6 parameters: m, n, C1, C2, C3 and C4. According to Ref. [26] C4 is a dimensional constant set to 1 m3/ kmol; some values of the other 5 parameters are listed in Table 2. Figs. 6 and 7 show, respectively, density calculations for n-octane and for water. The latter has been predicted using the volume shift correction. 4.5. Homogeneous two-phase speed of sound In thermodynamic equilibrium the formula proposed by Wallis [27] is commonly used for the calculation of the speed of sound in a two-phase flow

1

ra2

¼

a 1a þ rv a2v rl a2l

(35)

This expression strictly holds only for fixed void fraction. In the present work equation (35) is used to obtain a first reference value of the sound speed. Then the two-phase flow speed of sound is obtained by taking the derivative of equation (5) with respect to the pressure, as proposed by Catania et al. [28] for a vaporizing liquid

  a 1 1a 1 1 dY r $ ¼ þ   rv rl dp rv a2v ra2 rl a2l

(36)

800 700

Density [kg/m3]

Enthalpy [J/kg]

600

Ideal enthalpy PR-EOS enthalpy

500 400 Saturated vapour density n-octane Saturated liquid density n-octane

300 1 10 6 5 10

200

5

0 250

100

300

350

400

450

500

550

600

650

T[K] Fig. 5. Enthalpy calculation using the PengeRobinson EOS for water at P ¼ 1 bar.

0 200

250

300

350

400 450 T [K]

500

550 Tc

600

Fig. 6. Saturated densities calculated for n-octane without volume correction.

1466

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

1200

Density [kg/m3]

1000

800

600

400

200

0 200

Sat. vapour density water Sat. liquid density water corrected Sat. liquid density water

300

400

500 T [K]

600

700 Tc

Fig. 7. Saturated densities calculated for water considering the volume shift correction.

The third term in the right side of equation (36) takes into account the contribution of the evaporation process. In Le Métayer et al. [29] a similar derivation can also be found. The knowledge of the derivative dY/dp is required for the calculation of the speed of sound a from equation (36). At this purpose, it is assumed that during the vaporization process the local pressure is constant and is determined by the temperature (vapor pressure), thus rv and rl are constant, and the derivative of Y can be related to the mixture density variation as

rv rl dY dr ¼ dp ðrv  rl Þr2 dp

(37)

The term dr/dp is evaluated at constant vapor fraction and set equal to the inverse of the speed of sound calculated using equation (35). The third term of equation (36) is calculated using updated values of the densities. According to the present model, the phase change process, thus the variation of Y with p, makes the speed of sound tending to zero. This method has been considered as one of the simplest ways for capturing macroscopical effects of the vaporizing fluid. 4.6. Limits of stability of the liquid phase The PR-EOS can calculate pressures and densities for superheated liquid states as long as the mechanical stability criterion is respected, which is (dp/dr)T > 0. For the liquid phase this condition is met by all points above the liquid spinodal limit and below the saturation pressure, in the liquid branch (Fig. 3). Those states represent systems in meta-stable equilibrium. The PR-EOS was also found to be one of the few equations of state able to perform good results in spinodal limits calculation. Other conclusions can be found in Ref. [30]. 5. External flashing with surface evaporation These kinds of jets are usually composed of a liquid core surrounded by a two-phase flow and present compressible behavior associated with the liquid flashing outside the nozzle. The experimental investigation conducted by Vieira and Simoes-Moreira [8] demonstrated the absence of internal nucleation, showing photographic evidences of the regular and axisymmetrical shape of the liquid jet. According to them, flashing takes place on the surface of the liquid core through an ‘evaporation wave’ process.

The wave structure is composed by a rarefaction wave, propagating through the liquid that becomes metastable, followed by a subsonic phase transition front, propagating through the liquid as well, originating a high velocity liquidevapor mixture. SimoesMoreira and Shepherd [9] experimentally measured the velocity of the evaporation front in superheated dodecane, which was found to be in between 0.2 and 1.6 m/s. The velocity at which the mixture is ejected in the low pressure chamber is instead of the order of 100 m/s [29]. The present external-flashing model aims to find an approximated kinetic closure of the evaporation problem, which can be applied on the surface of the undisturbed superheated liquid jet, assuming a jet structure like the one described by Simoes-Moreira and Shepherd [9] and Vieira and SimoesMoreira [8]. As done in Refs. [8,9,31], the evaporation wave is seen as a discontinuity and is treated in the same way as a flame front in deflagrations, using a jump condition to link the upstream superheated liquid state to the vapor phase outside the nozzle. In this analogy the ‘educts’ of the reaction are represented by the superheated liquid in a meta-stable state and the ‘products’ are represented by the vapor phase. The evaporation rate is the ‘rate of reaction’ which functions as a kinetic closure of the problem. The reaction energy of the chemical bonds is here replaced by the energy stored in the meta-stable liquid state. In order to be consistent with the assumption of adiabatic phase transition the fluid must be retrograde [32]. The retrograde behavior of a fluid is determined by its molar complexity whose high degree of freedom is responsible for its large specific heat. Thus the higher the molar complexity, the higher is the chance to have a complete adiabatic evaporation of the liquid. To calculate the jump in thermodynamic properties across the meta-stable liquid, the conservations laws of mass, momentum and energy must be solved through a control volume that includes the evaporation wave. Referring to Fig. 8 it can be written (brackets refer to a ‘jump’ between two states)

½rw ¼ 0 h

(38)

i

rw2 þ P ¼ 0

(39)

w2 ¼ 0 hþ 2

(40)

w is the relative velocity between the front velocity u and liquid or vapor velocity, respectively u1 and u2, so then

w1 ¼ u þ u1

(41)

w2 ¼ u þ u2

(42)

Combining equations (38) and (39) we obtain the Rayleigh equation

u2-vapor

Vapour (products)

u1-liquid u

Superheated Liquid (reactants)

control volume Fig. 8. To calculate the jump in thermodynamic properties across the meta-stable liquid, the conservations laws of mass, momentum and energy must be solved through a control volume that includes the evaporation wave.

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

½P ðrwÞ2 ¼  ½v

(43)

where v express the specific volume, while combining equations (38) and (40) we obtain the Rankine-Hugoniot equation

½h ¼

v1 þ v2 ½P 2

(44)

where v1 and v2 are respectively the specific volumes of the liquid and of the vapor phase. Referring to Fig. 9, assuming that state 1 is known, the calculation of state 2 will be the target. The difficulty is due to the fact that ‘products’ and ‘educts’ of the reaction reside in disconnected part of the diagram and their thermodynamic states cannot be linked by an equation of state. Between the two states there is a reaction. A possible solution, as proposed by many authors, is to close the system with a kinetic relation in order to find the finite ‘rate of reaction’ of the transformation. This solution is known as the lower Chapman-Jouguet condition and determines a limit to the propagation speed of the discontinuity front. Equation (44) is the evaporation adiabatic curve locus of downstream states while equation (43) links the initial state (p1,v1), in the meta-stable liquid region, to one of the possible states of the RH curve. Equation (43) gives the slope of the Rayleigh line. The correct location of the downstream state 2 depends on the Rayleigh line slope, which is determined by the ‘lower’ CJ-point condition

dðrwÞ2 ¼ 0

v2 zvV >> vL The Rayleigh equation then becomes



w2 v2

2 z

½P v2

(46)

Furthermore at the CJ-point

w2 ¼ Cs ¼

P1 1þg

(48)

with the specific heats ratio g varying in the range from 1.0 to 1.4. Considering now a liquid volume issuing the nozzle with a velocity U, during an infinitesimal time interval dt, schematically represented in Fig. 10, it is reasonable to assume that this volume is cylindrical. In this work an axi-symmetric evaporation structure is assumed as well. Using this geometry we can locate the injection pressure P0 and the backpressure P3. After a first expansion process of the liquid phase, that takes place inside the nozzle [9], the liquid will be in a meta-stable state at the nozzle exit section (Fig. 9, point 1). Then an evaporation wave takes place on the surface of the unbroken liquid jet and propagates within the liquid phase, which undergoes a phase transition [9,29]. This state is calculated by imposing the CJ-point conditions. Point 2 can be either located inside or outside the saturation dome. The one- or two-phase flow originated becomes sonic and a shock wave can be formed to match the backpressure P3 [29]. Modeling the latter process is not the target of the work. The pressure P3 is not influent on the driving pressure drop. After steady conditions are reached, the surface evaporation rate outside the nozzle is driven by the pressures P1 and P2, before and after the evaporation front. These two pressures are the problem’s unknowns. 5.1. Meta-stable pressure and degree of metastability

(45)

This solution corresponds to the tangent point to the RH curve, is a sonic condition and gives the maximum slope of the Rayleigh line, thus it is also a maximum mass flux condition. A simplified formulation for the calculation of the CJ-point has been proposed by Simoes-Moreira and Shepherd [9] for ideal gas sufficiently far from the critical point. Assuming that the specific volume of the liquid phase is negligible compared to the specific volume of the vapor, we can write

ðrwÞ2 ¼

P2 ¼

1467

pffiffiffiffiffiffiffiffiffiffiffi gRT2

(47)

Thus substituting into (46)

Vieira and Simoes-Moreira [8] suggested an experimental way to find the meta-stable pressure P1 indirectly, measuring the mass flow rate and applying the orifice equation

P1 ¼ P0 

_ 2v m 2A2 CD2

(49)

where A is the orifice cross section and CD is the coefficient of discharge. Numerically the meta-stable pressure at the outlet section of the nozzle cannot be calculated using the nozzle governing equations (1)e(4), because the pressure is locked during choking conditions and the orifice cannot reach any imposed outlet pressure. The actual pressure on the outlet section is determined, but is not equal to, by the saturation pressure at the local superheated liquid temperature. The relation between the meta-stable pressure and the saturation pressure lies in the definition of the degree of metastability U

P ¼

ðPs  P1 Þ Ps

(50)

where Ps is the saturation pressure at the inlet liquid temperature and P1 is the meta-stable pressure at the nozzle exit section. The meta-stable pressure can be calculated as

P

P1 ¼ ð1  PÞPs

Rankine-Hugoniot curve

The knowledge of the degree of metastability P is necessary in order to calculate the meta-stable pressure at the outlet section. In

Rayleigh line (negative slope) Lower CJ-point

1

(51)

P2

nozzle

2

cylindrical expansion

liquid P0

V Fig. 9. A Rayleigh line, a Rankine-Hugoniot curve and a cubical equation of state (dashed line) cross the phase boundary curve. The tangent point 2 defines a deflagration-like process because the specific volume increases from 1 to 2. It is also a sonic outflow condition and corresponds to the maximum admissible flow rate.

P1

U liquid jet P3

Fig. 10. The evaporation process is modeled as an axi-symmetric jet expansion.

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

the present work experimental values of the degree of metastability from Vieira and Simoes-Moreira [8] have been used. Under the assumption of absence of in-flow nucleation, equation (7) can be adjusted to fit the assumption of external flashing, setting a low constant value for a (i.e. 0.01) on the outlet section

P ¼ 6:43$102 Q0:446

(52)

The weak dependency of P on the void fraction is neglected (constant a), thus the degree of metastability becomes a decreasing function of the relaxation time only.

12

10 Mass Flow Rate [g/s]

1468

8

Psat=1.98

6

4

5.2. Velocity of the evaporation front

Experimental

2

A method is here proposed to formulate an expression for the front velocity within the liquid. Equation (47) states that the evaporation front velocity with respect to the vapor flow velocity at CJ-point equals the local speed of sound. The evaporation front velocity with respect to the liquid is expected to be one to two orders of magnitude smaller than the velocity at which the vapor is ‘ejected’ from the liquid surface [29]. In the present formulation the variation with density of the pressure is considered the main driving factor and is expressed by differences, before and after the front. The velocity of the evaporation front determines the extinction jet length, which has an exponential dependency on the saturation to injection pressure ratio Ps/P0 that has been experimentally observed [8]. The degree of metastability is considered as an inverse measure of the ‘delay of evaporation’ (equation (52)) and it represents the ‘reactivity’ of the reaction. This dependency has been taken into account in the calculation of the evaporation front velocity as follow:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  u
P u met  PCJ  Vf ¼ xtP$  rl  rg

(53)

where Pmet ¼ P1 is the pressure before the evaporation front and PCJ ¼ P2 is the pressure after the evaporation front. x is a nondimensional constant that has been set equal to 0.025. The front velocity can be directly calculated once P2 is obtained.

4-Flash 0 0

P0

Vf

5

6

diameter is evaluated during a time step dt. Assuming that the evaporation front velocity Vf and the issuing velocity U are constant and perpendicular in direction, every liquid cylinder will undergo a reduction in diameter every dt. The result is a discretization of a conical shaped liquid jet using cylindrical elements (U*dt) long. Thus, summing up all cylindrical length contributions, the extinction length of the outside jet is calculated. 6. Results and discussion 6.1. Internal flashing In order to test the HRM models implemented in the 4-Flash code, the test case proposed by Rossmeissl and Wirth [11] has been considered. This comparison offers the possibility to check the models capability in capturing the mass flow rate choking and the influence of the superheating degree on flow choking.

j+1

j+n-1

P1

P3

10 Mass Flow Rate [g/s]

cylindrical evaporation P2 j

3 4 Outlet Pressure [bar]

12

The non-dimensional intact jet length ratio LE/D, where D is the nozzle diameter, is an indirect measurement of the evaporation rate which is taking place on the surface of the jet itself. Computationally, referring to Fig. 11, after P2 and P1 are calculated respectively from equations (48) and (51), the evaporation front velocity Vf is determined (53) and a decreasing in the liquid cylinder

nozzle

2

Fig. 12. Comparison between experimental from Rossmeissl and Wirth [11] and numerical mass flow at different discharge pressure obtained for a temperature of 413 K.

5.3. Non-dimensional extinction length

liquid

1

Psat=3.61 8

6

4 j+n U

liquid jet

Experimental 4-Flash

2

Vf Dx=Udt

0 Le Fig. 11. In modeling the evaporating liquid jet the conical liquid core is approximated with a number of cylinders. This approach allows a simplification in the calculation of the vapor front velocity since the direction is perpendicular to the liquid jet and the evaporation front velocity Vf is independent from the jet flow velocity U.

0

1

2

3 4 Outlet Pressure [bar]

5

6

Fig. 13. Comparison between experimental from Rossmeissl and Wirth [11] and numerical mass flow at different discharge pressure obtained for a temperature of 393 K.

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

40

40 35

4-Flash Vieira and Simoes-Moreira's experiment

35

4-Flash Vieira and Simoes-Moreira's experiment

30

30

25

25 Le/D

Le/D

1469

20

20

15

15

10

10

5

5 0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

Ps/Po

Ps/Po Fig. 14. Experimental from Vieira and Simoes-Moreira [8] and calculated (4-Flash code) liquid core extinction length for injection pressure Po ¼ 122 kPa.

According to our model, vaporization of the liquid within the nozzle is achieved if: - The local pressure falls below the vapor pressure at liquid temperature - The superheat of the liquid triggers a nucleation process - The relaxation time is lower than the residence time of the liquid within the nozzle The test consists of a pressurized liquid vessel containing deionized water, kept at constant pressure of 7 bar, and, at a determined temperature, the water was issued through a sharp-edged cylindrical nozzle having length-to-diameter ratio L/D ¼ 14 and D ¼ 0.7 mm. A uniform grid spacing of 0.1 mm has been adopted in simulations. Two different superheated levels have been considered, corresponding to the injection temperature 393 K and 413 K. The outlet pressure was varied between 5 and 1 bar and the mass flow rate was measured. Calculations have been compared to experimental results to establish model’s potential to detect choking conditions. The corresponding vapor pressures for water at the considered temperatures are, respectively, 1.98 bar and

Fig. 16. Experimental from Vieira and Simoes-Moreira [8] and calculated (4-Flash code) liquid core extinction length for injection pressure Po ¼ 500 kPa.

3.61 bar. Figs. 12 and 13 show the comparison between experimental and the predicted flow rates for different downstream pressures, considering the two different superheating levels. In both cases, the 4-Flash code more than reasonably reproduces the overall experimental mass flow trends. Calculations, as well as experiments, show that increasing the inlet temperature a reduction of the mass flow rate becomes evident. The higher vapor pressure related to the higher temperature affects the evaporation rate, which increases as well. The increasing of the evaporation rate drastically reduces the speed of sound in the two-phase flow, creating choking conditions within the nozzle. Once choking conditions are reached the mass flow rate becomes constant for further decreases in the outlet pressure. 6.2. External flashing An experimental investigation conducted by Vieira and SimoesMoreira [8] has been taken into consideration as a reference in order to test the numerical model versus external-flashing conditions. A set of simulations has been run using superheated iso-octane and working parameters proposed in their main experimental test. A conical converging nozzle with an exit section

40

50 Vieira and Simoes Moreira Experiment

4-Flash Vieira and Simoes-Moreira's experiment

35

4-Flash

40 30 30 Le/D

Le/D

25 20

20

15 10

10 5 0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Ps/Po

Ps/Po

Fig. 15. Experimental from Vieira and Simoes-Moreira [8] and calculated (4-Flash code) liquid core extinction length for injection pressure Po ¼ 250 kPa.

Fig. 17. Experimental from Vieira and Simoes-Moreira [8] and calculated (4-Flash code) liquid core extinction length for injection pressure Po ¼ 750 kPa.

1470

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

Table 3 Comparison between experimental [8] and calculated (4-Flash code) parameters for the external-flashing test. Po  1.8 (kPa)

Ps/Po (e)

Cd  0.018 (e)

P (e)

LE exp  0.5 (mm)

LE calc. (mm)

122.10 122.70 123.50 249.00 251.80 249.60 257.10 501.60 497.60 503.40 503.30 500.00 502.10 750.40 752.50 748.40 753.20 749.60 749.20 751.60

0.20 0.40 0.73 0.10 0.20 0.71 0.89 0.05 0.10 0.18 0.36 0.46 0.89 0.04 0.07 0.12 0.24 0.31 0.60 0.81

0.913 0.917 0.920 0.922 0.925 0.928 0.926 0.929 0.932 0.934 0.936 0.936 0.933 0.933 0.936 0.937 0.939 0.939 0.939 0.936

0.832 0.865 0.836 0.631 0.641 0.470 0.387 0.369 0.408 0.412 0.265 0.240 0.236 0.183 0.026 0.012 0.113 0.045 0.083 0.067

6.23 4.22 2.33 7.44 6.33 3.12 2.91 11.56 7.84 5.73 3.92 3.62 2.11 14.47 10.35 7.64 5.03 4.42 2.31 1.71

6.03 4.53 2.92 7.63 5.17 2.16 1.72 11.94 8.29 6.16 4.85 3.42 1.65 15.30 11.70* 9.06* 6.37 5.39* 3.22 2.14

of 0.31  0.01 mm, an inlet diameter of 3 mm, with rounded entrance, and length of 8 mm has been used in the experimental work and has been reproduced for calculation in the present 1D study. The roundness of the inlet geometry has been modeled using discharge coefficients experimentally provided by Ref. [8]. Numerically, in order to preserve the conservative laws in the convergent nozzle, the total variation diminishing flux scheme (TVD) was disabled, and the grid spacing was reduced to 0.016 mm. The purpose of the test was to measure the intact jet length when the jet undergoes flash evaporation outside the nozzle. Figs. 14e17 show the extinction liquid core calculated with 4-Flash code and compared with experimental data extracted from the work of Vieira and Simoes-Moreira [8]. The liquid core extinction length has been plotted versus the saturated-to-injection pressure ratio Ps/Po. Results have been collected in test series at constant pressure and different temperatures. The saturated pressure Ps is related to the injection temperature To.

60

LE calculated [mm]

40

30

20

10

0 0

10

20

7. Conclusion An analytical and numerical 1D-model for flashing liquid jets is presented. The model aims to provide boundary conditions for 3D CFD spray simulations, when the effects of the degree of superheat need to be considered in the mixture formation. Internal- and external-flashing mechanisms are modeled in the present work. The coupling of a HRM based evaporation model with a VTPR-EOS has been tested in internal flashing, detecting choking conditions inside the orifice. In external-flashing behavior, when the jet issues the orifice in a superheated unbroken liquid phase, thermodynamic analysis of the CJ-point has been used to get to an approximated kinetic closure of the system and a solution for the calculation of the extinction jet length has been proposed. The main goal of building a first step of a simple and versatile 1D two-phase homogeneous framework, suitable for many different fuel injection conditions, has been obtained with encouraging results. Acknowledgment The research is supported by University of Bologna. The first author thanks Graeme Supeene and Prof. C.P.T. Groth, UTIAS University of Toronto, for the invaluable exchange of material and Graeme’s useful suggestions for the implementation of the PR-EOS equation of state. Nomenclature

Po=122 kPa Po=250 kPa Po=500 kPa Po=750 kPa

50

Figs. 14e17 shows respectively test cases at Po ¼ 122 kPa, Po ¼ 250 kPa, Po ¼ 500 kPa and Po ¼ 750 kPa. The extinction jet length decreases with the increasing of the Ps/Po ratio. The numerical model follows the experimental data with good agreement, showing the exponential dependency of LE with the increasing Ps/Po. Table 3 displays test data: the injection pressure Po, saturated-to-injection pressure ratio Ps/Po, the coefficient of discharge Cd, the experimental degree of metastability P, the measured extinction length from the work of Vieira and SimoesMoreira [8], and, in the last column, the calculated extinction length. Runs with degree of metastability lower than 0.06 have been low bounded with a value set to 0.06, according to the arbitrary limit of equation (52). The corresponding intact jet lengths are reported with a (*) in Table 3. All numerical results and experimental data are plotted in Fig. 18.

30

40

50

60

LE measured [mm] Fig. 18. Comparison between experimental, from Vieira and Simoes-Moreira [8], and calculated liquid core extinction length.

a C(W) Cc Cd D eo f F G h H k L p P PC r R

speed of sound [m/s] source terms vector contraction coefficient discharge coefficient nozzle diameter [mm] total specific internal energy (e þ 0.5u2) fugacity [Pa] fluxes term vector Gibbs free energy [J/mol] molar enthalpy [J/mol] specific enthalpy per mass unit [J/kg] ratio of specific heat capacities nozzle length [mm] mixture pressure [Pa] pressure [Pa] critical pressure [Pa] fillet radius [m] universal gas constant

S. Negro, G.M. Bianchi / International Journal of Thermal Sciences 50 (2011) 1460e1471

t T Tr u v V W x Y

time [s] temperature [K] reduced temperature T/TC mixture velocity [m/s] specific volume [m3/kg] volume [m3] conserved variable vector axial coordinate [m] vapor mass fraction Y ¼ arv/r [e]

Greek letters r  rL void fraction [e], a ¼ rv  rL fugacity coefficient vapor mass fraction transport equation source term acentric factor relaxation time [s] superheating degree [K] mixture density [kg/m3]

a F G u Q Dq r

Subscript 0 e l s, sat v

initial value equilibrium liquid saturation condition vapor

References [1] G. Supeene, Numerical Simulation of Aviation Fuel Sprays at Elevated Temperatures, Report No. DEC3. UTIAS University of Toronto, 2009. [2] L. Ng, Modelling quasi one-dimensional flow of supercritical pre-heated jet fuel, BASc thesis, University of Toronto, 2006. [3] J. Lee, R. Madabhushi, C. Fotache, S. Gopalakrishnan, D. Schmidt, Flashing flow of superheated jet fuel, Proceedings of the Combustion Institute 32 (2) (2009) 3215e3222. [4] R.D. Reitz, A photographic study of flash-boiling atomization, Aerosol Science and Technology 12 (1990) 561e569. [5] B.S. Park, S.Y. Lee, An experimental investigation of the flash atomization mechanism, Atomization and Sprays 4 (1994) 159e179. [6] R.D. Oza, J.F. Sinnamon, An Experimental and Analytical Study of Flash-Boiling Fuel Injection (1984) SAE paper 830590. [7] D. Yildiz, et al., A Study on the Dynamic of a Flashing Jet Final Contract Research Report EAR0030/2002. von Karman Institute for Fluid Dynamics, 2002. [8] M.M. Vieira, J.R. Simoes-Moreira, Low-pressure flashing mechanisms in isooctane liquid jets, Journal of Fluid Mechanics 572 (2007) 121e144. [9] J.R. Simoes-Moreira, J.E. Shepherd, Evaporation waves in superheated dodecane, Journal of Fluid Mechanics 382 (1999) 63e86.

1471

[10] D.Y. Peng, D.B. Robinson, A new two-constant equation of state, Industrial and Engineering Chemistry Fundamentals 15 (1) (1976) 59e64. [11] M. Rossmeissl, K.E. Wirth, Critical mass-flow in orifice-nozzles at the disintegration of superheated liquids, in: Spray 2006 Workshop über Sprays, Erfassung von Sprühvor-gängen und Techniken der Fluidzerstäubungy. [12] P. Downar-Zapolski, Z. Bilicki, The non-equilibrium relaxation model for onedimensional flashing liquid flow, International Journal of Multiphase Flow 3 (22) (1996) 473e483. [13] G.M. Bianchi, C. Forte, S. Negro, P. Pelloni, A 1D model for the prediction of flash atomization in GDI multi-hole injectors: preliminary results, SAE International Journal of Engines 1 (1) (2008) 1278e1293. [14] G.M. Bianchi, S. Negro, C. Forte, G. Cazzoli, P. Pelloni, The Prediction of Flash Atomization in GDI Fuel Injectors (2009) SAE paper 2009-01-1501. [15] J. Senda, Y. Hojyo, H. Fujimoto, Modelling of Atomization Process in Flash Boiling Spray (1994) SAE Paper 941925. [16] D. Kawano, H. Ishii, H. Suzuki, Y. Goto, M. Odaka, J. Senda, Numerical study on flash-boiling spray of multicomponent fuel, Heat Transfer e Asian Research 35 (5) (2006) 369e385. [17] C. von Kuensberg Sarre, S.C. Kong, R.D. Reitz, Modeling the Effects of Injector Nozzle Geometry on Diesel Sprays (1999) SAE paper 1999-01-0912. [18] S.F. Davis, A simplified TVD finite difference scheme via artificial viscosity, SIAM Journal on Scientific and Statistical Computing 8 (1987) 1e18. [19] D.E. Winterbone, R.J. Pearson, Design Techniques for Engine Manifolds: Wave Action Methods for IC Engine. Wiley, 1999. [20] P. Ghosh, Prediction of vaporeliquid equilibria using PengeRobinson and Soave-Redlich-Kwong equations of state, Chemical Engineering & Technology 22 (5) (1999) 379e399. [21] J.O. Valderrama, The state of the cubic equations of state, Industrial & Engineering Chemistry Research 42 (2003) 1603e1618. [22] T.F. Irvine Jr., P.E. Liley, Steam and Gas Tables with Computer Equations. Academic Press, New York, 1984. [23] N.B. Vargaftik, Tables on the Thermophysical Properties of Liquids and Gases. Hemisphere Pub. Corp., Halsted Press, New York, 1975. [24] NIST, Scientific and Technical Databases, NSRDS, National Institute of Standards and Technology. [25] J. McBride, S. Gordon, M.A. Reno, Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species (1993) NASA Technical Memorandum 4513. [26] W.D. Monnery, W.Y. Svrcek, M.A. Satyro, Gaussian-like volume shifts for the PengeRobinson equation of state, Industrial & Engineering Chemistry Research 37 (1998) 1663e1672. [27] W. Wallis, One Dimensional Two-phase Flow. McGraw Hill, 1975. [28] A.E. Catania, A. Ferrari, M. Manno, E. Spessa, Comprehensive thermodynamic approach to acoustic cavitation simulation in high-pressure injection systems by a conservative homogeneous two-phase barotropic flow model, ASME Journal of Engineering for Gas Turbines and Power 128 (2) (2006) 434e445. [29] O. Le Métayer, J. Massoni, R. Saurel, Modelling evaporation fronts with reactive Riemann solvers, Journal of Computational Physics 205 (2005) 567e610. [30] J. Wisniak, K. Galon, Supersaturation, the spinodals, and their prediction using equations of state, Physics & Chemistry of Liquids 38 (2000) 643e661. [31] J.R. Simoes-Moreira, M.M. Vieira, E. Angelo, Highly expanded flashing liquid jets, Journal of Thermophysics and Heat Transfer 16 (3) (2002) 415e424. [32] T. Kurschat, H. Chaves, G.E.A. Meier, Complete adiabatic evaporation of highly superheated liquid jets, Journal of Fluid Mechanics 236 (1992) 43e59.