Evaluation of commonly used assumptions for isolated and cluster heptane drops in nitrogen at all pressures

Evaluation of Commonly Used Assumptions for Isolated and Cluster Heptane Drops in Nitrogen At All Pressures K. HARSTAD and J. BELLAN*

Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109-8099, USA A study is performed to assess commonly used assumptions in the modeling of drop behavior in moderate to high temperature surroundings and at all pressures. The model employed for this evaluation has been previously validated for isolated drops by using microgravity data, and is very general: it contains Soret and Dufour effects, does not assume mass transfer quasi-steadiness at the drop boundary, or necessarily the existence of a drop surface (i.e., phase discontinuity). Moreover, the numerical simulations are performed with accurate equations of state and transport properties over a wide range of thermodynamic variables. Consistent with low pressure conditions, the drop boundary is identified a posteriori of the calculations with the location of the largest density change. Simulations are here performed for isolated drops, and for monodisperse as well as binary size drop clusters. The results show that at locations arbitrarily near the boundary, the drop does not reach the mixture critical point within the wide range of conditions investigated (far-field temperatures of 470 –1000 K and pressures ranging from 0.1 to 5 MPa). However, the state arbitrarily near the boundary is closer to the critical condition for smaller drops in a cluster than for the larger drops. Evaluations of the effect of the relaxation time at the drop boundary show that quasi-steadiness of the mass transfer prevails for drops of radius as small as 2 ⫻ 10⫺3 cm. Finally, the diameter squared exhibits a linear time variation only at atmospheric pressure. At all other pressures investigated (1–5 MPa), the diameter squared displays a negative curvature with time which never becomes linear. In agreement with existing experimental data, the drop lifetime increases monotonically with pressure at low far field temperatures (470 K), but exhibits a maximum as a function of pressure at high temperatures (1000 K). On an appropriate scale, the slope of the diameter squared versus time is shown to be independent of the drop size at all pressures. © 2001 by The Combustion Institute

NOMENCLATURE A AJ Aq BJ Bq CJ Cp Cq d D Fems G h ជJ K m M n N Nu p

area coefficient in the transport matrix coefficient in the transport matrix coefficient in the transport matrix coefficient in the transport matrix coefficient in the transport matrix heat capacity at constant pressure coefficient in the transport matrix drop diameter mass diffusion coefficient mass emission rate Gibbs function molar enthalpy molar flux d(d 2 )/dt molar mass mass molar density number of drops Nusselt number pressure

*Corresponding author: E-mail: [email protected] COMBUSTION AND FLAME 127:1861–1879 (2001) © 2001 by The Combustion Institute Published by Elsevier Science Inc.

qជ r R Ru t T u v V X Y Z

heat flux radial coordinate specific radial position universal gas constant time temperature radial velocity molar volume volume molar fraction mass fraction compression factor

Greek Symbols

␣BK ␣D ␣h ␣IK ␣v ␥ ␩ ␬s

Bearman–Kirkwood form of the thermal diffusion factor mass diffusion factor defined by Eq. 13 Irwing–Kirkwood form of the thermal diffusion factor thermal expansion ratio activity coefficient viscosity isentropic compressibility 0010-2180/01/$–see front matter PII S0010-2180(01)00292-9

1862

␭ ␮ ␳ ␶ ␾ ⌽v ␸

K. HARSTAD AND J. BELLAN thermal conductivity chemical potential density stress tensor C7H16/N2 mass ratio viscous dissipation fugacity coefficient

Subscripts c C d e j r s ␣ 1

critical cluster drop external species reduced value; radial direction ‘optical’ boundary drop size class heptane

Superscripts 0 i si

initial condition, pure species interstitial sphere of influence

INTRODUCTION The study of drop behavior in high-pressure environments has been a topic of interest for many decades owing to its relevance to costly combustion devices such as gas turbine engines and liquid rocket engines. Because high pressure tests are expensive and sometimes dangerous, considerable effort has been devoted to the development of accurate models capable of portraying the physics of drop evolution at all pressures. A recent review by Bellan [1] tabulates current drop models published over more than a decade. Although most models are based on the transient conservation equations augmented by the real gas equation of state, and thermodynamic variable dependent transport properties, the overwhelming majority of these models contains simplifying assumptions meant to enable less computationally intensive simulations. Such common assumptions are the existence of a drop surface assumed to be at an existing critical locus (e.g., [2–5]) or assumed to be at an existing location of quasi-steadiness of the mass flux, sometimes called thermodynamic

equilibrium (e.g., [2, 5–12, 14]), and the neglect of the thermal diffusion effects represented by the Soret and Dufour term for the species and energy equations, respectively. Predictions from such models include, for example, the validity of the diameter squared, d 2 , linear variation with time even at high pressures (e.g., Yang et al. [3] for the LOX/H2 system); if such results were correct, it would enable a considerable simplification in the incorporation of drop models in complex Computational Fluid Dynamics (CFD) codes. Therefore, there is an incentive to evaluate the validity of assumptions commonly used in high pressure drop models so as to determine possible simplifications that may result in large reductions of computational costs. Harstad and Bellan [15] developed a very general fluid drop model, and validated it with the microgravity data of Nomura et al. [16] for the heptane/ nitrogen system before utilizing it to study the validity of some hypotheses listed above. For example, the recent study of Harstad and Bellan [17] addressed the correctness of the linear d 2 (t) variation, where t is the time, and the accuracy of the neglect of thermal diffusion effects in the context of the LOX/H2 system. Isolated drop results displayed a nearly linear d 2 (t) variation in the low (but higher than atmospheric) pressure subcritical regime, however, with increasing pressure the results exhibited considerable departures from linearity. Constant and linear fits of the d 2 (t) slope, which at subcritical pressures is the evaporation constant, showed that the assumption of it being constant may be in substantial error. The linear fits displayed both increasing and decreasing trends, indicating that d 2 (t) may have either positive or negative curvatures, according to the pressure. Because for isolated C7H16 drops in N2 Harstad and Bellan [15] found consistently a negative curvature in agreement with Nomura et al. [16], the conclusion based on these limited investigations is that the curvature is species dependent. In the context of drop clusters, estimates of characteristic gradient terms (describing relative flux contributions) for two values of the thermal diffusion factor showed that the thermal diffusion factor influences mostly the Soret contribution which, although not negligible, was dominated both by the mass gradient term (Fick’s or Dufour) and by the Fourier

EVALUATION OF HEPTANE DROPS IN NITROGEN term. An increase by a factor of 2.5 in the thermal diffusion factor yielded only modest changes in the magnitudes of the Dufour and Soret contributions. The present heptane/nitrogen study extends the previous assumption evaluations in several ways. First, we systematically examine the timewise variation of d 2 for isolated drops as a function of pressure for both moderate and elevated surrounding-fluid temperatures; the explored thermodynamic regime approximately corresponds to that over which the model has been validated in a previous study [15]. Then, we evaluate the assumptions imposed on the drop boundary, that is the existence of a discontinuity which may be a surface, the possibility that the critical mixture point is reached arbitrarily near it, or that it may be in a quasi-steady mass flux state (i.e., the value of the flux of molecules leaving the boundary is very nearly equal to that of molecules arriving at the boundary). Furthermore, we compare the evolution of isolated drops to drops in clusters having the same initial radius and exposed to otherwise the same far-field conditions; we consider both monodisperse clusters and binary size clusters. In a previous LOX in H2 study [17], the focus was in comparing the evolution of considerably more numerous and smaller drops in a binary size cluster with that of same size drops in a monodisperse cluster. The conclusion was that whereas the predicted small-drop lifetime was somewhat larger for the drops in the monodisperse compared to the binary size cluster, the lifetime of the cluster was substantially smaller in the monodisperse simulation. The present study is in this respect complementary to our previous investigation because it focuses on the larger, less numerous drops in a cluster. Finally, we also discuss comparisons between small and large drops belonging to the same cluster, in particular addressing the question of which size class reaches conditions closer to the mixture critical point. CONSERVATION EQUATIONS The model equations have been previously developed by Harstad and Bellan for isolated drops [15], for monodisperse drop clusters [18],

1863

