# Theoretical Analysis and Semianalytical Solutions for a Turbulent Buoyant Hydrogen-Air Jet

Authors: M. F. El-Amin; Shuyu Sun; Amgad Salama
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

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).


---

## Body

## 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 [1–6] studied the problem of hydrogen leakage in air. In [1–3], 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 [4–6], 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., [8–10]). On the other hand, CFD simulations of the problem have been done by many researcher such as Matssura and coauthors [11–14], 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(z,r) with the z-axis vertical, the source is located at z=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]:(2.1)∂(rρV)∂r+∂(rρU)∂z=0,(2.2)∂(rρUV)∂r+∂(rρU2)∂z+∂(rρuv¯)∂r=gr(ρ-ρ∞),(2.3)∂(rρVC)∂r+∂(rρUC)∂z+∂(rρvc¯)∂r=0,Figure 1
Schematic diagram of turbulent hydrogen-air jet.whereU 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, rdirections, 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]):(2.4)U(r,z)=Ucl(z)exp⁡(-r2b2(z)),(2.5)ρ∞-ρ(r,z)=(ρ∞-ρcl(z))exp⁡(-λ2r2b2(z)),(2.6)ρ(r,z)C(r,z)=ρcl(z)Ccl(z)exp⁡(-λ2r2b2(z)),(2.7)ρ=1([(1/ρ0)-(1/ρ∞)]C+(1/ρ∞)),whereU(r,z) and ρ(r,z) are the mean velocity and mean density at any point of the jet body; Ucl(z) and ρcl(z) are the centerline velocity and density. b(z)=cm(z-z0) is the jet/plume width which increases linearly with z,cm is the momentum spread rate of the jet. z0 is the virtual origin, which is the distance above/below the orifice where the flow appears to originate. The experimentally measured spread rate cm varies in the range 0.1–0.13. The buoyancy spreading factor λ=cm/cc 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 C=Cclexp(-0.693r2/b2). In general, the spread rate for the concentration cc is given in the formula C=Cclexp(-r2/cc2(z-z0)2). In the work of Schefer et al. [10], the momentum spread rate for the case of hydrogen jet was estimated as cm=0.103, from which one can find cc=0.124, and λ=0.832. It is well known that cC≠cm, that is, velocity and density spread at different rates.


### 2.2. Similarity Solutions

