Semianalytical solutions are developed for turbulent hydrogen-air plume. We derived analytical expressions for plume centerline variables (radius, velocity, and density deficit) in terms of a single universal function, called plume function. By combining the obtained analytical expressions of centerline variables with empirical Gaussian expressions of the mean variables, we obtain semianalytical expressions for mean quantities of hydrogen-air plume (velocity, density deficit, and mass fraction).

1. Introduction

One of the important safety issues of hydrogen energy is the hydrogen leakage into ambient air and the associated risk of fire or explosion. In fact, industry has already produced several prototype products using hydrogen as a fuel. Unfortunately, these products are not yet available for commercial use because of safety concerns related to hydrogen leakage. So studying hydrogen-air behavior is very important in order to estimate expected hazards from leakage as well as to propose recommendations when designing hydrogen-related facilities.

Recently, El-Amin and coauthors [16] studied the problem of hydrogen leakage in air. In [13], they introduced boundary layer theory approach to model the concentration layer adjacent to a ceiling wall at the impinging and far regions in both planar and axisymmetric cases for small-scale hydrogen leakage. While in [46], they studied the turbulent hydrogen-air jet/plume resulted from hydrogen leakage in open air. The laminar hydrogen jet is analyzed by Sánchez–Sanz et al. [7]. Also, experimental measurements for turbulent hydrogen jet have been performed by Schefer and coauthors (e.g., [810]). On the other hand, CFD simulations of the problem have been done by many researcher such as Matssura and coauthors [1114], Kikukawa [15], Agarant et al.[16], and Swain et al. [17, 18].

Hydrogen-air jet is an example of non-Boussinesq plume; since the initial fractional density difference is high. The initial fractional density difference is defined as Δ𝜌0/𝜌=(𝜌𝜌0)/𝜌, where 𝜌0 is the initial centerline density (density at the source) and 𝜌 is the ambient density. As an example, the initial fractional density differences for selected binary low-density gases at temperature 15°C are 0.93 for H2-Air, 0.86 for He-Air, 0.43 for CH4-Air, and 0.06 for C2H2-N2. Crapper and Baines [19] suggested that the upper bound of applicability of the Boussinesq approximation is that the initial fractional density difference Δ𝜌0/𝜌 does not exceed 0.05. In general, one can say that the Boussinesq approximation is valid for small initial fractional density difference, Δ𝜌0/𝜌1 (e.g., El-Amin and Kanayama [5]). This is correct only for the case of a plume produced by a positive source of buoyancy, that is, a plume composed of fluid less dense than the ambient. For the cases where this criterion is not met, Boussinesq approximation may not be used and a density equation needs to be incorporated. El-Amin [6] introduced a numerical investigation of a non-Boussinesq, low-density gas jet (hydrogen) leaking into a high-density ambient (air). The integral models of jet fluxes are obtained and transformed into a set of ordinary differential equations of the mean centerline quantities. Therefore, mean quantities are obtained in addition to cross-stream velocity, Reynolds stresses, and turbulent Schmidt number. Furthermore, the normalized jet-feed material density and momentum flux density are correlated.

It is worth mentioning that theoretical developments and analysis of jet/plume theory were studied by a number of authors since 1950s (see, e.g., Morton et al. [20]; Morton [21]; Morton and Middleton [22]; Delichatsios [23]; Rooney and Linden [24]; Hunt and Kaye [25, 26]; Carlotti and Hunt [27]). Recently, Michaux and Vauquelin [28] developed analytical solutions for centerlines quantities of turbulent plumes rising from circular sources of positive buoyancy in a quiescent environment of uniform density for both Boussinesq and non-Boussinesq cases.

In this paper, semianalytical solution and theoretical analysis are developed for round hydrogen jet leaking into air based on Michaux and Vauquelin [28]. It is assumed that the rate of entrainment is a function of the plume centerline velocity and the ratio of the mean plume and ambient densities. Analytical expressions for all plume variables (radius, velocity, and density deficit) in terms of plume function for a given source parameter are derived.

2. Mathematical Analysis and Similarity Solutions

2.1. Governing Equations

