Flow Characteristics of High-Pressure Hydrogen Gas
in the Critical Nozzle
Shigeru MATSUO*1, Soichiro KOYAMA*2, Junji NAGAO*3
and Toshiaki SETOGUCHI*4*1

Saga University, Dept. of Mechanical Engineering,
1 Honjo-machi, Saga-shi, Saga 840-8502, Japan

Hydrogen gas has been extensively receiving much interest as one of alternative energy sources infuture. Of the practical usages of hydrogen gas, a fuel cell car makes use of high-pressure hydrogen gas
nearly at atmospheric temperature. In such an application, metering technique of hydrogen gas is ofpractical importance for mileage and power output of the car. However, the precise measurement of flowrate of hydrogen gas is extremely difficult due to compressibility and real gas effects. For the purpose ofpractical use of high-pressure hydrogen gas, systematic research is required to clarify the critical nozzle
flow of high-pressure hydrogen gas. In the present study, a computational method has been applied topredict the critical nozzle flow of high-pressure hydrogen gas.
Key Words : Critical Nozzle, Compressible Flow, Hydrogen Gas, Real Gas Effects, Internal Flow1.

For instance, a hydrogen fuel cell, which isIntroduction

being received as the driving power system of

Critical nozzle is one of the flow metering

vehicles, yields high-pressure hydrogen gas

devices which are being extensively used innearly at atmospheric temperature. In such anindustrial area dealing with gases. It makes use

application, precise measurement of flow rateof the concept of flow choke that occurs at thenozzle

throat(1).

Underthechoked

flow

conditions, pressure variations in the flow

is of practical importance for mileage andpower output of the vehicle.
Only a few researches have been to datefield downstream of the nozzle have negligible

madeinfluence on the mass flow rate, and thehigh-pressure hydrogen gas through criticalcoefficient of discharge is easily obtained only

nozzle. Recently, Nakao (4 ) has conducted theby the flow properties measured upstream offlow rate measurement of hydrogen gas using athe nozzle.

critical(2 ), ( 3)onthemass

nozzle,and

flow

has

rateoffoundthatthethehave reported thedischarge coefficient of hydrogen gas exceeds

discharge coefficients of a variety of gases forunity in a specific Reynolds number regime.

a quite wide range of Reynolds number, using

No detailed explanation has been made for thisa computational fluid dynamics method. Theabnormal

flow characteristics through the critical nozzle

high-pressure hydrogen gas.

Recently, Kim et al.

discharge

coefficientof

both

The present study aims at investigating theconsiderably low and high Reynolds number

detailed flow of high-pressure hydrogen gas

regimes.

through a critical nozzle, with the help of a

have

beenwelldocumented

atOf many kinds of working gases employedcomputational fluid dynamics method. The

in industrial field, hydrogen gas is one of thecomputational

most promising gases as alternative energy.

some experimental data available.
2.

*1

佐賀大学理工学部 (〒840-8502 佐賀県佐賀市本庄町 1)
佐賀大学大学院 工学系研究科*3
佐賀大学大学院 工学系研究科．
*3
佐賀大学 海洋エネルギー研究センター
E-mail: matsuo@me.saga-u.ac.jp

results

werevalidated

with

Computational analysis

*2

The

high-pressurehydrogengasflow

through the critical nozzle is simulated using a
Computational Fluid Dynamics method. The

governingequations

aregivenbythec = v {C p /( R − C p )} /( ∂v / ∂p ) T

conservation forms of mass, momentum andenergy. The axisymmetric,time-dependent(3)

mass averaged,

Navier-Stokes

equations,

The dynamic viscosity is given as ( 7 ) ,which use a k- ε turbulent model, are employedin the present computations ( 1 ) , (2 ) .

µ (T ) = 6.3 ×10 −7

Hydrogen gas has little attraction between

M w0.5 pc0.6666  Tr1.5 
Tc0.1666  Tr + 0.8 (4)molecules so that it is inert in character until

where T r is the reduced temperature (=T/T c )the molecules are disrupted. From the point ofand M w is the molecular weight of real gas.

view

statistical

Using this viscosity, the thermal conductivitygas,

is obtained from Eucken formula ( 8 ) .