Integrating the continuity (2.1) radially gives(2.8)ddz∫0∞rU(r,z)ρ(r,z)dr=-rρV(r,z)|r=∞=-rV(r,z)|r=∞ρ∞.Since U(r,z) is negligible for r>b, then integrating (2.1) for b<r<∞ gives(2.9)∫b∞∂∂r(rV(r,z)ρ(r,z))dr=0.This implies that(2.10)-rV(r,z)|r=∞=bVe,whereVe denotes the inflow velocity at the plume edge which is known as the entrainment velocity. Therefore, we have(2.11)ddz∫0∞rU(r,z)ρ(r,z)dr=bVeρ∞.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 boundaryb. Batchelor [33] concluded that a vigorous entrainment of the ambient will be obtained as the density ratio tends to unity, ρcl/ρ∞→1. While as the density ratio tends to zero, ρcl/ρ∞→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(2.12)Ve=α(ρclρ∞)1/2Ucl,Ve(ρclρ∞⟶1)=αUcl,Ve(ρclρ∞⟶0)=0,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 ρcl/ρ∞ and Reynolds stresses which have a magnitude proportional to ρclUcl2. Therefore, the local entrainment velocity may be obtained as α(ρcl/ρ∞)1/2Ucl, which has also been suggested by Thomas [36], Steward [37], and Townsend [38]. Therefore, (2.11) can be written in the form(2.13)ddz∫0∞2πrU(r,z)ρ(r,z)dr=2πbρ∞α(ρcl(z)ρ∞)1/2Ucl(z).For calculating the momentum flux, let us integrate (2.2) with respect to r, from r=0 to r=∞, noting that |rρUV|0∞=0and |rρuv¯|0∞=0, we get(2.14)ddz∫0∞2πrU2(r,z)ρ(r,z)dr=∫0∞2πrg(ρ∞-ρ(r,z))dr.Similarly, for concentration flux, integrating (2.3) with respect to r from r=0 to r=∞, noting also that |rρCV|0∞=0and |rρuc¯|0∞=0, we get(2.15)ddz∫0∞2πrU(r,z)ρ(r,z)C(r,z)dr=0.This equation may be equivalent to the buoyancy flux equation which can be written in the following form [19]:(2.16)ddz∫0∞2πrU(r,z)(ρ∞-ρ(r,z))dr=0.
Substituting (2.4)–(2.6) into (2.16), one obtains(2.17)ddz(b2Ucl(z)(1-ρcl(z)ρ∞))=0.
Substituting (2.5), (2.6), and (2.17) into (2.13), (2.14), and (2.15) gives(2.18)ddz(b2(z)Ucl(z)ρcl(z)ρ∞)=2bVe,(2.19)ddz(b2(z)Ucl2(z)ρcl(z)ρ∞)=ρ∞-ρcl(z)ρ∞gb2(z).In homogeneous surroundings of densityρ∞, the density deficit flux, (2.17), is equivalent to(2.20)B=πgb2(z)Ucl(z)(1-ρcl(z)ρ∞),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]:(2.21)β(z)=(ρcl(z)ρ∞)1/2b(z),η(z)=ρ∞-ρcl(z)ρcl(z).Therefore, (2.17)–(2.19) can be rewritten in the following form:(2.22)ddz(β2ηUcl)=0,ddz(β2Ucl)=2αβUclddz(β2Ucl2)=ηgβ2.,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 β(z)=Cβzm,Ucl(z)=CUclzn, and η(z)=Cηzp, the constants Cβ,CUcl,and Cη; exponents m, n, p can be determined to obtain similarity solutions for β,Ucl, and η as(2.23)β(z)=6α5z,w(z)=(34)1/3(6α5)-2/3(Bπ)1/3z-1/3,η(z)=1g(34)-1/3(6α5)-4/3(Bπ)2/3z-5/3,z in these expressions may be replaced by z-z0 to adapt solutions at near-source region and z0 is virtual origin.


### 2.3. Plume Function and Source Parameter