Consider a vertical axisymmetric hydrogen-air buoyant jet resulting from a small-scale hydrogen leakage in the air with a finite circular source. Using cylindrical polar coordinates (𝑧,𝑟) with the 𝑧-axis vertical, the source is located at 𝑧=0. The continuity, momentum, and concentration equations in cylindrical coordinate system (Figure 1) for the steady vertical axisymmetric buoyant free jet can be written as [29]:𝜕(𝑟𝜌𝑉)+𝜕𝑟𝜕(𝑟𝜌𝑈)𝜕𝑧=0,(2.1)𝜕(𝑟𝜌𝑈𝑉)+𝜕𝜕𝑟𝑟𝜌𝑈2+𝜕𝜕𝑧𝑟𝜌𝑢𝑣𝜕𝑟=𝑔𝑟𝜌𝜌,(2.2)𝜕(𝑟𝜌𝑉𝐶)+𝜕𝑟𝜕(𝑟𝜌𝑈𝐶)+𝜕𝜕𝑧𝑟𝜌𝑣𝑐𝜕𝑟=0,(2.3)

where U is the mean streamwise velocity, and V is the mean cross-stream velocity, and C is the hydrogen concentration (mass fraction). The overbar denotes the time-averaged quantities, u, v are the components of velocity fluctuations in z, r directions, respectively, c is the concentration fluctuation, and 𝜌 is the mixture density.

On the other hand, from the experimental observations, the equations for the vertical velocity, density deficiency, and mass fraction profiles, assuming that the hydrogen-air mixture behaves as an ideal gas, are as follows (Fisher et al. [30], Hussein et al. [31], Shabbir and George [32], and Schefer et al. [9, 10]):𝑈(𝑟,𝑧)=𝑈𝑐𝑙𝑟(𝑧)exp2𝑏2𝜌(𝑧),(2.4)𝜌𝜌(𝑟,𝑧)=𝜌𝑐𝑙(𝑧)exp𝜆2𝑟2𝑏2(𝑧),(2.5)𝜌(𝑟,𝑧)𝐶(𝑟,𝑧)=𝜌𝑐𝑙(𝑧)𝐶𝑐𝑙(𝑧)exp𝜆2𝑟2𝑏21(𝑧),(2.6)𝜌=1/𝜌01/𝜌𝐶+1/𝜌,(2.7)

where 𝑈(𝑟,𝑧) and 𝜌(𝑟,𝑧) are the mean velocity and mean density at any point of the jet body; 𝑈𝑐𝑙(𝑧) and 𝜌𝑐𝑙(𝑧) are the centerline velocity and density. 𝑏(𝑧)=𝑐𝑚(𝑧𝑧0) is the jet/plume width which increases linearly with z, 𝑐𝑚 is the momentum spread rate of the jet. 𝑧0 is the virtual origin, which is the distance above/below the orifice where the flow appears to originate. The experimentally measured spread rate 𝑐𝑚 varies in the range 0.1–0.13. The buoyancy spreading factor 𝜆=𝑐𝑚/𝑐𝑐 expresses the ratio of spreading rates between the velocity and density deficiency profiles. The corresponding streamwise concentration for the axisymmetric hydrogen-air, free jet as detected experimentally by Schefer et al. [10] is given as 𝐶=𝐶𝑐𝑙exp(0.693𝑟2/𝑏2). In general, the spread rate for the concentration 𝑐𝑐 is given in the formula 𝐶=𝐶𝑐𝑙exp(𝑟2/𝑐𝑐2(𝑧𝑧0)2). In the work of Schefer et al. [10], the momentum spread rate for the case of hydrogen jet was estimated as 𝑐𝑚=0.103, from which one can find 𝑐𝑐=0.124, and 𝜆=0.832. It is well known that 𝑐𝐶𝑐𝑚, that is, velocity and density spread at different rates.

2.2. Similarity Solutions

Integrating the continuity (2.1) radially gives𝑑𝑑𝑧0||𝑟𝑈(𝑟,𝑧)𝜌(𝑟,𝑧)𝑑𝑟=𝑟𝜌𝑉(𝑟,𝑧)𝑟=||=𝑟𝑉(𝑟,𝑧)𝑟=𝜌.(2.8)Since 𝑈(𝑟,𝑧) is negligible for 𝑟>𝑏, then integrating (2.1) for 𝑏<𝑟< gives𝑏𝜕𝜕𝑟(𝑟𝑉(𝑟,𝑧)𝜌(𝑟,𝑧))𝑑𝑟=0.(2.9)This implies that||𝑟𝑉(𝑟,𝑧)𝑟==𝑏𝑉𝑒,(2.10)

where 𝑉𝑒 denotes the inflow velocity at the plume edge which is known as the entrainment velocity. Therefore, we have𝑑𝑑𝑧0𝑟𝑈(𝑟,𝑧)𝜌(𝑟,𝑧)𝑑𝑟=𝑏𝑉𝑒𝜌.(2.11)