ofkinetic

thermodynamics

theoryoforhydrogenthecompressibility factor (Z) is frequently given

In the present study, the thermodynamic

as a polynomial function of specific volume v,

properties obtained by the above equations arepv
B (T ) C (T ) D (T )Z== 1++ 2 + 3 + ⋅⋅⋅⋅RTvvv

incorporated(1)

where B(T), C(T) and D(T) are the virialcoefficients that are dependent only ontemperature.
Similar equation of state of real gas was also
reported by Benedict, Webb and Rubin (BWRequation of state) ( 5) , and later on, Lee &Redlich and Kwong have improved the van
Waals

equationcompressibility

factortomoregoverning

equation

system and numerically solved to assess thereal gas effects of high-pressure hydrogen gas.
Inthepresent

study,

Redlich-Kwong’s

equation of state is selected to simulate thehigh-pressure hydrogen gas through a criticalnozzle, since it bears a better agreement withthe virial equation of state.

obtaintheaccurately,astemperature

andpressure

onthecompressibility factor of hydrogen gas, Figure
1 shows the calculation results using thefollows,

2.0RTa (T )~−(+ b0 )vv
v−b(2)

where,nT a(T ) = a0  c  ,T ~
b = b0 − c0 ,R 2Tc2a0 = 0.42747,
pc

c0 = RTc

RT
b0 = 0.08664 c ,pc

a0 + b0 − vc
 pc +vc (vc + b0 ) In case of hydrogen gas, the n value in theCompressibility factor [Z]p=theIn order to scrutinize the dependence ofKesler ( 6 ) modified BWR equation of state.
Derinto(Eq.103))
Redlich-Kwong's equation of state (Eq.

1.560
T = 50 K

7080

125102050

100 200

pa , atm

Fig. 1 Dependence of pressure on thecompressibility factor ofhydrogen gas
R = 2.0 D
Pressure
inlet

3oAdiabatic
no-slip wall
4ovolume, respectively, based on the critical

Pressure
outletFlow
D/2
direction

point.
The speed of sound of real gas is determined

160 200

minimum point

0.5

function of a(T) is given by n=0.31. p c , T c andv c are the pressure, temperature and specific100

1.0-1.0 D

5.0 D

D = 0.6 mm

from the following thermodynamic relation,

Fig. 2 Schematic diagram of criticalnozzle flow field

500

1.0

Redlich-Kwong’s equation of state. As thepressure is quite low, the compressibility
factor approaches 1.0, being the same to thatas

thecompressibility

pressure

factorincreases,

becomes

athestrong

ps / poof ideal gas, regardless of the temperature.
However,

Ideal gas
Real gas0.80.6
po = 1.01 MPa

0.4function of temperature: At low temperature,
the compressibility factor decreases and thenincreases

withpressure,

afterreaching

35.46 MPa
1.01 MPa35.46 MPa

0.2aminimum value. It seems that the minimum
value is dependent on the gas temperature.
However, at high temperature, such a trend is

0.0
-1.0

0.0

1.0

2.0

3.0

4.0

5.0

x/D

no longer found, and the compressibility factor

Fig. 3 Predicted static pressure distributions
along the nozzle axis

is a simply increasing function of the pressure.

properties at the inlet of nozzle. It is foundThis is because

the molecularweightofthat the coefficient of discharge is a stronghydrogenisgas

intermolecularverylowandthefunction

attraction

force

isless,predicted coefficients of discharge were incompared with the other gases.

ofReynolds

number,andthegood agreement with the experimental results.

The critical nozzle employed in the present

It is believed that the present computation

study is of a typical conical type ( 1 ), ( 2) and its

method predicts well the gas flow through thediameter D at the throat is D=0.6 mm. A

critical nozzle.

convergent part with a radius of curvature of

Figure 3 describes the static pressure (p s )2.0D is given upstream of the nozzle throat.

distributions of hydrogen gas flow along theThe straight divergent part has a half angle ofnozzle axis. The pressure ratio (p a /p o ) is fixed

4 degrees and its axial length is 5.0D. The

at 0.5. The static pressure decreases with thecomputational domain and boundary conditions

distance and then suddenly increases when theused in the present study are illustrated inflow meets the shock wave. The shock wave