and more recently for polydisperse drop clusters [17]. To avoid undue repetition of published material, only the highlights and final conservation equations are given herein, and the reader is referred to those publications for details on these models. The fundamental base of the single drop model is Keizer’s fluctuation-dissipation (FD) theory [19] which formally accounts for nonequilibrium processes and furthermore, naturally relates fluxes and forces for a general fluid by providing the form of the transport matrix. In contrast, the kinetic theory of rarefied gases is frequently extended non-rigorously to obtain the form of the transport matrix for more general gases. According to FD theory, the general flux equations are given by additive forms whereby both the molar fluxes, Jជ␣, with ␣ denoting the species, and the heat flux, qជ , are the sums of coefficients multiplied by ⵜX ␤ with ␤ 僆 [1, N sp ], X being the molar fraction and N sp being the total number of species, a coefficient multiplied by ⵜT, where T is the temperature, and a coefficient multiplied by ⵜp, where p is the pressure. These matrix coefficients are functions of thermodynamic quantities and transport coefficients, as given by Harstad and Bellan [15]. Noteworthy, the general form of the transport matrix includes, additional to the well known Fick mass diffusion and Fourier heat conduction terms, the Soret term in the molar flux and the Dufour term in the heat flux. For a binary set of species, these new terms are proportional to a single new transport coefficient, the thermal diffusion factor. For a multispecies system, a set of thermal diffusion factors is involved. One issue in exercising this fundamental model is the knowledge of the value of the thermal diffusion coefficient(s) which are generally unknown. Experimentally determined values are scarce (see Hirshfelder et al. [20] and Chapman and Cowling [21]) even for low pressures, and measurements or calculations at high pressures are very scant. The simulations performed in this study all employ the thermal diffusion factor value determined by Harstad and Bellan [15] from comparisons with part of the data of Nomura et al. [16], and further validated by comparisons with the remaining part of that data.

1864

K. HARSTAD AND J. BELLAN

Governing Equations for an Isolated Drop The configuration modeled is that of a spherical, isolated drop in quiescent surroundings having the initial far-field boundary at a specified location, R si,0 . The far-field boundary location is a function of time and is followed in a Lagrangian manner. Conservation Equations According to Harstad and Bellan [15], the isolated drop conservation equations are: ●

Continuity

⭸ ␳ 1 ⭸共r ␳ u兲 ⫹ ⫽ 0. ⭸t r 2 ⭸r

(1)

⭸共 ␳ u兲 1 ⭸共r 2␳ uu兲 ⭸p ⭸ ␶ rr 3 ␶ rr ⫹ 2 ⫹ ⫽ ⫹ . ⭸t r ⭸r ⭸r ⭸r r ●

Species conservation

␳

DY j ⫽ ⫺m jⵜ 䡠 ជJ j. Dt

●

Energy equation

r

q

BJ Aq

CJ Bq

冊

冢冣

(5)

where ⫺J 1r and ⫺q r are the radial fluxes and A J ⫽ 共m/m 1兲nD ␣ D

(6)

B J ⫽ 共m 2/m兲nDX 1X 2共 ␣ IK ⫺ ␣ h兲/T

(7)

C J ⫽ 共m 2/m兲nD共m 1m 2X 1X 2/m兲共v 1/m 1

l

(9) (10)

(11) where m ⫽ X 1 m 1 ⫹ X 2 m 2 is the mixture molar mass, R u is the universal gas constant, ␭ is the thermal conductivity, ␣ IK is the Irwing–Kirkwood form of the thermal diffusion factor (see Harstad and Bellan [15]), and D is the mixture diffusivity. The binary mass diffusion factor can be calculated from thermodynamics

␣ D ⫽ 1 ⫹ X 1⭸共ln ␥ 1兲/⭸X 1

Nsp

l

A q ⫽ ␭ ⫹ 共 ␣ IK ⫺ ␣ h兲 ␣ IKR unDX 1X 2

(2)

(3)

冘 共h ⵜ 䡠 ជJ 兲.

(8)

B q ⫽ nD ␣ IK共m 1m 2X 1X 2/m兲共v 1/m 1 ⫺ v 2/m 2兲

Dp DT ⫽ ␣ vT ⫺ ⵜ 䡠 qជ ⫹ ⌽ v Dt Dt ⫹

J

C q ⫽ 关m 2/共m 1m 2兲兴nD ␣ D␣ IKR uT

Momentum conservation

nC p

1r

⫺ v 2/m 2兲/共R uT兲

2

●

冉 ⫺J⫺q 冊 ⫽ 冉 CA

⭸Y 1 ⭸r ⭸T ⭸r ⭸p ⭸r

(4)

l⫽1

In these equations, ␳ is the mass density, u is the radial velocity, r is the radial coordinate, ␶ rr ⫽ (4/3) ␩ [⭸u/⭸r ⫺ u/r] is the stress tensor, ␩ is the mixture viscosity, Y j is the mass fraction of species j, m j is the molar mass, n is the molar density, C p is the mixture molar heat capacity at constant pressure, h j ⫽ (⭸h/⭸X j ) p,T,X l⫽j is the partial molar enthalpy where h is the total molar enthalpy, ␣ v ⫽ [(⭸v/⭸T) p,X l]/v is the thermal expansion ratio with v being the molar volume, ⌽ v ⫽ (4/3) ␩ [⭸u/⭸r ⫺ u/r] 2 is the viscous dissipation and D/Dt ⬅ ⭸/⭸t ⫹ u(⭸/⭸r). For a binary species mixture (species subscripts 1 and 2 with m 1 ⬎ m 2 ), the fluxes can be written in compact form as related to the transport matrix

(12)

␸ 1 / ␸ o1

where ␥ 1 ⬅ is the activity coefficient and ␸ is the fugacity coefficient with the superscript o denoting the pure (X 1 ⫽ 1) limit. According to Harstad and Bellan [15], the Bearman–Kirkwood form of the thermal diffusion factor, ␣ BK , is related to ␣ IK by

␣ BK ⫽ ␣ IK ⫺ ␣ h,

␣ h ⬅ 共m 1m 2/m兲

䡠 共h 1/m 1 ⫺ h 2/m 2兲/共R uT兲.

(13)

These two factors correspond to the two forms of the heat flux vector (see Sarman and Evans [22]). In the above derivation (Eq. 9), it has been carefully ensured that the expression for the thermal conductivity converges to the kinetic value in the low pressure limit [15], and thus, that its value is indeed that tabulated or calculated through corresponding states procedures [23].

EVALUATION OF HEPTANE DROPS IN NITROGEN To close the system of equations, the real gas equation of state (EOS) presented by Harstad et al. [25] is coupled to the above system of differential equations. Boundary Conditions The above conservation equations pertain to the entire field, 0 ⱕ r ⱕ R si , where R si is the instantaneous location of the far field boundary which is initially specified for single drop calculations and is further calculated as part of the solution (see below). The far-field and drop boundary conditions are those derived in Harstad and Bellan [15], and the reader is referred to that publication for the detailed derivation. Only the highlights of the model are discussed below. Several aspects of these boundary conditions are noteworthy. As stated above, the value of R si is calculated as a function of time by following the initial far-field boundary in a Lagrangian way. At r ⫽ R si , the values of the thermodynamic variables are specified (they can be either constant or functions of time). At the drop boundary, the conservation statements are stated as jump conditions across the boundary, ⫹ relating locations within (R ⫺ d ) and outside (R d ) the drop; that is, the drop boundary is treated as a mathematical discontinuity. Within the context of this study relying on thermodynamic concepts, this is the only physically, and therefore mathematically sound approach since in ⫹ the layer between R ⫺ d and R d , usual thermodynamic concepts do not apply because the molecular forces in the layer are not homogeneous as the molecular distribution is inhomogeneous. The layer thickness is of the order of the molecular scale, and thus it is approximately four orders of magnitude smaller than that of the drop sizes used in the calculations presented below. This emphasizes the fact that one cannot assume the state of the drop boundary, unless corroborated by molecular dynamics calculations resolving the details of the layer, whereby ⫹ the conditions at R ⫺ d and R d are imposed, and the molecular dynamics obtained solution determines the structure of the boundary layer at R d . We distinguish between the Lagrangian location of the initial drop boundary, R 0d , which is followed in time, R d (t), and the location of

1865