This equation indicates that the increase in plume volume flux is supplied by a radial influx from the far field which in turn implies a flow across the plume boundary 𝑏. Batchelor [33] concluded that a vigorous entrainment of the ambient will be obtained as the density ratio tends to unity, 𝜌𝑐𝑙/𝜌1. While as the density ratio tends to zero, 𝜌𝑐𝑙/𝜌0, the entrainment falls to zero. Between these two limits, there will be a smooth transition of entrainment pattern. The experiments by Ricou and Spalding [34] suggest that the entrainment velocity may be obtained using the following formula𝑉𝑒𝜌=𝛼𝑐𝑙𝜌1/2𝑈𝑐𝑙,𝑉𝑒𝜌𝑐𝑙𝜌1=𝛼𝑈𝑐𝑙,𝑉𝑒𝜌𝑐𝑙𝜌0=0,(2.12)where 𝛼 is the entrainment coefficient.

Also, Morton [35] assumed that the rate of entrainment into a strongly buoyant plume is a function of both density ratio 𝜌𝑐𝑙/𝜌 and Reynolds stresses which have a magnitude proportional to 𝜌𝑐𝑙𝑈2𝑐𝑙. Therefore, the local entrainment velocity may be obtained as 𝛼(𝜌𝑐𝑙/𝜌)1/2𝑈𝑐𝑙, which has also been suggested by Thomas [36], Steward [37], and Townsend [38]. Therefore, (2.11) can be written in the form𝑑𝑑𝑧02𝜋𝑟𝑈(𝑟,𝑧)𝜌(𝑟,𝑧)𝑑𝑟=2𝜋𝑏𝜌𝛼𝜌𝑐𝑙(𝑧)𝜌1/2𝑈𝑐𝑙(𝑧).(2.13)

For calculating the momentum flux, let us integrate (2.2) with respect to 𝑟, from 𝑟=0 to 𝑟=, noting that |𝑟𝜌𝑈𝑉|0=0and |𝑟𝜌𝑢𝑣|0=0, we get𝑑𝑑𝑧02𝜋𝑟𝑈2(𝑟,𝑧)𝜌(𝑟,𝑧)𝑑𝑟=0𝜌2𝜋𝑟𝑔𝜌(𝑟,𝑧)𝑑𝑟.(2.14)

Similarly, for concentration flux, integrating (2.3) with respect to 𝑟 from 𝑟=0 to 𝑟=, noting also that |𝑟𝜌𝐶𝑉|0=0and |𝑟𝜌𝑢𝑐|0=0, we get𝑑𝑑𝑧02𝜋𝑟𝑈(𝑟,𝑧)𝜌(𝑟,𝑧)𝐶(𝑟,𝑧)𝑑𝑟=0.(2.15)

This equation may be equivalent to the buoyancy flux equation which can be written in the following form [19]:𝑑𝑑𝑧0𝜌2𝜋𝑟𝑈(𝑟,𝑧)𝜌(𝑟,𝑧)𝑑𝑟=0.(2.16) Substituting (2.4)–(2.6) into (2.16), one obtains𝑑𝑏𝑑𝑧2𝑈𝑐𝑙𝜌(𝑧)1𝑐𝑙(𝑧)𝜌=0.(2.17) Substituting (2.5), (2.6), and (2.17) into (2.13), (2.14), and (2.15) gives𝑑𝑏𝑑𝑧2(𝑧)𝑈𝑐𝑙𝜌(𝑧)𝑐𝑙(𝑧)𝜌=2𝑏𝑉𝑒,(2.18)𝑑𝑏𝑑𝑧2(𝑧)𝑈2𝑐𝑙𝜌(𝑧)𝑐𝑙(𝑧)𝜌=𝜌𝜌𝑐𝑙(𝑧)𝜌𝑔𝑏2(𝑧).(2.19)

In homogeneous surroundings of density 𝜌, the density deficit flux, (2.17), is equivalent to𝐵=𝜋𝑔𝑏2(𝑧)𝑈𝑐𝑙𝜌(𝑧)1𝑐𝑙(𝑧)𝜌,(2.20)

which has the dimension of a buoyancy flux.

Recently, Michaux and Vauquelin [28] developed analytical solutions for centerlines quantities of turbulent plumes. Now, let us introduce a modified radius 𝛽 and a dimensionless density deficit 𝜂 [28]:𝜌𝛽(𝑧)=𝑐𝑙(𝑧)𝜌1/2𝜌𝑏(𝑧),𝜂(𝑧)=𝜌𝑐𝑙(𝑧)𝜌𝑐𝑙.(𝑧)(2.21)