Now, let us introduce the plume function as [28]:(2.24)Γ(z)=5g8αηβUcl2.At source(z=0), Γ(0)=Γ0, corresponding to the source parameter initially introduced by Morton [21] and defined as [26]:(2.25)Γ0=5Q2Bf4αM5/2,whereQ, M, and Bf are the initial values of specific mass flux, specific momentum flux, and specific buoyancy flux, respectively, defined as(2.26)Q=14πd2u0,M=14πd2u02,Bf=gQΔρ∞ρ∞,d is the inlet diameter, U0 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Γ0≪1 (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:(2.27)dβdz=4α5(52-Γ),(2.28)dUcldz=-8α5(Uclβ)(54-Γ),(2.29)dηdz=-16α25g(Uclβ)2Γ.Using plume functionΓ, (2.24), and (2.27)–(2.29), we may write(2.30)dΓdz=4αΓβ(1-Γ).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(2.31)dββ=12dΓΓ+310dΓ1-Γ.
Integrating this equation subject to the source conditions, one obtains(2.32)ββ0=(ΓΓ0)1/2(1-Γ01-Γ)3/10,
therefore,(2.33)bb0=(ρ0ρcl)1/2(ΓΓ0)1/2(1-Γ01-Γ)3/10.
Similarly, we can find that(2.34)dUclUcl=-12dΓΓ-110dΓ1-Γ,
therefore,(2.35)UclU0=(Γ0Γ)1/2(1-Γ1-Γ0)1/10.
Finally, using (2.24), (2.32), and (2.34), we get(2.36)ηη0=(Γ0(1-Γ)Γ(1-Γ0))1/2,
or(2.37)ρ0(ρ∞-ρcl)(ρ∞-ρ0)ρcl=(Γ0(1-Γ)Γ(1-Γ0))1/2.
Now, Γ is a function of z to relate each plume variable to z. Substituting (2.32) into (2.30), one gets(2.38)dΓdz=1Λ0Γ1/2(1-Γ)13/10,
where Λ0=(β0/4α)(|1-Γ0|3/10/Γ01/2), and β0=(ρ0/ρ∞)1/2b0 is the characteristic length defined from initial plume condition. For the case under consideration, we find Λ0=0.17.Integrating (2.38), one obtains(2.39)zΛ0=∫Γ0Γγ-1/2(1-γ)-13/10dγ=I(Γ)-I(Γ0).The above integration function has no explicit form, so Michaux and Vauquelin [28] computed and tabulated the integral function ℑ(X) for several values of X. Unfortunately, they did not provide very small values as we find in this investigation when the source parameter is Γ0≪1 (of order 10−4). Alternatively, we can write the integration function as(2.40)I(Γ)=∫Γ0Γγ-1/2(1-γ)-13/10dγ=34Γ1/2F21(310,12;32;Γ)-103Γ1/2Γ-1(1-Γ)7/10,whereF21 is the hypergeometric function defined by(2.41)F21(a,b;c;Γ)=∑k=0∞(a)k(b)k(c)kΓkk!.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.42)F21(a,b;c;Γ)=1+abcΓ.This may lead to(2.43)I(Γ)=34Γ1/2(1+Γ10)+103Γ1/2(1-Γ)3/10.According to this relation,ℑ(Γ0)=0.039294, therefore,(2.44)I(Γ)=0.039294+z0.17.From (2.43) and (2.44), one gets(2.45)z0.17=34Γ1/2(1+Γ10)+103Γ1/2(1-Γ)3/10-0.039294.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 z=0, which is equal to Γ0. Therefore, using (2.36), one can find Ucl/U0=1, and the maximum velocity is Ucl=U0.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:(2.46)U(r,z)=U0(Γ0Γ)1/2(1-Γ1-Γ0)1/10×exp⁡(-r2(Γ0/Γ)((1-Γ)/(1-Γ0))3/5b02ρ∞/ρ0((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1)),(2.47)ρ(r,z)=ρ∞-(ρ∞-ρ∞((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1))×exp⁡(-λ2r2(Γ0/Γ)((1-Γ)/(1-Γ0))3/5b02ρ∞ρ0((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1)),(2.48)C(r,z)=-ρ0ρ∞-ρ0+ρ0ρ∞ρ(r,z)(ρ∞-ρ0).


## 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(z,r) with the z-axis vertical, the source is located at z=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]:(2.1)∂(rρV)∂r+∂(rρU)∂z=0,(2.2)∂(rρUV)∂r+∂(rρU2)∂z+∂(rρuv¯)∂r=gr(ρ-ρ∞),(2.3)∂(rρVC)∂r+∂(rρUC)∂z+∂(rρvc¯)∂r=0,Figure 1
Schematic diagram of turbulent hydrogen-air jet.whereU 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, rdirections, 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]):(2.4)U(r,z)=Ucl(z)exp⁡(-r2b2(z)),(2.5)ρ∞-ρ(r,z)=(ρ∞-ρcl(z))exp⁡(-λ2r2b2(z)),(2.6)ρ(r,z)C(r,z)=ρcl(z)Ccl(z)exp⁡(-λ2r2b2(z)),(2.7)ρ=1([(1/ρ0)-(1/ρ∞)]C+(1/ρ∞)),whereU(r,z) and ρ(r,z) are the mean velocity and mean density at any point of the jet body; Ucl(z) and ρcl(z) are the centerline velocity and density. b(z)=cm(z-z0) is the jet/plume width which increases linearly with z,cm is the momentum spread rate of the jet. z0 is the virtual origin, which is the distance above/below the orifice where the flow appears to originate. The experimentally measured spread rate cm varies in the range 0.1–0.13. The buoyancy spreading factor λ=cm/cc 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 C=Cclexp(-0.693r2/b2). In general, the spread rate for the concentration cc is given in the formula C=Cclexp(-r2/cc2(z-z0)2). In the work of Schefer et al. [10], the momentum spread rate for the case of hydrogen jet was estimated as cm=0.103, from which one can find cc=0.124, and λ=0.832. It is well known that cC≠cm, that is, velocity and density spread at different rates.


## 2.2. Similarity Solutions

Integrating the continuity (2.1) radially gives(2.8)ddz∫0∞rU(r,z)ρ(r,z)dr=-rρV(r,z)|r=∞=-rV(r,z)|r=∞ρ∞.Since U(r,z) is negligible for r>b, then integrating (2.1) for b<r<∞ gives(2.9)∫b∞∂∂r(rV(r,z)ρ(r,z))dr=0.This implies that(2.10)-rV(r,z)|r=∞=bVe,whereVe denotes the inflow velocity at the plume edge which is known as the entrainment velocity. Therefore, we have(2.11)ddz∫0∞rU(r,z)ρ(r,z)dr=bVeρ∞.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 boundaryb. Batchelor [33] concluded that a vigorous entrainment of the ambient will be obtained as the density ratio tends to unity, ρcl/ρ∞→1. While as the density ratio tends to zero, ρcl/ρ∞→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(2.12)Ve=α(ρclρ∞)1/2Ucl,Ve(ρclρ∞⟶1)=αUcl,Ve(ρclρ∞⟶0)=0,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 ρcl/ρ∞ and Reynolds stresses which have a magnitude proportional to ρclUcl2. Therefore, the local entrainment velocity may be obtained as α(ρcl/ρ∞)1/2Ucl, which has also been suggested by Thomas [36], Steward [37], and Townsend [38]. Therefore, (2.11) can be written in the form(2.13)ddz∫0∞2πrU(r,z)ρ(r,z)dr=2πbρ∞α(ρcl(z)ρ∞)1/2Ucl(z).For calculating the momentum flux, let us integrate (2.2) with respect to r, from r=0 to r=∞, noting that |rρUV|0∞=0and |rρuv¯|0∞=0, we get(2.14)ddz∫0∞2πrU2(r,z)ρ(r,z)dr=∫0∞2πrg(ρ∞-ρ(r,z))dr.Similarly, for concentration flux, integrating (2.3) with respect to r from r=0 to r=∞, noting also that |rρCV|0∞=0and |rρuc¯|0∞=0, we get(2.15)ddz∫0∞2πrU(r,z)ρ(r,z)C(r,z)dr=0.This equation may be equivalent to the buoyancy flux equation which can be written in the following form [19]:(2.16)ddz∫0∞2πrU(r,z)(ρ∞-ρ(r,z))dr=0.
Substituting (2.4)–(2.6) into (2.16), one obtains(2.17)ddz(b2Ucl(z)(1-ρcl(z)ρ∞))=0.
Substituting (2.5), (2.6), and (2.17) into (2.13), (2.14), and (2.15) gives(2.18)ddz(b2(z)Ucl(z)ρcl(z)ρ∞)=2bVe,(2.19)ddz(b2(z)Ucl2(z)ρcl(z)ρ∞)=ρ∞-ρcl(z)ρ∞gb2(z).In homogeneous surroundings of densityρ∞, the density deficit flux, (2.17), is equivalent to(2.20)B=πgb2(z)Ucl(z)(1-ρcl(z)ρ∞),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]:(2.21)β(z)=(ρcl(z)ρ∞)1/2b(z),η(z)=ρ∞-ρcl(z)ρcl(z).Therefore, (2.17)–(2.19) can be rewritten in the following form:(2.22)ddz(β2ηUcl)=0,ddz(β2Ucl)=2αβUclddz(β2Ucl2)=ηgβ2.,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 β(z)=Cβzm,Ucl(z)=CUclzn, and η(z)=Cηzp, the constants Cβ,CUcl,and Cη; exponents m, n, p can be determined to obtain similarity solutions for β,Ucl, and η as(2.23)β(z)=6α5z,w(z)=(34)1/3(6α5)-2/3(Bπ)1/3z-1/3,η(z)=1g(34)-1/3(6α5)-4/3(Bπ)2/3z-5/3,z in these expressions may be replaced by z-z0 to adapt solutions at near-source region and z0 is virtual origin.