the ‘optical’ drop boundary (not necessarily an interface), R s (t), which is also followed in time being calculated a posteriori of the solution as the location of the largest density change. The largest density change corresponds to either the location of the largest density gradient, or a density jump as it is an infinite gradient which is not calculable. This definition is relevant to optical observations which identify changes in the index of refraction. Of importance is the fact that we do not assume the existence of a material interface in the region 0 ⱕ r ⱕ R si . If such an interface exists at any time, it is totally consistent with our definition whereby the drop boundary is identified to be at the location of the largest density change, r ⫽ R s . Just as important, consistent with the discussion above noting that no assumption can be made regarding the state of the drop boundary, here we do not assume phase equilibrium at the drop boundary, or anywhere else in the computational domain, and we do not necessarily relate the critical point of the mixture to any of the radial locations in the flow. This very general formulation makes the model ideally suited for evaluating assumptions about the state of the drop arbitrarily near the boundary, that is about ⫹ the states at R ⫺ d and R d because within the formalism of the mathematical discontinuity, only the location of R d is known and the conditions at R d are not defined. Governing Equations for a Polydisperse Cluster The equations derived elsewhere [17] are simply recalled, with highlights, below. Cluster Conservation Equations A cluster of drops is defined by an ensemble of N drops partitioned into K size classes denoted by ␣, with N ␣ drops in each size class where N ⫽ ⫽K ¥ ␣␣ ⫽1 N ␣ . The cluster is assumed spherical, of volume V C , boundary area A C and radius R C . Consistent with the cluster sphericity, we assume that the interstitial region between drops is uniform and quiescent with respect to the cluster, therefore, addressing only average interstitial fluid properties. As in Harstad and Bellan [17], each cluster drop is surrounded by a

1866

K. HARSTAD AND J. BELLAN

fictitious ‘sphere of influence’ centered at the drop center and having a radius R si (superscript “si” refers here to the sphere of influence). The value of R si is one half of the mean distance between adjacent drops, or less, and is found from the condition that the spheres of influence must be tightly packed; in contrast with the isolated drop simulations where R si,0 is specified, here it is calculated from the specified initial mass loading. Each individual drop obeys the conservation equations listed above, but the far-field boundary conditions at the edge of each sphere of influence can no longer be specified, and instead are calculated from the solution of the conservation equations for the ⫽K interstitial cluster volume, V i ⫽ V C ⫺ ¥ ␣␣ ⫽1 N ␣ V ␣ (the superscript ‘i’ identifies the interstitial region), where V ␣ is the volume corresponding to the Lagrangian value of R si for each size class, R ␣ , and contains the drop and its surrounding fluid having by definition a fixed mass. The final form of these equations is for a total interstitial mass M i (see for details [17]): V i共t兲 ⫽ M i/ ␳ i

冉

dY ji m j ⫽ i J jeA C ⫺ dt M

冋

冘NJ

␣⫽K

␣ j,␣A ␣

␣⫽1

(14)

冊

(15)

冘 N 共q ␣

␣

册

˜ J 1,␣兲 A ␣ ⫺ m 1h

␣⫽1

p e dp e ␳ i dt

(16)

˜ V ⬅ C pi/m i ⫺ ␣ vip e/ ␳ i, C

(17)

˜ Vm i␣ viT i/共 p eC pi兲, ␬ ˜ s ⬅ ␬ is ⫺ C

(18)

⫺␬ ˜s

J 1e ⫽ A Ji共⭸Y 1/⭸r兲 r⫽e ⫹ B Ji共⭸T/⭸r兲 r⫽e ⫹ C Ji共⭸p/⭸r兲 r⫽e

with

˜ ⬅ 共h i1/m 1 ⫺ h i2/m 2兲 ⫺ p e共v i1/m 1 ⫺ v i2/m 2兲, h (19) where C pi , ␣ vi , ␬ is , h ji , and v ji are the molar heat capacity, the thermal expansion ratio, the isentropic compressibility, the partial molar enthalpy and the partial molar volume of species j,

(20)

q e ⫽ A qi共⭸T/⭸r兲 r⫽e ⫹ C qi共⭸Y 1/⭸r兲 r⫽e ⫹ B qi共⭸p/⭸r兲 r⫽e

(21)

where j ⫽ 1 refers to heptane, r is the direction normal to the cluster boundary enclosing V C and the coefficients A Ji , B Ji , C Ji , A qi , C qi , and B qi are calculated similarly to those given by Eqs. 6 through 11, except that the values of all variables are those in the interstitial region. Boundary Conditions Additional to the boundary conditions at the drop boundary and at the edge of the sphere of influence, here there is a set of boundary conditions which must be satisfied at the cluster boundary. These latter are calculated as 共⭸Y 1/⭸r兲 r⫽e ⫽ 共Y 1e ⫺ Y i1兲/r e, ⫽ 共T e ⫺ T i兲/r e

1 dT i ˜V ˜ J 1e兲 A C ⫺ C ⫽ i 共q e ⫺ m 1h dt M ␣⫽K

all in V i ; and the subscript ‘e’ identifies values of the variables at a far-field location external to the cluster. The values of J 1, ␣ and q ␣ are calculated at the edge of the sphere of influence for each size class and fluxes J 1e and q e are calculated at the cluster boundary:

共⭸T/⭸r兲 r⫽e (22)

where r e is an external scale which may be related to what is equivalent to a Nusselt number, Nu C , through r e ⫽ R C /Nu C , because it is the characteristic distance over which transport occurs at the cluster boundary. A parametric study of Nu C was presented by Harstad and Bellan [18] and [17] for monodisperse and polydisperse clusters, respectively, and it was shown to be a weak parameter in this diffusion dominated situation. TRANSPORT PROPERTIES The calculation of transport properties generally follows the procedure outlined in Harstad and Bellan [26]. However, following accepted practice [23, 27], the excess of thermal conductivity and viscosity of vapor-like mixture components over their temperature-dependent lowpressure behavior is calculated as functions of

EVALUATION OF HEPTANE DROPS IN NITROGEN reduced density alone instead of using the prior correlations in p and T. Also, binary diffusion coefficients for vapor-like mixture components are calculated using semi-empirical collision integral ratios [28]. NUMERICAL METHOD In all simulations presented below, the pressure gradient terms have been neglected in Eqs. 20 and 21 because the Mach number ⬍⬍ 1; a posteriori evaluations of these terms showed that they are indeed negligible. The numerical method for obtaining the solution for the variables in each sphere of influence was explained in detail elsewhere [15], and thus will not be repeated. In the interstitial region, the unknowns are T i , Y i1 , and ␳ i determined from Eqs. 15, 16, and the EOS. Once ␳ i is known, V i and V C are calculated. During each time step, the total differential Eqs. 15 and 16 are solved by using a second order predictorcorrector method.

1867

attaining the binary mixture critical point is that the mass diffusion factors, ␣ Dij , be null because for a multicomponent system

␣ Dij ⬅ 关1/共R uT兲兴关X i⭸ ␮ i/⭸X j兴 ⫽ 关1/共R uT兲兴 䡠 关X i共⭸ 2G/⭸X i⭸X j兲 p,T,Xk兴

(23)

where ␮ is the chemical potential. According to the Gibbs–Duhem relationship, for a binary mixture there is a single value of the mass diffusion factor, ␣ D ⫽ ␣ D11 ⫽ ␣ D22 ⫽ ⫺ ␣ D12 ⫽ ⫺ ␣ D21 . This observation is of significance because for a binary species system the single value of ␣ D can serve as a ‘diagnostic’ to determine whether there is a possibility that the critical point has been reached in the computational domain. This information can further support the evidence provided by the values of the reduced pressure, p r ⫽ p/p c , and reduced temperature, T r ⫽ T/T c , where the subscript ‘c’ denotes the critical point.

RESULTS MIXTURE CRITICAL POINT According to Peng and Robinson [29], the critical point occurs when both the determinant of ញ 兩 (i.e., U ញ ij ⫽ ⭸ 2 G/ the second derivatives, 兩U ⭸X i ⭸X j ), of the Gibbs function, G, and that of the second derivatives of G combined with the ញ 兩 are null. Therefore, the first derivatives of 兩U accuracy of mixture critical point calculations is highly dependent of the accuracy of the EOSs. To verify the reliability of the present mixture critical locus, comparisons of critical locus calculations were performed based upon both the Peng–Robinson definition [29] and the free energy method of Reid et al. [23], and they were further compared to data for two binary species systems. Thus, calculations for decane/CO2 were compared with the data of Reamer and Sage [30], and calculations for ethane/heptane were compared with the well known data of Kay [31]. For both binary mixtures, the computational results using the two methods agreed and duplicated the data with very good accuracy (not shown). A consequence of the definition of the critical point is the fact that a necessary condition for