Fig.2. Inlet total pressures and back pressure

moves downstream as the inlet total pressureof nozzle are p o and p a , respectively. The

increases. It is noted that for the same p o , thesymmetric conditions are assumed at the axis

real gas effect causes the shock wave to beof

located further downstream, compared with thecriticalnozzle.

The

adiabatic,no-slip

conditions are applied to the solid walls.

ideal gas. This is qualitatively the same as

In the present study, the nozzle pressure

what the inlet total pressure influences theratio is defined as p a / p o and its value

shock locations. In this figure, it should bemaintains constant at 0.5, but the inlet total

noted that the real gas effect reduces the staticpressure is varied in the range from 2.0 bar topressure at the nozzle throat, compared with

350.0 bar at a fixed total temperature, T 0 =288

the ideal gas. This implies that the mass flow

K. A structured grid system with about 35,000

through a critical nozzle can be different due

grid points was employed

to the real gas effect.

in the present

computations.
3.

Meanwhile, the coefficient of discharge C d isResults and discussion

The present computations were validatedwith the experimental results(9)usually defined as,

C d = m& m& theo(5)

using hydrogengas, where the Reynolds number is based on
the diameter of nozzle throat and the total

where m& is the theoretical mass flow rate
through a critical nozzle and m& theo means themass flow rate calculated by one-dimensional
gasdynamics theory.
Figure 4 shows the comparison of thepredictedandexperimental ( 4 )discharge
coefficients. At present, the only experimental
data for high-pressure hydrogen gas areavailable in Ref.(4). In 1965, Johnson ( 1 0 )calculated the ratio of specific heats using theVirial equation of state and tabulated thethermodynamic properties of the hydrogen gas
flow through a critical nozzle. It is interesting
to note that the experimented data show a

Discharge coefficient [Cd ]1.04
Present CFD (Ideal gas)
Present CFD (Real gas)Ref. ( (4)14 )Ref. ( (10)
10 )

1.02

higher coefficient of discharge than unity, inthe range of Reynolds number below 6.0 × 10 5 .
At present, this unreasonable trend is not wellunderstood, whether it comes from the real gas
effect or from the experimental error. The
present computation taking account for the realgas effect is qualitatively similar to Ref. (10)and (4), but still fails in predicting thecoefficient of discharge higher than unity.
Unfortunately,aclearandpersuasive
reasoning for this is, at present, not known.
More study is needed to elucidate thisambiguous problem.
In order to, in more detail, investigate thereal gas effect of hydrogen gas flow, thepredicted compressibility factor and mass flux

1.00

along the nozzle axis are presented in Fig.5

0.98

iso-density contours. Upstream of the shock(p 0 =20.27MPa), together with the computedwave,
0.96

thecompressibility

factorgreatly

decreases with the distance, and at the shock
wave location, it sharply increases due to the

0.94
0.0

4.0

8.0

12.0

16.0

20.0

Re (10× -5)

compression effect of shock wave. There isalso

Fig. 4 Variation of the predicted andexperimental dischargecoefficients with Reynolds
numbersome

difference

inthemassflux

distributions upstream of the shock wave. Thus,
it is concluded that at high pressure conditions,the compressibility factor of the real gas canbe one of the reasons for the coefficients ofDensity level

discharge that were discussed in Fig.4.

( kg / m3 )Figure 6 shows the axial distributions of theIdeal gas

4.90

15.90

computed values of the ratio of specific heats

4.90
9.90

6.90

11.90
9.90

5.90

13.90

Real gas

3.90

3.90

8.90

6.90
8.90

of hydrogen gas. The ratio of specific heats for
8.90

7.90

ideal gas is constant at 1.41. The specific heatthroat1.0000.960

mass flux
compressibility factor

Specific heat ratio [γ]1.040

Normalized mass flux [ρu / ρutheo ]Compressibility factor [Z]

1.080

1.43

1.20

1.120

1.00

0.80

Ideal gas

0.60

-1.0

Ideal gas
( Z=1.0 )
0.0

1.0

1.42

2.02
1.01
0.5

po = 0.2 MPa
5.06
10.13

1.41

20.27

Ideal gas (γ = 1.41)

