- Email: [email protected]
A homogeneous non-equilibrium two-phase critical flow model
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
A homogeneous non-equilibrium two-phase critical flow model J.R. Travis a, D. Piccioni Koch b,*, W. Breitung c a
Engineering & Scientific Software, Inc., 2128 S. Ensenada Circle, Rio Rancho, NM 87124, USA Steinbuch Centre for Computing (SCC), Karlsruher Institut fu¨r Technologie (KIT), Postfach 3640, 76021 Karlsruhe, Germany c simaps GmbH, Rietburgweg 5, 76751 Jockgrim, Germany b
article info
abstract
Article history:
A non-equilibrium two-phase single-component critical (choked) flow model for cryogenic
Received 16 January 2012
fluids is developed from first principle thermodynamics. Modern equations-of-state (EOS)
Received in revised form
based upon the Helmholtz free energy concepts are incorporated into the methodology.
15 June 2012
Extensive validation of the model is provided with the NASA cryogenic data tabulated for
Accepted 18 July 2012
hydrogen, methane, nitrogen, and oxygen critical flow experiments performed with four
Available online 11 August 2012
different nozzles. The model is used to develop a hydrogen critical flow map for stagnation states in the liquid and supercritical regions.
Keywords:
Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights
Cryogenic fluids
reserved.
Critical flows Critical flow model Helmholtz free energy
1.
Introduction
The work described in the present paper was initiated in the frame of the icefuel (Integrated Cable Energy system for FUEL and power) project [1], which aimed to create an Infrastructure System for Cryogenic Hydrogen Storage, Distribution and Decentral Reconversion. The concept is based on the production and liquefaction of hydrogen to buffer large quantities of excess renewable electricity, mainly produced in facilities such as wind and solar power plants. A schematic view of the icefuel system is shown in Fig. 1. The cryogenic hydrogen (<40 K) involved in this project was distributed through a grid of superinsulated pipes. Safety related consequences of cryogenic hydrogen releases from icefuel cables were investigated both theoretically and experimentally at the Institute for Nuclear and Energy Technology (IKET) of the Karlsruhe Institute of Technology (KIT) (Germany). In order to investigate the safety
consequences of small leaks or holes in the cables and estimate consequent discharges, dedicated simulation codes were developed and applied. In this paper, the development of a homogeneous, nonequilibrium, two-phase, critical flow model for cryogenic fluids based upon first principal thermodynamicsis described. The model can be used to accurately calculate discharge mass flow rates from high pressure reservoirs.
2. Critical flow and two-phase critical flow models The critical flow, also referred to as chocking or chocked flow, is defined in terms of Mach number, M ¼ U/c, where U is the fluid flow speed and c is the speed of sound. For critical flows, the Mach number is unitary, increases for supersonic flows (M > 1) and decreases for subcritical ones (M < 1).
* Corresponding author. E-mail addresses: [email protected] (J.R. Travis), [email protected] (D. Piccioni Koch), [email protected] (W. Breitung). 0360-3199/$ e see front matter Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijhydene.2012.07.077
17374
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
Fig. 1 e icefuel schematic [1].
Starting from 1947 with the work of J.G. Burnell [2] many theories have been developed for the calculation of the twophase critical flow and several classification criteria have been used to order such theories. A first classification distinguishes between homogeneous and non-homogeneous models. Homogeneous models assume that the liquid and vapor are mixed together and can be treated as a mixture, while in the non-homogeneous models, the liquid and vapor exist as separated phases. Another relevant criterion of classification is the thermodynamic equilibrium. Some models assume Thermodynamic Equilibrium, i.e. both phases coexist at the same saturation conditions (Thermal Equilibrium), are well-mixed, with equal velocity (Dynamic Equilibrium) and with phase densities constant during expansion (Chemical Equilibrium); the non-equilibrium models assume no thermodynamic equilibrium. Some relevant two-phase critical flow models are reviewed in [3e5] and are not further discussed in this paper. Hereafter, the development of our homogeneous, nonequilibrium two-phase critical flow model is described.
3.
Modern equations of state
Modern equations-of-state [6] are often formulated using the Helmholtz energy as the fundamental property with independent variables of temperature and density, aðT; rÞ ¼ a0 ðT; rÞ þ ar ðT; rÞ;
(1)
where a is the Helmholtz energy, a0(T,r) is the ideal gas contribution to the Helmholtz energy, and ar(T,r) is the residual Helmholtz energy, which corresponds to the influence of intermolecular forces in real gases. Thermodynamics properties can be calculated as derivatives of the Helmholtz energy. For example, the pressure can be expressed as va : (2) p ¼ r2 vp T In practical applications, the functional form is explicit in the dimensionless Helmholtz energy, a, using independent variables of dimensionless density and temperature. The form of this equation is
aðT; rÞ ¼ aðs; dÞ ¼ a0 ðs; dÞ þ ar ðs; dÞ RT
(3)
where s ¼ Tc/T, the inverse reduced temperature, d ¼ r/rc, the reduced density and R is the universal gas constant (8.314510 J/(mol K)). The ideal gas Helmholtz energy is often represented in the computational convenient parameterized form
a0 ðs; dÞ ¼ lnd þ a0 lns þ a1 þ a2 s þ
N X
ak ln½1 expðbk sÞ;
(4)
k¼3
and the residual contribution to the Helmholtz free energy takes the form ar ðs;dÞ¼
l X i¼1
Ni ddi sti þ
m X i¼lþ1
Ni ddi sti expðdpi Þþ
i þ4i ðdDi Þ2 þbi ðsgi Þ2 ;
n X
h Ni ddi sti exp
i¼mþ1
(5)
where the parameters and coefficients in these expressions are given for hydrogen (normal, parahydrogen and orthohydrogen) [7], oxygen [8], nitrogen [9], methane [10], and water [11]. The advantages of this explicit formulation in the Helmholtz free energy become apparent for the calculation of enthalpy, entropy, and sound speed, respectively: r r 0 va va va þ þ1 ; þd hðT; rÞ ¼ RT s vs d vs d vd s
(6)
0 r va va sðT; rÞ ¼ R s a0 ar ; and þ vs d vs d
(7)
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 2 r 2 3 r u va v a u ds 1þd 2 r uRT6 7 var vd s vdvs u 2 v a 6 þd wðT;rÞ¼u 41þ2d 2 0 2 r 7 5: vd s tM vd2 s v a v a s2 þ vs2 d vs2 d (8) Other fluid properties can be found in references [6,12]. The saturation line can be described by the ancillary equation [6] for the saturated vapor-pressure, psat as
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
k n psat Tc X T i ln Ni 1 ¼ Tc pc T i¼1
(9)
where pc is the critical pressure. The derivative of the vaporpressure, which shall be used later in this report, reads as follows: " k 1 # n X dpsat psat ps T i ¼ ki Ni 1 þ (10) ln Tc dT T pc i¼1 Table 1 lists the critical constants and molecular weights for each of the substances that are addressed in this report.
4. Critical discharge analysis from a high pressure reservoir 4.1.
Single-phase choking of a pure substance
17375
which is the definition of the classical critical flow or choked condition, or should the maximum occur at the lowest pressure in the system, the flow is considered subcritical. This approach is referred to as the “Homogeneous Direct Integration” (HDI) method [13,14]. With the use of the above equation of state, it is straight forward to generate a table of paired pressureedensity values from the stagnation state along an isentrope to a pressure less than the choked pressure (for example discharging into the atmosphere at 0.1 MPa) and perform the direct integration of Eq. (15). The maximum value of the integration is found to be the critical discharge mass flux. An equivalent, but more rigorous methodology, which we name the “Homogeneous Direct Evaluation” (HDE) method, that directly exploits the equations-of-state discussed above, is to consider the energy Eq. (12) while neglecting the upstream velocity, i.e., U0 ¼ 0,
The development starts with the differential form of the first law of thermodynamics
h0 ¼ ht þ
dh ¼ Tds þ vdp;
which is arranged to the convenient general mass flux form
(11)
and the control volume form of the conservation of energy 1 1 h0 þ U20 ¼ h þ U2 : 2 2
(12)
Assuming that the process is reversible and adiabatic (an isentropic process with ds ¼ 0) then Eq. (11) can be integrated and combined with Eq. (12) to obtain the famous compressible Bernoulli Equation 1 2 1 2 U U0 ¼ 2 2
Zp0 v$dp
(13)
p
The upsteam reservoir variables, the stagnation state, where the velocity is often assumed zero at location “0”, at any instant are considered in a quasi-steady-state, and as such, the velocity at the choked or critical location “t” can be expressed in terms of the integral along a streamline outside the boundary layer flow to yield 312 2 Zp0 1 7 6 dp5 (14) Ut ¼ 42 r pt
U2t ; 2
1
G ¼ r½2ðh0 hÞ2 :
As above, critical flow requires a local maximum of Eq. (16), or ðvG=vpÞt ¼ 0. When this condition is applied to Eq. (16) one obtains vh 2ðh0 hÞ vp ¼ s : (17) vv v vp s From the first law of thermodynamics Eq. (1) vh ¼ v; vp s
2 6 Gt ¼ rt Ut ¼ rt 42
Zp0
312 1 7 dp5 : r
(15)
(18)
and noting from the definition of sound speed squared w2 ¼
vp v2 ¼ ; vv vr s vp s
(19)
one can combine Eqs. (17)e(19) to get 2ðh0 ht Þ ¼ w2t
At location “t”, the critical discharge mass flux is then
(16)
(20)
It’s not surprising that the maximum velocity of a critical flow condition is the sound speed. The interesting fact is that by directly solving the coupled isentropic and critical flow conditions
pt
The task is to find the maximum of this function, that is, to find the pressure, pt, such that the mass flux is maximum,
Table 1 e Critical properties and molecular weights for hydrogen, methane, nitrogen, and oxygen. Fluid Hydrogen Methane Nitrogen Oxygen
Pc (MPa)
Tc (K)
rc (kg/m3)
M (kg/kmol)
1.2964 4.5992 3.3958 5.043
33.145 190.564 126.192 154.581
31.263 162.66 313.300 436.1
2.01588 16.0428 28.01348 31.9988
s0 ðT0 ; r0 Þ ¼ sðTt ; rt Þ ; 2½h0 ðT0 ; r0 Þ hðTt ; rt Þ ¼ wðTt ; rt Þ2
(21)
for Tt and rt, the exact critical flow state for the given stagnation condition is obtained. The advantages of the HDE method over the HDI method are two: 1. A table of paired densityepressures need not be created, and 2. The local maximum mass flux need not be found using a search technique. The HDE method, by solving system (21), directly determines the exact critical mass flux condition. In the impressive work at NASA by Simoneau and Hendricks [15], four different nozzles were used to investigate choked flow for a number of cryogenic fluids (hydrogen, nitrogen and
17376
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
methane). Gaseous nitrogen was used to calibrate the four nozzles, and a Table was presented with the results [15]. The same Table is presented here for completeness as well as to compare, in the last two columns, the HDE method, the solution of system (21), with the original Table results. The agreement is excellent between the HDE model and the NASA experiments for both critical mass flux and ratio of throat to stagnation pressure. Note that the mass flux ratio is the effective discharge coefficient for each of the individual nozzles. This discharge coefficient shall be applied in the twophase analysis described below.
4.2.
Two-phase choking of a pure substance
The development of the two-phase methodology follows directly from the single-phase approach. There are a number of assumptions that should be noted: 1. The stagnation condition is a pure substance at saturated liquid, subcooled liquid, or supercritical such that an isentropic expansion from the stagnation state to the saturation line, the saturation locus, occurs through the compressed liquid region and not the superheated vapor side of the critical point. This assumption is not very critical to the final results provided the stagnation state is supercritical and not superheated, but since we’re mostly interested in the liquid side, we state this condition, 2. The two-phase flow is homogeneous, 3. The twophase flow is in mechanical equilibrium; that is, the phases have equal velocities. The methodology could be extended into regimes with slip or relative velocity between the phase, but in our direct application (see below), the choked vapor volume fraction is usually less than 10%, so the mechanical coupling between the phases is large; and therefore, relative velocities small, 4. The vapor phase is at saturation, 5. The liquid phase may be in a metastable state (superheated state). 6. The phases share a common pressure (the vapor saturation pressure), 7. The mixture flow is adiabatic and frictionless; and therefore isentropic, and 8. The discharge location has a short L/D ratio. For example, one can imagine a rupture or puncture of a high pressure reservoir wall, or a short nozzle. The HDE method applied to two-phase critical conditions requires two steps: 1. Expand from the stagnation conditions to the liquid saturation line, the saturation locus, (ss ¼ s0, hs, Ts, ps, and rs), and 2. Expand from the liquid saturation locus into the two-phase coexistence region. If the stagnation state is saturated liquid, step 1 is omitted. It remains implicit in this two-step procedure that the maximum liquid superheat allowed for any vapor temperature less than TS is that the liquid temperature Tl ¼ TS. The general pure substance two-phase relationships between various fluid properties and the quality, x, are introduced
s2p ¼ x$sðTv ; rv Þ þ 1 x $sðTl ; rl Þ ¼ x$sv þ 1 x $sl v2p ¼ x$vðTv ; rv Þ þ 1 x $vðTl ; rl Þ ¼ x$vv þ 1 x $vl : h2p ¼ x$hðTv ; rv Þ þ 1 x $hðTl ; rl Þ ¼ x$hv þ 1 x $hl
(22)
If one determines the temperature and density for each phase at the choke plane, then the problem is solved. An analysis is presented below to determine those 4 properties. The two-phase mass flux equation, derived from the conservation of energy, is a direct extension of Eq. (16) with the enthalpy and specific volume states replace with two-phase conditions
G¼
1 1 2 h0 h2p 2 v2p
(23)
The two-phase extension of system (21) governing the twophase critical flow requires four equations in the four unknowns Tl, rl, Tv, and rv. This system is: (1) the conservation of energy, (2) the conservation of entropy, ds ¼ 0, (3) the vapor component is saturated, and (4) both phases share the same pressure. This system is written
2 h0 h2p ¼ w22p s2p ¼ s0 ; psat ¼ pv pl ¼ pv ¼ p
(24)
where the squared two-phase sound speed can be written w22p
vp vr2p
¼
! s
v22p ¼ vv2p vp s
v2 2p : ¼ vvv vvl vx x þ 1x þ vv vl vp s vp s vp s
(25)
Making use of Eq. (19) for each phase results in v22p ; ð1 xÞ vx þ ðvv vl Þ 2 2 vp s ðrv wv Þ ðrl wl Þ
w22p ¼
(26)
x
or in terms of the mass flux, G ¼ rw, G22p ¼
1 : ð1 xÞ vx þ ðv v Þ v l vp s ðrv wv Þ2 ðrl wl Þ2
(27)
x
Should ðvx=vpÞs ¼ 0, then the so called “frozen” mass flux or sound speed is defined, which for a homogeneous two-phase mixture in mechanical and thermal equilibrium is the maximum sound speed of the system. The task now is to find the derivative of the quality with respect to pressure holding the system entropy constant. This is accomplished by solving the system entropy in Eq. (22) for the quality x¼
s0 sðTl ; rl Þ ; sðTv ; rv Þ sðTl ; rl Þ
(28)
and performing the required differentiation yields
vx vp
¼
" # vsl vsv vsl ðsv sl Þ þ ðs0 sl Þ vp sl vp sv vp sl
s
ðsv sl Þ2
:
(29)
Taking sv and sl to be sv ¼ sv(Tv,p) and sl ¼ sl(Tl,p), respectively, and then writing the total differentials gives vsv vsv vp vsv vTv ¼ þ vp vp v vp sv vT vp v p sv : sv T vsl vsl vp vsl vTl ¼ þ vp sl vp Tl vp sl vTl p vp sl
(30)
Upon recognizing Maxwell’s fourth relationship
vs vp
T
vv b ¼ ¼ vT p r
(31)
and the definition of the specific heat at constant pressure
17377
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
30000
1484 Elliptical Nozzle + - 10%
27500 2
one can write Eq. (30) cpv vsv b vTv ¼ vþ vp sv rv Tv vp sv : cpl vsl b vTl ¼ lþ vp sl rl Tl p vp sl
NASA TP 1484 Hydrogen Critical Flow Data and the HDE Model
(32)
Calculated Mass Flux [kg/(m *s)]
vs cp ¼ T ; vT p
(33)
Note that the isobaric heat capacity, cp, and the volume expansivity, b, can be directly calculated from the EOS in section 2, respectively, as 2 r 2 r va v a ds 1þd vd s vdvs cp ðT; rÞ ¼ cv þ R r 2 r and va v a þ d2 1 þ 2d vd s vd2 s 1 vv 1 vp vr ¼ bðT; rÞ ¼ v vT p r vT r vp T where
25000 22500 20000 17500 15000 12500 10000 7500 7500
10000 12500 15000 17500 20000 22500 25000 27500 30000 2
Measured Mass Flux [kg/(m *s)]
Fig. 2 e HDE calculated critical mass fluxes and the NASA hydrogen data [15].
r 2 r vp va v a ¼ RT 1 þ 2d þ d2 and vr T vd s vd2 s NASA TP 1484 Methane Critical Flow Data and the HDE Model
vp vT
r
r 2 r va v a : ¼ Rr 1 þ d ds vd s vdvs
80000
Relating the liquid and vapor temperatures, by defining a “non-equilibrium” parameter, h, can be defined in the following manner: Tl ¼ h$TS þ ð1 hÞ$Tv :
(34)
As stated above, the maximum thermal non-equilibrium liquid superheat allowed is Tl ¼ TS (the saturation locus from step 1 where the fluid is expanded to the liquid saturation line) when h ¼ 1, and least superheat is Tl ¼ Tv when h ¼ 0. The latter case, h ¼ 0, defaults to the well know Homogeneous Equilibrium Model (HEM), that is with the mixture in both thermal and mechanical equilibrium. This analysis provides all degrees of liquid superheat, from none, the HEM, to liquid temperatures at the saturation locus. Principally because of Eq. (34) we’ve restricted, by assumption 1 of the model, to the liquid side of the critical point; otherwise, an assumption concerning metastable vapor, supercooled vapor, would be necessary, and where vapor volume fractions become greater than 0.5, mechanical equilibrium may not be valid as the vapor can accelerate more quickly than the liquid droplets. Differentiating Eq. (34) where the pressure is only a function of temperature on the saturation line, Eq. (10), Eq. (33) for the two distinct phases becomes cp ðTv ; rv Þ vsv bðTv ; rv Þ Tv ¼ þ dps ðTv Þ rv vp sv dTv (35) : cp ðTl ; rl Þ Tl vsl bðTl ; rl Þ p ¼ þ ð1 hÞ dps ðTv Þ rl vp sl dTv The sound speed (26), or the mass flux based sound speed (27), can be computed knowing the four phasic unknowns of temperatures and densities along with Eqs. (28), (29), (34) and
70000
1484 Elliptical Nozzle + - 10%
2
Calculated Mass Flux [kg/(m *s)]
60000 50000 40000 30000 20000 10000 10000
20000
30000
40000
50000
60000
70000
80000
2
Measured Mass Flux [kg/(m *s)]
Fig. 3 e HDE calculated critical mass fluxes and the NASA methane data [15].
Fig. 4 e HDE calculated critical mass fluxes and the NASA nitrogen data [15,16].
17378
HDE model validation
The NASA cryogenic critical flow data [15,16] was used to validate the homogeneous, non-equilibrium, two-phase, critical flow model described by the system (36). The results are shown for hydrogen [15], methane [15], nitrogen [15,16], and oxygen [16], respectively, in Figs. 2e5. The calculated values have been corrected with the discharge coefficient, the mass flux ratio, given in Table 2. In each Figs.2e6, a TeS diagram insert is included to display the analyzed stagnation conditions. The computed results appear to be consistently greater than the measured mass fluxes. This is likely due to fluid dynamic area change and frictional effects, which can become relevat for particular length-to-diameter ratios, pipe geometries and fluid properties. These effects are not addressed in this paper. However, in the overall, the solution of system (36) provides very good agreement with the experimental data.
6. HDE model calculated hydrogen critical mass fluxes The HDE model was used to develop a critical flow map for liquid and supercritical hydrogen. Stagnation conditions are
0.523 0.523 0.523 0.524 873 853 821 837 0.970 0.962 0.963 0.982
Mass flux ratio Gmeas/Gcalc Maximum calculated mass flux, Gcalc (g/cm2 s)
872 852 820 835 846 820 790 820 0.524 0.524 0.524 0.524 0.522 0.537 0.565 0.495
5.
356 351 343 313
Note that the last two expressions are only used for convenience since they are not independent relationships; and therefore are already in terms of the four unknown variables.
272 276.5 284 233
Ts Tv Tl Tv
7 Concial 3.5 Concial 2D Elliptical
h¼
Maximum measured mass flux, Gmeas (g/cm2 s)
s0 sl ðTl ; rl Þ sv ðTv ; rv Þ sl ðTl ; rl Þ
(36)
Ratio of throat to stagnation pressure (calculated)
x¼
:
Ratio of throat to stagnation pressure (measured)
psat ðTv Þ ¼ pv ðTv ; rv Þ pl ðTl ; rl Þ ¼ pv ðTv ; rv Þ
Stagnation pressure (N/cm2)
s2p ðx; Tl ; rl ; Tv ; rv Þ ¼ s0
Stagnation temperature (K)
(35). The system of Eq. (24) is closed, and the details are reviewed here 2 h0 h2p ðx; Tl ; rl ; Tv ; rv Þ ¼ w2p ðh; x; Tl ; rl ; Tv ; rv Þ2
Nozzle
Fig. 5 e HDE calculated critical mass fluxes and the NASA oxygen data [16].
Table 2 e NASA table [15] for gaseous nitrogen with the proposed HDE model results in bold font added in the last two columns.
Proposed HDE model mass flux (g/cm2 s)
Proposed HDE model throat to stagnation pressure
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
references
HDE Model Calculated Critical Mass Fluxes 30000
Critical Point Temperature (33.145 K)
B
2
Mass Flux [kg/(m *s)]
25000
20000
C
15000
10000
D
5000
CP
A
0 0
1
2
3
4
5
6
7
Stagnation Pressure [MPa]
Fig. 6 e HDE calculated critical mass fluxes for hydrogen with stagnation states in the liquid and supercritical regions.
shown in the inserted hydrogen TeS diagram (Fig. 6), where the stagnation temperature, 26 K T0 40 K, and pressure, P0 6 MPa, states are always in the single phase region with entropy, S0 Scritical. After determining the mass flux from the critical flow map in Fig. 6, one should correct it with the relevant discharge coefficient.
7.
17379
Conclusions
A homogeneous non-equilibrium, two-phase, critical flow model, the homogeneous direct evaluation model (HDE), has been developed from first principal thermodynamics and modern equation-of-state formulations. The model has been validated with extensive cryogenic data involving liquid and supercritical hydrogen, methane, nitrogen, and oxygen. A critical discharge flow map for hydrogen is presented that allows the reader a straight forward procedure to determine critical mass fluxes for a range of stagnation conditions.
Acknowledgments This work was performed within the icefuel project which was co-funded by the German Federal Ministry of Education and Research and the Institute for Nuclear and Energy Technology (IKET) of the Karlsruhe Institute of Technology (KIT). The authors are very grateful to the program manager Dr. G. Markowz for continuous support and fruitful discussions throughout this research.
[1] Markowz G, Dylla A, Elliger T. icefuel e an infrastructure system for cryogenic hydrogen storage, distribution and decentral use. In: proceedings of the 18th world hydrogen energy conference (WHEC) 2010. Essen: May 16e21 2010. [2] Burnell JG. Flow of boiling water through nozzles, orifices and pipes. Engineering 1947;164:572e6. [3] D’Auria F, Vigni P. Two-phase critical flow model. CNSI Report No. 49; May 1980. [4] Wallis GB. Critical two-phase flow. Int J Multiphase Flow 1980;6:97e112. Pergamon/Elsevier. [5] Sokolowski L, Kozlowski T. Assessment of two-phase critical flow models performance in RELAP5 and TRACE against Marviken critical flow tests. International Agreement Report, NUREG/IA-0401; February 2012. [6] Jacobsen RT, Penoncello SG, Lemmon EW. Thermodynamic properties of cryogenic fluids. New York: Plenum Press; 1997. [7] Leachman JW, Jacobsen RT, Penoncello SG, Lemmon EW. Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen. J Phys Chem Ref Data 2009; 38(No. 3):721e48. [8] Span R, Lemmon EW, Jacobsen RT, Wagner W, Yokozeki A. A reference equation of state for the thermodynamic properties of nitrogen for temperatures from 63.151 to 1000 K and pressures to 2200 MPa. J Phys Chem Ref Data 2000;29(No. 6):1361e433. [9] Schmidt R, Wagner W. A new form of the equation of state for pure substances and its application to oxygen. Fluid Phase Equilibria 1985;19:175e200. [10] Setzmann U, Wagner W. A new equation of state and tables of thermodynamic properties for methane covering the range from the melting line to 625 K at pressures upto 1000 MPa. J Phys Chem Ref Data 1991; 20(No. 6):1061e151. [11] Wagner W, Pruss A. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J Phys Chem Ref Data 2002;31(No. 2):387e535. [12] Lemmon EW, Jacobsen RT, Penoncello SG, Beyerlein SW. Computer programs for the calculation of thermodynamic properties of cryogens and other fluids. Adv Cryog Eng 1994; 39:1891e7. [13] Darby R, Self FE, Edwards VH. Properly size pressure-relief valves for two-phase flow. Chem Eng 2002:68e74. [14] Darby R. On two-phase frozen and flashing flows in safety relief valves. Recommended calculation method and the proper use of the discharge coefficient. J Loss Prevent Proc 2004;17:255e559. [15] Simoneau RJ, Hendricks RC. Two-phase choked flow of cryogenic fluids in converging-diverging nozzles. NASA Technical Paper 1484; July 1979. [16] Hendricks RC, Simoneau RJ, Barrows RF. Two-phase choked flow of subcooled oxygen and nitrogen. NASA Technical Note TN-8169; February 1976.
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
A homogeneous non-equilibrium two-phase critical flow model J.R. Travis a, D. Piccioni Koch b,*, W. Breitung c a
Engineering & Scientific Software, Inc., 2128 S. Ensenada Circle, Rio Rancho, NM 87124, USA Steinbuch Centre for Computing (SCC), Karlsruher Institut fu¨r Technologie (KIT), Postfach 3640, 76021 Karlsruhe, Germany c simaps GmbH, Rietburgweg 5, 76751 Jockgrim, Germany b
article info
abstract
Article history:
A non-equilibrium two-phase single-component critical (choked) flow model for cryogenic
Received 16 January 2012
fluids is developed from first principle thermodynamics. Modern equations-of-state (EOS)
Received in revised form
based upon the Helmholtz free energy concepts are incorporated into the methodology.
15 June 2012
Extensive validation of the model is provided with the NASA cryogenic data tabulated for
Accepted 18 July 2012
hydrogen, methane, nitrogen, and oxygen critical flow experiments performed with four
Available online 11 August 2012
different nozzles. The model is used to develop a hydrogen critical flow map for stagnation states in the liquid and supercritical regions.
Keywords:
Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights
Cryogenic fluids
reserved.
Critical flows Critical flow model Helmholtz free energy
1.
Introduction
The work described in the present paper was initiated in the frame of the icefuel (Integrated Cable Energy system for FUEL and power) project [1], which aimed to create an Infrastructure System for Cryogenic Hydrogen Storage, Distribution and Decentral Reconversion. The concept is based on the production and liquefaction of hydrogen to buffer large quantities of excess renewable electricity, mainly produced in facilities such as wind and solar power plants. A schematic view of the icefuel system is shown in Fig. 1. The cryogenic hydrogen (<40 K) involved in this project was distributed through a grid of superinsulated pipes. Safety related consequences of cryogenic hydrogen releases from icefuel cables were investigated both theoretically and experimentally at the Institute for Nuclear and Energy Technology (IKET) of the Karlsruhe Institute of Technology (KIT) (Germany). In order to investigate the safety
consequences of small leaks or holes in the cables and estimate consequent discharges, dedicated simulation codes were developed and applied. In this paper, the development of a homogeneous, nonequilibrium, two-phase, critical flow model for cryogenic fluids based upon first principal thermodynamicsis described. The model can be used to accurately calculate discharge mass flow rates from high pressure reservoirs.
2. Critical flow and two-phase critical flow models The critical flow, also referred to as chocking or chocked flow, is defined in terms of Mach number, M ¼ U/c, where U is the fluid flow speed and c is the speed of sound. For critical flows, the Mach number is unitary, increases for supersonic flows (M > 1) and decreases for subcritical ones (M < 1).
* Corresponding author. E-mail addresses: [email protected] (J.R. Travis), [email protected] (D. Piccioni Koch), [email protected] (W. Breitung). 0360-3199/$ e see front matter Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijhydene.2012.07.077
17374
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
Fig. 1 e icefuel schematic [1].
Starting from 1947 with the work of J.G. Burnell [2] many theories have been developed for the calculation of the twophase critical flow and several classification criteria have been used to order such theories. A first classification distinguishes between homogeneous and non-homogeneous models. Homogeneous models assume that the liquid and vapor are mixed together and can be treated as a mixture, while in the non-homogeneous models, the liquid and vapor exist as separated phases. Another relevant criterion of classification is the thermodynamic equilibrium. Some models assume Thermodynamic Equilibrium, i.e. both phases coexist at the same saturation conditions (Thermal Equilibrium), are well-mixed, with equal velocity (Dynamic Equilibrium) and with phase densities constant during expansion (Chemical Equilibrium); the non-equilibrium models assume no thermodynamic equilibrium. Some relevant two-phase critical flow models are reviewed in [3e5] and are not further discussed in this paper. Hereafter, the development of our homogeneous, nonequilibrium two-phase critical flow model is described.
3.
Modern equations of state
Modern equations-of-state [6] are often formulated using the Helmholtz energy as the fundamental property with independent variables of temperature and density, aðT; rÞ ¼ a0 ðT; rÞ þ ar ðT; rÞ;
(1)
where a is the Helmholtz energy, a0(T,r) is the ideal gas contribution to the Helmholtz energy, and ar(T,r) is the residual Helmholtz energy, which corresponds to the influence of intermolecular forces in real gases. Thermodynamics properties can be calculated as derivatives of the Helmholtz energy. For example, the pressure can be expressed as va : (2) p ¼ r2 vp T In practical applications, the functional form is explicit in the dimensionless Helmholtz energy, a, using independent variables of dimensionless density and temperature. The form of this equation is
aðT; rÞ ¼ aðs; dÞ ¼ a0 ðs; dÞ þ ar ðs; dÞ RT
(3)
where s ¼ Tc/T, the inverse reduced temperature, d ¼ r/rc, the reduced density and R is the universal gas constant (8.314510 J/(mol K)). The ideal gas Helmholtz energy is often represented in the computational convenient parameterized form
a0 ðs; dÞ ¼ lnd þ a0 lns þ a1 þ a2 s þ
N X
ak ln½1 expðbk sÞ;
(4)
k¼3
and the residual contribution to the Helmholtz free energy takes the form ar ðs;dÞ¼
l X i¼1
Ni ddi sti þ
m X i¼lþ1
Ni ddi sti expðdpi Þþ
i þ4i ðdDi Þ2 þbi ðsgi Þ2 ;
n X
h Ni ddi sti exp
i¼mþ1
(5)
where the parameters and coefficients in these expressions are given for hydrogen (normal, parahydrogen and orthohydrogen) [7], oxygen [8], nitrogen [9], methane [10], and water [11]. The advantages of this explicit formulation in the Helmholtz free energy become apparent for the calculation of enthalpy, entropy, and sound speed, respectively: r r 0 va va va þ þ1 ; þd hðT; rÞ ¼ RT s vs d vs d vd s
(6)
0 r va va sðT; rÞ ¼ R s a0 ar ; and þ vs d vs d
(7)
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 2 r 2 3 r u va v a u ds 1þd 2 r uRT6 7 var vd s vdvs u 2 v a 6 þd wðT;rÞ¼u 41þ2d 2 0 2 r 7 5: vd s tM vd2 s v a v a s2 þ vs2 d vs2 d (8) Other fluid properties can be found in references [6,12]. The saturation line can be described by the ancillary equation [6] for the saturated vapor-pressure, psat as
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
k n psat Tc X T i ln Ni 1 ¼ Tc pc T i¼1
(9)
where pc is the critical pressure. The derivative of the vaporpressure, which shall be used later in this report, reads as follows: " k 1 # n X dpsat psat ps T i ¼ ki Ni 1 þ (10) ln Tc dT T pc i¼1 Table 1 lists the critical constants and molecular weights for each of the substances that are addressed in this report.
4. Critical discharge analysis from a high pressure reservoir 4.1.
Single-phase choking of a pure substance
17375
which is the definition of the classical critical flow or choked condition, or should the maximum occur at the lowest pressure in the system, the flow is considered subcritical. This approach is referred to as the “Homogeneous Direct Integration” (HDI) method [13,14]. With the use of the above equation of state, it is straight forward to generate a table of paired pressureedensity values from the stagnation state along an isentrope to a pressure less than the choked pressure (for example discharging into the atmosphere at 0.1 MPa) and perform the direct integration of Eq. (15). The maximum value of the integration is found to be the critical discharge mass flux. An equivalent, but more rigorous methodology, which we name the “Homogeneous Direct Evaluation” (HDE) method, that directly exploits the equations-of-state discussed above, is to consider the energy Eq. (12) while neglecting the upstream velocity, i.e., U0 ¼ 0,
The development starts with the differential form of the first law of thermodynamics
h0 ¼ ht þ
dh ¼ Tds þ vdp;
which is arranged to the convenient general mass flux form
(11)
and the control volume form of the conservation of energy 1 1 h0 þ U20 ¼ h þ U2 : 2 2
(12)
Assuming that the process is reversible and adiabatic (an isentropic process with ds ¼ 0) then Eq. (11) can be integrated and combined with Eq. (12) to obtain the famous compressible Bernoulli Equation 1 2 1 2 U U0 ¼ 2 2
Zp0 v$dp
(13)
p
The upsteam reservoir variables, the stagnation state, where the velocity is often assumed zero at location “0”, at any instant are considered in a quasi-steady-state, and as such, the velocity at the choked or critical location “t” can be expressed in terms of the integral along a streamline outside the boundary layer flow to yield 312 2 Zp0 1 7 6 dp5 (14) Ut ¼ 42 r pt
U2t ; 2
1
G ¼ r½2ðh0 hÞ2 :
As above, critical flow requires a local maximum of Eq. (16), or ðvG=vpÞt ¼ 0. When this condition is applied to Eq. (16) one obtains vh 2ðh0 hÞ vp ¼ s : (17) vv v vp s From the first law of thermodynamics Eq. (1) vh ¼ v; vp s
2 6 Gt ¼ rt Ut ¼ rt 42
Zp0
312 1 7 dp5 : r
(15)
(18)
and noting from the definition of sound speed squared w2 ¼
vp v2 ¼ ; vv vr s vp s
(19)
one can combine Eqs. (17)e(19) to get 2ðh0 ht Þ ¼ w2t
At location “t”, the critical discharge mass flux is then
(16)
(20)
It’s not surprising that the maximum velocity of a critical flow condition is the sound speed. The interesting fact is that by directly solving the coupled isentropic and critical flow conditions
pt
The task is to find the maximum of this function, that is, to find the pressure, pt, such that the mass flux is maximum,
Table 1 e Critical properties and molecular weights for hydrogen, methane, nitrogen, and oxygen. Fluid Hydrogen Methane Nitrogen Oxygen
Pc (MPa)
Tc (K)
rc (kg/m3)
M (kg/kmol)
1.2964 4.5992 3.3958 5.043
33.145 190.564 126.192 154.581
31.263 162.66 313.300 436.1
2.01588 16.0428 28.01348 31.9988
s0 ðT0 ; r0 Þ ¼ sðTt ; rt Þ ; 2½h0 ðT0 ; r0 Þ hðTt ; rt Þ ¼ wðTt ; rt Þ2
(21)
for Tt and rt, the exact critical flow state for the given stagnation condition is obtained. The advantages of the HDE method over the HDI method are two: 1. A table of paired densityepressures need not be created, and 2. The local maximum mass flux need not be found using a search technique. The HDE method, by solving system (21), directly determines the exact critical mass flux condition. In the impressive work at NASA by Simoneau and Hendricks [15], four different nozzles were used to investigate choked flow for a number of cryogenic fluids (hydrogen, nitrogen and
17376
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
methane). Gaseous nitrogen was used to calibrate the four nozzles, and a Table was presented with the results [15]. The same Table is presented here for completeness as well as to compare, in the last two columns, the HDE method, the solution of system (21), with the original Table results. The agreement is excellent between the HDE model and the NASA experiments for both critical mass flux and ratio of throat to stagnation pressure. Note that the mass flux ratio is the effective discharge coefficient for each of the individual nozzles. This discharge coefficient shall be applied in the twophase analysis described below.
4.2.
Two-phase choking of a pure substance
The development of the two-phase methodology follows directly from the single-phase approach. There are a number of assumptions that should be noted: 1. The stagnation condition is a pure substance at saturated liquid, subcooled liquid, or supercritical such that an isentropic expansion from the stagnation state to the saturation line, the saturation locus, occurs through the compressed liquid region and not the superheated vapor side of the critical point. This assumption is not very critical to the final results provided the stagnation state is supercritical and not superheated, but since we’re mostly interested in the liquid side, we state this condition, 2. The two-phase flow is homogeneous, 3. The twophase flow is in mechanical equilibrium; that is, the phases have equal velocities. The methodology could be extended into regimes with slip or relative velocity between the phase, but in our direct application (see below), the choked vapor volume fraction is usually less than 10%, so the mechanical coupling between the phases is large; and therefore, relative velocities small, 4. The vapor phase is at saturation, 5. The liquid phase may be in a metastable state (superheated state). 6. The phases share a common pressure (the vapor saturation pressure), 7. The mixture flow is adiabatic and frictionless; and therefore isentropic, and 8. The discharge location has a short L/D ratio. For example, one can imagine a rupture or puncture of a high pressure reservoir wall, or a short nozzle. The HDE method applied to two-phase critical conditions requires two steps: 1. Expand from the stagnation conditions to the liquid saturation line, the saturation locus, (ss ¼ s0, hs, Ts, ps, and rs), and 2. Expand from the liquid saturation locus into the two-phase coexistence region. If the stagnation state is saturated liquid, step 1 is omitted. It remains implicit in this two-step procedure that the maximum liquid superheat allowed for any vapor temperature less than TS is that the liquid temperature Tl ¼ TS. The general pure substance two-phase relationships between various fluid properties and the quality, x, are introduced
s2p ¼ x$sðTv ; rv Þ þ 1 x $sðTl ; rl Þ ¼ x$sv þ 1 x $sl v2p ¼ x$vðTv ; rv Þ þ 1 x $vðTl ; rl Þ ¼ x$vv þ 1 x $vl : h2p ¼ x$hðTv ; rv Þ þ 1 x $hðTl ; rl Þ ¼ x$hv þ 1 x $hl
(22)
If one determines the temperature and density for each phase at the choke plane, then the problem is solved. An analysis is presented below to determine those 4 properties. The two-phase mass flux equation, derived from the conservation of energy, is a direct extension of Eq. (16) with the enthalpy and specific volume states replace with two-phase conditions
G¼
1 1 2 h0 h2p 2 v2p
(23)
The two-phase extension of system (21) governing the twophase critical flow requires four equations in the four unknowns Tl, rl, Tv, and rv. This system is: (1) the conservation of energy, (2) the conservation of entropy, ds ¼ 0, (3) the vapor component is saturated, and (4) both phases share the same pressure. This system is written
2 h0 h2p ¼ w22p s2p ¼ s0 ; psat ¼ pv pl ¼ pv ¼ p
(24)
where the squared two-phase sound speed can be written w22p
vp vr2p
¼
! s
v22p ¼ vv2p vp s
v2 2p : ¼ vvv vvl vx x þ 1x þ vv vl vp s vp s vp s
(25)
Making use of Eq. (19) for each phase results in v22p ; ð1 xÞ vx þ ðvv vl Þ 2 2 vp s ðrv wv Þ ðrl wl Þ
w22p ¼
(26)
x
or in terms of the mass flux, G ¼ rw, G22p ¼
1 : ð1 xÞ vx þ ðv v Þ v l vp s ðrv wv Þ2 ðrl wl Þ2
(27)
x
Should ðvx=vpÞs ¼ 0, then the so called “frozen” mass flux or sound speed is defined, which for a homogeneous two-phase mixture in mechanical and thermal equilibrium is the maximum sound speed of the system. The task now is to find the derivative of the quality with respect to pressure holding the system entropy constant. This is accomplished by solving the system entropy in Eq. (22) for the quality x¼
s0 sðTl ; rl Þ ; sðTv ; rv Þ sðTl ; rl Þ
(28)
and performing the required differentiation yields
vx vp
¼
" # vsl vsv vsl ðsv sl Þ þ ðs0 sl Þ vp sl vp sv vp sl
s
ðsv sl Þ2
:
(29)
Taking sv and sl to be sv ¼ sv(Tv,p) and sl ¼ sl(Tl,p), respectively, and then writing the total differentials gives vsv vsv vp vsv vTv ¼ þ vp vp v vp sv vT vp v p sv : sv T vsl vsl vp vsl vTl ¼ þ vp sl vp Tl vp sl vTl p vp sl
(30)
Upon recognizing Maxwell’s fourth relationship
vs vp
T
vv b ¼ ¼ vT p r
(31)
and the definition of the specific heat at constant pressure
17377
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
30000
1484 Elliptical Nozzle + - 10%
27500 2
one can write Eq. (30) cpv vsv b vTv ¼ vþ vp sv rv Tv vp sv : cpl vsl b vTl ¼ lþ vp sl rl Tl p vp sl
NASA TP 1484 Hydrogen Critical Flow Data and the HDE Model
(32)
Calculated Mass Flux [kg/(m *s)]
vs cp ¼ T ; vT p
(33)
Note that the isobaric heat capacity, cp, and the volume expansivity, b, can be directly calculated from the EOS in section 2, respectively, as 2 r 2 r va v a ds 1þd vd s vdvs cp ðT; rÞ ¼ cv þ R r 2 r and va v a þ d2 1 þ 2d vd s vd2 s 1 vv 1 vp vr ¼ bðT; rÞ ¼ v vT p r vT r vp T where
25000 22500 20000 17500 15000 12500 10000 7500 7500
10000 12500 15000 17500 20000 22500 25000 27500 30000 2
Measured Mass Flux [kg/(m *s)]
Fig. 2 e HDE calculated critical mass fluxes and the NASA hydrogen data [15].
r 2 r vp va v a ¼ RT 1 þ 2d þ d2 and vr T vd s vd2 s NASA TP 1484 Methane Critical Flow Data and the HDE Model
vp vT
r
r 2 r va v a : ¼ Rr 1 þ d ds vd s vdvs
80000
Relating the liquid and vapor temperatures, by defining a “non-equilibrium” parameter, h, can be defined in the following manner: Tl ¼ h$TS þ ð1 hÞ$Tv :
(34)
As stated above, the maximum thermal non-equilibrium liquid superheat allowed is Tl ¼ TS (the saturation locus from step 1 where the fluid is expanded to the liquid saturation line) when h ¼ 1, and least superheat is Tl ¼ Tv when h ¼ 0. The latter case, h ¼ 0, defaults to the well know Homogeneous Equilibrium Model (HEM), that is with the mixture in both thermal and mechanical equilibrium. This analysis provides all degrees of liquid superheat, from none, the HEM, to liquid temperatures at the saturation locus. Principally because of Eq. (34) we’ve restricted, by assumption 1 of the model, to the liquid side of the critical point; otherwise, an assumption concerning metastable vapor, supercooled vapor, would be necessary, and where vapor volume fractions become greater than 0.5, mechanical equilibrium may not be valid as the vapor can accelerate more quickly than the liquid droplets. Differentiating Eq. (34) where the pressure is only a function of temperature on the saturation line, Eq. (10), Eq. (33) for the two distinct phases becomes cp ðTv ; rv Þ vsv bðTv ; rv Þ Tv ¼ þ dps ðTv Þ rv vp sv dTv (35) : cp ðTl ; rl Þ Tl vsl bðTl ; rl Þ p ¼ þ ð1 hÞ dps ðTv Þ rl vp sl dTv The sound speed (26), or the mass flux based sound speed (27), can be computed knowing the four phasic unknowns of temperatures and densities along with Eqs. (28), (29), (34) and
70000
1484 Elliptical Nozzle + - 10%
2
Calculated Mass Flux [kg/(m *s)]
60000 50000 40000 30000 20000 10000 10000
20000
30000
40000
50000
60000
70000
80000
2
Measured Mass Flux [kg/(m *s)]
Fig. 3 e HDE calculated critical mass fluxes and the NASA methane data [15].
Fig. 4 e HDE calculated critical mass fluxes and the NASA nitrogen data [15,16].
17378
HDE model validation
The NASA cryogenic critical flow data [15,16] was used to validate the homogeneous, non-equilibrium, two-phase, critical flow model described by the system (36). The results are shown for hydrogen [15], methane [15], nitrogen [15,16], and oxygen [16], respectively, in Figs. 2e5. The calculated values have been corrected with the discharge coefficient, the mass flux ratio, given in Table 2. In each Figs.2e6, a TeS diagram insert is included to display the analyzed stagnation conditions. The computed results appear to be consistently greater than the measured mass fluxes. This is likely due to fluid dynamic area change and frictional effects, which can become relevat for particular length-to-diameter ratios, pipe geometries and fluid properties. These effects are not addressed in this paper. However, in the overall, the solution of system (36) provides very good agreement with the experimental data.
6. HDE model calculated hydrogen critical mass fluxes The HDE model was used to develop a critical flow map for liquid and supercritical hydrogen. Stagnation conditions are
0.523 0.523 0.523 0.524 873 853 821 837 0.970 0.962 0.963 0.982
Mass flux ratio Gmeas/Gcalc Maximum calculated mass flux, Gcalc (g/cm2 s)
872 852 820 835 846 820 790 820 0.524 0.524 0.524 0.524 0.522 0.537 0.565 0.495
5.
356 351 343 313
Note that the last two expressions are only used for convenience since they are not independent relationships; and therefore are already in terms of the four unknown variables.
272 276.5 284 233
Ts Tv Tl Tv
7 Concial 3.5 Concial 2D Elliptical
h¼
Maximum measured mass flux, Gmeas (g/cm2 s)
s0 sl ðTl ; rl Þ sv ðTv ; rv Þ sl ðTl ; rl Þ
(36)
Ratio of throat to stagnation pressure (calculated)
x¼
:
Ratio of throat to stagnation pressure (measured)
psat ðTv Þ ¼ pv ðTv ; rv Þ pl ðTl ; rl Þ ¼ pv ðTv ; rv Þ
Stagnation pressure (N/cm2)
s2p ðx; Tl ; rl ; Tv ; rv Þ ¼ s0
Stagnation temperature (K)
(35). The system of Eq. (24) is closed, and the details are reviewed here 2 h0 h2p ðx; Tl ; rl ; Tv ; rv Þ ¼ w2p ðh; x; Tl ; rl ; Tv ; rv Þ2
Nozzle
Fig. 5 e HDE calculated critical mass fluxes and the NASA oxygen data [16].
Table 2 e NASA table [15] for gaseous nitrogen with the proposed HDE model results in bold font added in the last two columns.
Proposed HDE model mass flux (g/cm2 s)
Proposed HDE model throat to stagnation pressure
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 1 7 3 7 3 e1 7 3 7 9
references
HDE Model Calculated Critical Mass Fluxes 30000
Critical Point Temperature (33.145 K)
B
2
Mass Flux [kg/(m *s)]
25000
20000
C
15000
10000
D
5000
CP
A
0 0
1
2
3
4
5
6
7
Stagnation Pressure [MPa]
Fig. 6 e HDE calculated critical mass fluxes for hydrogen with stagnation states in the liquid and supercritical regions.
shown in the inserted hydrogen TeS diagram (Fig. 6), where the stagnation temperature, 26 K T0 40 K, and pressure, P0 6 MPa, states are always in the single phase region with entropy, S0 Scritical. After determining the mass flux from the critical flow map in Fig. 6, one should correct it with the relevant discharge coefficient.
7.
17379
Conclusions
A homogeneous non-equilibrium, two-phase, critical flow model, the homogeneous direct evaluation model (HDE), has been developed from first principal thermodynamics and modern equation-of-state formulations. The model has been validated with extensive cryogenic data involving liquid and supercritical hydrogen, methane, nitrogen, and oxygen. A critical discharge flow map for hydrogen is presented that allows the reader a straight forward procedure to determine critical mass fluxes for a range of stagnation conditions.
Acknowledgments This work was performed within the icefuel project which was co-funded by the German Federal Ministry of Education and Research and the Institute for Nuclear and Energy Technology (IKET) of the Karlsruhe Institute of Technology (KIT). The authors are very grateful to the program manager Dr. G. Markowz for continuous support and fruitful discussions throughout this research.
[1] Markowz G, Dylla A, Elliger T. icefuel e an infrastructure system for cryogenic hydrogen storage, distribution and decentral use. In: proceedings of the 18th world hydrogen energy conference (WHEC) 2010. Essen: May 16e21 2010. [2] Burnell JG. Flow of boiling water through nozzles, orifices and pipes. Engineering 1947;164:572e6. [3] D’Auria F, Vigni P. Two-phase critical flow model. CNSI Report No. 49; May 1980. [4] Wallis GB. Critical two-phase flow. Int J Multiphase Flow 1980;6:97e112. Pergamon/Elsevier. [5] Sokolowski L, Kozlowski T. Assessment of two-phase critical flow models performance in RELAP5 and TRACE against Marviken critical flow tests. International Agreement Report, NUREG/IA-0401; February 2012. [6] Jacobsen RT, Penoncello SG, Lemmon EW. Thermodynamic properties of cryogenic fluids. New York: Plenum Press; 1997. [7] Leachman JW, Jacobsen RT, Penoncello SG, Lemmon EW. Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen. J Phys Chem Ref Data 2009; 38(No. 3):721e48. [8] Span R, Lemmon EW, Jacobsen RT, Wagner W, Yokozeki A. A reference equation of state for the thermodynamic properties of nitrogen for temperatures from 63.151 to 1000 K and pressures to 2200 MPa. J Phys Chem Ref Data 2000;29(No. 6):1361e433. [9] Schmidt R, Wagner W. A new form of the equation of state for pure substances and its application to oxygen. Fluid Phase Equilibria 1985;19:175e200. [10] Setzmann U, Wagner W. A new equation of state and tables of thermodynamic properties for methane covering the range from the melting line to 625 K at pressures upto 1000 MPa. J Phys Chem Ref Data 1991; 20(No. 6):1061e151. [11] Wagner W, Pruss A. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J Phys Chem Ref Data 2002;31(No. 2):387e535. [12] Lemmon EW, Jacobsen RT, Penoncello SG, Beyerlein SW. Computer programs for the calculation of thermodynamic properties of cryogens and other fluids. Adv Cryog Eng 1994; 39:1891e7. [13] Darby R, Self FE, Edwards VH. Properly size pressure-relief valves for two-phase flow. Chem Eng 2002:68e74. [14] Darby R. On two-phase frozen and flashing flows in safety relief valves. Recommended calculation method and the proper use of the discharge coefficient. J Loss Prevent Proc 2004;17:255e559. [15] Simoneau RJ, Hendricks RC. Two-phase choked flow of cryogenic fluids in converging-diverging nozzles. NASA Technical Paper 1484; July 1979. [16] Hendricks RC, Simoneau RJ, Barrows RF. Two-phase choked flow of subcooled oxygen and nitrogen. NASA Technical Note TN-8169; February 1976.