Simulations are here performed for isolated drops and drops in clusters, for both monodisperse and binary size drops. All simulations were conducted for heptane (T c ⫽ 540.3 K, p c ⫽ 2.76 MPa) drops in nitrogen (T c ⫽ 126.2 K, p c ⫽ 3.39 MPa), and the initial and/or far-field boundary conditions are displayed in Tables 1 through 3. For isolated drops, the radius of the sphere of influence is initially specified; for monodisperse clusters, R si,0 is calculated from the heptane to nitrogen specified initial mass ratio, ␾0 and the given initial drop size; for binary size clusters, this radius is initially calculated from the heptane to nitrogen specified initial mass ratio, ␾0, the given initial drop sizes, and the imposed value of the ratio of the number of drops for the two drop sizes. The value of Nu C is also specified. In Table 1, the coefficients ␣ a1 and ␣ a2 represent the accommodation coefficients in the kinetic law governing the molecular motion to and from the drop boundary: F ems ⫽

冘 关␣ m u aj

j⫽1,2

G j Tj共n j,equil

⫺ n jG兲兴

(24)

1868

K. HARSTAD AND J. BELLAN TABLE 1 Simulation Conditions for Isolated Drops R d0 , R si,0 , cm

p e , MPa

T e, K

␣ a1

␣ a2

0.002; 0.017

0.006; 0.05

0.035; 0.3

0.1 0.5 1 2 4 5 0.1 0.5 1 2 4 5 5 5 5 5

470 470 470 470 470 470 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000

1 1 1 1 1 1 1 1 1 1 1 1 0.1 3 10 10

1 1 1 1 1 1 1 1 1 1 1 1 0.1 3 3 10

— — — — — — ISTH-0.1 — — ISTH-2 — ISTH-5 — — — IS10-10

IMTL-0.1 IMTL-0.5 IMTL-1 IMTL-2 IMTL-4 IMTL-5 IMTH-0.1 IMTH-0.5 IMTH-1 IMTH-2 IMTH-4 IMTH-5 IM0.1-0.1 IM3-3 IM10-3 IM10-10

— — — — — — ILTH-0.1 — — ILTH-2 — ILTH-5 — — — —

In all simulations T 0d ⫽ 325 K.

where F ems ⫽ ⫺(1/A d )dM/dt (M is the drop mass and A d is the drop boundary area) is the mass ‘emission’ rate from the drop boundary corresponding to the evaporation rate at atmoG spheric pressure, n j,equil ’s are calculated (superscript G indicates the initial nitrogen side of the drop boundary) from thermodynamic relationships ([32]) and u Tj is the mean molecular velocity crossing a plane in one direction; details are given in Harstad and Bellan [15]. Throughout this study, a thermodynamic state will be called ‘supercritical’ if either the pressure or the temperature exceed their critical value. We note that this definition of the supercritical state is not without controversy as sometimes the supercritical state is defined to be in the region where both the pressure and temper-

ature are above their critical values [33]; however, we also note that there are many subtleties to the definition of the supercritical region. For a single species fluid, p c , and T c are the maximum pressure and temperature at which two phases (liquid and/or vapor) may occur [20, 32], and therefore, the single phase region is well delineated ( p r ⬎ 1 or T r ⬎ 1). The possible complications of retrograde condensation near the critical point [20, 32] are not taken into account in this characterization. However, a binary mixture may have one or several critical loci, as shown by Prausnitz et al. [32]; in fact, at least six types or classes of binary mixtures exist with five having multiple critical loci. To avoid TABLE 3 Simulation Conditions for Binary Size Drop Clusters

TABLE 2 Far Field Pressure for the Monodisperse Cluster Simulations p e , MPa 0.1 1 2 5

M-0.1 M-1 M-2 M-5

In all simulations the accommodation coefficients are equal to unity and T e ⫽ 1000 K, T 0d ⫽ 325 K, R 0d ⫽ 6 ⫻ 0 10 ⫺3 cm, ␾0 ⫽ 4, and R C ⫽ 2 cm.

p e , MPa

␾0

0.1 1 2 5 5 5

4 4 4 4 1 2

B4-0.1 B4-1 B4-2 B4-5 B1-5 B2-5

In all simulations the accommodation coefficients are 0 equal to unity, T e ⫽ 1000 K, T 0d ⫽ 325 K, R d1 ⫽ 4.5 ⫻ ⫺3 0 ⫺3 3 10 cm, R d2 ⫽ 6 ⫻ 10 cm, Nu C ⫽ 10 , there are three 0 times more smaller than larger drops, and R C ⫽ 2 cm.

EVALUATION OF HEPTANE DROPS IN NITROGEN

1869

Fig. 1. Variation of an isolated drop boundary area, A s , and K, with time for various pressures, at subcritical (Figs. 1a and 1b), and supercritical (Figs. 1c and 1d) ambient temperatures. For all simulations R 0d ⫽ 6 ⫻ 10 ⫺3 cm, R si,0 ⫽ 5 ⫻ 10 ⫺2 cm, and T 0d ⫽ 325 K. For Figs. 1a and 1b, T e ⫽ 470 K. For Figs. 1c and 1d, T e ⫽ 1000 K. In all figures, the legend is as follows: 0.1 MPa (——), 0.5 MPa (- - - -), 1 MPa (– 䡠 – 䡠 –), 2 MPa ( 䡠 䡠 䡠 䡠 ), 4 MPa (– – –), and 5 MPa (– 䡠 䡠 – 䡠 䡠 –).

the confusion that may be introduced for such mixtures having several critical loci, we do not consider these types of mixtures; instead we focus on simple mixtures (i.e., those having a single critical locus). We define the supercritical region as that for which only a single phase is possible and one in which the transport properties are generally pressure as well as temperature dependent. We keep the criterion of p r ⬎ 1 or T r ⬎ 1, as one for which only a single phase is possible; this corresponds to our definition of the supercritical region. This definition is consistent with the use of a Peng–Robinson type equation of state and the critical locus discussed above. Isolated Drop Behavior as a Function of Temperature and Pressure Nomura et al. [16] measured the drop lifetime, t d , under microgravity conditions and found that once corrected by subtracting the drop heat up time, the dependence of t d on the ambient

pressure changed according to the ambient temperature; these results were reproduced by the simulations of Harstad and Bellan [15]. Here, we utilize the results obtained from simulations listed in Table 1 to complement the existing information by exploring the dependence of t d on the pressure at two different ambient temperatures, one being slightly more elevated than the largest temperature used in the experiments of [16]. Moreover, we also investigate the thermodynamic state arbitrarily near the drop boundary. Time History of d2 One of the most fundamental quantities depicting the evolution of a drop embedded into surroundings at a higher temperature is the time history of d 2 . Under atmospheric conditions, the variation of d 2 with t is linear, and the slope, K, of this line is termed the evaporation constant. Illustrated in Fig. 1 are d 2 and K at an ambient temperature of 470 K (Figs. 1a and 1b),

1870 and of 1000 K (Fig. 1c and 1d); the results used in these plots are those from simulations IMTL0.1 to IMTL-5, and IMTH-0.1 to IMTH-5, respectively. At low ambient temperature, the heat up time increases with pressure (see Fig. 1a), and for the larger pressure, the drop even increases in size because of the decrease in the drop density (see more on this below). Past their initial heat up time whose duration increases with pressure, the results show a nearly linear d 2 (t) variation and a monotonic increase in the drop lifetime with increasing pressure. The global non-linear variation of d 2 (t) at larger than atmospheric pressure is documented in Fig. 1b. In contrast with the variation at atmospheric pressure, where the initial heat up time is short with respect to the drop lifetime and K is globally constant, for larger pressures K is an increasing function of time and a decreasing function of pressure. Parallel results at high temperature (Fig. 1c) display a different variation. Whereas the linearity of d 2 (t) at atmospheric pressure and the non-linear behavior at larger pressure is still observed, the heat up time and the drop lifetime have a non-monotonic variation with pressure, with a maximum lifetime reached in the interval 2– 4 MPa, encompassing the critical pressure of heptane. Figure 1c qualitatively duplicates the results in Fig. 5f of Nomura et al. [16] which were obtained for the same species and in the same range of pressures, but for a larger drop and at somewhat lower temperatures (648 – 661 K). The value of K corresponding to the results of Fig. 1c is displayed in Fig. 1d and, as expected, it exhibits a larger magnitude than at the lower temperature. At atmospheric pressure, K is nearly constant, however, for the other pressures K is a much stronger function of time than at lower temperature (Fig. 1b). The change in the magnitude of K (Fig. 1d) for p ⬎ 0.1 MPa leads to a maximum lifetime (Fig. 1c) at an intermediate pressure value. At large times, K is an increasing function of pressure, whereas at the subcritical temperature with respect to the drop species it was a decreasing function of pressure. Reduced Pressure and Temperature As discussed in the Introduction, a common assumption of many models is that the drop