Therefore, (2.17)–(2.19) can be rewritten in the following form:𝑑𝛽𝑑𝑧2𝜂𝑈𝑐𝑙𝑑=0,𝛽𝑑𝑧2𝑈𝑐𝑙=2𝛼𝛽𝑈𝑐𝑙𝑑𝛽𝑑𝑧2𝑈2𝑐𝑙=𝜂𝑔𝛽2.,(2.22)

In the previous work of El-Amin and Kanayama [5] and El-Amin [6], they developed the similarity formulation and numerical solutions of the centerline quantities such as velocity and concentration.

In the current work, we follow the work of Michaux and Vauquelin [28] to obtain analytical/semianalytical solutions for centerline plume quantities. Using 𝛽(𝑧)=𝐶𝛽𝑧𝑚,𝑈𝑐𝑙(𝑧)=𝐶𝑈𝑐𝑙𝑧𝑛, and 𝜂(𝑧)=𝐶𝜂𝑧𝑝, the constants 𝐶𝛽,𝐶𝑈𝑐𝑙,and 𝐶𝜂; exponents m, n, p can be determined to obtain similarity solutions for 𝛽,𝑈𝑐𝑙, and 𝜂 as𝛽(𝑧)=6𝛼5𝑧,𝑤3(𝑧)=41/36𝛼52/3𝐵𝜋1/3𝑧1/3,1𝜂(𝑧)=𝑔341/36𝛼54/3𝐵𝜋2/3𝑧5/3,(2.23)

𝑧 in these expressions may be replaced by 𝑧𝑧0 to adapt solutions at near-source region and 𝑧0 is virtual origin.

2.3. Plume Function and Source Parameter

Now, let us introduce the plume function as [28]:Γ(𝑧)=5𝑔8𝛼𝜂𝛽𝑈2𝑐𝑙.(2.24)

At source (𝑧=0), Γ(0)=Γ0, corresponding to the source parameter initially introduced by Morton [21] and defined as [26]:Γ0=5𝑄2𝐵𝑓4𝛼𝑀5/2,(2.25)

where Q, M, and 𝐵𝑓 are the initial values of specific mass flux, specific momentum flux, and specific buoyancy flux, respectively, defined as1𝑄=4𝜋𝑑2𝑢0,1𝑀=4𝜋𝑑2𝑢20,𝐵𝑓=𝑔𝑄Δ𝜌𝜌,(2.26)

𝑑 is the inlet diameter, 𝑈0 is velocity at source, and Δ𝜌 is the difference in density between the receiving fluid and the fluid being discharged. Based on source parameter value Morton and Middleton [22] have categorized plumes with positive buoyancy as simple (pure) plume (Γ0=1), forced plume (Γ0<1), and lazy plume (Γ0>1). Other possibilities (Hunt and Kaye [25]) for Γ0=0, flow is pure jet without buoyancy, and, for Γ0<0, flow is weak fountains (negative buoyancy).

For hydrogen-air plume the source parameter is Γ01 (of order 10−4), so, based on Morton et al. [20] classification, it is a forced plume.

Equation (2.22) can be rewritten in terms of Γ as follow:𝑑𝛽=𝑑𝑧4𝛼552Γ,(2.27)𝑑𝑈𝑐𝑙𝑑𝑧=8𝛼5𝑈𝑐𝑙𝛽54Γ,(2.28)𝑑𝜂𝑑𝑧=16𝛼2𝑈5𝑔𝑐𝑙𝛽2Γ.(2.29)

Using plume function Γ, (2.24), and (2.27)–(2.29), we may write𝑑Γ=𝑑𝑧4𝛼Γ𝛽(1Γ).(2.30)

One can deduce that for Γ0<1, Γ increases monotonically with height and tends asymptotically toward unity.

Using (2.27) and (2.30), one may get𝑑𝛽𝛽=12𝑑ΓΓ+310𝑑Γ1Γ.(2.31) Integrating this equation subject to the source conditions, one obtains𝛽𝛽0=ΓΓ01/21Γ01Γ3/10,(2.32) therefore,𝑏𝑏0=𝜌0𝜌𝑐𝑙1/2ΓΓ01/21Γ01Γ3/10.(2.33) Similarly, we can find that𝑑𝑈𝑐𝑙𝑈𝑐𝑙1=2𝑑ΓΓ110𝑑Γ1Γ,(2.34) therefore,𝑈𝑐𝑙𝑈0=Γ0Γ1/21Γ1Γ01/10.(2.35) Finally, using (2.24), (2.32), and (2.34), we get𝜂𝜂0=Γ0(1Γ)Γ1Γ01/2,(2.36) or𝜌0𝜌𝜌𝑐𝑙𝜌𝜌0𝜌𝑐𝑙=Γ0(1Γ)Γ1Γ01/2.(2.37) Now, Γ is a function of z to relate each plume variable to z. Substituting (2.32) into (2.30), one gets𝑑Γ=1𝑑𝑧Λ0Γ1/2(1Γ)13/10,(2.38) where Λ0=(𝛽0/4𝛼)(|1Γ0|3/10/Γ01/2), and 𝛽0=(𝜌0/𝜌)1/2𝑏0 is the characteristic length defined from initial plume condition. For the case under consideration, we find Λ0=0.17.