## 2.3. Plume Function and Source Parameter

Now, let us introduce the plume function as [28]:(2.24)Γ(z)=5g8αηβUcl2.At source(z=0), Γ(0)=Γ0, corresponding to the source parameter initially introduced by Morton [21] and defined as [26]:(2.25)Γ0=5Q2Bf4αM5/2,whereQ, M, and Bf are the initial values of specific mass flux, specific momentum flux, and specific buoyancy flux, respectively, defined as(2.26)Q=14πd2u0,M=14πd2u02,Bf=gQΔρ∞ρ∞,d is the inlet diameter, U0 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Γ0≪1 (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:(2.27)dβdz=4α5(52-Γ),(2.28)dUcldz=-8α5(Uclβ)(54-Γ),(2.29)dηdz=-16α25g(Uclβ)2Γ.Using plume functionΓ, (2.24), and (2.27)–(2.29), we may write(2.30)dΓdz=4αΓβ(1-Γ).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(2.31)dββ=12dΓΓ+310dΓ1-Γ.
Integrating this equation subject to the source conditions, one obtains(2.32)ββ0=(ΓΓ0)1/2(1-Γ01-Γ)3/10,
therefore,(2.33)bb0=(ρ0ρcl)1/2(ΓΓ0)1/2(1-Γ01-Γ)3/10.
Similarly, we can find that(2.34)dUclUcl=-12dΓΓ-110dΓ1-Γ,
therefore,(2.35)UclU0=(Γ0Γ)1/2(1-Γ1-Γ0)1/10.
Finally, using (2.24), (2.32), and (2.34), we get(2.36)ηη0=(Γ0(1-Γ)Γ(1-Γ0))1/2,
or(2.37)ρ0(ρ∞-ρcl)(ρ∞-ρ0)ρcl=(Γ0(1-Γ)Γ(1-Γ0))1/2.
Now, Γ is a function of z to relate each plume variable to z. Substituting (2.32) into (2.30), one gets(2.38)dΓdz=1Λ0Γ1/2(1-Γ)13/10,
where Λ0=(β0/4α)(|1-Γ0|3/10/Γ01/2), and β0=(ρ0/ρ∞)1/2b0 is the characteristic length defined from initial plume condition. For the case under consideration, we find Λ0=0.17.Integrating (2.38), one obtains(2.39)zΛ0=∫Γ0Γγ-1/2(1-γ)-13/10dγ=I(Γ)-I(Γ0).The above integration function has no explicit form, so Michaux and Vauquelin [28] computed and tabulated the integral function ℑ(X) for several values of X. Unfortunately, they did not provide very small values as we find in this investigation when the source parameter is Γ0≪1 (of order 10−4). Alternatively, we can write the integration function as(2.40)I(Γ)=∫Γ0Γγ-1/2(1-γ)-13/10dγ=34Γ1/2F21(310,12;32;Γ)-103Γ1/2Γ-1(1-Γ)7/10,whereF21 is the hypergeometric function defined by(2.41)F21(a,b;c;Γ)=∑k=0∞(a)k(b)k(c)kΓkk!.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.42)F21(a,b;c;Γ)=1+abcΓ.This may lead to(2.43)I(Γ)=34Γ1/2(1+Γ10)+103Γ1/2(1-Γ)3/10.According to this relation,ℑ(Γ0)=0.039294, therefore,(2.44)I(Γ)=0.039294+z0.17.From (2.43) and (2.44), one gets(2.45)z0.17=34Γ1/2(1+Γ10)+103Γ1/2(1-Γ)3/10-0.039294.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 z=0, which is equal to Γ0. Therefore, using (2.36), one can find Ucl/U0=1, and the maximum velocity is Ucl=U0.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:(2.46)U(r,z)=U0(Γ0Γ)1/2(1-Γ1-Γ0)1/10×exp⁡(-r2(Γ0/Γ)((1-Γ)/(1-Γ0))3/5b02ρ∞/ρ0((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1)),(2.47)ρ(r,z)=ρ∞-(ρ∞-ρ∞((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1))×exp⁡(-λ2r2(Γ0/Γ)((1-Γ)/(1-Γ0))3/5b02ρ∞ρ0((ρ∞/ρ0-1)(Γ0/Γ)1/2((1-Γ)/(1-Γ0))1/2+1)),(2.48)C(r,z)=-ρ0ρ∞-ρ0+ρ0ρ∞ρ(r,z)(ρ∞-ρ0).


## 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.



---

*Source: 710172-2012-02-20.xml*