K. HARSTAD AND J. BELLAN boundary coincides with the mixture critical point. This statement includes actually two assumptions: (1) that the critical point is reached at a location in the mixture, and (2) that the location where the critical point is reached coincides with the drop boundary. The set of simulations discussed above offers an opportunity to evaluate this assumption because it encompasses supercritical and subcritical pressures and temperatures with respect to the pure species. Shown in Fig. 2 is the low (subcritical with respect to heptane) temperature behavior at 0.1 MPa. The plots represent the spatial variation of ␳ (Fig. 2a), p r (Fig. 2b), T r (Fig. 2c), and ␣ D (Fig. 2d) at several time stations. The density profiles show a fluid drop whose size is continuously reduced, but whose density remains constant and consistent with that of a liquid, as it should be at atmospheric pressure; outside the drop, the density magnitude is that of a gas, and therefore, does not appear on the graph scale which has been chosen purposefully to be the same for Figs. 2a, 3a, 4a, and 5a (see below). The pressure is everywhere subcritical, whereas the temperature is supercritical for the pure nitrogen (both air and nitrogen at atmospheric conditions are in the supercritical regime) and remains consistently subcritical in the drop as its temperature does not exceed the boiling point of heptane. Because the temperature is the thermodynamic variable which relaxes first, due to the very large effective Lewis number (see Harstad and Bellan [24]), the extent of the radial domain for Fig. 2c is enlarged compared to Figs. 2a and 2b, to display the full profiles. The mixture of heptane and nitrogen outside of the drop is nearly ideal ( ␣ D ⯝ 1), as expected at atmospheric pressure, and ␣ D reaches a minimum at the drop boundary, a feature prevailing in all of the results presented below. Therefore, the physical picture corresponds to the well known process of evaporation at low pressure and temperature. When the ambient temperature is increased, but the pressure remains atmospheric, the basic characteristics of the situation are unchanged with the exception of the time scale which becomes much shorter. The illustrations in Figs. 3a–3d parallel those of Figs. 2a–2d. However, the physical picture changes appre-

EVALUATION OF HEPTANE DROPS IN NITROGEN

1871

Fig. 2. Characteristics of isolated drops at 0.1 MPa and T e ⫽ 470 K. The other initial conditions are: R 0d ⫽ 6 ⫻ 10 ⫺3 cm, R si,0 ⫽ 5 ⫻ 10 ⫺2 cm, and T 0d ⫽ 325 K. Plotted are ␳ (Fig. 2a), p r (Fig. 2b), T r (Fig. 2c), and ␣ D (Fig. 2d) at several time stations: 0 s (——), 2 ⫻ 10⫺2 s (- - - -), 4 ⫻ 10⫺2 s (– 䡠 – 䡠 –), 6 ⫻ 10⫺2 s ( 䡠 䡠 䡠 䡠 ), 8.5 ⫻ 10⫺2 s (– – –), and 9.5 ⫻ 10⫺2 s (– 䡠 䡠 – 䡠 䡠 –).

ciably once the pressure is increased. Displayed in Figs. 4a– 4d are equivalent plots to those of Figs. 2a–2d, but at a pressure of 5 MPa. Instead of the constant density, we witness here a continuously decreasing density with time (Fig. 4a) resulting from the nitrogen dissolving in the heptane drop and from the drop heating. Noteworthy, although the density decreases inside the drop, it remains nearly uniform at all times. Throughout the drop lifetime, the pressure is supercritical inside the drop and subcritical (except for the initial condition) outside the drop (Fig. 4b), with the opposite occurring for the temperature (Fig. 4c). This is to say that in both sides of the drop boundary, the mixture is supercritical (according to our definition), but there is a different supercritical state. The remarkable fact is that although there is no phase change, the fluid properties on each side of the boundary are very distinct, and the characteristics of the internal liquid-like and external gaslike fluids are maintained over the drop lifetime. This is in contrast to LOX drops in H2

(see Harstad and Bellan [18]) where profiles of the primitive variables are smooth, with no jumps. The indication that the critical point is not reached anywhere in the field is emphasized by the mass diffusion factor shown in Fig. 4d which does not reach a null value (a necessary condition for being at the critical point; see Eq. 23), although the mixture exhibits considerable departures from ideality. That is, 0 ⬍ ␣ D ⬍ 1, with the null value corresponding to a possible critical point and the unity value corresponding to the ideal mixture state. At far-field supercritical temperature and pressure (Figs. 5a–5d), the features encountered at subcritical (with respect to the fuel) ambient temperature are maintained, with the difference that now the density not only decreases with time, but exhibits substantial nonuniformities inside the drop (Fig. 5a). Again, the mixture on the two sides of the drop boundary is supercritical (Figs. 5b and 5c), a conclusion reinforced by the plots of ␣ D in Fig. 5d which show greater departures from non-ideal-

1872

K. HARSTAD AND J. BELLAN

Fig. 3. Characteristics of isolated drops at 0.1 MPa and T e ⫽ 1000 K. The other initial conditions are: R d0 ⫽ 6 ⫻ 10 ⫺3 cm, R si,0 ⫽ 5 ⫻ 10 ⫺2 cm, and T 0d ⫽ 325 K. Plotted are ␳ (Fig. 3a), p r (Fig. 3b), T r (Fig. 3c) and ␣ D (Fig. 3d) at several time stations: 0 s (——), 7.5 ⫻ 10⫺3 s (- - - -), 1.5 ⫻ 10⫺2 s (– 䡠 – 䡠 –), 2 ⫻ 10⫺2 s ( 䡠 䡠 䡠 䡠 ), 2.25 ⫻ 10⫺2 s (– – –), and 2.5 ⫻ 10⫺2 s (– 䡠 䡠 – 䡠 䡠 –).

ity than at far field subcritical temperatures, but remain far from the null value. Compared to the LOX/H2 system, the heptane/nitrogen mixture stays further away from the critical point, because values of ␣ D as low as 0.25 were encountered [17] for LOX/H2 at an ambient pressure of 6 MPa and temperature of 1000 K (for LOX, p c ⫽ 5.043 MPa and T c ⫽ 154.6 K). Nevertheless, even at these ambient conditions close to the critical pressure, the LOX drop boundary never reached the critical state. With increasing ambient pressure, the state arbitrarily near the drop boundary moved further away from the critical point. Therefore, the studies of these two systems of species do not support the assumption that the critical point is reached anywhere between the drop center and the far field, or that the drop boundary coincides with the critical mixture condition. When either the far-field pressure or the farfield temperature is supercritical, the mixture on both sides of the drop boundary is supercritical, but, depending on the solubility of the

species, may be at a different supercritical state corresponding to the fluid properties of that mixture. And, as discussed above, in the context of thermodynamic modeling it is only appropri⫹ ate to discuss the conditions at R ⫺ d and R d because the conditions at R d cannot be readily defined if there is a jump at R d . Variation of d2 as a Function of Drop Size and Pressure One of the important results of the atmospheric drop evaporation theory is that K is constant and independent of drop size, being a function mainly of the fuel properties. This result is based on the quasi-steady state theory yielding the well-known d 2 law. Figures 1b and 1c show that K is no longer constant at high pressures. The question then arises as to the dependence of K on the drop size at high pressures. To investigate this aspect, K is plotted in Fig. 6 for three disparate drop sizes at 0.1, 2 and 5 MPa (complete conditions are available in Table 1)

EVALUATION OF HEPTANE DROPS IN NITROGEN

1873

