Real-Gas Effects and Phase Separation in Underexpanded Jets at Engine-Relevant Conditions
RReal-Gas Effects and Phase Separation inUnderexpanded Jets at Engine-Relevant Conditions (cid:73)
Christoph Traxinger a, ∗ , Matthias Banholzer a , Michael Pfitzner a a Institute for Thermodynamics, Bundeswehr University Munich,Werner-Heisenberg-Weg 39, 85577 Neubiberg, Germany
Abstract
A numerical framework implemented in the open-source tool OpenFOAM is presented in this workcombining a hybrid, pressure-based solver with a vapor-liquid equilibrium model based on the cubicequation of state. This framework is used in the present work to investigate underexpanded jets atengine-relevant conditions where real-gas effects and mixture induced phase separation are probable tooccur. A thorough validation and discussion of the applied vapor-liquid equilibrium model is conductedby means of general thermodynamic relations and measurement data available in the literature. Engine-relevant simulation cases for two different fuels were defined. Analyses of the flow field show that the usedfuel has a first order effect on the occurrence of phase separation. In the case of phase separation twodifferent effects could be revealed causing the single-phase instability, namely the strong expansion andthe mixing of the fuel with the chamber gas. A comparison of single-phase and two-phase jets disclosedthat the phase separation leads to a completely different penetration depth in contrast to single-phaseinjection and therefore commonly used analytical approaches fail to predict the penetration depth.
1. Introduction
Many of today’s and future transportation and power generation systems are based on the combustionof fossil fuels. In recent years, the concerns of environmental protection and global warming increasedsignificantly. Thus, the need of the reduction of emissions and fuel consumption of engines are the majordrivers for new innovations and technical improvements. One of the major trends in all types of enginesis the steady increase in operating pressure to fulfill the aforementioned goals. This leads to situationswhere the mixture and the combustion of fuel and oxidizer takes place at supercritical pressures ( p > p c )with respect to the pure components value. At supercritical state, the thermodynamic properties arenon-linear functions of temperature and pressure and the widely used ideal gas law is not valid anymore,especially at low/cryogenic and moderate temperatures. This aspect is crucial during the injectionprocess and many researchers have therefore focused on understanding and optimizing the high-pressureinjection prior to combustion, mostly using computational fluid dynamics (CFD).In the field of liquid or liquid-like fuel injection into a gaseous environment at supercritical pressuretypically encountered in gasoline, diesel and rocket engines, it is state-of-the-art to conduct large-eddysimulations (LES) based on real-gas thermodynamics. Pioneer work was done by Oefelein and Yang [1]as well as Zong et al. [2]. Many research groups, e.g., Schmitt et al. [3] and M¨uller et al. [4], have followedtheir method and used the dense-gas approach to do LES at rocket engine relevant conditions. Morerecently, Matheis and Hickel [5] and Traxinger et al. [6] applied a vapor-liquid equilibrium (VLE) modeland showed that phase separation due to non-linear mixing phenomena is likely to occur at typical in-jection conditions of diesel and rocket engines. In contrast to these thoroughly investigated applications, (cid:73) Preprint submitted to AIAA Scitech 2018, Kissimmee, Florida ∗ Corresponding author: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J a n he injection of gaseous fuel, especially at conditions typical for piston engines, is a topic which has notreceived as much attention as the aforementioned applications. This is changing considerably since thelast few years as environmental and global warming aspects are getting more and more important andnew fuel concepts have to be considered also for piston engines. Challenges resulting from the directinjection of gaseous fuels at high pressure ratios are the gas dynamics of the near-nozzle flow structure,the gaseous jet penetration and the fuel-oxidizer mixing. These fluid dynamic aspects are accompaniedby strong pressure and temperature changes leading to non-negligible real-gas effects, see, e.g., Khaksar-fard et al. [7] and Bonelli et al. [8]. In Fig. 1 the schematic of a highly underexpanded jet [9] forminga series of shock and expansion structures is illustrated like it occurs during the injection process whenhigh nozzle pressure ratios are present as it will be the case in this study. Mach diskSlip lineReflected shockCompressionwaves InterceptingshockJet boundaryExpansionfanMa = 1 Ma > (cid:29) < Figure 1: Schematic of an underexpanded jet showing a series of expansions and shock structures [9].
Recently, Banholzer et al. [10] presented first results of the injection of gaseous hydrogen into airwith the focus on the jet penetration depth as well as the Mach disk position. They showed an excellentagreement between experimental Schlieren images and their simulations. The present study is expandingthese results focusing on the comparison of two different fuels, namely a hydrogen and a methane basedone, and on the thermodynamics itself because they are a crucial aspect at engine-relevant conditions.We will extend the real-gas framework presented by Banholzer et al. [10] with the VLE-model used inTraxinger et al. [6] and show that phase separation is likely to occur at engine-relevant conditions. Tothe authors knowledge this is the first time that simulations of underexpanded jets at such conditionsare presented. As experimental investigations in terms of underexpanded jets at real-gas conditions arelacking we will strongly focus on a thorough discussion and validation of the VLE-model by using availableexperimental data from VLE measurements of the appropriate mixtures. The paper is structured asfollows. In the sections 2 and 3 the numerical approach and the VLE model are discussed in detail.Section 4 presents the validation of the numerical framework (solver and thermodynamic model) usingexperimental and numerical findings recently presented by Traxinger et al. [6]. In section 5 the engine-relevant simulation cases for two different fuels are defined and a-priori analyses are carried out concerningphase separation probability and choked nozzle conditions. The main results are discussed in section 6while the conclusions and some aspects concerning future work will be outlined in section 7.
2. Numerical Method
Numerical simulations were carried out based on the inert and compressible conservation equationsfor mass, momentum and energy and the transport equation for species k∂ρ∂t + ∂ ( ρu i ) ∂ζ i = 0 , (1) ∂ ( ρu i ) ∂t + ∂ ( ρu i u j ) ∂ζ j = − ∂p∂ζ i + ∂σ ij ∂ζ j , (2) ∂ρe t ∂t + ∂ ( ρu i e t ) ∂ζ i = − ∂ ( u i p ) ∂ζ i + ∂ ( u i σ ij ) ∂ζ j − ∂q i ∂ζ i , (3) ∂ ( ρY k ) ∂t + ∂ ( ρu i Y k ) ∂ζ i = − ∂D k,i ∂ζ i , (4)2herein t is the time, ζ is the vector of Cartesian coordinates, ρ is the density, u is the velocity vector, p is the static pressure, σ ij is a component of the viscous stress tensor, q is the heat flux vector and e t is thetotal energy given as e t = e + u /
2, where e is the internal energy. Furthermore, Y is the vector of themass fractions and D is the mass flux vector. The heat flux vector and the mass flux vector are modeledusing Fourier’s and Fick’s law, respectively. In the numerical simulations the Reynolds-Averaged Navier-Stokes (RANS) equations are solved. The turbulence closure was achieved with the widely-used SSTturbulence model. The final closure problem in terms of primitive and derived thermodynamic propertieswas solved by applying real-gas thermodynamics, which will be thoroughly discussed in section 3. In order to solve the system of equations a hybrid, pressure-based solver is used which was originallyderived by Kraposhin et al. [11]. It combines the Pressure Implicit with Splitting of Operator (PISO)algorithm proposed by Issa [12] and the Kurganov-Tadmor [13] (KT) scheme. While the PISO algorithmis suitable for fluid flows with small Mach numbers, it lacks of accuracy for high speed flows due tothe occurrence of numerical oscillations in regions with discontinuities. The KT scheme however is ahigh-resolution central scheme for trans- and supersonic flows. A blending function based on the localMach number and the Courant-Friedrichs-Lewy (CFL) criterion switches between the incompressible andcompressible flux formulations. Wherever the flow approaches trans- and supersonic flow regimes theblending function manages the switch to the high-resolution central scheme of Kurganov-Tadmor. Thisscheme uses more precise information of local characteristic propagation speeds at the cell boundaries,which are projected from the cell centers onto the cell faces. This leads to a separate treatment of smoothand non-smooth regions. Non-smooth parts of the computed approximations are averaged over smallercells than smooth parts. The introduced numerical diffusion is independent of ∆ t . As the high-resolutionscheme relies on the correct local characteristic propagation speeds at the cell faces, which are dependenton the local speed of sound, the exact evaluation of those is a necessity, for single-phase regions andespecially for two-phase regions where the speed of sound drops significantly at the phase boundariesdue to inhomogeneities in the fluid and discontinuities become present (see chapter 3, section 3.3).Recently, Kraposhin et al. [14] extended the source code of the solver based on the work of Jarczykand Pfitzner [15] and M¨uller et al. [16] such that it consistently accounts for real-gas thermodynamics.Additionally, the multi-species transport equations were implemented by Banholzer et al. [10]. In flowswith real-gas mixing the thermodynamic properties are much more sensitive to the transported fields,i.e., species, enthalpy and pressure, than in ideal-gas flows. Therefore, it is important to include thespecies equation, the enthalpy equation as well as the evaluation of the thermodynamic properties intothe PISO loop [15]. In addition, the modeling of the density ρ in the transient term ∂ρ/∂t in the pressureequation has to be done differently compared to the ideal-gas approach [15]. As the density is not alinear function of the pressure in the real-gas regime a Taylor expansion for the evaluation of the densityinside the pressure equation is used ρ = ρ + ∂ρ∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( p − p ) (5)where the index 0 refers to the base point of the Taylor expansion, i.e., the last time or iteration step.Due to the fact that the transport properties are solved sequentially only the values of the respectivetransported field change when it is solved. This in turn implies that the enthalpy and the species fields areconstant during the evaluation of the pressure equation and therefore ∂ρ/∂p is evaluated at isenthalpicconditions ( h = const.) and at constant compositions ( Y k = const.), i.e., ∂ρ/∂p (cid:12)(cid:12) h,Y k . Therefore, we needto use the isenthalpic compressibility ψ h as the basis for the evaluation due to the following mathematicalrelationship: ∂ρ∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h,Y k = − v ∂v∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h,Y k = 1 v ψ h . (6)Here, v is the specific volume. The evaluation of the compressibility in single and two-phase systemswill be discussed in the next section. For a detailed derivation of the pressure equation and a thoroughvalidation of the solver see Kraposhin et al. [14]. 3 . Thermodynamic Modeling At engine-relevant injection conditions characterized by supercritical pressures and moderate tem-peratures real-gas effects are a prominent feature [17, 18] and have to be taken into account within thenumerical framework to accurately model the fluid behavior. Commonly accepted in computational fluiddynamics (CFD) due to their efficiency and acceptable accuracy are cubic equation of states (EoS) basedon the corresponding states principle, which can be written in the following general, pressure-explicitform [19]: p = R Tv − b − a ( T ) v + ubv + wb = R Tv − b − a c α ( T ) v + ubv + wb . (7)Here, R is the gas constant, T is the static temperature and v = 1 /ρ is the specific volume. Theparameters a ( T ) = a c α ( T ) and b account for the intermolecular attractive and repulsive forces, respec-tively, and u and w are model constants. The most popular cubic EoS-models are the ones by Pengand Robinson [20] (PR-EoS) and Soave, Redlich and Kwong [21] (SRK-EoS) where u and w are (2,-1)and (1,0), respectively. Due to the different model constants, a c and b are different among the cubicEoS-models which results in turn in different EoS-specific critical compressibility factors Z c, PR = 0 . Z c, SRK = 0 .
333 for the SRK-EoS. As the critical compressibility factor is a fluidspecific parameter, see Tab. 2, neither the PR-EoS nor the SRK-EoS can predict the p , v , T -behavior ofall real fluids accurately and an appropriate EoS-model has to be chosen based on the particular probleminvestigated. As an example, Fig. 2 shows that the PR-EoS as well as the SRK-EoS can only predict alimited number of alkanes from the homologue series within an acceptable error [22, 23]. Generally, thedensity/compressibility factor prediction of the SRK-EoS is better for simple fluids like methane CH and hydrogen H , while the PR-EoS gives superior results for longer alkanes like n-hexane C H . Thisfinding can be further emphasized by comparing the isotherms of these different fluids over a wide pres-sure range, see Fig. 2 b). For hydrogen the PR-EoS gives a mean deviation from the reference data [24]of about -3.81% and the SRK-EoS of about -0.13%. In contrast, for n-hexane the mean deviation of thePR-EoS is -1.38% and for the SRK-EoS it is 10.46%. In this study, we will investigate three differentinjectants, namely a hydrogen fuel mixture, compressed natural gas (consisting mostly of methane) andn-hexane. In order to get the best fluid modeling possible, we will use the PR-EoS for n-hexane and theSRK-EoS for the other two fuels containing hydrogen and methane. Table 1: Parameters of the general cubic equation of state.
Parameter PR-EoS [20] SRK-EoS [21] u w -1 0 a c . R T c p c . R T c p c α ( T ) (cid:104) κ (cid:16) − (cid:113) TT c (cid:17)(cid:105) κ . . ω − . ω .
480 + 1 . ω − . ω b . R T c p c . R T c p c Z c .
307 0 . a and b independent of the chosen EoS: a = N c (cid:88) i N c (cid:88) j ξ i ξ j a ij and b = N c (cid:88) i ξ i b i . (8)Here, ξ i is the mole fraction of species i , whereby in the following we denote the overall mole fractionby z = { z , ..., z N c } and the liquid and vapor mole fractions by x = { x , ..., x N c } and y = { y , ..., y N c } ,respectively. The variables a ij and b i in Eq. (8) are calculated using the corresponding state principleand can therefore be evaluated dependent on the EoS, see Tab. 1. For the calculation of the diagonal4 2 3 4 5 6 7 8 9 10400500600 Number of carbons ρ [ k g/ m ] CoolPropPR-EoSSRK-EoS
Methanen-HexaneHydrogen p [bar] Z [ - ] CoolPropPR-EoSSRK-EoS a) b)
Figure 2: Comparison of PR-EoS and SRK-EoS with reference data from CoolProp [24]: a) Density ρ for different alkanesfrom the homologue series at p r = pp c = 1 . T r = TT c = 0 .
75; b) Compressibility factor Z = pv R T of methane, hydrogenand n-hexane at T = 300 K. elements of a ij the respective critical parameters of the pure components are used, see Tab. 2. Theoff-diagonal elements of a ij are estimated using the pseudo-critical combination rules [19]: ω ij = 0 . ω i + ω j ) , v c,ij = (cid:16) v / c,i + v / c,j (cid:17) , Z c,ij = 0 . Z c,i + Z c,j ) ,T c,ij = (cid:112) T c,i T c,j (1 − k ij ) and p c,ij = Z c,ij R T c,ij /v c,ij . (9)The binary interaction parameter k ij in Eq. (9) is usually used to fit the mixture to available measurementdata. For this study we are setting k ij to zero and will discuss later how this assumption compares tothe experimental data of the appropriate multi-component mixtures investigated.For the calculation of the caloric properties, like enthalpy and specific heat, the departure functionformalism is used, see, e.g., Poling et al. [25]. The reference condition is determined using the seven-coefficient NASA polynomials proposed by Goos et al. [26]. The viscosity and the thermal conductivityare modeled with the empirical correlation proposed by Chung et al. [27]. Cubic equations of state are generally suitable to describe the complete p , v , T -behavior of single-component as well as multi-component fluids. Therefore, they are also able to predict the phase sepa-ration in fluids which might occur when a fluid enters the subcritical regime (temperature and pressurebelow the critical values). In contrast to single component fluids, which have a distinct, fluid specificcritical point, see black dots in Fig. 3, multi-component mixtures show a critical locus, see gray line inFig. 3, reaching over a finite temperature and pressure range. For binary mixtures with type I phasebehavior the critical locus is spanning a line from the low volatile component to the high volatile compo-nent [28, 29], see Fig. 3. The mixture critical temperature is limited by the pure components values whilethe critical pressure of the mixture often by far exceeds the critical pressure of the pure components [30].Thus, there is a clear upper border in temperature above which no phase separation will occur, namelythe temperature of the high volatile component (Fig. 3 b): T c , n-C H = 469 . T c , CH =190 .
56 K). In terms of critical pressure, Fig. 3 shows that the specific form of the critical locus and there-fore the maximum critical pressure is strongly dependent on the pure components forming the mixture.Compared to the critical pressure of the high volatile component, the maximum mixture critical pressurefor the binary mixture ethane + n-heptane is approximately a factor of 3 . . able 2: Critical properties ( p c , T c and Z c ) and acentric factor ω of hydrogen H , methane CH , ethane C H , n-hexaneC H and nitrogen N taken from CoolProp [24]. Specie p c T c Z c ω [MPa] [K] [-] [-]Hydrogen (H ) 1.2964 33.15 0.303 -0.2190Methane (CH ) 4.5992 190.56 0.286 0.0114Ethane (C H ) 4.8722 305.32 0.280 0.099n-Hexane (C H ) 3.0340 507.82 0.266 0.299Nitrogen (N ) 3.3958 126.19 0.289 0.0372due to this separation compared to the single phase state or the state with lesser phases, see, e.g.,Michelsen and Mollerup [32]. As we will only consider a maximum of two phases within this study,a phase separation and therefore a reduction in Gibbs energy will only be executed from the single-phase state to a two-phase state. To take into account this possible phase separation process within thepresent study, we are using a vapor-liquid equilibrium (VLE) model, see, e.g., Matheis and Hickel [5]and Traxinger et al. [6]. In this model, the stability check is done with the help of the tangent planedistance (TPD) method of Michelsen [33] T P D ( w ) = (cid:88) i w i [ln w i + ln ϕ i ( w ) − ln z i − ln ϕ i ( z )] , (10)where the Gibbs energy of a trial phase composition w = { w , ..., w N c } is compared to the Gibbsenergy of the respective feed mixture z . The Gibbs energy is expressed based on the fugacity coefficientof component i ϕ i which is defined as f li / ( x i p ) for the liquid and as f vi / ( y i p ) for the vapor phase,respectively. The fugacity of component i in the respective phase is denoted by f i . The natural logarithmof ϕ i can be calculated directly from the applied cubic EoS:ln ( ϕ i ) = ∂ (cid:0) F − F ig (cid:1) / ( RT ) ∂n i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T,V − ln ( Z ) . (11)Here, ∂ ( F − F ig ) / ( RT ) ∂n i (cid:12)(cid:12) T,V is the partial derivative of the departure function of the normalized FreeHelmholtz energy F with respect to the mole number of the i -th component n i . This partial derivativereads − ln (1 − bρ ) + b i b ( Z − − a √ bRT (cid:18) (cid:80) j x j a ij a − b i b (cid:19) ln (cid:34) bρ (cid:0) √ (cid:1) bρ (cid:0) − √ (cid:1) (cid:35) for the PR-EoS and − ln (1 − bρ ) + b i b ( Z − − abRT (cid:18) (cid:80) j x j a ij a − b i b (cid:19) ln (cid:18) bv (cid:19) for the SRK-EoS. If the TPD-analysis of Eq. (10) leads to a negative value for any of the trial phasecompositions w , the mixture is unstable. In this investigation, a separation in two phases, namely avapor ( v ) and a liquid ( l ) phase, is done yielding a decrease in Gibbs energy. It is important to noticethat this phase separation is assumed to occur instantaneously, so no non-equilibrium effects are takeninto account, which might be present especially at high Mach numbers as it is for example the case in laststages of low-pressure steam turbines [34]. After the detection of an unstable mixture an iso-energeticand isobaric flash is solved. Due to the fact that we are using an enthalpy based energy conservationequation the flash can be considered a hpn -flash which is constructed as a nested loop whereby the outerloop is updating the temperature in order to meet the energy-criterion and in the inner loop a T pn -flashfor the appropriate temperature is solved. The solution of the
T pn -flash is characterized by the equalityof the fugacities of each component in the considered phases [30]: f li ( p, T, x ) = f vi ( p, T, y ) . (12)6or solving the TPD-analysis and the T pn -flash we strongly followed the suggestions of Michelsen andMollerup [32] and implemented three different methods, a successive substitution method (SSM), adominant eigenvalue method (DEM) and Newton’s method (NM). Usually, we are starting with thesuccessive substitution for solving both problems and after a certain number of evaluation steps weare switch to the higher order methods (DEM or NM) to speed-up the iteration process. This is alsoimportant for near-critical calculations as the SSM might not converge or only with a very large numberof iterations [32]. . . . . x C2H6 , y C2H6 , z C2H6 [-] p [ b a r ] T = 338 .
71 K T = 366 .
48 K T = 394 .
26 K T = 422 .
04 K T = 449 .
82 K 250 300 350 400 450 500 550 600 65020406080100 C H n-C H T [K] p [ b a r ] z = 0 . z = 0 . z = 0 . z = 0 . z = 0 . . . . . x N2 , y N2 , z N2 [-] p [ b a r ] T = 150 K T = 160 K T = 170 K T = 180 K 100 120 140 160 180 20010254055 N CH T [K] p [ b a r ] z = 0 . z = 0 . z = 0 . z = 0 . a) b)c) d) Figure 3: VLEs for binary mixtures of type I compared to experimental data [35, 36, 37, 38, 39, 40, 41] (shown as dots): a) pxy -diagram for ethane + n-heptane, b) pT -diagram for ethane + n-heptane, c) pxy -diagram for nitrogen + methane andd) pT -diagram for nitrogen + methane. The black line shows the saturation curve of the pure component and the blackdot marks the pure component’s critical point. The gray line is the critical locus of the binary mixture calculated basedon the approach of Heidemann and Khalil [42]. The PR-EoS was used to calculate the ethane + n-heptane mixture. TheSRK-EoS was used to calculate the nitrogen + methane mixture. Compressibility effects play an important role in real-gas flows and flows with high Mach numbers.Basically, three different definitions are used to define compressibilities in thermodynamics, the isother-mal compressibility ψ T , the isentropic compressibility ψ s and the isenthalpic compressibility ψ h , seeTab. 3. These three compressibilities can be related to each other by the Gr¨uneisen parameter [43] seeTab. 3, φ = v ∂p∂e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v = vc v ∂p∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v = ρT ∂T∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s (13)7 able 3: Definition of thermodynamic parameters. Name Symbol EquationIsothermal compressibility ψ T − v ∂v∂p (cid:12)(cid:12) T = ψ s (1 + T α p φ )Isentropic compressibility ψ s − v ∂v∂p (cid:12)(cid:12) s Isenthalpic compressibility ψ h − v ∂v∂p (cid:12)(cid:12) h = ψ s (1 + φ )Thermal expansivity α p v ∂v∂T (cid:12)(cid:12) p which can also be expressed by different general thermodynamic response functions [44] as: φ = vα p c v ψ T = α p a s c p . (14)In Equations (13) and (14) e is the internal energy, c v is the isochoric heat capacity and c p is the isobaricheat capacity. In fluids without hydrogen bondings, as it is the case for this study, the Gr¨unsteinparameter is always greater than zero because this type of fluids are having no density anomalies like,e.g., water has [44]. As a result, one can deduce from the equations in Tab. 3 that ψ T > ψ h > ψ s . (15)It is very important to distinguish between these three compressibilities and to choose the suitableone depending on the application and numerical framework. In our case we are using the isenthalpiccompressibility in the pressure equation as the energy conservation equation is solved in terms of anenthalpy equation. The second compressibility being highly important within this study is the isentropiccompressibility because the speed of sound a s , which is an important ingredient of the KT scheme, is afunction of the inverse of the isentropic compressibility ψ s and the density ρ : a s = (cid:118)(cid:117)(cid:117)(cid:116) ∂p∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s,z i = (cid:118)(cid:117)(cid:117)(cid:116) − v ∂p∂v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s,z i = 1 √ ρ ψ s . (16)In single phase systems the calculation of the speed of sound and the isenthalpic compressibility isalmost straightforward, although it implies the calculation of some partial derivatives, which can be solvedanalytically for the cubic EoS used in this investigation, see, e.g., M¨uller et al. [16] for the PR-EoS. Inmulti-phase systems this is not true anymore and the evaluation gets way more complicated. This is alsotrue for measurements and therefore many researchers focus on the numerical investigation of the speed ofsound in multi-component, multi-phase systems, e.g., Picard and Bishnoi [45], Firoozabadi and Pan [46],Nichita et al. [47] and Castier [48]. A simplified and very popular method for calculating the speedof sound in two-phase systems is the correlation of Wood [49]. This correlation is based on a volume-weighted approach of the bulk modulus respectively its inverse, the isentropic compressibility ψ s = (1 − β v ) ψ s, l ( x , p, T ) + β v ψ s, v ( y , p, T ) = (1 − β v ) ρ l a s, l + β v ρ g a s, g (17)where β v denotes the vapor volume fraction. Applying Eq. (17) to Eq. (16) yields Wood’s correlationfor the speed of sound in a two-phase mixture: a s = (cid:34) ρ (cid:32) (1 − β v ) ρ l a s, l + β v ρ g a s, g (cid:33)(cid:35) ( − / . (18)This correlation is suitable to capture the basic trends of the variation of the speed of sound in the two-phase region but misses the abrupt/discontinues changes at the phase boundaries, see Nichita et al. [47] fora more thorough discussion. Different approaches can be found in literature to numerically/analyticallycalculate the thermodynamic speed of sound in two-phase mixtures [45, 46, 47, 48]. In this study we8ollow the approach devised by Nichita et al. [47] as it is a general and efficient numerical scheme [48].Here we are able to use the already implemented T pn -flash as a basis to calculate the thermodynamicspeed of sound.In a single-phase fluid the isothermal and the isentropic compressibility are related by the Gr¨uneisenparameter, see Tab. 3, but can also be related by the ratio of the specific heat at constant pressure toconstant volume κ : ψ s = ψ T κ . (19)By applying some basic mathematical and thermodynamical relations Eq. (19) can be rewritten as [47] ψ s = ψ T − α p T vc p , (20)where α p is the thermal expansivity, see Tab. 3. This basic derivation holds also for two-phase systemsif one takes into account the phase split when evaluating the different properties and partial derivatives,for more details see Nichita et al. [47]. From the Eq. (20) and Tab. 3 it is obvious that Eq. (20) onlycontains partial derivatives at constant temperature and pressure, respectively. Therefore, we followedthe suggestion of Nichita et al. [47] and calculate this partially derivatives numerically by applying T pn -flashes at ( p , T ± (cid:15) ) and ( p ± (cid:15) , T ). These flashes converge very fast within 1-2 Newton iterations, sincean excellent initial guess is available from the hpn -flash performed during every iteration/time step (ifthe mixture was characterized unstable by the TPD-analysis).For calculating the isenthalpic compressibility ψ h we follow the same idea and rewrite the basicdefinition in such a way that it only contains partial derivatives with respect to temperature and pressure.This can be done by using the chain rules ∂v∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h ∂p∂h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v ∂h∂v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p = − ∂T∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v ∂p∂v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T ∂v∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p = − , (22)the total derivative of the enthalpy with respect to temperature and pressured h = ∂h∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p d T + ∂h∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T d p , (23)the Gibbs equation d h = T d s + v d p (24)as well as the Maxwell identity ∂s∂p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T = − ∂v∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p . (25)Applying Eqs. (21)-(25) to the definition of the isenthalpic compressibility, see Tab. 3, yields ψ h = 1 v − c p ∂v∂p (cid:12)(cid:12) T / ∂v∂T (cid:12)(cid:12) p + v − T ∂v∂T (cid:12)(cid:12) p c p / ∂v∂T (cid:12)(cid:12) p (26)and therefore an equation containing only partial derivatives with respect to temperature and pressure.As a result, the speed of sound a s and the isenthalpic compressibility ψ h in a multi-component, two-phasesystem can be calculated numerically by solving the T pn -flashes at ( p , T ± (cid:15) ) and ( p ± (cid:15) , T ) after the hpn -flash-problem has been solved successfully. This makes it numerically efficient compared to iso-energy,iso-volume or iso-entropy flashes because those flashes always require a nested loop or a further objectivefunction in addition to the equality of fugacities in the T pn -flash.9
40 160 180 200 220 240 260102030405060708090 CP T [K] p [ b a r ]
130 150 170 190 210 230 250 2702004006008001 , , T [K] a s [ m / s ]
130 150 170 190 210 230 250 2700 . . . . T [K] ψ s [ b a r ]
130 150 170 190 210 230 250 2700 . . . . T [K] ψ h [ b a r ]
130 150 170 190 210 230 250 2700 . . . . T [K] ψ T [ b a r ]
130 150 170 190 210 230 250 2700 . . T [K] φ [ - ] a) b) c)d) e) f)
10 bar 30 bar 50 bar 70 bar Thermodynamics (PR-EoS) Wood
Figure 4: Investigation of the 14 component Prudhoe Bay mixture [45] at four different pressure levels (10 bar, 30 bar,50 bar and 70 bar): a) Vapor-liquid equilibrium, b) Speed of sound (solid: thermodynamic calculation, dashed: Wood’scorrelation [49]), c) Isentropic compressibility, d) Isenthalpic compressibility, e) Isothermal compressibility and f) Gr¨uneisenparameter. For all diagrams the PR-EoS was used.
Different research groups [45, 47, 48] have used the Prudhoe Bay gas mixture to study the vari-ation of the speed of sound inside the two-phase region by means of numerics and by applying thePR-EoS to model the real-gas behavior. We will also follow their example in this study but will gobeyond this and also discuss additional thermodynamic characteristics like the isenthalpic compress-ibility and the Gr¨uneisen parameter, which relates the different compressibilities to each other. ThePrudhoe Bay mixture contains 14 components with the following mole fractions listed in parenthe-ses [45]: methane (83.3310), ethane (9.6155), propane (3.5998), iso-butane (0.3417), n-butane (0.4585),iso-pentane (0.0403), n-pentane (0.0342), n-hexane (0.0046), n-heptane (0.0003), n-octane (0.0001),toluene (0.0002), nitrogen (1.4992), oxygen (0.0008) and carbon-dioxide (1.0738). In Fig. 4 a) the VLEof this 14 component mixture is shown together with the four pressure levels used for the more thor-ough investigation of this mixture. The shape of the VLE is in very good agreement with the results ofPicard and Bishnoi [45]. However, the critical point is estimated at a different location (Picard and Bish-noi [45]: 206 K, 5.5 MPa; present study: 226.3 K, 7.48 MPa). The speed of sound shown in Fig. 4 b) is invery good agreement with the results of other research groups [45, 47, 48]. At the bubble-point line (lowtemperatures) a drastic, discontinuous jump is observed, which reduces with increasing pressure. Forthe lowest pressure (10 bar) the speed of sound is almost dropping by a factor of 30 as the fluid crossesthe coexistence-line and first vapor bubbles are forming. Afterwards, as the fluid is heating up and thevapor fraction is steadily rising, the speed of sound is increasing continuously. A small discontinuity isalso observed at the dew-point line (large temperatures). This jump gets larger with increasing pressurewhich is contrary to the behavior of a s at the bubble-point line. Outside of the two-phase region the fluidbehaves liquid-like to the left of the VLE resulting in an increase of a s with reducing temperature. Onthe right side of the VLE the opposite is true as the fluid behaves gas-like. The comparison of Wood’s10orrelation [49] with the thermodynamic speed of sound in Fig. 4 b) underlines the statement that thecorrelation is able to capture the basic trend within the two-phase region but misses the significant dropsof the speed of sound at the phase boundaries, compare Nichita et al. [47] for similar findings. This ob-servation holds also for the isentropic compressibility as it behaves indirectly proportional to the speedof sound, see Fig. 4 c) and Eq. (16).In order to do a further evaluation of the numerical approach described above, we also show the differentcompressibilities ( ψ s , ψ h and ψ T ) in Fig. 4. The comparison of the plots shows that the basic relationbetween the three compressibilities ψ T > ψ h > ψ s , see Eq.(15), holds. Furthermore, Figs. 4 c)-e) showthat the isenthalpic and the isentropic compressibility are very similar within the VLE. In contrast, theisothermal compressibility is much larger than the other two ones (note the scaling of the y -axis by afactor of two). The Gr¨uneisen parameter in Fig 4 f) shows a very similar pattern compared to the speedof sound outside as well as inside of the VLE. The identical pattern of a s and φ was also reported byArp et al. [43], who found this kind of behavior for single component fluids. As an overall result of thisdiscussion we conclude that the thermodynamic framework is able to calculate both the speed of soundas well as the isenthalpic compressibility inside the two-phase region.
4. Validation of the VLE-model in OpenFOAM
Recently, experiments and numerical simulations focusing on the process of phase separation dueto multi-component mixing at high-pressure conditions were conducted by Traxinger et al. [6]. In theexperiments n-hexane was injected into a pressurized chamber with pure nitrogen ( p ch = 50 bar, T ch =293 K) at rest. The injected fluid was heated to three different temperatures such that the test cases coverregimes in which thermodynamic non-idealities, namely phase-separation induced by multi-componentmixing at supercritical pressure with respect to the pure components values, are significant. To capturethe flow structure as well as the phase separation in the experiments, shadowgraphy and elastic lightscattering were performed simultaneously. For additional information about the experimental setup seeLamanna et al. [50] and Baab et al. [51]. In addition to the experiments, large eddy simulations werecarried out applying a VLE-model. The results were in qualitatively very good agreement with theexperimental findings in terms of phase separation and overall flow structure. In order to validate thethermodynamic modeling coupled with the hybrid solver, numerical simulations with the same boundaryconditions were performed, see Tab. 4. As the boundary conditions for the experimental setup wereonly given as static inlet properties, the corresponding total plenum conditions were calculated underthe constant isentropic and isenthalpic nozzle flow assumptions [10]: s = s ( p , T ) = s ( p, T ) h = h ( p , T ) = h ( p, T ) + 12 u . (27)For the three test cases T600, T560 and T480, named after the total temperature T t,inlet of theinjected n-hexane, the resulting total pressures are 56 .
40 bar, 56 .
57 bar and 54 .
56 bar. Case T600 withthe highest injection temperature of 600 K is likely to be a dense-gas jet and is therefore showing nophase separation effects. For a decreased temperature of 560 K minor two-phase effects can be expectedwhile a temperature decrease to 480 K will show strong two-phase effects in form of a spray-like jet [6].
Table 4: Overview of the numerical boundary conditions used in the validation cases.
Case EoS p t,fuel p fuel T t,fuel T fuel p ch T ch [bar] [bar] [K] [K] [bar] [K]T600 PR 56.40 50.0 600.0 595.0 50.0 293.0T560 56.57 50.0 560.0 554.8 50.0 293.0T480 54.56 50.0 480.0 479.3 50.0 293.0The numerical domain used for the validation case is a two degree rotationally symmetric sectionof a single hole injector with a diameter of D = 0 .
236 mm connecting the high pressure side (here:n-hexane) with the low pressure side (here: nitrogen), see Fig. 5. The low pressure side was modeledusing a constant volume chamber with L x /D = 150 and L r /D = 100. A structured mesh is used andthe coordinate system is located at the nozzle exit aligned with the flow direction, starting at x = 0.11nitially, both sides (low and high pressure) are separated by a membrane located at x = − . D whichbursts at simulation start resulting in the development of an expanded/underexpanded jet (dependingon the prescribed pressure ratio Π = p t,fuel /p ch ) into the positive x -direction. As this validation caseis almost isobaric [6], a low Mach number flow and a weakly expanded jet can be expected. Hence, thefocus will be on the overall flow structure and the two-phase effects in the far field and a resolution of∆ r/D = 20 in the radial nozzle direction was therefore chosen. The geometric specifications as well asthe radial mesh resolution are summarized in Tab. 5. As in Traxinger et al. [6], we are using the PR-EoSfor the simulation of these three test cases. This is in good agreement with our findings in section 3 as theprediction of the PR-EoS are superior for high number alkanes like n-hexane and n-heptane compared tothe SRK-EoS which gives better predictions for short alkanes, see Fig. 2. Traxinger et al. [6] furthermoreshow that the prediction of the PR-EoS is also excellent for the binary VLE of n-hexane and nitrogenat the pressure considered for these test cases and that not fitting of the binary interaction parameterin Eq. (9) is necessary. x -axis x =0 D /2 L x L r Dp high -chamber p low -chamberInit. Figure 5: Schematic of the numerical domain used both in the validation case and in the test cases. The low and highpressure chamber are initially separated.
Figure 6 shows the results of the experimental studies [6] on the left compared to snapshots of thenumerical simulations performed in this work for the three cases T600, T560 and T480 on the right. Eachframe in the left column represents a shadowgraphy image superimposed by the resulting elastic lightscattering image (bottom half). In the right column the frame represents the temperature contour plotsuperimposed by the vapor mass fraction β in the bottom half. For all three test cases the agreementof the numerical results with the experimental findings is remarkably good as it was also observed byTraxinger et al. [6] for the conducted LESs. The experimental results for the temperature T t,fuel = 600 Kshow no significant scattering signal and therefore Traxinger et al. [6] deduced a single-phase state for thisinjection conditions. The same can be concluded from the numerical simulation on the right, where noregion with phase separation is visible (indicated by a vapor mass fraction β between zero and one) whichin turn leads to the conclusion that the TPD-analysis was always positive for every single-phase statetested. Decreasing the total temperature to 560 K leads to a very dark jet in the shadowgraphy imageindicating the formation of dense droplet clouds. The increase in scattering intensity over some ordersof magnitude confirms this statement as this measurement technique is very sensitive to particle sizesrespectively (small) droplets. On the right the URANS simulation shows also a two-phase region with avapor mass fraction of 0 .
992 to 1. As also reported by Traxinger et al. [6], the area of phase separationin the numerical results coincides very well with regions of high scattering signal in the experiments.The third test case at T t,fuel = 480 K shows a spray-like characteristic and therefore strong two-phaseeffects in the shadowgram as well as in the scattering. The scattering for this test case is comprisingseveral orders of magnitude and occupies almost the entire jet domain. In addition, the shadowgramshows again a very dark jet but this time with distinct droplets visible in the outer jet region underliningthe spray-like character. The numerical simulation confirms this finding. The mixture-induced phaseseparation is almost covering the complete jet and β is ranging from 0 in the inner jet region to 1 in theouter jet areas. Table 5: Geometry specifications for the validation case and the numerical study.
Case
D L x /D L r /D ∆ r/D [mm] [-] [-] [-]Validation 0.236 150 100 20High-pressure injection 1.5 80 70 4012ll the findings for the three different test cases are in very good agreement with the LESs conductedby Traxinger et al. [6]. We therefore conclude that our numerical framework consisting basically of thehybrid, pressure-based solver and the VLE-model is able to capture the main two-phase effects in weaklyexpanded and isobaric jets. Hence, we can use this framework to simulate highly underexpanded jetswith probable phase separation.Experiment01020 Simulation0 30 60 90-20-10 0 30 60 9001020 0 30 60 90-20-10 0 30 60 9001020 0 30 60 90-20-10 x/D [-] 0 30 60 90 x/D [-] r / D [ - ] r / D [ - ] r / D [ - ] a ) C a s e T b ) C a s e T c ) C a s e T T [K]293 595 T [K]293 554 . β [ − ]0 .
992 1 T [K]293 479 . β [ − ]0 1˜ I [Cnts / mJ]10 . ˜ I [Cnts / mJ]10 . ˜ I [Cnts / mJ]10 . Figure 6: Comparison of experimental (left, figures taken from Traxinger et al. [6]) and numerical results (right, presentstudy) for the different test cases T600, T560 and T480.
5. Engine-relevant Test Cases
The numerical investigation of the injection process of gaseous fuels at engine-relevant conditionsare the focus of this study. As detailed experiments and/or simulations are lacking in literature, ownsetups and operating conditions were defined as follows. Based on the validation case with a chamberpressure of 50 bar and the recent investigations of Banholzer et al. [10], who simulated high-pressureinjections of hydrogen ( p t,fuel = 500 bar) into air at rest ( p ch = 100 bar), the test cases were chosensuch that application-relevant conditions and penetration depths are met. As the purity of hydrogenis not given for industrial applications, nitrogen with a volume fraction of 0 . . ), ethane (C H ) and nitrogen (N ) [52, 53, 54]. Concluding the CNG-compositions available in publications a generic composition was defined, see Tab. 6. Furthermore, thepressure levels were adapted to engine-relevant conditions. Typical chamber pressures range from 5 bar13o 80 bar, depending on the ICE-type and the start of injection (SOI). As the validation case was carriedout with a chamber pressure of 50 bar the same pressure value was chosen for the test cases. Similarconsiderations lead to a maximum total fuel pressure of 600 bar [55, 56, 57, 58, 59]. For a comparisonof different pressure ratios additional simulations were performed with a total fuel pressure of 300 bar,resulting in pressure ratios Π = p t,fuel /p ch of 12 and 6. A total of four numerical simulations werecarried out, see Tab. 7. Table 6: Volume fractions [%] of CHG- and CNG-compositions.
Fuel Methane (CH ) Hydrogen (H ) Ethane (C H ) Nitrogen (N )CNG 95.8 0.0 2.6 1.6CHG 0.0 99.9 0.0 0.1The schematic of the numerical domain used for the test cases is shown in Fig. 5 and is similarto the one for the previous validation case, except for the dimensions. For the diameter D a value of1 . L x /D = 80 and L r /D = 70.The computational setup of Hamzehloo et al. [60, 61, 62] and a grid independence study performed byBanholzer et al. [10] lead to a resolution of ∆ x = ∆ r = D/
40 = 0 . Table 7: Overview of the boundary conditions for the four test cases.
Case EoS Fluid Fluid p t,fuel p ch Π T t,fuel T ch p fuel -chamber p ch -chamber [bar] [bar] [-] [K] [K]CNG-p600 SRK CNG Nitrogen (N ) 600 50 12 300 300CNG-p300 300 50 6CHG-p600 SRK CHG Nitrogen (N ) 600 50 12 300 300CHG-p300 300 50 6 At their rest conditions both multi-component fuels CHG and CNG are in a single-phase state. As weare expecting strongly underexpanded jets due to the large prescribed pressure ratios, two different waysare feasible how the formation of a vapor-liquid equilibrium of the appropriate mixture could take place:Firstly, by an strong expansion lowering pressure and temperature at constant overall feed compositionand secondly, by nitrogen dilution of the fuel and a simultaneous expansion.In the CHG test cases two fluids are forming the multi-component mixture, hydrogen and nitrogen,see Tab. 6, whereby hydrogen is the low volatile and nitrogen the high volatile fluid. Hence, nitrogenis defining the upper bound of the multi-component critical locus in terms of temperature at its criticaltemperature 126 .
19 K. Above this temperature no phase separation is possible and therefore a strongexpansion and reduction of temperature is necessary to enter the multi-component VLE. In Fig. 7 someselected vapor-liquid equilibria at relevant temperatures (100 . . . , respectively, has to mix into the CHG to en-ter the respective VLE. Hence, the probability of phase separation can be assessed low for the CHG cases.14 0 . . . .
50 bar x H , y H , z H [-] p [ b a r ] T = 100 . T = 107 . T = 113 . Figure 7: Vapor-liquid equilibria of hydrogen and nitrogen at three different temperatures. The solid lines are calculatedusing the SRK-EoS. The dots are representing experimental data [63, 64, 65].
For the CNG the VLE and the probability of phase separation is significantly different. Firstly, CNGis defined as a ternary mixture consisting of methane, ethane and nitrogen, see Tab. 6, and thereforethe mixture space is increased by one dimension compared to binary mixtures. Secondly, nitrogen isnow the low volatile component because methane and ethane are having a significantly higher criticaltemperature, see Tab. 2. Therefore, the critical locus has in the case of CNG an upper bound at 305.32 Kwhich is the critical temperature of ethane. This temperature in turn is almost three times larger thanthe critical temperature of nitrogen which is the high volatile component in the CHG mixture. Thehigher maximum critical temperature makes phase separation during the CNG injection significantlymore probable than in the CHG-case because the injection temperature and the critical temperature ofethane are almost equal (Recall: For the CHG-case an expansion below 126.19 K is at least necessary toenter the binary VLE.). Figure 8 shows VLEs of the ternary mixture at three different pressures (40 bar,20 bar and 8 bar) and two different temperatures (200 K and 160 K). The dew-point lines are shown inred and the bubble-point lines in blue while the experiments are represented by dots and the solid linesare calculated using the SRK-EoS. The feed composition of CNG is marked by the white dot in each plot.In order to give a more detailed glance onto these ternary VLE-diagrams, we also plotted the relatedbinary mixtures into pxy -diagrams, see Fig. 3 c for nitrogen + methane, Fig. 9 a for nitrogen + ethaneand Fig. 9 b for methane + ethane. In the appropriate diagrams the three investigated pressure levels(8 bar, 20 bar and 40 bar) are marked with dashed lines and the chamber pressure of 50 bar is markedwith a dashed dotted line. Especially in Fig. 9 a, the deviations between the experimental data and thecubic EoS at the bubble-point line become very obvious, where cubic EoS are always facing considerableinaccuracies. With decreasing temperature the prediction error is rising steadily. Similar findings canbe made for the binary mixture nitrogen + methane, but compared to the nitrogen-ethane mixture theoverall prediction accuracy is acceptable over the investigated pressure and temperature range, comparealso the very good prediction of the critical locus in Fig. 3 d. An excellent prediction is evident for thebinary mixture methane + ethane over the relevant pressure and temperature range. These findings forthe different binary mixtures transfer directly to the VLEs of the ternary mixture CNG. Apart from the40 bar and 160 K diagram, see Fig. 8 b, the comparison to experimental data is very good, see Figs. 8 a, cand d. Especially the dew-point line is predicted excellently by the SRK-EoS and the already mentionedproblems in terms of prediction accuracy are obvious at the bubble-point line. This in turn is not abig issue in this investigation as most of our two-phase states will occupy the left corner of the ternaryVLE-diagram, where the prediction of the SRK-EoS is very good. This can be seen in detail in Fig. 8 awhere the feed composition of the CNG and the dilution process of CNG with nitrogen are shown. Asin our case the VLE will be entered across the dew-point line due to the strong expansion process, wecan expect an excellent prediction concerning the onset of phase separation.1520406080100 0 20 40 60 80 100 020406080100
Dilution of CNGwith nitrogen. D e w - p o i n t li n e B u b b l e - p o i n t l i n e CNG a) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 40 bar, T = 200 K 020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B ubb l e - p o i n t li n e CNG b) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 40 bar, T = 160 K020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B u b b l e - p o i n t l i n e CNG c) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 20 bar, T = 200 K 020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B ubb l e - p o i n t li n e CNG d) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 20 bar, T = 160 K020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B ubb l e - p o i n t li n e CNG e) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 8 bar, T = 200 K 020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B u b b l e - p o i n t l i n e CNG f) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 8 bar, T = 160 K Figure 8: VLEs for the ternary mixture consisting of methane, ethane and nitrogen at two different temperatures (160 Kand 200 K) and three different pressures (40 bar, 20 bar and 8 bar). The CNG mixture as defined in Tab. 6 is marked witha white dot. The solid lines are calculated using the SRK-EoS. The dots are representing experimental data [66]. . . . . x N2 , y N2 , z N2 [-] p [ b a r ] T = 149 .
82 K T = 172 .
04 K T = 194 .
26 K 0 0 . . . . x CH4 , y CH4 , z CH4 [-] p [ b a r ] T = 144 .
26 K T = 158 .
15 K T = 172 .
04 K T = 186 .
09 K T = 199 .
93 K a) b)
Figure 9: VLEs for binary mixtures: a) nitrogen + ethane and b) methane + ethane. The solid lines are calculated usingthe SRK-EoS. The dots represent experimental data [40, 67].
Using the one-dimensional flow analysis presented in Banholzer et al. [10] the choked flow propertiesof the two fuels for the different operating conditions can be calculated, see Tab. 8. The larger criticalpressure p ∗ of CHG leads to an earlier choking of the nozzle compared to the CNG cases. As we areusing a constant chamber pressure of 50 bar for all four test cases the flow in the nozzle is always choked.For this type of flow the maximum mass flow rate through the orifice can be evaluated as˙ m ∗ n = ρ ∗ a ∗ s A (28)with A being the cross-sectional area of the nozzle. Reduction of the cross section of the vena contractadue to detachments caused by the shape of the nozzle and friction effects inside the nozzle can decreasethe mass flow by up to 20 % [68]. The reduction of the achievable mass flow is known as dischargecoefficient and is defined as C d = ˙ m/ ˙ m ∗ n . For the lower fuel pressure of 300 bar the maximum mass flowof CNG is 3 . ρ ∗ of CNG compared to CHG. Doubling the pressure to 600 bar leads to an increase of the critical massflow rate of nearly a factor of 2 for both fuels. Table 8: Choked nozzle flow properties for CHG and CNG ( T t,fuel = 300 K, SRK-EoS). Fuel p t,fuel p ∗ ρ ∗ a ∗ s ˙ m ∗ n ˙ M ∗ n [bar] [bar] [kg m − ] [m s − ] [g s − ] [kg m s − ]CNG 600 191.37 220.55 556.42 216.86 120.67300 124.66 148.87 439.36 115.59 50.78CHG 600 279.61 23.87 1469.18 62.02 91.11300 148.09 13.65 1330.11 32.09 42.68The main influencing parameter of the penetration depth of the jet is the momentum flux beingpresent at the nozzle exit [69, 70, 71]. The momentum flux in a choked flow is defined as:˙ M ∗ n = ρ ∗ a ∗ s A . (29)For the lower fuel pressure of 300 bar the value evaluated for the CNG-case is higher than for the CHG-case. The difference of the two fluids is not as big as in case of the mass flux because the momentum17ux is proportional to the squared speed of sound. For the fuel pressure of 600 bar the momentum fluxincreases by approximately a factor of 2. Comparing the values for both fuels and operating points itis presumed that the CNG-jet with 600 bar (CNG-p600) penetrates the fastest, followed by CHG-p600,CNG-p300 and CHG-p300.
6. Results
Fig. 10 shows snapshots of the temporal evolution of the high-pressure injection of CNG with a totalpressure of 600 bar and a total temperature of 300 K into the pressurized chamber filled with nitrogenat rest ( p ch = 50 bar , T ch = 300 K). For the detailed discussion we are only showing the first 20 %of the simulated chamber which is in our opinion sufficient to discuss all the main fluid-dynamic andthermodynamic effects happening during the transient part of the injection process. The left columnrepresents the temperature distribution superimposed by the vapor mass fraction β . In the right columnthe pressure field is superimposed by the Mach number, in which regions with Ma < t = 5 µ s r / D [ − ] expansion fan t = 10 µ s traveling spherical vortex r / D [ − ] t = 20 µ s single-phase state r / D [ − ] t = 40 µ s r / D [ − ] intercepting shock/sharp boundaryFig. 11 t = 60 µ s r / D [ − ] t = 80 µ s r / D [ − ] t = 250 µ s x/D [ − ] r / D [ − ] x/D [ − ] T [K]150 350 β [ − ]0 . . p [bar]10 50 Ma [ − ]1 . . Figure 10: Snapshots of the temporal evolution of the underexpanded jet for the test case CNG-p600 ( p t,fuel = 600 bar).In the left column the temperature and the vapor mass fraction β are superimposed to visualize regions where the jet entersthe VLE. In the right column the superposition of pressure and Mach number is shown to get a more detailed glance ofthe supersonic flow accompanying the jet. µ s) the fluid enters themixing chamber and a Prandtl-Mayer expansion fan is formed at the nozzle edge. The flow acceleratesacross the expansion fan and therefore the Mach number increases while the static pressure and tem-perature drop inside the jet. As the jet penetrates further into the chamber a traveling spherical vortexforms in front of the jet. This spherical vortex is growing steadily and breaks up at approximately 60 µ s.The fluid inside this vortex gets pushed into the mixing chamber by the subsequent steady-state jetregion which is apparent in the Mach number plots at t = 60 µ s and t = 80 µ s. Apart from this dominantvortex structure shock phenomena can be observed during the injection process which can be seen forthe first time at 40 µ s. The expansion fans being present at the outer part of the jet interact with thefree jet boundary leading to weak compression waves and an intercepting shock, which distinguishes thejet from the surrounding nitrogen. This intercepting shock in turn ends in a strong curved shock, theMach disk being a very prominent feature in such underexpanded jets. In our case, this concise shock israther weak/oblique and therefore the flow remains at supersonic conditions (Ma > µ s wherea total of four shock barrels is visible and all flow properties are only fluctuating around their mean values.From the left column in Fig. 10 it is obvious that larger parts of the jet have entered the VLE andtherefore the underexpanded jet features pronounced two-phase effects which can be identified from thevapor mass fraction β ranging from 0 . µ s a first small region with phase separation has formedvery close to the nozzle trailing edge resulting from the decrease of the static pressure and temperatureinside the Prandtl-Mayer expansion fan. Related to this is a discontinuous drop in the speed of sound andin turn a non-linear increase in the Mach number from approximately unity to Mach numbers rangingfrom 2 up to 3. The subsequent expansion leads to a further decrease of pressure and temperature andat 10 µ s a large part of the jet has entered the VLE. This trend continues until the formation of a firstshock barrel sets in. The start of the formation of this shock structure can be seen in terms of a smallsingle-phase region close to the centerline at t = 20 µ s. At t = 40 µ s the first shock barrel is already fullydeveloped and can be seen at x/D <
4. Due to the shock and the accompanying increase in temperatureand pressure the jet gets back to a single-phase state upstream of the shock barrel. As the flow isaccelerated after this shock, the continuous expansion leads again to a phase instability and therefore toa phase separation. This phenomenon happens several times during the propagation of the jet into thechamber and leads to four distinct shock barrels containing one large region of two-phase flow each at t = 250 µ s. Outside of this region with strong expansion large areas of phase separation are also visiblein the traveling spherical vortex. Compared to the jet center, the vapor mass fraction is larger and moreinhomogeneous, i.e. covering a wider range. This fact can be explained by taking a deeper look in thephase separation process itself and how it is triggered in these two different regions of the jet. This willbe done with results from the snapshot at t = 40 µ s, where we will focus on a small region at the jettip marked by the red box in Fig. 10 where both phase separation regions are present (jet center andtraveling spherical vortex). In Fig. 11 the detailed view of this region is visualized in terms of the staticpressure p , the static temperature T , the methane mass fraction Y CH , the vapor mass fraction β , thespeed of sound a s , the Mach number Ma and the partial derivative ∂ρ/∂p (cid:12)(cid:12) h .Beginning on the left, the first snapshot reveals the large pressure gradient being present at the tipof the underexpanded jet. The minimum pressure of the first expansion wave is approximately 8 barwhile the compression process in front of the jet increases the fluid-pressure up to 80 bar being 1.6 timeslarger than the prescribed chamber pressure of 50 bar. Due to this strong expansion and compression,the static temperature ranges approximately from 150 K to 350 K (second frame) showing a patternsimilar to the static pressure. In the third frame the mass fraction of methane is shown which is themain component of CNG, see Tab. 6. At 40 µ s almost no nitrogen has mixed into the jet core yetwhich is therefore almost at the feed composition. In the region of the traveling spherical vortex theopposite is true as the convection inside of this vortex is forcing the dilution of CNG with the surroundingnitrogen. As a consequence, two different phase separation regions can be distinguished during the CNG19 p [bar] 80150 T [K] 350 0 Y CH4 [ − ]0 . . β [ − ] 1 . a s [m s − ]410 0 Ma [ − ] 4 8 e − ∂ρ/∂p (cid:12)(cid:12) h [s / m ]4 e − Figure 11: Detailed view of different flow properties for the CNG-p600 test case at 40 µ s after injection. Snapshots aretaken at the tip of the jet including both the jet center as well as the traveling spherical vortex, for more details see Fig. 10red box. The yellow and orange boxes mark the sampling points for Fig. 12. injection process, a phase separation solely caused by expansion and a phase separation triggered byexpansion and nitrogen dilution. In order to emphasize this finding we have sampled two appropriateareas of data points at the positions marked by the yellow and orange box in Fig. 11 and plotted themin ternary mixture diagrams, see Fig. 12. The constant temperature and pressure in these diagrams020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B ubb l e - p o i n t li n e CNG Expansion inducedphase separation a) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 8.9 bar, T = 147.3 K 020406080100 0 20 40 60 80 100 020406080100 D e w - p o i n t li n e B ubb l e - p o i n t li n e CNG Mixture inducedphase separation b) x N , y N , z N [ % ] x C H , y C H , z C H [ % ] x C H , y C H , z C H [%] p = 28.5 bar, T = 167.2 K Figure 12: Ternary phase diagrams for the CNG test-case CNG-p600. The colored points mark the sampling data high-lighted in Fig. 11 by the appropriate colored boxes. The solid lines as well as the sampling data for the numerical simulationsat 40 µ s are both calculated with the SRK-EoS. were chosen such that a large number of meaningful points can be plotted. In Fig. 12 a) the scatteringdata of the expansion induced phase separation (yellow) is plotted together with the VLE. Almost nomixing with nitrogen has taken place as already discussed and therefore all data points are groupingaround the feed composition. Due to the strong expansion occurring at the jet center the VLE is almostoccupying the complete mixture space and not surprisingly all points in the CFD have entered theregion of single-phase instability. For the sample data taken from the traveling spherical vortex the soleexpansion process would not have been sufficient to lead to a phase separation. From Fig. 12 b) it getsobvious that a dilution of the CNG with the surrounding nitrogen is necessary to enter the VLE at theappropriate temperature and pressure. Due to the convective mixing inside the vortex the individualmixtures contain more than 20 mole-% nitrogen and are therefore able to enter the VLE. Going backto Fig. 11 one can now reconstruct the snapshot for the vapor mass fraction β shown in the fourthframe. In the jet center the variation in temperature, pressure and composition is low and therefore analmost constant β -value can be observed. Inside the vortex a wider range of β -values can be seen due tothe steady dilution/mixing with the surrounding nitrogen and also the wider pressure and temperaturerange. As already discussed earlier, the speed of sound a s (fifth frame), the Mach number Ma (sixthframe) and the partial derivative ∂ρ/∂p (cid:12)(cid:12) h exhibit a discontinuous jump as the mixture enters the VLE.The speed of sound drops significantly and the Mach number defined as Ma = u/a s increases to valuesup to 4. 20 = 20 µ s 600 bar t = 60 µ s 600 bar300 bar -3-2-1
300 bar t = 40 µ s 600 bar t = 80 µ s 600 bar300 bar x/D [ − ]
300 bar x/D [ − ] T [K]150 350 β [ − ]0 . . r / D [ − ] r / D [ − ] Figure 13: Comparison of the near-nozzle flow structure between the two CNG test-cases CNG-p600 (upper half) andCNG-p300 (lower half).
In Fig. 13 the near-nozzle flow structures of the direct injection of CNG for both total fuel pressures600 bar and 300 bar are compared in terms of snapshots at different times after injection. The snapshotsare a superposition of the temperature field with the corresponding vapor mass fraction. First of all,the overall flow-structure is very similar for both cases concerning fluid-dynamic and thermodynamiceffects. For the case CNG-p300 a weaker expansion due to the lower pressure ratio of Π = 6 can be seencompared to the case CNG-p600. Nevertheless, the expansion process is still sufficient to lead to theinstability of the single-phase and to enter the VLE. The weaker expansion manifests itself directly ina larger vapor mass fraction compared to the CNG-p600 case. In addition, the comparison shows thatthe traveling spherical vortex is also visible and shows phase separation effects, but in contrast to thehigh pressure jet the break-up is weakened in its intensity and spatial expansion. Due to the lower totalfuel pressure and therefore a lower density and speed of sound at the nozzle exit, the penetration of theCNG-p300 jet is slower. Furthermore, the length of the stationary shock barrels is smaller by a factor ofalmost two, which correlates well with the predictions of Velikorody and Kudriakov [72].
In the following we will compare the two different fuels, namely CNG and CHG, to each other. In thispart we will focus on the 300 bar cases because the comparison of the 600 bar cases is almost identical.Figure 14 shows the effects of the fuel change from CNG to CHG for a total pressure of 300 bar in terms ofthe temperature field superimposed by the vapor mass fraction, whereby the temperature ranges between150 K and 350 K for the CNG-p300 case and between 100 K and 350 K for the CHG-p300 case. Bothjets show the typical flow structures accompanying underexpanded jets but for the hydrogen case theshock and vortex structures are slightly more pronounced and show a larger spatial extent whereas thelength of the shock barrels is a little bit shorter. Although the expansion of the CHG is much stronger(manifesting in a lower static temperature), the mixture remains in a single-phase state through out thecomplete jet. For the jet center this gets more clear when reconsidering the VLEs shown in Fig. 7. Asnearly no dilution occurs in the jet core and the CHG mixture is almost pure hydrogen, see Tab 6, anexpansion to 100 K is by far not sufficient to enter the VLE of the binary mixture. In the region of thetraveling spherical vortex the dilution of the fuel with the surrounding nitrogen would be sufficient buthere in turn the expansion is too weak to enter the VLE.21 = 20 µ s CNG t = 60 µ s CNGCHG -3-2-1 CHG t = 40 µ s CNG t = 80 µ s CNGCHG x/D [ − ] CHG x/D [ − ] T [K]150 350 β [ − ]0 . . T [K]100 350 r / D [ − ] r / D [ − ] Figure 14: Comparison of the near-nozzle flow structure in terms of snapshots at four different times for CNG (upper half)and CHG (lower half) for a fuel pressure of 300 bar.
In Fig. 15 the normalized jet penetration depth Z tip /D tracked with a mass fraction of Y fuel = 0 . . J ∗ = a ∗ s A is dominant [69] and the CHG-p300-jet penetrates faster thanthe CNG-p300-jet with the same total pressure, see also Fig. 14. In the beginning of the injection thepenetration curve of the CHG-jet shows a linear incline up to 0 .
08 ms. Comparing snapshots of the startand end time of the linear profile the reason for this behavior becomes obvious, being a detachment ofthe jet from the centerline, see Fig. 14 on the right side. A blocking cone is formed in front of the jetconsisting of a higher pressure and therefore higher density region which redirects the main flow awayfrom the centerline to the sides, an effect which was also noticed for high pressure gas injections by otherresearchers, see, e.g., Hamzehloo et al. [61]. After the cone vanishes the jet reattaches to the centerlineand the injection profile is in a very good agreement with the correlation of Hill and Ouellette [70] withthe fitting parameter of Γ = 2 .
968 obtained in Banholzer et al. [10], see dashed lines in Fig. 15. A firstconclusion can be drawn that although the injection is performed at elevated pressures and therefore real-gas effect need to be modeled, the penetration depth for gas-like, single-phase jets can still be estimatedby a rather simple correlation.0 0 . linear incline Z ∝ √ t t [ms] Z t i p / D [ − ] CNG-p300CHG-p300Corr. [ ? ] . . t [ms] Z t i p / D [ − ] Figure 15: Left: Normalized jet penetration depth. Right: Detailed view of the injection profile up to 0 . .
02 ms and 0 .
08 ms, respectively. Thecontour plots are colored using the mass fraction Y CH ranging from 0 (blue) to 1 (red). .
7. Conclusion and Outlook
In order to simulate and discuss real-gas effects and phase separation processes occurring in under-expanded jets at engine-relevant conditions two major numerical challenges were addressed in this work.On the one hand, it is indispensable to derive numerical schemes which are capable of solving highlycompressible flows including accurate predictions of shock- and flow-discontinuities. On the other hand,a consistent, appropriate and efficient modeling of the thermodynamic state considering real-gas effectsand phase instability is needed. Hence, a numerical framework implemented in the open-source toolOpenFOAM is presented in this work combining a hybrid, pressure-based solver with a vapor-liquidequilibrium (VLE) model based on the cubic equation of state (EoS). This framework is used to investi-gate underexpanded jets at engine-relevant conditions. To the authors knowledge such an investigationis carried out for the first time in literature and therefore both numerical as well as experimental studiesare lacking in this field. As a consequence, we strongly focused on a thorough validation of both parts ofthe numerical framework (solver and thermodynamics), whereby the validation of the solver was recentlycarried out by Kraposhin et al. [14] by means of one-dimensional flow problems. Therefore, in the presentstudy the focus was mainly on the validation of the thermodynamics employing general thermodynamicrelations and measurement data available in the literature. In this context, the range of applicabilityof two different EoSs was discussed and the most important parts of the VLE-model were described.Additionally, compressibility effects and the variation of the speed of sound in the VLE-region was ex-amined thoroughly. A large number of measurement data was used to discuss the prediction accuracy ofthe VLE-model. Besides, a recent (experimental and numerical) investigation [6] about mixture inducedphase separation for high-pressure injection in the low Mach number regime was simulated to validatethe numerical framework set up in OpenFOAM. In this sense, three simulations were carried out wheren-hexane is injected into nitrogen at elevated pressure. Based on all these successful validations, engine-relevant simulation cases for two different fuels (CNG and CHG) were defined with total pressures up to600 bar and pressure ratios up to 12. While CHG showed no phase separation effects for the consideredhigh-pressure injections, two-phase regions for the CNG-jets evolved. Analyses revealed that the phaseseparation are caused by two different effects. Firstly, by the strong expansion due to the large pressureratio and secondly, by the mixing of the fuel with the chamber gas. A comparison of the single-phasewith the two-phase jets disclosed that the phase separation leads to a completely different penetrationdepth in contrast to single-phase injection and therefore commonly used analytical approaches fail topredict the penetration depth.In future studies several different topics have to be addressed. Different fluids and injection conditionsadditional to the ones presented in this work are important to be investigated in order to reveal furtherconditions where phase separation may occur. In this context, the difference in penetration depth betweensingle- and two-phase jets has to be further examined. Possible non-equilibrium effects in terms of phaseseparation have to be inquired as they can delay phase separation especially in high Mach number jetsas it is for instance the case in last stages of low-pressure steam turbines.
Acknowledgements
The authors would like to thank Munich Aerospace ( ) and The ResearchAssociation for Combustion Engines eV (FVV, Frankfurt, ) for the funding. Furthermore,we acknowledge the Dortmund Data Bank Software & Separation Technology (DDBST) for their supportin finding appropriate measurement data for the investigated binary and ternary mixtures.23 eferences [1] J. C. Oefelein and V. Yang. Modeling high-pressure mixing and combustion processes in liquid rocketengines.
Journal of Propulsion and Power , 14(5):843–857, 1998.[2] N. Zong, H. Meng, S. Hsieh, and V. Yang. A numerical study of cryogenic fluid injection and mixing undersupercritical conditions.
Physics of Fluids , 16(12):4248–4261, 2004.[3] T. Schmitt, L. Selle, B. Cuenot, and T. Poinsot. Large-eddy simulation of transcritical flows.
ComptesRendus M´ecanique , 337(6-7):528–538, 2009.[4] H. M¨uller, M. Pfitzner, J. Matheis, and S. Hickel. Large-eddy simulation of coaxial ln2/gh2 injection attrans- and supercritical conditions.
Journal of Propulsion and Power , 32(1):46–56, 2015.[5] J. Matheis and S. Hickel. Multi-component vapor-liquid equilibrium model for les and application to ecnspray a. In
Proceedings of the 2016 Summer Program, Center for Turbulence Research, Stanford University,arXiv:1609.08533 , pages 25–34, 2016.[6] C. Traxinger, H. M¨uller, M. Pfitzner, S. Baab, G. Lamanna, B. Weigand, J. Matheis, C. Stemmer, N. A.Adams, and S. Hickel. Experimental and numerical investigation of phase separation due to multi-componentmixing at high-pressure conditions. In
Proceedings of the 2017 ILASS in Valencia , 2017.[7] R. Khaksarfard, M. R. Kameshki, and M. Paraschivoiu. Numerical simulation of high pressure release anddispersion of hydrogen into air with real gas model.
Shock Waves , 20(3):205–216, 2010.[8] F. Bonelli, A. Viggiano, and V. Magi. A numerical analysis of hydrogen underexpanded jets under real gasassumption.
Journal of Fluids Engineering , 135(12):121101, 2013.[9] S. Crist, D. R. Glass, and P. M. Sherman. Study of the highly underexpanded sonic jet.
AIAA Journal ,4(1):68–71, 1966.[10] M. Banholzer, H. M¨uller, and M. Pfitzner. Numerical investigation of the flow structure of underexpandedjets in quiescent air using real-gas thermodynamics. ,2017.[11] M. Kraposhin, A. Bovtrikova, and S. Strijhak. Adaptation of kurganov-tadmor numerical scheme for applyingin combination with the piso method in numerical simulation of flows in a wide range of mach numbers.
Procedia Computer Science , 66:43–52, 2015.[12] R. I. Issa. Solution of the implicitly discretised fluid flow equations by operator-splitting.
Journal ofcomputational physics , 62(1):40–65, 1986.[13] A. Kurganov and E. Tadmor. New high-resolution central schemes for nonlinear conservation laws andconvection–diffusion equations.
Journal of Computational Physics , 160(1):241–282, 2000.[14] M. V. Kraposhin, M. Banholzer, M. Pfitzner, and I. K. Marchevsky. A hybrid pressure-based solver fornon-ideal single-phase fluid flows at all speeds.
International Journal for Numerical Methods in Fluids,Manuscript in preparation. , 2017.[15] M. Jarczyk and M. Pfitzner. Large eddy simulation of supercritical nitrogen jets. In , page 1270, 2012.[16] H. M¨uller, C. A. Niedermeier, J. Matheis, M. Pfitzner, and S. Hickel. Large-eddy simulation of nitrogeninjection at trans-and supercritical conditions.
Physics of Fluids , 28(1):015102–1–015102–28, 2016.[17] B. Chehroudi. Recent experimental efforts on high-pressure supercritical injection for liquid rockets andtheir implications.
International Journal of Aerospace Engineering , 2012:1–31, 2012.[18] M. Oschwald, J. J. Smith, R. Branam, J. Hussong, A. Schik, B. Chehroudi, and D. Talley. Injection of fluidsinto supercritical environments.
Combustion Science and Technology , 178(1-3):49–100, 2006.[19] R. C. Reid, J. M. Prausnitz, and B. E. Poling. The properties of liquids and gases.
McGraw-Hill, New York ,1987.[20] D. Peng and D. B. Robinson. A new two-constant equation of state.
Industrial & Engineering ChemistryFundamentals , 15(1):59–64, 1976.[21] G. Soave. Equilibrium constants from a modified redlich-kwong equation of state.
Chemical EngineeringScience , 27(6):1197–1203, 1972.[22] M. Cismondi and J. M. Mollerup. Development and application of a three-parameter rk–pr equation ofstate.
Fluid Phase Equilibria , 232(1):74–89, 2005.[23] Seong-Ku Kim, Hwan-Seok Choi, and Yongmo Kim. Thermodynamic modeling based on a generalized cubicequation of state for kerosene/lox rocket combustion.
Combustion and Flame , 159(3):1351–1365, 2012.[24] I.H. Bell, J. Wronski, S. Quoilin, and V. Lemort. Pure and pseudo-pure fluid thermophysical property eval-uation and the open-source thermophysical property library coolprop.
Industrial & Engineering ChemistryResearch , 53(6):2498–2508, 2014.[25] B. E. Poling, J. M. Prausnitz, and J. P. O’Connell. The properties of gases and liquids.
McGraw-Hill, NewYork , 2001.[26] E. Goos, A. Burcat, and B. Ruscic. Report ANL 05/20 TAE 960. Technical report, 2005.http://burcat.technion.ac.il/dir.[27] T. Chung, M. Ajlan, L. Lee, and K. E. Starling. Generalized multiparameter correlation for nonpolar andpolar fluid transport properties.
Ind. Eng. Chem. Res , 27(4):671–679, 1988.
28] L. Qiu and R. D. Reitz. An investigation of thermodynamic states during high-pressure fuel injection usingequilibrium thermodynamics.
International Journal of Multiphase Flow , 72:24–38, 2015.[29] P. H. Van Konynenburg and R. L. Scott. Critical lines and phase equilibria in binary van der waals mix-tures.
Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and EngineeringSciences , 298(1442):495–540, 1980.[30] J. R. Elliott and C. T. Lira.
Introductory chemical engineering thermodynamics . Prentice Hall, 2012.[31] D. L¨udecke and C. L¨udecke.
Thermodynamik: Physikalisch-chemische Grundlagen der thermischen Ver-fahrenstechnik . Springer-Verlag, 2013.[32] M. L. Michelsen and J. M. Mollerup. Thermodynamic models: Fundamentals & computational aspects.
Tie-Line Publications, Holte , 2007.[33] M. L. Michelsen. The isothermal flash problem. part i. stability.
Fluid Phase Equilibria , 9(1):1–19, 1982.[34] G. Gyarmathy.
Grundlagen einer Theorie der Nassdampfturbine . PhD thesis, 1962.[35] W. B. Kay. Liquid-vapor phase equilibrium relations in the ethane-n-heptane system.
Industrial & Engi-neering Chemistry , 30(4):459–465, 1938.[36] V. S. Mehra and G. Thodos. Vapor-liquid equilibrium in the ethane-n-heptane system.
Journal of Chemicaland Engineering Data , 10(3):211–214, 1965.[37] A. J. Kidnay, R. C. Miller, W. R. Parrish, and M. J. Hiza. Liquid-vapour phase equilibria in the n2 ch4system from 130 to 180 k.
Cryogenics , 15(9):531–540, 1975.[38] S. Chang and B. C. Y. Lu. Vapor-liquid equilibriums in the nitrogen-methane-ethane system. In
Chem.Eng. Prog., Symp. Ser.;(United States) , volume 63. Univ. of Ottawa, Ontario, 1967.[39] O. T. Bloomer and J. D. Parent. Liquid-vapor phase behavior of the methane-nitrogen system. In
ChemicalEngineering Progress Symposium Series , volume 49, pages 11–24. AMER INST CHEMICAL ENGINEERS345 E 47TH ST, NEW YORK, NY 10017, 1953.[40] R. Stryjek, P. S. Chappelear, and R. Kobayashi. Low-temperature vapor-liquid equilibriums of nitrogen-methane system.
Journal of Chemical and Engineering Data , 19(4):334–339, 1974.[41] M. R. Cines, J. T. Roach, R. J. Hogan, and C. H. Roland. Nitrogen-methane vapor-liquid equilibria.In
Chemical Engineering Progress Symposium Series , volume 49, pages 1–10. AMER INST CHEMICALENGINEERS 345 E 47TH ST, NEW YORK, NY 10017, 1953.[42] R. A. Heidemann and A. M. Khalil. The calculation of critical points.
AIChE journal , 26(5):769–779, 1980.[43] V. Arp, J. M. Persichetti, and C. Guo-bang. The gruneisen parameter in fluids.
Journal of fluids engineering ,106(2):193–200, 1984.[44] P. Mausbach, A. K¨oster, G. Rutkai, M. Thol, and J. Vrabec. Comparative study of the gr¨uneisen parameterfor 28 pure fluids.
The Journal of chemical physics , 144(24):244505, 2016.[45] D. J. Picard and P. R. Bishnoi. Calculation of the thermodynamic sound velocity in two-phase multicom-ponent fluids.
International journal of multiphase flow , 13(3):295–308, 1987.[46] A. Firoozabadi and H. Pan. Two-phase isentropic compressibility and two-phase somic velocity formulticomponent-hydrocarbon mixtures. In
SPE Reservoir Eval. & Eng. , 2000.[47] D. V. Nichita, P. Khalid, and D. Broseta. Calculation of isentropic compressibility and sound velocity intwo-phase fluids.
Fluid Phase Equilibria , 291(1):95–102, 2010.[48] M. Castier. Thermodynamic speed of sound in multiphase systems.
Fluid Phase Equilibria , 306(2):204–211,2011.[49] A. B. Wood.
A textbook of sound . The Macmillan Company, 1930.[50] G. Lamanna, E. Oldenhof, S. Baab, I. Stotz, and B. Weigand. Disintegration regimes near the critical point.In , page5914, 2012.[51] S. Baab, F. J. F¨orster, G. Lamanna, and B. Weigand. Combined elastic light scattering and two-scaleshadowgraphy of near critical fuel jets. In ,2014.[52] B. H. Min, J. T. Chung, H. Y. Kim, and S. Park. Effects of gas composition on the performance andemissions of compressed natural gas engines.
Journal of Mechanical Science and Technology , 16(2):219–226,2002.[53] K. Kim, H. Kim, B. Kim, and K. Lee. Effect of natural gas composition on the performance of a cng engine.
Oil & Gas Science and Technology-Revue de l’IFP , 64(2):199–206, 2009.[54] A. Demirbas.
Methane Gas Hydrate . Green Energy and Technology. Springer London, 2010.[55] V. V. N. Bhaskar, R. H. Prakash, and B. D. Prasad. Hydrogen fuelled ic engine–an overview.
IJITR ,1(1):046–053, 2013.[56] M. Tuner. Combustion of alternative vehicle fuels in internal combustion engines.
Report within project Apre-study to prepare for interdisciplinary research on future alternative transportation fuels, financed by TheSwedish Energy Agency , 3(2), 2016.[57] M. Tuner. Review and benchmarking of alternative fuels in conventional and advanced engine concepts withemphasis on efficiency, co 2, and regulated emissions. Technical report, SAE Technical Paper, 2016.[58] E. Hu, Z. Huang, B. Liu, J. Zheng, X. Gu, and B. Huang. Experimental investigation on performance and missions of a spark-ignition engine fuelled with natural gas–hydrogen blends combined with egr. Interna-tional journal of hydrogen energy , 34(1):528–539, 2009.[59] M. Sierra-Aznar, D. I. Pineda, B. S. Cage, X. Shi, J. P. Corvello, J. Chen, and R. W. Dibble. Work-ing fluid replacement in gaseous direct-injection internal combustion engines: A fundamental and appliedexperimental investigation. 2017.[60] A. Hamzehloo and P. G. Aleiferis. Large eddy simulation of highly turbulent under-expanded hydrogen andmethane jets for gaseous-fuelled internal combustion engines.
International Journal of Hydrogen Energy ,39(36):21275–21296, 2014.[61] A. Hamzehloo and P. G. Aleiferis. Gas dynamics and flow characteristics of highly turbulent under-expandedhydrogen and methane jets under various nozzle pressure ratios and ambient pressures.
International Journalof Hydrogen Energy , 41(15):6544–6566, 2016.[62] A. Hamzehloo and P. G. Aleiferis. Numerical modelling of transient under-expanded jets under differentambient thermodynamic conditions with adaptive mesh refinement.
International Journal of Heat and FluidFlow , 2016.[63] W. B. Streett and J. C. G. Calado. Liquid-vapour equilibrium for hydrogen+ nitrogen at temperatures from63 to 110 k and pressures to 57 mpa.
The Journal of Chemical Thermodynamics , 10(11):1089–1100, 1978.[64] J. Xiao and K. Liu.
Huaxue-gongcheng , 18(2):8–12, 1990.[65] F.A. Shtekkel and N.M. Tsinn.
Zh.Khim.Prom. , 16(8):24–28, 1939.[66] G. Trappehl and H. Knapp. Vapour-liquid equilibria in the ternary mixtures n2 ch4 c2h6 and n2 c2h6 c3h8.
Cryogenics , 27(12):696–716, 1987.[67] R. Wichterle, I.and Kobayashi. Vapor-liquid equilibrium of methane-ethane system at low temperatures andhigh pressures.
Journal of Chemical and Engineering Data , 17(1):9–12, 1972.[68] S. P. Tang and J. B. Fenn. Experimental determination of the discharge coefficients for critical flow throughan axisymmetric nozzle.
AIAA Journal , 16(1):41–46, 1978.[69] J. Abraham. Entrapment characteristics of transient gas jets.
Numerical Heat Transfer, Part A Applications ,30(4):347–364, 1996.[70] P. G. Hill and P. Ouellette. Transient turbulent gaseous fuel jets for diesel engines.
Journal of fluidsengineering , 121(1):93–101, 1999.[71] J. Gerold, P. Vogl, and M. Pfitzner. New correlation of subsonic, supersonic and cryo gas jets validated byhighly accurate schlieren measurements.
Experiments in fluids , 54(6):1542, 2013.[72] A. Velikorodny and S. Kudriakov. Numerical study of the near-field of highly underexpanded turbulent gasjets.
International journal of hydrogen energy , 37(22):17390–17399, 2012., 37(22):17390–17399, 2012.