35.46

1.40
1.39

Real gas

2.0

3.0

4.0

x/D

5.0

1.38
-1.0

0.0

1.0

2.0

3.0

4.0

x/D
Fig. 5 Distributions of predicted
compressibility factors and mass
fluxes along the nozzle axis (po =20.27 MPa)

Fig. 6 Influence of inlet total pressure on
the specific heat ratio of hydrogen
gas

5.0

ratio for real gas increases as the flow isaccelerated through the nozzle throat, and it
sharply decreases at the location of the shock
wave.

Sucha

tendency

appears

similar,regardless of P 0 , but the sudden drop in theratioofspecific

heats

seems

tomoresignificant when P 0 increases. It is interesting
to note that at the nozzle throat, the ratio ofspecific heats remains nearly constant at 1.41,
whenisp0lessthen1.01MPa,

but

it

significantly decreases as p 0 increases. At thenozzle throat, this variation in the ratio ofspecific heats with p 0 has an appreciable
importance on the mass flow rate of hydrogengas through a critical nozzle. Therefore, thereal gas effect should be included in evaluating
theperformanceandaccuracy

ofcriticalReferences(1)

Kim, H. D., Kim, J. H., Park, K. A., Setoguchi,T., and Matsuo, S., Computational Study of the Gas
Flow through a Critical Nozzle, IMechE J
Mechanical Eng Science , 2003, 217(10), pp.
1179-1189.
(2) Kim, H. D., Kim, J. H., Park, K. A., Setoguchi,T., and Matsuo, S., Study of the Effects of UnsteadyDownstream Conditions on the Gas Flow through a
Critical Nozzle, IMechE J Mechanical Eng
Science , 2004, 218, pp. 1163-1173.
(3) Kim, J. H., Kim, H. D., Park, K. A., Setoguchi,T., and Matsuo, S., A Fundamental Study of a
Variable Critical Nozzle Flow, Experiments inFluids , 2006, 40(1), pp. 127-134.
(4) Nakao, S., Development of critical nozzlefor High pressure Hyd rogen Gas Flow
Measurements, Proceedings of JSME , Fluid
Dynamics Section, 2005.
(5) Giorgio, S. S., An Effective Modificationof the Benedict-Webb-Rubin Equation ofState, Fluid Phase Equilibria , 1999, 164, pp.

nozzle as flow metering device.
4.
SeveralConcluding remarks
importantandmeaningfulconclusions obtained from the present study
are summarized; Redlich-Kwong’s equation ofstatepredicts

therealgaseffectsofhigh-pressure hydrogen gas comparatively well.
However, unlike the coefficient of discharge ofideal gases which have been obtained to date,the coefficient of discharge of real gas throughcritical nozzle decreases with an increase in
Reynolds number, as Reynolds number exceeds
a certain value. It is believed that this mainly
results from the thermodynamic properties ofreal gas, such as the compressibility factor andthe ratio of specific heats, which appears more
remarkable as the pressure of hydrogen gas
increases.
Acknowledgement
This study was partially supported by Harada
Memorial Foundation.

157-172.
Jyisy, Y., Peter R, G., Anthon y, R. H. G.,
Comparison of Methods for Calculating
ThermodynamicPropertiesofBinary
Mixtures in the Sub and Super Critical State:
Lee-Kesler and Cubic Equations of State forBinary Mixtures Containing either CO2 or
H2S, J. Chemical Thermodynamics , 2003, 35,
pp. 1521-1539.
(7)
Cheremisinoff, N. P., Fluid Flow Pocket
Handbook , Gulf Publishing Co., Houston,TX., 1984.
(8)
Eckert, ERG., Drake, RM., Analysis ofHeat and Mass Transfer , McGraw-Hill Co.,
1972.
(9) Nakao, S., Irayama, T. and Takamoto, M.,(6)

Relations
betweenthedischargecoefficients of the sonic venture nozzle
and a kind of gas, JSME J MechanicalEngineering , Ser. B, 2000, 66(642), pp.
438-444.
(10) Johnson

R. C., Real-Gas Effects in
Critical-Flow-throughNozzlesandTabulated
ThermodynamicProperties,
NASA TN D-2565, 1965.