Fig. 4. Characteristics of isolated drops at 5 MPa and T e ⫽ 470 K. The other initial conditions are: R d0 ⫽ 6 ⫻ 10 ⫺3 cm, R si,0 ⫽ 5 ⫻ 10 ⫺2 cm, and T 0d ⫽ 325 K. Plotted are ␳ (Fig. 4a), p r (Fig. 4b), T r (Fig. 4c), and ␣ D (Fig. 4d) at several time stations: 0 s (——), 4 ⫻ 10⫺2 s (- - - -), 6 ⫻ 10⫺2 s (– 䡠 – 䡠 –), 1.2 ⫻ 10⫺1 s ( 䡠 䡠 䡠 䡠 ), 1.5 ⫻ 10⫺1 s (– – –), and 1.8 ⫻ 10⫺1 s (– 䡠 䡠 – 䡠 䡠 –).

as a function of t/(d 0 ) 2 . At atmospheric pressure, we recover the classical result of independence of drop size and near constancy for K. At 2 and 5 MPa, the values of K do not coincide during the initial heat up time, however, this result may also be influenced by the initial conditions which although nominally the same, may be numerically somewhat different (see Harstad and Bellan [15] for a discussion on this topic). The remarkable result is the coincidence of K after this initial transient, at a time when each of the drops has lost only 9% of their initial mass (not shown). This high pressure result is not predictable from simple mathematical manipulations of the equations, and must await a thorough explanation. From these limited simulations, the indication is that on the scale t/(d 0 ) 2 , the solution is self-similar in a parameter proportional to the drop radius squared. That is, although the quasi-steady assumption that is one of the bases of the d 2 law does not hold, rendering K variable (here increasing) with time instead of constant, in this range of

drop sizes covering more than one order of magnitude, and in this range of pressures encompassing both subcritical and supercritical pressures with respect to the drop species, K is only a function of pressure at fixed temperature. Although large single drops studied under microgravity conditions provide valuable scientific information, the lack of convective flow past the drops and of drop interactions indicates that the relevance of that information to practical combustion systems remains to be established. The topic of the differences between isolated and interacting drops in quiescent surroundings is addressed in sections to follow. Effect of the Relaxation Time Scale at Supercritical Pressure The relaxation time scale at the drop boundary measures how fast molecules leave the boundary or come to it from the surrounding field. Under atmospheric, evaporative conditions, the relaxation time measures the evaporation time

1874

K. HARSTAD AND J. BELLAN

Fig. 5. Characteristics of isolated drops at 5 MPa and T e ⫽ 1000 K. The other initial conditions are: R 0d ⫽ 6 ⫻ 10 ⫺3 cm, R si,0 ⫽ 5 ⫻ 10 ⫺2 cm and T 0d ⫽ 325 K. Plotted are ␳ (Fig. 5a), p r (Fig. 5b), T r (Fig. 5c), and ␣ D (Fig. 5d) at several time stations: 0 s (——), 7.5 ⫻ 10⫺3 s (- - - -), 1.5 ⫻ 10⫺2 s (– 䡠 – 䡠 –), 2 ⫻ 10⫺2 s ( 䡠 䡠 䡠 䡠 ), 2.25 ⫻ 10⫺2 s (– – –), and 2.5 ⫻ 10⫺2 s (– 䡠 䡠 – 䡠 䡠 –).

scale. Because at supercritical pressures there is no longer a surface, and the latent heat is null, the concept of evaporation no longer applies. However, there is still a mass flux from the drop boundary and it can be calculated in a general

Fig. 6. Dependence of K on the drop size and pressure. Three drop sizes and pressures are considered, as follows: R 0d ⫽ 2 ⫻ 10 ⫺3 cm: 0.1 MPa (——), 2 MPa (- - - -), 5 MPa (– 䡠 – 䡠 –); R 0d ⫽ 6 ⫻ 10 ⫺3 cm: 0.1 MPa ( 䡠 䡠 䡠 䡠 ), 2 MPa (– – –), 5 MPa (– 䡠 䡠 – 䡠 䡠 –); R 0d ⫽ 3.5 ⫻ 10 ⫺2 cm: 0.1 MPa (●), 2 MPa (■), 5 MPa (Œ). Other properties, including the values of R si,0 , are available in Table 1.

manner employing Eq. 24; Harstad and Bellan [15] have called this mass flux ‘the emission rate’ to enable the utilization of a terminology including both subcritical and supercritical situations. Equation 24 states that the net mass flux is the difference between the flux leaving the boundary and that coming to the boundary, and the time scale of the process is determined by the values of the accommodation coefficients, ␣ aj . Under subcritical conditions, a common assumption is that this relaxation time is very small, leading to a quasi-steadiness of this process with respect to all other time scales. This small relaxation time renders the set of equations very stiff, and thus, difficult to solve. This difficulty, together with the uncertainties associated with the values of the accommodation coefficients provide the incentive to neglect Eq. 24 and find a simpler thermodynamic statement, such as the equality of the fugacities at the boundary, as a replacement. The assumption of very small relaxation time was investigated under the atmospheric condi-

EVALUATION OF HEPTANE DROPS IN NITROGEN tion by Bellan and Summerfield [34], and more recently by Miller et al. [35], and it was found to hold reasonably well for drops of diameter larger than 50 ␮m; for smaller drops, the assumption proved to be unreliable. Because as discussed in the Introduction, the hypothesis of the small relaxation time is used in some models of high pressure drops, we examine here its validity by varying the accommodation coefficients as shown in Table 1. For the large drops, R 0d ⫽ 6 ⫻ 10 ⫺3 cm, the plotted d 2 and K coincide for all five conditions representing simulations IMTH-5, and IM0.1-0.1 to IM10-10 (not shown), indicating that indeed the relaxation time is small. For small drops of R 0d ⫽ 2 ⫻ 10 ⫺3 cm, equivalent plots for the smallest and largest values of the accommodation coefficients (see Table 1) obtained from simulations ISTH-5 and IS10-10, show again complete coincidence. When plotted versus t/(d 0 ) 2 , the K curves coincide (not shown) for the two drop sizes and the different accommodation coefficients (total of seven simulations) over the entire drop lifetime with the exception of the short drop-heat-up duration, indicating that on that scale the process of drop disappearance is only pressure dependent. Single Drops Versus Monodisperse Drop Clusters Most currently existing experiments featuring high pressure drops consist of observations of single drops. Conclusions based on these observations are interpreted not only to derive scientific information, but also sometimes to obtain information relevant to high pressure energy producing devices. To understand the applicability of single drop data to situations where drops congregate in clusters (see Snyder et al. [36]), illustrated in Fig. 7a is the value of K for simulations IMTH-0.1, IMTH-1 to IMTH-5 (see Table 1), and M-01 to M-5 (see Table 2) comparing single drops to monodisperse drop clusters. At all pressures, but dramatically more at larger than atmospheric, the lifetime of drops in clusters is considerably longer due to drops competing for the available heat. This statement is supported by the plot of the interstitial temperature in the cluster, T i , displayed in Fig. 7b which shows an initial reduction due to drop

1875

Fig. 7. Comparison of K for isolated and monodisperse cluster drops (Fig. 8a) at different pressures. The interstitial temperature in the monodisperse cluster during it lifetime (Fig. 8b) at different pressures. For both isolated and monodisperse drops R 0d ⫽ 6 ⫻ 10 ⫺3 cm, T e ⫽ 1000 K and T 0d ⫽ 325 K. For the isolated drops R si,0 ⫽ 5 ⫻ 10 ⫺2 cm, whereas for the monodisperse cluster ␾0 ⫽ 4 and R C ⫽ 2 cm. The legend is as follows for isolated drops: 0.1 MPa (——), 1 MPa (- - - -), 2 MPa (– 䡠 – 䡠 –), 5 MPa ( 䡠 䡠 䡠 䡠 ). For monodisperse clusters, 0.1 MPa (—●—), 1 MPa (- - ● - -), 2 MPa (– 䡠 –●– 䡠 –), 5 MPa ( 䡠 䡠 ● 䡠 䡠 ).

heating. The further recovery of the interstitial temperature is governed by the value of the Nusselt number (see [17]) and corresponds to the increasing portion of K in Fig. 7a; however, for the ␾0 at which these simulations were conducted (the influence of ␾0 is discussed below), T i remains consistently below the initial far-field temperature. The value of T i is determined by the characteristic time associated with drop heating, tending to reduce it, and by the characteristic time associated with Nu C , tending to increase it. Because in this diffusion dominated situation, Nu C is a weak parameter [18], the drop heating process has the overwhelming influence on the lifetime. Also noteworthy is the fact that, for single drops K is a much weaker function of the ambient pressure than it is for drops in mono-