Integrating (2.38), one obtains𝑧Λ0=ΓΓ0𝛾1/2(1𝛾)13/10Γ𝑑𝛾=(Γ)0.(2.39)

The above integration function has no explicit form, so Michaux and Vauquelin [28] computed and tabulated the integral function (𝑋) for several values of X. Unfortunately, they did not provide very small values as we find in this investigation when the source parameter is Γ01 (of order 10−4). Alternatively, we can write the integration function as(Γ)=ΓΓ0𝛾1/2(1𝛾)13/10=3𝑑𝛾4Γ21/2𝐹13,1102;32;Γ103Γ1/2(Γ11Γ)7/10,(2.40)

where 2𝐹1 is the hypergeometric function defined by2𝐹1(𝑎,𝑏;𝑐;Γ)=𝑘=0(𝑎)𝑘(𝑏)𝑘(𝑐)𝑘Γ𝑘𝑘!.(2.41)

But Γ1 is very small (of order 10−4) for the case of hydrogen-air plume, so two terms of the above series can approximate the function. Thus,2𝐹1(𝑎,𝑏;𝑐;Γ)=1+𝑎𝑏𝑐Γ.(2.42)

This may lead to3(Γ)=4Γ1/2Γ1++10103Γ1/2(1Γ)3/10.(2.43)

According to this relation, (Γ0)=0.039294, therefore,𝑧(Γ)=0.039294+0.17.(2.44)

From (2.43) and (2.44), one gets𝑧=30.174Γ1/2Γ1++10103Γ1/2(1Γ)3/100.039294.(2.45)

Again from (2.43), we can determine values of Γ from (Γ) given by (2.44). It is clear that Γ decreases as (Γ) increases and (Γ) increases as z increases. Thus, Γ decreases as z increases and the maximum value of Γ is located at 𝑧=0, which is equal to Γ0. Therefore, using (2.36), one can find 𝑈𝑐𝑙/𝑈0=1, and the maximum velocity is 𝑈𝑐𝑙=𝑈0.

Finally, by combining (2.4)–(2.6) with (2.33), (2.36), and (2.37), one can obtain analytical expressions for vertical velocity, density deficiency, and mass fraction for hydrogen-air mixture in terms of the universal variable Γ as follows:𝑈(𝑟,𝑧)=𝑈0Γ0Γ1/21Γ1Γ01/10×exp𝑟2Γ0/Γ(1Γ)/(1Γ0)3/5𝑏20𝜌/𝜌0𝜌/𝜌0Γ10/Γ1/2(1Γ)/(1Γ0)1/2,+1(2.46)𝜌(𝑟,𝑧)=𝜌𝜌𝜌𝜌/𝜌0Γ10/Γ1/2(1Γ)/(1Γ0)1/2+1×exp𝜆2𝑟2Γ0/Γ(1Γ)/(1Γ0)3/5𝑏20𝜌𝜌0𝜌/𝜌0Γ10/Γ1/2(1Γ)/(1Γ0)1/2,+1(2.47)𝐶𝜌(𝑟,𝑧)=0𝜌𝜌0+𝜌0𝜌𝜌𝜌(𝑟,𝑧)𝜌0.(2.48)

3. Conclusion

This paper introduces the reader to a set of features of hydrogen-air plume, which is very important to assess the potential hazard resulting from hydrogen sources upon leakage into the ambient atmosphere. Throughout this work, we derived profiles of the mean quantities for a turbulent hydrogen-air plume. These mean quantities, such as plume radius, velocity, and density deficit, are expressed in terms of the plume function for a given source parameter. These quantities are determined by integral relations and by analysis using similarity variables. Therefore, mean quantities are expressed solely in terms of the plume function and the source parameter. The plume function is valid for a range of small values of Γ as required for hydrogen-air plume. The hypergeometric function is exploited, and the profiles of the mean quantities are obtained. These results may be generalized and extended in a future work to cover more complex flow types found in a prober accident of hydrogen leaks such as leakage in hydrogen station.