1876

K. HARSTAD AND J. BELLAN

disperse clusters. Therefore, it is not only the quantitative information that may be in error when basing spray information on isolated drop studies, but also is the qualitative aspect embodied by the parametric variation that become questionable as well. Thus, high pressure, isolated drop simulations and experiments should be interpreted cautiously when inferring results for spray applications. Monodisperse and Binary Size Drop Clusters Harstad and Bellan [17] have shown that if the smaller drops are considerably more numerous in a cluster, at far-field supercritical conditions with respect to the pure drop species, their lifetime is overestimated by that of drops in a monodisperse cluster of same ␾0. The present study is complementary to that of [17] in that the emphasis is here on the large drops. The comparison is made between the simulations listed in Tables 2 and 3. Plotted in Fig. 8a is K for drops of same size and belonging to clusters having the same initial and far-field conditions, except that in one cluster the drops are of same size whereas in the other there is a binary drop size with the larger drops being less numerous; these larger drops have the same initial radius as that of the drops in the monodisperse cluster. Because the small drops heat up faster, they use the available heat in the interstitial region at a higher rate, thus decreasing the K of drops in binary clusters with respect to those in monodisperse clusters. This effect is clearly seen in Fig. 8b where T i is consistently smaller for the binary cluster until the small drops disappear. The kink in the curves obtained from binary size simulations corresponds to the time of the small drops disappearance, and is followed by a larger rate of increase in T i . Although K is larger for drops in monodisperse clusters than for drops of same initial size in a cluster where they are the less numerous large drops, dK/dt appears to be similar. Thus, the present results contrast those from the previous study [17] that focused on the smaller drops in a binary cluster. Whereas small, more numerous drops in binary clusters disappear faster than similar size drops in a monodisperse cluster, the larger, less numerous drops disappear slower than same initial size

Fig. 8. (a) Comparison of K for monodisperse and binary cluster drops for drops of same initial size; these drops are the larger, less numerous drops in the binary cluster. For both clusters T e ⫽ 1000 K, T 0d ⫽ 325 K, ␾ 0 ⫽ 4, Nu C ⫽ 10 3 and R C ⫽ 2 cm. For the monodisperse cluster R d0 ⫽ 6 ⫻ 0 10 ⫺3 cm, and for the binary size cluster R d1 ⫽ 4.5 ⫻ 10 ⫺3 0 cm, R d2 ⫽ 6 ⫻ 10 ⫺3 cm and there are three times more smaller than larger drops. (b) Comparison of T i for the monodisperse and binary size cluster drops as described in (a). In both figures the legend is as follows: p e ⫽ 0.1 MPa, 0 0 R d2 (- - - -), R 0d (——), R 0d (—●—); p e ⫽ 1 MPa, R d2 0 (- - ● - -); p e ⫽ 2 MPa, R d2 (– 䡠 – 䡠 –), R 0d (– 䡠 –●– 䡠 –); p e ⫽ 5 0 MPa, R d2 ( 䡠 䡠 䡠 䡠 ), R 0d ( 䡠 䡠 ● 䡠 䡠 ).

drops in monodisperse clusters. This conclusion indicates that approximation of a polydisperse high pressure fluid spray by a collection of monodisperse drops could lead to errors whose magnitude is difficult to determine a priori. Because the small and large drops in a binary cluster have different heat up histories, it is of interest to inquire which size class is closer to the critical point during its lifetime. Because ␣ D is a measure of the departure from ideality ( ␣ D ⫽ 1) and from the critical point ( ␣ D ⫽ 0), depicted in Fig. 9 is a typical plot for both drop sizes at 5 ⫻ 10⫺2 s for pressures 1, 2, and 5 MPa; the drops from the simulation conducted at 0.1 MPa are expected to be far from the critical point and have disappeared before that time

EVALUATION OF HEPTANE DROPS IN NITROGEN

Fig. 9. The mass diffusion factor for small and large size drops belonging to the same, binary size cluster, as a function of pressure. All initial conditions are those stated in the caption of Fig. 8. The plots correspond to t ⫽ 5 ⫻ 0 10 ⫺2 s, and the legend is as follows: p e ⫽ 1 MPa, R d1 0 0 0 (- - - -), R d2 (- - ● - -); p e ⫽ 2 MPa, R d1 (– 䡠 – 䡠 –), R d2 0 0 (– 䡠 –●– 䡠 –); p e ⫽ 5 MPa, R d1 ( 䡠 䡠 䡠 䡠 ), R d2 ( 䡠 䡠 ● 䡠 䡠 ).

station. Clearly, the smaller drops are somewhat closer to the critical point, although drops of both sizes are considerably far from it. Figure 10 portrays the residual mass of the two drop sizes in the cluster and shows the increasing difference with increasing pressure. Influence of the Drop Number Density at Supercritical Pressure To assess the influence of the heptane/nitrogen mass ratio on the above results, simulations were conducted with several values of ␾0 (see Table 3) and the results are illustrated in Fig. 11

1877

Fig. 11. Mass diffusion factor at drop half mass time as a function of the radial coordinate for several values of ␾0; 1, 2, and 4. All other initial conditions are those stated in the 0 caption of Fig. 8. The legend in the figure is as follows: R d1 : 0 ␾ 0 ⫽ 1 (——), ␾0 ⫽ 2 (- - - -), ␾0 ⫽ 4 (– 䡠 – 䡠 –); R d2 : ␾0 ⫽ 1 ( 䡠 䡠 䡠 䡠 ), ␾0 ⫽ 2 (– – –), and ␾0 ⫽ 4 (– 䡠 䡠 – 䡠 䡠 –).

depicting ␣ D at the drop half-mass times for each size class. With increasing ␾0 the mixture is increasingly non-ideal both inside and outside the drop. In the drop, this is because of the decrease in the reduced pressure towards unity (not shown), whereas outside the drop this is because of the decrease in the reduced temperature towards unity (not shown). Indeed, as the drop number density is increased, T i declines faster as there is an increased heat transfer because of the larger number of drops. For the smaller drop size, the conditions at the drop boundary move towards ideality and away from the critical point with increasing ␾0, whereas for the larger size drops it is exactly the opposite. However, even for these larger drops, the state of the drop at the boundary is never close to the critical point. Therefore, the conclusion that the drop boundary does not reach the critical point is further supported by these results. CONCLUSIONS

Fig. 10. Residual mass of drops in a binary cluster as a function of time. All initial conditions are those stated in the caption of Fig. 8. The legend is as follows: p e ⫽ 0.1 MPa, 0 0 0 0 R d1 (- - - -), R d2 (——), R d2 (—●—); p e ⫽ 1 MPa, R d1 0 0 (- - ● - -); p e ⫽ 2 MPa, R d1 (– 䡠 – 䡠 –), R d2 (– 䡠 –●– 䡠 –); p e ⫽ 0 0 5 MPa, R d1 ( 䡠 䡠 䡠 䡠 ), R d2 ( 䡠 䡠 ● 䡠 䡠 ).

A numerical study has been conducted to assess some of the assumptions commonly employed in studies of high pressure drops in quiescent surroundings. The study was based on a model of isolated drop behavior at all pressures which has been previously validated with microgravity experimental data. This model is based on fluctuation-dissipation theory which is consistent with non-equilibrium thermodynamics, and

1878 leads to the most general form of the flux equations, including the Soret and Dufour effects. To close the system, conservation equations are completed by the equation of state, valid for all pressures. The traditional transport coefficients used in the simulations were calculated as functions of the thermodynamic variables; the thermal diffusion factor, which is the new transport coefficient introduced by the Soret and Dufour effects, was chosen according to the value determined during the previous validation study. The drop boundary conditions were based on the conservation of mass, energy and species across the boundary, with no assumption made about the occurrence of a quasisteadiness in mass transfer to/from the boundary. Calculated after obtaining the solution, the drop boundary was identified with the location of the maximum density change because of the relevance to optical measurements. Results from this isolated drop model were also compared with those from monodisperse and binary size clusters which were also based upon the same fundamental fluid dynamic model and high pressure thermodynamics; for clusters, the isolated drop equations were augmented by the cluster conservation equations. The results show that the far-field subcritical (with respect to the pure drop species) temperature drop behavior is different from that at supercritical (also with respect to the pure drop species) far-field temperature. Whereas the drop lifetime increases monotonically with pressure at subcritical temperatures, at supercritical temperatures a maximum is reached. This result is confirmed by existing microgravity data. Moreover, it was also shown that at supercritical pressure, the critical point is never reached arbitrarily close to the drop boundary. Instead, two supercritical states, each portraying the specific mixture, coexist on the two sides of the drop boundary. Therefore, these results do not support previous models based on the assumption that the critical point is attained in the domain, and that it is reached at the drop boundary. An examination of the assumption of linear dependency of the drop diameter squared with time shows that this functional behavior is approximately encountered only at atmospheric pressure. As the pressure increases, the slope of

K. HARSTAD AND J. BELLAN the diameter squared, K, becomes an increasing function of time; this slope progressively increases with pressure at subcritical far-field temperatures, however, at supercritical temperatures a maximum is reached. Results obtained for different drop diameters spanning more than an order of magnitude display the traditional independence of K on the drop size at atmospheric pressure when time is scaled as t/(d 0 ) 2 . As the pressure is increased, the early part of K(t) [corresponding both to the initial drop heating and to possible transients arising from initial conditions that are nominally the same but perhaps numerically somewhat different (see explanation in [15])] is drop size dependent, however, during the majority of the drop lifetime corresponding to more that 90% of the mass, the values of K are again size independent. This indicates that the solution is selfsimilar in a variable proportional to the drop radius squared; this fact cannot be predicted from simple scrutiny of the equations and should be further investigated. Similar to results previously obtained at atmospheric pressure, the effect of the relaxation time scale at the drop boundary was found to be negligible for large drops. For atmospheric pressure, it has been previously shown that mass flux unsteadiness (usually called somewhat misleadingly thermodynamic non-equilibrium effects) in the drop evaporation rate, become important for drops of diameter smaller than 50 ␮m. Comparably small size drops at supercritical pressure were, however, shown to display a quasi-steady mass flux behavior. Therefore, the assumption of quasi-steadiness in the mass fluxes to/from the drop surface appears to be valid within the range of drop sizes investigated. From the practical viewpoint, this finding could be important, if of general validity, and used judiciously could simplify CFD calculations. Comparisons of isolated drops and drops belonging to monodisperse clusters at supercritical pressure lead to the conclusion that not only the drop lifetime is shorter because of the lack of drop interaction, but also the pressure dependency of K is considerably weaker for isolated drops than for drops in clusters, indicating that results from isolated drop parametric studies must be cautiously interpreted in the context of spray studies. Furthermore, it was

EVALUATION OF HEPTANE DROPS IN NITROGEN also found that the less numerous, larger size drops in binary size clusters disappear slower than same size drops in monodisperse clusters of the same fuel-to-surrounding-fluid mass ratio, because of the larger heat transfer to the small drops. The boundary of the smaller drops in a binary size cluster was shown to be at a state closer to the critical point than that of larger drops, although in both cases the state was far from the critical point. Finally, the influence of the fuel-to-surrounding-fluid mass ratio, that is, of the drop number density, upon these results was examined for binary size clusters, and it was shown that the conclusions focussing on the thermodynamic state of the boundary hold independently of the drop number density. This research was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under sponsorship from the National Aeronautics and Space Administration, the Glenn Research Center at Lewis Field, with Dr. Daniel L. Bulzan as technical contract monitor. His continuing interest and support are greatly appreciated. REFERENCES 1. 2. 3. 4.

5. 6. 7. 8.

9. 10. 11.

Bellan, J., Prog. Energy Combust. Sci. 26(4 – 6):329 –366 (2000). Shuen, J.-S., Yang, V., and Hsiao, G. C., Combust. Flame 89:299 –319 (1992). Yang, V., Lin, N. N., and Shuen, J.-S., Comb. Sci. and Tech. 97:247–270 (1994). Hsiao, G. C., Yang, V., and Shuen, J. S., Supercritical vaporization and dynamics of liquid oxygen (LOX) droplet in hydrogen stream, AIAA 95-0383, 33rd Aerospace Sciences Meeting, Reno, NV, 1995. Haldenwang, P., Nicoli, C., and Daou, J., Int. J. Heat Mass Transfer 39(16):3453–3464 (1996). Umemura, A., Proc. of the Comb. Inst. 21:463– 471 (1986). Curtis, E. W., and Farrell, P. V., Acta Astronautica 17(11/12):1189 –1193 (1988). Chiang, C. H., and Sirignano, W. A., Axisymmetric vaporizing oxygen droplet computations, AIAA 910281, 29th Aerospace Sciences Meeting, Reno, NV, 1991. Curtis, E. W., and Farrell, P. V., Combust. Flame 90:85–102 (1992). Jia, H., and Gogos, G., J. Thermophysics and Heat Transfer 6(4):738 –745 (1992). Jia, H., and Gogos, G., Int. J. Heat Mass Transfer 36(18):4419 – 4431 (1993).

12. 13. 14. 15. 16. 17. 18. 19. 20.

21.

22. 23.

24. 25. 26. 27.

28. 29. 30. 31. 32.

33. 34. 35. 36.

1879

Jiang, T. L., and Chiang, W.-T., Comb. Flame 97:17–34 (1994). Jiang, T. L., and Chiang, W.-T., Int. J. Heat Mass Transfer 39(5):1023–1031 (1996). Umemura, A., and Shimada, Y., Proc. of the Comb. Inst. 27:2659 –2665 (1998). Harstad, K., and Bellan, J., Int. J. of Multiphase Flow 26(10):1675–1706 (2000). Nomura, H., Ujiie, Y., Rath, H. J., Sato, J., and Kono, M., Proc. of the Comb. Inst. 26:1267–1273 (1996). Harstad, K., and Bellan, J., Combustion and Flame, 124:535–550 (2001). Harstad, K., and Bellan, J., Int. J. Heat and Mass Transfer 41:3551–3558 (1998). Keizer, J., Statistical Thermodynamics of Nonequilibrium Processes, Springer–Verlag, New York, 1987. Hirshfelder, J. O., Curtis, C. F., and Bird, R. B., Molecular Theory of Gases and Liquids, John Wiley and Sons, Inc., New York, NY, 1964. Chapman, S., and Cowling, T. G., The Mathematical Theory of Nonuniform Gases, Cambridge University Press, Cambridge, 1970. Sarman, S., and Evans, D. J., Phys. Rev. A45(4):2370 – 2379 (1992). Reid, R. C., Prausnitz, J. M., and Polling, B. E., The Properties of Gases and Liquids, 4th Edition, McGraw– Hill Book Company, New York, NY, 1987. Harstad, K., and Bellan, J., Int. J. Heat and Mass Transfer 42:961–970 (1999). Harstad, K., Miller, R. S., and Bellan, J., AIChE J. 43(6):1605–1610 (1997). Harstad, K., and Bellan, J., Int. J. Heat and Mass Transfer 41:3537–3550 (1998). Rabinovich, V. A., and Abdulagatov, I. M., Viscosity and Thermal Conductivity of Individual Substances in the Critical Region, Begell House, New York, NY, 1996. Bzowski, J., Kestin, J., Mason, E. A., and Uribe, F. J., J. Phys. Chem. Ref. Data 19(5):1179 –1232 (1990). Peng, D.-Y., and Robinson, D. B., AIChE J. 23(2):137– 144 (1977). Reamer, H. H., and Sage, B. H., J. Chem. Eng. Data 8(4):508 –513 (1963). Kay, W. B., Ind. Eng. Chem. 30:459 – 465 (1938). Prausnitz, J., Lichtenthaler, R., and de Azevedo, E., Molecular Thermodynamics for Fluid-Phase Equilibrium, 2nd Edition, Prentice–Hall, Inc., Englewood Cliffs, NJ, 1986. Clifford, T., Fundamentals of Supercritical Fluids, Oxford University Press, Oxford, UK, 1999. Bellan, J., and Summerfield, M., Combust. Flame 33:107–122 (1978). Miller, R. S., Harstad, K., and Bellan, J., Int. J. of Multiphase Flow 24:1025–1055 (1998). Snyder, R., Herding, G., Rolon, J. C., and Candel, S., Combust. Sci. and Tech. 124:331–370 (1997).

Received 27 November 2000; revised 28 June 2001; accepted 8 July 2001