- Email: [email protected]
Characteristics of flammable, buoyant hydrogen plumes rising from open vertical containers
international journal of hydrogen energy 34 (2009) 6568–6579
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
Characteristics of flammable, buoyant hydrogen plumes rising from open vertical containers S. Fardisi, Ghazi A. Karim* University of Calgary, Mechanical and Manufacturing Engineering Department, 2500 University Drive, NW, Calgary, Alberta, Canada T2N 1N4
article info
abstract
Article history:
The dynamics of the dispersion of a fixed mass of the highly buoyant hydrogen when
Received 24 February 2009
exposed to overlaying atmosphere with a negligible pressure difference from open vertical
Received in revised form
cylindrical enclosures are examined. Features of the rapid formation and dispersion of
6 May 2009
flammable mixtures both inside and immediate outside of the enclosure and their corre-
Accepted 10 May 2009
sponding propagation rates were examined using a 3-D CFD model. For the cases consid-
Available online 27 June 2009
ered, the puffs of the fuel–air mixture appear to produce lean flammable boundaries that move mainly at a near constant rate for much of the time. A similar simulation that used
Keywords:
an axis-symmetrical 2-D model tended to under-predict the dynamics of the lean and rich
Hydrogen safety
mixture boundaries. Hydrogen plume characteristics were compared with that of the less
Flammable mixtures
buoyant methane and helium release. Unlike methane, helium propagation rate was found
Convective dispersion
fairly close to that of hydrogen.
Buoyancy
ª 2009 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
The special thermo-chemical characteristics of hydrogen make it a very unique safety challenge. Due to its low density, low viscosity and small molecule, hydrogen is a strong leaking, fast moving gas that can be ignited much easier than most other conventional fuels. The ignition energy of the stoichiometric hydrogen/oxygen mixture at ambient conditions is about 20 mJ which is around 1/10th of that for most other conventional fuels such as methane, 0.29 mJ, and propane, 0.26 mJ [1]. In addition, the flammability limits of hydrogen in air at atmospheric conditions, 5.0–75.0% by volume, are much wider than for most other conventional fuels under the same conditions. Consequently, hydrogen can rapidly produce dangerously large flammable zones following its release from fuel installations. Accordingly, any operation that involves hydrogen has the potential to represent significant fire and
explosion hazards. Moreover it is projected that the use of hydrogen will be increased as progress is made towards the ‘‘hydrogen economy’’ [1]. Such changes are feasible only if hydrogen’s safety issues are thoroughly investigated and made comparable with those of other common fuels. Among these issues is the improvement to our understanding of the mechanism and processes of hydrogen dispersion and developing effective predictive models which should help in developing better guidelines to detect and mitigate hydrogen hazards and facilitate hydrogen acceptance by the public [2]. The isothermal and steady state dispersion of a fluid has been widely investigated, both experimentally and numerically. Karim and Tsang [3] examined the temporally changing concentration gradients of flammable mixtures of methane and air by monitoring experimentally the consequent flame propagation patterns following ignition. Meroney [4] investigated the validity of CFD techniques for modelling the
* Corresponding author. Fax: þ1 403 282 8406. E-mail address: [email protected] (G.A. Karim). 0360-3199/$ – see front matter ª 2009 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2009.05.123
international journal of hydrogen energy 34 (2009) 6568–6579
Nomenclature Cr D d g Gr h H L M R P Re Ra Sc t u0
RMS Courant number diffusion coefficient diameter gravitational acceleration Grashof number penetration height height of the fuel in the cylinder cylinder length molecular weight cylinder radius pressure reference Reynolds number Rayleigh number Schmitt number time reference velocity
dispersion of a toxic cloud. He found that the modeled clouds are qualifiedly similar in appearance to those experimentally observed and travel with similar rates. Ilegbusi and Mat [5] simulated the mixing of two fluids and compared the accuracy of different numerical approaches for modelling buoyancy driven flows. They also considered the perturbations introduced by a valve opening and concluded that when the Grashof number is relatively large, this effect is not negligible. Alhajraf et al. [6] explored the effects of local atmospheric conditions on the dispersion of hazardous chemicals from a single source. They indicated that the weather stability condition can play an important role on the dispersion processes depending on the stability level of the atmosphere. Bunama and Karim [7] examined the formation of flammable atmospheres within vertical enclosures containing a vaporizing liquid fuel and noted on different factors such as fuel type, the exposed liquid surface area and the size of the vent area. Many of the studies on hydrogen dispersion and safety have dealt with hydrogen flame characteristics [8] and hydrogen diffusion for different large-scale leakage scenarios [9]. Wilkening and Baraldi [10] compared the hydrogen release from a pipeline with that of methane. They concluded that the lower density of hydrogen is partially compensated by its higher heating value with the overall energy release is greater for hydrogen. Thus the consequences of a possible hydrogen explosion might be worse compared to that of methane. Kikukawa [11] modeled hydrogen leakage from fueling stations and visualized the effects of barriers on hydrogen diffusion. Liu and Schreiber [12] looked into hydrogen leakage in a hydrogen fueled vehicle and concluded that the ventilation system can be modified to greatly reduce the combustible hydrogen in the car interior. Venetsanos et al. [13] considered hydrogen release from an automobile in open atmosphere and in a tunnel. They concluded that a similar release of hydrogen is more severe in a tunnel. Most researchers have considered continuous and largescale hydrogen releases. Very few studies were concerned with minor, transient leakage scenarios which are often
6569
! V velocity vector w mass fraction a A constant m static viscosity y dynamic viscosity r density %Dev percent deviation CFD computational fluid dynamics RTI Rayleigh Taylor instability SIP strongly implicit procedure Subscripts f fuel a air Max maximum N infinity 3-D three-dimensional 2-D two-dimensional
technically impossible to eliminate from pipes, valves, pinholes or joints. In addition the characteristics of minor and major hydrogen releases do not appear to be similar. The previous work by Cisse and Karim [14] considered the formation and dispersion of a hydrogen cloud using a 2-D axis-symmetrical model. They compared the resulting hydrogen plume with that of the nearly buoyancy neutral ethylene, and the heavier than air propane, and noted the overwhelming importance of the fuel density relative to that of air on the dynamics of the resulting plumes. The current contribution examines the transient dispersion of a fixed mass of gaseous hydrogen which is suddenly permitted to spread with a negligible pressure difference into air at atmospheric conditions within vertical cylindrical enclosures that are open at the top to the atmosphere using the 3-D Computational Fluid Dynamics (CFD) simulation. Calculated results appear to be in good agreement with a corresponding set of experimental results and those using other numerical simulations. CFD modelling was chosen since experimenting with hydrogen is not easy due to its rapid mixing and ignition properties and the associated hazards especially within closed enclosures. Also the poor reproducibility of the tests and high costs make experimentation less practical. The associated complex temporal changes of the spatial concentrations fields and the flammable zones formed are described and compared to those of other buoyant gases especially methane and helium. Comparing the behavior of helium with that of hydrogen is important since helium is often employed to simulate the behavior of hydrogen [15] to avoid experimental difficulties and associated hazards. It is also useful to compare hydrogen with methane which is a well developed and widely used fuel resource with known safety problems. Since, there is a big difference in terms of the computational complexity and associated costs between two- and three-dimensional simulation models, the present contribution also examines whether the results of a corresponding 2-D simulation would be sufficiently accurate in comparison with
6570
international journal of hydrogen energy 34 (2009) 6568–6579
the results of the more complex and time consuming 3-D simulation for modelling hydrogen dispersion.
in which Dr ¼ ra rf, u0 is the reference velocity defined by: u20 ¼ gDr=ra L [16] and Re and Gr are the reference Reynolds and Grashof numbers, respectively.
2.
Theory and calculations
2.2.
2.1.
Problem description
A fluid with a low density cannot stably rest below a fluid with a larger density in a uniform gravitational field. This unstable situation results in the production of a flow field when the plane interface is perturbed. The mixing process is mainly driven by natural convection with molecular diffusion playing a smaller role. The resulting concentration field depends strongly on both the radial and axial directions and the bulk batches of fuel rising and forming the plume. The structure of the flow under such a condition is dominated by a special type of instabilities known as Rayleigh–Taylor Instability (RTI), described in detail in the relevant literature, e.g. [17 and 18]. Fig. 2 shows the modeled stages of the 3-D hydrogen fingering into air as a result of RTI instabilities with Gr0 ¼ 3.6 106, Re0 ¼ 1500. Some wrinkles first appear on the interface between the fuel and air. These instabilities grow and may neck at some point to release a batch of fuel to the overlaying atmosphere. Then a plume begins to form which later develops into a ‘puff’ of mixture rich in fuel that heads to the top of the cylinder. Fig. 3 depicts similar results for higher Reynolds number of Gr0 ¼ 3.6 106, Re0 ¼ 6000 with the larger diameter of d ¼ 40 mm. Obviously much more complex velocity and concentration fields will develop while similar stages of fingering process are observed. Whenever there is a continuous leakage or release of a fuel, it takes a relatively long time to have a stable concentration field established in the compartment. This does not mean that
As shown schematically in Fig. 1, the lower section of the open to atmosphere circular cylindrical vertical vessel of diameter d and length L is assumed to initially contain a fixed mass of hydrogen while the remaining upper part is filled with air. It is assumed that at a certain instant of time, assigned time zero, hydrogen is exposed to the air to commence spreading and dispersing into the overlaying atmosphere under isothermal conditions. Transient velocity and concentration fields begin to develop through the coupled transport processes of mass and momentum. Mass is transferred by the combined effects of molecular diffusion and natural convection. Molecular diffusion is driven by the local concentration gradients while, for lighter than air gases such as hydrogen (M ¼ 2 kg/kmol), helium (M ¼ 4 kg/kmol) or methane (M ¼ 16 kg/kmol) the convective mass transfer which is the bulk flow of the moving gaseous mixture due to the buoyancy effects is dominant. The variation of the concentration field of the different species within the cylinder at any time instance and location is to be determined so as to establish the development of the flammable regions. These are represented by the mixture portions where the fuel concentration lies between those corresponding to the values of the lean and rich limits, under the prevailing conditions. Particular interest is to be focused on the early stages of the evolution of the resulting plume. The key dimensionless numbers involved for this configuration were defined as follows: Re ¼
Gr ¼
ra u0 D ma gDr L3 ra y2a
Fig. 1 – Schematic representation of the set up.
Mixing stages
(1)
(2)
Fig. 2 – 3-D Transient development of the lean limit concentration surface, Gr0 [ 3.6 3 106, Re0 [ 1500, d [ 10 mm, L [ 50 mm, and H [ 20 mm.
international journal of hydrogen energy 34 (2009) 6568–6579
6571
fraction contours, was shown by Yang et al. [17] to grow quadratically with time when assuming a self-similar fluid flow and constant density difference across the advancing fronts. Thus, the penetration height, h, is given by: h¼a
Fig. 3 – 3-D Transient development of the lean surface, Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
there is no fire hazard during the transition since each puff of fuel produces a transient flammable zone. This phenomenon is extremely important when considering safety issues of buoyant fuel releases. Fig. 4 demonstrates the concentration fields at a particular moment for different diameters. Clearly, the resolved concentration fields become much more complex as the diameter increases due to different interactions between the developing fingerings. For example a rising finger might deviate from normal direction, become unstable and collapse onto another finger. Two fingers may unite or, a single finger might bifurcate into several smaller sub-fingers. These complexities are associated with multiple simultaneous plumes and their interactions, [19]. As it can be seen for methane fingering process in Fig. 4, the mean penetration of the fingers, unlike the fingering structure, seems to be independent of the diameter of the container. To reduce the complexities associated with multiple plumes, typical representative results are to be shown for the case of a vertical cylinder of 40 mm diameter where the free evolution of the plume is not affected significantly by the side walls and being confined to a single fuel finger plume. The fuel penetration height, which is defined as the maximum distance between the 0.01 and 0.99 fuel molar
ra rf 2 gt ra þ rf
(3)
where ‘‘a’’ is a constant that was shown experimentally to be approximately 0.06 and independent of the density ratios [17] and g is the gravitational acceleration. Since the density difference is only constant between immiscible fluids, even though the flow may not reach a selfsimilar condition, the quadratic growth relationship is just used for RTI with immiscible fluid. Miscible gaseous fluids involve mixing in the molecular scale that renders the buoyancy force gradually weaker. After some relatively long time, the concentration gradients become so small that the buoyancy force can no longer establish a velocity field. From this stage on, it is only molecular diffusion that controls the mixing processes and consequently the penetration length varies approximately with the square root of time. For miscible fluids such as gases, the penetration height grows approximately quadratically with time at the early stages when the diffusion time scale is much smaller than the momentum time scale and diffusion effects are negligible compared to buoyancy effects. After the lapse of some relatively long period of time buoyancy effects almost vanish and the problem becomes essentially diffusion driven with the penetration height tending to grow with the square root of time. Within the transition stage the behaviour is seen to be changing gradually from a quadratic to a square root dependence on time. Experimental and numerical studies suggest that the miscible penetration rate approaches a limiting velocity which remains relatively constant during most of the transition stage [18]. Thus, fingers of the fuel rich mixture travel at an essentially constant speed leaving a chaotic mixing zone full of vortices behind. In the present work, the evolution of the flammable mixture plumes is examined so as to improve the understanding of the development of the fire and explosion hazard in related situations. Interest is focused on the vertical propagation velocity of the lean and rich limit flammable boundaries, which are the key safety parameters, for hydrogen as compared to methane and helium. The relative accuracy of the results of the 2-D and 3-D axis-symmetric modelling was evaluated by comparing their predictions of the limit vertical velocity with the experimentally derived values.
2.3.
Modelling
The transient equations of the conservation of mass, momentum and energy are the governing equations which are essentially the transient compressible Navier–Stokes equations coupled with a species transport equation.
Fig. 4 – Methane fingering process for top diameter of one meter, twenty centimeters and four centimeters.
vr / ! þ V $r V ¼ 0 vt
(4)
! / v V ! /! ! m // ! þ V $ V V ¼ V PþmV2 V þ V V $ V þr! g r vt 3
(5)
6572
international journal of hydrogen energy 34 (2009) 6568–6579
/ / vw ! / r þ ð V $ V Þw ¼ V $ðrd V wÞ vt
(6)
The thermodynamic and transport properties are taken to change continuously all over the domain of calculation according to the transient local values of concentration and temperature [20]. The flow of the Newtonian fluid is assumed to be isothermal and non reactive for the conditions examined in this contribution. Two distinct numerical codes were employed to solve the equations to ensure that the results are sufficiently independent of the computational method. A FORTRAN code developed by Cisse [14] was used in which the governing equations were discretized using a staggered non-uniform grid. The advection terms are approximated using the QUICK scheme combined with the flux limiter ULTRA-SHARP to eliminate the numerical diffusion errors and the non-physical oscillations [21]. The enhanced SIMPLE algorithm [22] is used to couple the momentum and continuity equations. The algebraic equations resulting from the discretized equations are solved iteratively using the Strongly Implicit Procedure (SIP). In addition, Ansys-CFX Commercial software, which uses a coupled scheme and a number of modifying subroutines are used for both the 3-D and the 2-D axis-symmetrical simulation. Advection terms are approximated using a second order QUICK scheme. Also, a second order discrimination is employed to deal with the transient terms. The transition value of the Rayleigh number for natural convection in a cavity was found to be around 107–109 [23]. This number corresponds to the release of hydrogen in a vessel of 20 mm diameter at atmospheric conditions. This supports the fact that most buoyancy driven problems that involve hydrogen may be turbulent in nature. Thus, for 3-D cases a large eddy simulation (LES) with Smagorinsky subgrid scale model [10] was employed to account for turbulence effects, in parallel with a buoyancy-modified K–3 model [10] for 2-D cases that do not allow LES. A special treatment extended the calculation domain beyond the physical limits of the cylinder typically to five diameters in the radial direction and five cylinder lengths in the vertical direction outside, as depicted in Fig. 5. This was made so as to examine the
consequences of the emergence of hydrogen into the immediate vicinity of the cylinder which has important practical safety implications. A structured non-uniform grid is used for both the 2-D and 3-D simulations. In order to check the accuracy of the results, a series of mesh refinement tests were carried out. It was observed that due to the pronounced effects of the walls, the small diameter cases were stiffer to solve and needed a more robust numerical treatment. Thus, the number of iterations per time step for 3-D cases was set to 100 for small diameter cases in comparison to 50 for other larger diameter problems. This way the same convergence behaviour was established for all problems. The number of iterations per time step for 2-D cases was set to 200 for small diameters in comparison to 100 for other larger diameter problems. The convergence criterion was set to be 105 for both 3-D and 2-D simulations. In order to ensure that the time step is small enough to minimize numerical errors, time step independence tests were carried out. Simulations were performed by choosing an initial time step and then running a small period of the simulation in physical time. This initial time step was estimated [22] according to the RMS Courant number defined by: Cr ¼
Vdt dX
(7)
where, V is a characteristic velocity, dt is the time step and dX is the dimension of the grid cell. As an example for hydrogen diffusing into air, a maximum velocity of roughly 200 mm/s was expected. Accordingly, the time step was chosen to be 0.2 ms for the grid step used to ensure a Currant number below unity and a stable solution.
3.
Results
3.1. Characteristics of flammable hydrogen plume, 3-D versus 2-D simulation Cases of a finite quantity of the pure hydrogen diffusing upwards into an overlaying atmosphere of air within the
Fig. 5 – The 2-D and 3-D geometrics showing the extended control volume used in the simulation.
international journal of hydrogen energy 34 (2009) 6568–6579
cylinder of Fig. 1 were considered. The fuel fingering processes were resolved in 3-D simulations of upward motion of the gaseous fuel. In addition, 2-D axis-symmetric simulations were performed using both the FORTRAN code and the AnsysCFX code which were in very good agreement with each other. Figs. 6 and 7 show the 3-D and 2-D simulation results of the dispersion of a certain mass of hydrogen into air, respectively and illustrate the stages of mixture development from interface wrinkling till developing an effective constant vertical velocity. It can be seen that the agreement is very poor and also that the resolved 3-D concentration fields are more complex than those of the 2-D simulations. This would indicate that imposing symmetry via the 2-D simulations is unjustified in general for such a basically 3-D phenomenon when a more authentic representation of all aspects of the phenomenon is aimed at. Figs. 8 and 9 show the flammable zones for the corresponding cases of Figs. 6 and 7, which are bounded outwardly by the lean limit concentration and inwardly by the rich limit. The flammable regions grow from being initially very thin and confined to a narrow band along the interface between the hydrogen and air, to thicken later on while moving upwards. The formation and decay of the flammable mixture zones are strongly dependant on the dynamics of the fuel fingering process. It is evident that the evolution of the flammable zone is dominated by the growth and mobility of the rising fuel and the fingering structure. The flammable region which does not necessarily always extend over the whole cross-sectional area of the cylinder is of irregular shape and develops both inside and outside the cylinder. Moreover, the flammable zone can become fragmented. This would have an important implication about the characteristic spread of the fire hazard within and outside of these confined configurations. It was observed that the 2-D simulation tends to be adequately accurate in situations when the 3-D results display near symmetrical behaviour and the concentration field is not so complex. As an example, the dispersion of the fuel is effectively resolved by 2-D axis-symmetric modelling under low gravitational acceleration or for a small density difference
6573
between the fuel and air. Figs. 10 and 11 show the 3-D and 2-D simulation results of the dispersion of a certain mass of hydrogen into air in an atmosphere with a reduced gravitational acceleration of g ¼ 1.0 m/s2 with their corresponding flammable zone depicted in Figs. 12 and 13, respectively. Obviously, the flow field is much simpler compared to when g ¼ 10 m/s2. In addition, there is an acceptable level of agreement between the two sets of solutions which indicates that when the buoyancy force is weaker and the symmetry assumption is more valid, 2-D simulation gives an accurate solution of the problem. It was noticed that the rich flammability limit rapidly disappears after the dispersion, as shown in Figs. 8 and 12. This phenomenon was not observed in methane dispersion results and may be due to very large rich flammability limit, 0.75, as compared to that of methane, 0.15, and the relatively faster plume motion of hydrogen. This shows that while a methane flammable region is limited by the rich flammability limit boundary, the hydrogen flammable region extends down to the base of the vessel and covers a much larger volume. It was also observed that the life time of the flammable volume is over predicted by the 2-D model by up to 10% for both Hydrogen and Methane. This is since the 2-D model underestimates the propagation rate and it is reasonable to expect that it would overestimate the dissipation time. The rapid disappearance of the rich limit and the time required after the spill before the region is safe and noncombustible depends strongly on the initial amount of fuel to be dispersed from the vessel, the molecular weight and viscosity of the fuel and its diffusion coefficient. Further research is needed while considering a wider variety of situations to establish in more detail the behavior of the flammable region within and outside containers.
3.2. Comparison of hydrogen plume dynamics with that of methane and helium Apart from the structure of the fluid flow, it is the velocity of the fuel-first-arrival and the lean limit boundary that are
Fig. 6 – 3-D simulation results of hydrogen dispersing into air, (a) Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1 and 0.5, respectively, Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
6574
international journal of hydrogen energy 34 (2009) 6568–6579
Fig. 7 – Corresponding 2-D Simulation results of hydrogen dispersing into air, Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5 fuel molar fractions, respectively Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
particularly important for safety analysis. The growth and spread pattern of the fuel-first-arrival-boundary taken as the 1% fuel molar fraction contours was monitored and shown in Fig. 14 for both 3-D and 2-D simulations. This boundary appears to accelerate at the beginning of its motion to reach effectively a limiting vertical velocity. The 2-D axis-symmetric model appears to underestimate this velocity by up to 22% and 17% for hydrogen and methane, respectively for the cases examined. This would indicate that the 2-D models may be of insufficient accuracy compared to the 3-D counterpart in predicting the vertical dispersion rates and patterns, for both hydrogen and methane. The observed slower mixture evolution in the 2-D simulation may be due to the known tendency
of large-scale features to grow more rapidly in a 3-D domain than in a 2-D one. Moreover, a fuel concentration front can be considered to be similar to a moving object in the flow and it is known that 2-D models tend to overestimate the drag force on submerged objects confining the RTI and forcing them to move slower. It is also interesting to note in Fig. 14 that the helium first-arrival-boundary is fairly close to that of hydrogen. This supports the view that helium may be adequately used instead of hydrogen in safety experiments to reduce the difficulties and hazards as has been emphasized by others [15]. The observed velocity of the fuel-first-arrival-boundary obtained from the 3-D simulation along with 3-D deviation is
Fig. 8 – Development of the corresponding 3-D flammable zones.
international journal of hydrogen energy 34 (2009) 6568–6579
6575
Fig. 9 – Development of the corresponding 2-D flammable zones.
plotted in Fig. 15 against the analytical-experimental results given in reference [18]. The experimental value of the limit mixture velocity appears to be sufficiently close to the 3-D constant velocity growth pattern predicted. This confirms the accuracy of the numerical solution for solving this problem. Unlike the fuel-first-arrival-boundary, i.e. the 1% fuel molar fraction front there is little experimental and numerical
data on the lean flammable boundary, i.e. 5% fuel molar front for hydrogen in air, the boundary which is more important in terms of safety issues. Accordingly, focus is made on the dynamics of this specific concentration surface. The lean flammable boundary showed a similar behaviour to that of the fuel-first-arrival-boundary (1% fuel molar front) and had a vertical velocity which was lower than the velocity of the
Fig. 10 – Corresponding 3-D simulation results of hydrogen dispersing into air, (a) Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5, respectively, Gr0 [ 1.7 3 104, Re0 [ 330, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
6576
international journal of hydrogen energy 34 (2009) 6568–6579
Fig. 11 – 2-D Simulation results of hydrogen dispersing into air, Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5 fuel molar fractions, respectively, Gr0 [ 1.7 3 104, Re0 [ 330, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
first-arrival-boundary by 5–10%. This velocity, as stated in Section 2.2 does not remain constant but after some time subsides towards zero. However, for the consideration of the early stages of the potential spread of the fire hazard, this constant value of the lean flammable front needs to be considered.
Figs. 16 and 17 display the 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of gravitational acceleration and the molecular weight of the fuel, respectively. Fig. 16 was produced using 6 simulation points at M ¼ 2(hydrogen), 4(helium), 16(methane), 8, 22 and 27. M ¼ 8, 22 and 27 was considered to have been the result of
Fig. 12 – Development of the corresponding 3-D flammable zones.
international journal of hydrogen energy 34 (2009) 6568–6579
6577
Fig. 13 – Development of the corresponding 2-D flammable zones.
homogenous mixing of proper portions of hydrogen and air and values of the corresponding fluid and flow properties such as diffusion coefficient, viscosity, etc. were obtained accordingly. Thus, the 5% molar fraction front represents the lean limit for only methane and not for other gases considered here. Fig. 17 was produced using 4 simulation points at g ¼ 1, 5, 10, and 20 m/s2. The 5% molar fraction limit velocity increases as the gravitational acceleration is increased or the molecular weight is decreased (density difference is increased) but the variation is not linear. Different formulations were tried to correlate the 5% molar fraction front velocity of the gas with density difference to the density of the air, Dr=r and the gravitational acceleration to the reference gravity ( gref ¼ 10 m/s2 in atmospheric conditions), g=gref , for the cases considered. It was
28
Elevation(cm)
24
3D 2D
observed that the 3-D and 2-D 5% molar fraction front velocities can be best represented by the following formulations that are shown to agree well with the calculated results in Figs. 16 and 17.
V3D ¼ a3D
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !0:4814 0:4960 Dr g Dr g za3D r gref r gref
i:e: V3D za3D
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dr g r gref
(9)
H2
He
20 16 12 CH4
8 4 0 0.0
0.5
1.0
1.5
2.0
Time(s) Fig. 14 – Elevation of the fuel-first-arrival-boundary as the 1% fuel molar concentration as a function of time.
Fig. 15 – Velocity of the fuel-first-arrival-boundary taken as the 1% fuel concentration as predicted by the 3-D simulation relative to the analytical-experimental correlation given in reference [18], as a function of molecular weight of the fuel.
6578
international journal of hydrogen energy 34 (2009) 6568–6579
Equ. 9
Plume velocity (cm/s)
25 3D
20 15 2D
10
Equ. 10
5 0
0
5
10
15
20
25
30
M(kg/kmol) Fig. 16 – 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of molecular weight of the fuel, g [ 10 m/s2.
V2D ¼ a2D
!0:4093 0:4316 Dr g r gref
(10)
Where, a3D ¼ 225.9 mm/s and a2D ¼ 182.1 mm/s. Thus for the cases considered, the 3-D 5% molar fraction front velocity appears to vary effectively with pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðDr=rÞ ðg=gref Þ. Fig. 18 displays the percent deviation observed in the 5% molar fraction front velocity as predicted by the 2-D compared pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to that of the 3-D simulation as a function of ðDr=rÞ ðg=gref Þ when the molecular weight of the fuel, M, and the gravitational acceleration, g, are changed. It can be observed that the deviation in the 5% molar fraction front velocity increases as g is increased or M is decreased, but the amount of variation is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi almost the same for equal changes in the ðDr=rÞ ðg=gref Þ. This was further illustrated by considering two separate cases in which both M and g were changed at the same time, as shown in Fig. 18. Again it can be seen that the deviation between the 2-D and 3-D models in the 5% molar fraction front pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi velocity is directly a function of ðDr=rÞ ðg=gref Þ and not g and M individually. Different formulations were tested to correlate the deviation observed in the 5% molar fraction front velocity as predicted by the 2-D compared to that of the 3-D simulation as a function of ðDr=rÞ ðg=gref Þ. For the cases considered Eq. (11)
Plume velocity (cm/s)
25 Equ. 9
3D
20 15 10 2D
5 0
Equ. 10 0
5
10
15
20
25
g(m/s2) Fig. 17 – 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of gravitational acceleration, M [ 16.
Fig. 18 – Deviation in the 5% molar fraction front velocity prediction by the 2-D and 3-D simulations.
fits best with the data presented in Fig. 17 which is a near linear correlation. Dr g %Dev ¼ 23:778 r gref
!1:013
Dr g z23:778 r gref
(11)
Accordingly, if 5% deviation is considered tolerable in modelling, then ðDr=rÞ ðg=gref Þ should be around 0.21. Thus for atmospheric conditions, g ¼ 10 m/s2, the evolution of the 5% molar fraction front of gases with M ¼ 23 kg/kmol or higher can be predicted with less than 5% error using the 2-D model considered in this paper. Similarly, this 2-D modelling would be sufficient to predict, the evolution of the 5% molar fraction front of hydrogen, helium and methane in effective gravity fields of g ¼ 2.3 m/s2, g ¼ 2.5 m/s2 and g ¼ 4.8 m/s2, respectively with the same 5% deviation from the 3-D predictions.
4.
Conclusions
The predicted results of the transient spatial concentration fields when a fixed mass of hydrogen is released with a negligible pressure difference into the atmosphere of an open vertical cylindrical enclosure showed that: - The temporal changes of the fuel concentration field and hence the formation and decay of the flammable mixture zones are strongly dependant on the dynamics of the fuel fingering produced. - The boundary of the first arrival of the fuel, when considered to be represented by the 1% fuel molar concentration contours, travel for the initial stages of the dispersion at approximately constant rate. The prediction of such a rate is better made with 3-D modelling. The 2-D axis-symmetric model may under-predict the value by up to 15% for methane in the cases considered and did not show the complex structure of the flow observed with the corresponding 3-D simulations particularly in the mixing zone behind the fuel-first-arrival-boundary. - The 3-D simulations displayed more complex velocity and concentration fields particularly in the mixing zone
international journal of hydrogen energy 34 (2009) 6568–6579
behind the fuel-first-arrival-boundary compared to the 2D ones - The lean flammable boundary (5% fuel molar front for hydrogen in air) also appeared to travel in the vertical direction approximately with a constant velocity which was around 0.20 m/s and 0.16 m/s for hydrogen and methane spreading into air, respectively, at ambient conditions. This velocity was lower than what was observed for the first-arrival-boundary by 5–10% for the cases considered. - The velocity of the 5% molar fraction boundary was correlated and seen to be best presented by: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V3D ¼ 225:9 ðDr=rÞ ðg=gref Þ mm/s. - The deviation in the 5% molar fraction front velocity prediction by the 2-D versus 3-D modelling increases as the gravitational acceleration is increased or the molecular weight is decreased, but the amount of variation is the same for equal changes in ðDr=rÞ ðg=gref Þ. This deviation was correlated according to: Dr g %Dev ¼ 23:78 r gref
!
- The evolution of the 5% molar fraction front of gases with M ¼ 23 kg/kmol or higher at atmospheric conditions can be predicted with less than 5% deviation using the 2-D model. - Helium propagation rates were found to be sufficiently close to those of hydrogen for it to be used adequately instead of hydrogen in experimentation.
Acknowledgements The financial assistance of Natural Sciences and Engineering Research Council of Canada (NSERC) and Canadian Natural Resources Limited (CNRL), are greatly acknowledged.
references
[1] Karim GA. Hydrogen as a spark ignition engine fuel. Int J Hydrogen Energy 2003;28(5):569–77. [2] MacIntyrea I, Tchouvelevb AV, Hayc DR, Wongd J, Grante J, Benardf P. Canadian hydrogen safety program. Int J Hydrogen Energy 2007;32:2134–43. [3] Karim GA, Tsang PW. Flame propagation through atmospheres involving concentration gradients formed by mass transfer phenomena. J Fluid Eng 1975:615–7. [4] Meroney RN. Validation of fluid modeling techniques for assessing hazards of dense gas cloud dispersion. J Hazard Mater 1987;13:377–417.
6579
[5] Ilegbusi OJ, Mat MD. A comparison of prediction and measurement of kinematical mixing of two fluids. Appl Math Model 2000;24:199–233. [6] Alhajraf S, Al-Awadhi L, Al-Fadala S, Al-Khubaizi A, Khan A, Baby S. Real-time response system for the prediction of the atmospheric transport of hazardous materials. J Loss Prev Process Indust 2005;18:520–7. [7] Bunama R, Karim GA. Investigation of the formation and venting of flammable mixtures formed within liquid fuel vessels. J Eng Gas Turbines Power 1999;121:68–73. [8] Schefer RW, Houf WG, Williams TC, Bourne B, Colton J. Characterization of high pressure, under-expanded hydrogen-jet flames. Int J Hydrogen Energy 2007;32(12): 2081–93. [9] Proust CH, Lacome JM, Jamois D, Perrette L. Process of the formation of large unconfined clouds following a massive spillage of liquid hydrogen on the ground. Spain. In: International conference on hydrogen safety; 2007. [10] Wilkening H, Baraldi DCFD. Modeling of accidental hydrogen release from pipelines. Int J Hydrogen Energy 2007;32: 2206–13. [11] Kikukawa S. Consequence analysis and safety verification of hydrogen fueling stations using CFD simulation. Int J Hydrogen Energy 2008;33:1425–34. [12] Liu H, Schreiber W. The effect of ventilation system design on hydrogen dispersion in a sedan. Int J Hydrogen Energy 2008;33:5113–9. [13] Venetsanos AG, Baraldi D, Adams P, Heggem PS, Wilkening HCFD. Modeling of hydrogen release, dispersion and combustion for automotive scenarios. J Loss Prev Process Indust 2008;21:162–84. [14] Cisse P, Karim GA. The rapid formation and dispersion of flammable zones within cylindrical vertical enclosures following the release of a fixed mass of hydrogen and other gaseous fuels into air. Int J Hydrogen Energy 2007;32:630–6. [15] Swain MR, Filoso P, Grilliot ES, Swain MN. Hydrogen leakage into simple geometric enclosures. Int J Hydrogen Energy 2003;28:229–48. [16] Incropera FP, DeWitt DP, Bergman TL, Lavine AS. Fundamentals of heat and mass transfer. New York, NY: John Wiley; 2006. [17] Young YN, Tufo H, Dubey A, Rosner R. On the miscible Rayleigh–Taylor instability. J Fluid Mech 2001;447:377–408. [18] Goncharov VN. Analytical model of nonlinear, single-mode, classical Rayleigh–Taylor instability at arbitrary Atwood numbers. Phys Rev Lett 2002;88:134502/1–134502/4. [19] Fardisi S, Karim GA. Dynamics of the flammable plumes resulting from the convective dispersion of a fixed mass of the buoyant gaseous fuel, methane, into air. J Hazard Mater 2009;169:852–9. [20] Karim GA, Wierzba I, Boon S. Some considerations of the lean flammability limits of mixtures involving hydrogen. Int J Hydrogen Energy 1985;10:117–23. [21] Chen Y, Falconer RA. Advection–diffusion modeling using the modified QUICK scheme. Int J Numer Meth Fluid 1992;15: 1171–96. [22] Barton I. Comparison of SIMPLE and PISO-type algorithms for transient flows. Int J Numer Meth Fluid 1998;26:459–83. [23] Gebhart B. Instability, transition, and turbulence in buoyancy-induced flows. Annu Rev Fluid Mech 1973;5: 213–46.
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
Characteristics of flammable, buoyant hydrogen plumes rising from open vertical containers S. Fardisi, Ghazi A. Karim* University of Calgary, Mechanical and Manufacturing Engineering Department, 2500 University Drive, NW, Calgary, Alberta, Canada T2N 1N4
article info
abstract
Article history:
The dynamics of the dispersion of a fixed mass of the highly buoyant hydrogen when
Received 24 February 2009
exposed to overlaying atmosphere with a negligible pressure difference from open vertical
Received in revised form
cylindrical enclosures are examined. Features of the rapid formation and dispersion of
6 May 2009
flammable mixtures both inside and immediate outside of the enclosure and their corre-
Accepted 10 May 2009
sponding propagation rates were examined using a 3-D CFD model. For the cases consid-
Available online 27 June 2009
ered, the puffs of the fuel–air mixture appear to produce lean flammable boundaries that move mainly at a near constant rate for much of the time. A similar simulation that used
Keywords:
an axis-symmetrical 2-D model tended to under-predict the dynamics of the lean and rich
Hydrogen safety
mixture boundaries. Hydrogen plume characteristics were compared with that of the less
Flammable mixtures
buoyant methane and helium release. Unlike methane, helium propagation rate was found
Convective dispersion
fairly close to that of hydrogen.
Buoyancy
ª 2009 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
The special thermo-chemical characteristics of hydrogen make it a very unique safety challenge. Due to its low density, low viscosity and small molecule, hydrogen is a strong leaking, fast moving gas that can be ignited much easier than most other conventional fuels. The ignition energy of the stoichiometric hydrogen/oxygen mixture at ambient conditions is about 20 mJ which is around 1/10th of that for most other conventional fuels such as methane, 0.29 mJ, and propane, 0.26 mJ [1]. In addition, the flammability limits of hydrogen in air at atmospheric conditions, 5.0–75.0% by volume, are much wider than for most other conventional fuels under the same conditions. Consequently, hydrogen can rapidly produce dangerously large flammable zones following its release from fuel installations. Accordingly, any operation that involves hydrogen has the potential to represent significant fire and
explosion hazards. Moreover it is projected that the use of hydrogen will be increased as progress is made towards the ‘‘hydrogen economy’’ [1]. Such changes are feasible only if hydrogen’s safety issues are thoroughly investigated and made comparable with those of other common fuels. Among these issues is the improvement to our understanding of the mechanism and processes of hydrogen dispersion and developing effective predictive models which should help in developing better guidelines to detect and mitigate hydrogen hazards and facilitate hydrogen acceptance by the public [2]. The isothermal and steady state dispersion of a fluid has been widely investigated, both experimentally and numerically. Karim and Tsang [3] examined the temporally changing concentration gradients of flammable mixtures of methane and air by monitoring experimentally the consequent flame propagation patterns following ignition. Meroney [4] investigated the validity of CFD techniques for modelling the
* Corresponding author. Fax: þ1 403 282 8406. E-mail address: [email protected] (G.A. Karim). 0360-3199/$ – see front matter ª 2009 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2009.05.123
international journal of hydrogen energy 34 (2009) 6568–6579
Nomenclature Cr D d g Gr h H L M R P Re Ra Sc t u0
RMS Courant number diffusion coefficient diameter gravitational acceleration Grashof number penetration height height of the fuel in the cylinder cylinder length molecular weight cylinder radius pressure reference Reynolds number Rayleigh number Schmitt number time reference velocity
dispersion of a toxic cloud. He found that the modeled clouds are qualifiedly similar in appearance to those experimentally observed and travel with similar rates. Ilegbusi and Mat [5] simulated the mixing of two fluids and compared the accuracy of different numerical approaches for modelling buoyancy driven flows. They also considered the perturbations introduced by a valve opening and concluded that when the Grashof number is relatively large, this effect is not negligible. Alhajraf et al. [6] explored the effects of local atmospheric conditions on the dispersion of hazardous chemicals from a single source. They indicated that the weather stability condition can play an important role on the dispersion processes depending on the stability level of the atmosphere. Bunama and Karim [7] examined the formation of flammable atmospheres within vertical enclosures containing a vaporizing liquid fuel and noted on different factors such as fuel type, the exposed liquid surface area and the size of the vent area. Many of the studies on hydrogen dispersion and safety have dealt with hydrogen flame characteristics [8] and hydrogen diffusion for different large-scale leakage scenarios [9]. Wilkening and Baraldi [10] compared the hydrogen release from a pipeline with that of methane. They concluded that the lower density of hydrogen is partially compensated by its higher heating value with the overall energy release is greater for hydrogen. Thus the consequences of a possible hydrogen explosion might be worse compared to that of methane. Kikukawa [11] modeled hydrogen leakage from fueling stations and visualized the effects of barriers on hydrogen diffusion. Liu and Schreiber [12] looked into hydrogen leakage in a hydrogen fueled vehicle and concluded that the ventilation system can be modified to greatly reduce the combustible hydrogen in the car interior. Venetsanos et al. [13] considered hydrogen release from an automobile in open atmosphere and in a tunnel. They concluded that a similar release of hydrogen is more severe in a tunnel. Most researchers have considered continuous and largescale hydrogen releases. Very few studies were concerned with minor, transient leakage scenarios which are often
6569
! V velocity vector w mass fraction a A constant m static viscosity y dynamic viscosity r density %Dev percent deviation CFD computational fluid dynamics RTI Rayleigh Taylor instability SIP strongly implicit procedure Subscripts f fuel a air Max maximum N infinity 3-D three-dimensional 2-D two-dimensional
technically impossible to eliminate from pipes, valves, pinholes or joints. In addition the characteristics of minor and major hydrogen releases do not appear to be similar. The previous work by Cisse and Karim [14] considered the formation and dispersion of a hydrogen cloud using a 2-D axis-symmetrical model. They compared the resulting hydrogen plume with that of the nearly buoyancy neutral ethylene, and the heavier than air propane, and noted the overwhelming importance of the fuel density relative to that of air on the dynamics of the resulting plumes. The current contribution examines the transient dispersion of a fixed mass of gaseous hydrogen which is suddenly permitted to spread with a negligible pressure difference into air at atmospheric conditions within vertical cylindrical enclosures that are open at the top to the atmosphere using the 3-D Computational Fluid Dynamics (CFD) simulation. Calculated results appear to be in good agreement with a corresponding set of experimental results and those using other numerical simulations. CFD modelling was chosen since experimenting with hydrogen is not easy due to its rapid mixing and ignition properties and the associated hazards especially within closed enclosures. Also the poor reproducibility of the tests and high costs make experimentation less practical. The associated complex temporal changes of the spatial concentrations fields and the flammable zones formed are described and compared to those of other buoyant gases especially methane and helium. Comparing the behavior of helium with that of hydrogen is important since helium is often employed to simulate the behavior of hydrogen [15] to avoid experimental difficulties and associated hazards. It is also useful to compare hydrogen with methane which is a well developed and widely used fuel resource with known safety problems. Since, there is a big difference in terms of the computational complexity and associated costs between two- and three-dimensional simulation models, the present contribution also examines whether the results of a corresponding 2-D simulation would be sufficiently accurate in comparison with
6570
international journal of hydrogen energy 34 (2009) 6568–6579
the results of the more complex and time consuming 3-D simulation for modelling hydrogen dispersion.
in which Dr ¼ ra rf, u0 is the reference velocity defined by: u20 ¼ gDr=ra L [16] and Re and Gr are the reference Reynolds and Grashof numbers, respectively.
2.
Theory and calculations
2.2.
2.1.
Problem description
A fluid with a low density cannot stably rest below a fluid with a larger density in a uniform gravitational field. This unstable situation results in the production of a flow field when the plane interface is perturbed. The mixing process is mainly driven by natural convection with molecular diffusion playing a smaller role. The resulting concentration field depends strongly on both the radial and axial directions and the bulk batches of fuel rising and forming the plume. The structure of the flow under such a condition is dominated by a special type of instabilities known as Rayleigh–Taylor Instability (RTI), described in detail in the relevant literature, e.g. [17 and 18]. Fig. 2 shows the modeled stages of the 3-D hydrogen fingering into air as a result of RTI instabilities with Gr0 ¼ 3.6 106, Re0 ¼ 1500. Some wrinkles first appear on the interface between the fuel and air. These instabilities grow and may neck at some point to release a batch of fuel to the overlaying atmosphere. Then a plume begins to form which later develops into a ‘puff’ of mixture rich in fuel that heads to the top of the cylinder. Fig. 3 depicts similar results for higher Reynolds number of Gr0 ¼ 3.6 106, Re0 ¼ 6000 with the larger diameter of d ¼ 40 mm. Obviously much more complex velocity and concentration fields will develop while similar stages of fingering process are observed. Whenever there is a continuous leakage or release of a fuel, it takes a relatively long time to have a stable concentration field established in the compartment. This does not mean that
As shown schematically in Fig. 1, the lower section of the open to atmosphere circular cylindrical vertical vessel of diameter d and length L is assumed to initially contain a fixed mass of hydrogen while the remaining upper part is filled with air. It is assumed that at a certain instant of time, assigned time zero, hydrogen is exposed to the air to commence spreading and dispersing into the overlaying atmosphere under isothermal conditions. Transient velocity and concentration fields begin to develop through the coupled transport processes of mass and momentum. Mass is transferred by the combined effects of molecular diffusion and natural convection. Molecular diffusion is driven by the local concentration gradients while, for lighter than air gases such as hydrogen (M ¼ 2 kg/kmol), helium (M ¼ 4 kg/kmol) or methane (M ¼ 16 kg/kmol) the convective mass transfer which is the bulk flow of the moving gaseous mixture due to the buoyancy effects is dominant. The variation of the concentration field of the different species within the cylinder at any time instance and location is to be determined so as to establish the development of the flammable regions. These are represented by the mixture portions where the fuel concentration lies between those corresponding to the values of the lean and rich limits, under the prevailing conditions. Particular interest is to be focused on the early stages of the evolution of the resulting plume. The key dimensionless numbers involved for this configuration were defined as follows: Re ¼
Gr ¼
ra u0 D ma gDr L3 ra y2a
Fig. 1 – Schematic representation of the set up.
Mixing stages
(1)
(2)
Fig. 2 – 3-D Transient development of the lean limit concentration surface, Gr0 [ 3.6 3 106, Re0 [ 1500, d [ 10 mm, L [ 50 mm, and H [ 20 mm.
international journal of hydrogen energy 34 (2009) 6568–6579
6571
fraction contours, was shown by Yang et al. [17] to grow quadratically with time when assuming a self-similar fluid flow and constant density difference across the advancing fronts. Thus, the penetration height, h, is given by: h¼a
Fig. 3 – 3-D Transient development of the lean surface, Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
there is no fire hazard during the transition since each puff of fuel produces a transient flammable zone. This phenomenon is extremely important when considering safety issues of buoyant fuel releases. Fig. 4 demonstrates the concentration fields at a particular moment for different diameters. Clearly, the resolved concentration fields become much more complex as the diameter increases due to different interactions between the developing fingerings. For example a rising finger might deviate from normal direction, become unstable and collapse onto another finger. Two fingers may unite or, a single finger might bifurcate into several smaller sub-fingers. These complexities are associated with multiple simultaneous plumes and their interactions, [19]. As it can be seen for methane fingering process in Fig. 4, the mean penetration of the fingers, unlike the fingering structure, seems to be independent of the diameter of the container. To reduce the complexities associated with multiple plumes, typical representative results are to be shown for the case of a vertical cylinder of 40 mm diameter where the free evolution of the plume is not affected significantly by the side walls and being confined to a single fuel finger plume. The fuel penetration height, which is defined as the maximum distance between the 0.01 and 0.99 fuel molar
ra rf 2 gt ra þ rf
(3)
where ‘‘a’’ is a constant that was shown experimentally to be approximately 0.06 and independent of the density ratios [17] and g is the gravitational acceleration. Since the density difference is only constant between immiscible fluids, even though the flow may not reach a selfsimilar condition, the quadratic growth relationship is just used for RTI with immiscible fluid. Miscible gaseous fluids involve mixing in the molecular scale that renders the buoyancy force gradually weaker. After some relatively long time, the concentration gradients become so small that the buoyancy force can no longer establish a velocity field. From this stage on, it is only molecular diffusion that controls the mixing processes and consequently the penetration length varies approximately with the square root of time. For miscible fluids such as gases, the penetration height grows approximately quadratically with time at the early stages when the diffusion time scale is much smaller than the momentum time scale and diffusion effects are negligible compared to buoyancy effects. After the lapse of some relatively long period of time buoyancy effects almost vanish and the problem becomes essentially diffusion driven with the penetration height tending to grow with the square root of time. Within the transition stage the behaviour is seen to be changing gradually from a quadratic to a square root dependence on time. Experimental and numerical studies suggest that the miscible penetration rate approaches a limiting velocity which remains relatively constant during most of the transition stage [18]. Thus, fingers of the fuel rich mixture travel at an essentially constant speed leaving a chaotic mixing zone full of vortices behind. In the present work, the evolution of the flammable mixture plumes is examined so as to improve the understanding of the development of the fire and explosion hazard in related situations. Interest is focused on the vertical propagation velocity of the lean and rich limit flammable boundaries, which are the key safety parameters, for hydrogen as compared to methane and helium. The relative accuracy of the results of the 2-D and 3-D axis-symmetric modelling was evaluated by comparing their predictions of the limit vertical velocity with the experimentally derived values.
2.3.
Modelling
The transient equations of the conservation of mass, momentum and energy are the governing equations which are essentially the transient compressible Navier–Stokes equations coupled with a species transport equation.
Fig. 4 – Methane fingering process for top diameter of one meter, twenty centimeters and four centimeters.
vr / ! þ V $r V ¼ 0 vt
(4)
! / v V ! /! ! m // ! þ V $ V V ¼ V PþmV2 V þ V V $ V þr! g r vt 3
(5)
6572
international journal of hydrogen energy 34 (2009) 6568–6579
/ / vw ! / r þ ð V $ V Þw ¼ V $ðrd V wÞ vt
(6)
The thermodynamic and transport properties are taken to change continuously all over the domain of calculation according to the transient local values of concentration and temperature [20]. The flow of the Newtonian fluid is assumed to be isothermal and non reactive for the conditions examined in this contribution. Two distinct numerical codes were employed to solve the equations to ensure that the results are sufficiently independent of the computational method. A FORTRAN code developed by Cisse [14] was used in which the governing equations were discretized using a staggered non-uniform grid. The advection terms are approximated using the QUICK scheme combined with the flux limiter ULTRA-SHARP to eliminate the numerical diffusion errors and the non-physical oscillations [21]. The enhanced SIMPLE algorithm [22] is used to couple the momentum and continuity equations. The algebraic equations resulting from the discretized equations are solved iteratively using the Strongly Implicit Procedure (SIP). In addition, Ansys-CFX Commercial software, which uses a coupled scheme and a number of modifying subroutines are used for both the 3-D and the 2-D axis-symmetrical simulation. Advection terms are approximated using a second order QUICK scheme. Also, a second order discrimination is employed to deal with the transient terms. The transition value of the Rayleigh number for natural convection in a cavity was found to be around 107–109 [23]. This number corresponds to the release of hydrogen in a vessel of 20 mm diameter at atmospheric conditions. This supports the fact that most buoyancy driven problems that involve hydrogen may be turbulent in nature. Thus, for 3-D cases a large eddy simulation (LES) with Smagorinsky subgrid scale model [10] was employed to account for turbulence effects, in parallel with a buoyancy-modified K–3 model [10] for 2-D cases that do not allow LES. A special treatment extended the calculation domain beyond the physical limits of the cylinder typically to five diameters in the radial direction and five cylinder lengths in the vertical direction outside, as depicted in Fig. 5. This was made so as to examine the
consequences of the emergence of hydrogen into the immediate vicinity of the cylinder which has important practical safety implications. A structured non-uniform grid is used for both the 2-D and 3-D simulations. In order to check the accuracy of the results, a series of mesh refinement tests were carried out. It was observed that due to the pronounced effects of the walls, the small diameter cases were stiffer to solve and needed a more robust numerical treatment. Thus, the number of iterations per time step for 3-D cases was set to 100 for small diameter cases in comparison to 50 for other larger diameter problems. This way the same convergence behaviour was established for all problems. The number of iterations per time step for 2-D cases was set to 200 for small diameters in comparison to 100 for other larger diameter problems. The convergence criterion was set to be 105 for both 3-D and 2-D simulations. In order to ensure that the time step is small enough to minimize numerical errors, time step independence tests were carried out. Simulations were performed by choosing an initial time step and then running a small period of the simulation in physical time. This initial time step was estimated [22] according to the RMS Courant number defined by: Cr ¼
Vdt dX
(7)
where, V is a characteristic velocity, dt is the time step and dX is the dimension of the grid cell. As an example for hydrogen diffusing into air, a maximum velocity of roughly 200 mm/s was expected. Accordingly, the time step was chosen to be 0.2 ms for the grid step used to ensure a Currant number below unity and a stable solution.
3.
Results
3.1. Characteristics of flammable hydrogen plume, 3-D versus 2-D simulation Cases of a finite quantity of the pure hydrogen diffusing upwards into an overlaying atmosphere of air within the
Fig. 5 – The 2-D and 3-D geometrics showing the extended control volume used in the simulation.
international journal of hydrogen energy 34 (2009) 6568–6579
cylinder of Fig. 1 were considered. The fuel fingering processes were resolved in 3-D simulations of upward motion of the gaseous fuel. In addition, 2-D axis-symmetric simulations were performed using both the FORTRAN code and the AnsysCFX code which were in very good agreement with each other. Figs. 6 and 7 show the 3-D and 2-D simulation results of the dispersion of a certain mass of hydrogen into air, respectively and illustrate the stages of mixture development from interface wrinkling till developing an effective constant vertical velocity. It can be seen that the agreement is very poor and also that the resolved 3-D concentration fields are more complex than those of the 2-D simulations. This would indicate that imposing symmetry via the 2-D simulations is unjustified in general for such a basically 3-D phenomenon when a more authentic representation of all aspects of the phenomenon is aimed at. Figs. 8 and 9 show the flammable zones for the corresponding cases of Figs. 6 and 7, which are bounded outwardly by the lean limit concentration and inwardly by the rich limit. The flammable regions grow from being initially very thin and confined to a narrow band along the interface between the hydrogen and air, to thicken later on while moving upwards. The formation and decay of the flammable mixture zones are strongly dependant on the dynamics of the fuel fingering process. It is evident that the evolution of the flammable zone is dominated by the growth and mobility of the rising fuel and the fingering structure. The flammable region which does not necessarily always extend over the whole cross-sectional area of the cylinder is of irregular shape and develops both inside and outside the cylinder. Moreover, the flammable zone can become fragmented. This would have an important implication about the characteristic spread of the fire hazard within and outside of these confined configurations. It was observed that the 2-D simulation tends to be adequately accurate in situations when the 3-D results display near symmetrical behaviour and the concentration field is not so complex. As an example, the dispersion of the fuel is effectively resolved by 2-D axis-symmetric modelling under low gravitational acceleration or for a small density difference
6573
between the fuel and air. Figs. 10 and 11 show the 3-D and 2-D simulation results of the dispersion of a certain mass of hydrogen into air in an atmosphere with a reduced gravitational acceleration of g ¼ 1.0 m/s2 with their corresponding flammable zone depicted in Figs. 12 and 13, respectively. Obviously, the flow field is much simpler compared to when g ¼ 10 m/s2. In addition, there is an acceptable level of agreement between the two sets of solutions which indicates that when the buoyancy force is weaker and the symmetry assumption is more valid, 2-D simulation gives an accurate solution of the problem. It was noticed that the rich flammability limit rapidly disappears after the dispersion, as shown in Figs. 8 and 12. This phenomenon was not observed in methane dispersion results and may be due to very large rich flammability limit, 0.75, as compared to that of methane, 0.15, and the relatively faster plume motion of hydrogen. This shows that while a methane flammable region is limited by the rich flammability limit boundary, the hydrogen flammable region extends down to the base of the vessel and covers a much larger volume. It was also observed that the life time of the flammable volume is over predicted by the 2-D model by up to 10% for both Hydrogen and Methane. This is since the 2-D model underestimates the propagation rate and it is reasonable to expect that it would overestimate the dissipation time. The rapid disappearance of the rich limit and the time required after the spill before the region is safe and noncombustible depends strongly on the initial amount of fuel to be dispersed from the vessel, the molecular weight and viscosity of the fuel and its diffusion coefficient. Further research is needed while considering a wider variety of situations to establish in more detail the behavior of the flammable region within and outside containers.
3.2. Comparison of hydrogen plume dynamics with that of methane and helium Apart from the structure of the fluid flow, it is the velocity of the fuel-first-arrival and the lean limit boundary that are
Fig. 6 – 3-D simulation results of hydrogen dispersing into air, (a) Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1 and 0.5, respectively, Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
6574
international journal of hydrogen energy 34 (2009) 6568–6579
Fig. 7 – Corresponding 2-D Simulation results of hydrogen dispersing into air, Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5 fuel molar fractions, respectively Gr0 [ 3.6 3 106, Re0 [ 6000, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
particularly important for safety analysis. The growth and spread pattern of the fuel-first-arrival-boundary taken as the 1% fuel molar fraction contours was monitored and shown in Fig. 14 for both 3-D and 2-D simulations. This boundary appears to accelerate at the beginning of its motion to reach effectively a limiting vertical velocity. The 2-D axis-symmetric model appears to underestimate this velocity by up to 22% and 17% for hydrogen and methane, respectively for the cases examined. This would indicate that the 2-D models may be of insufficient accuracy compared to the 3-D counterpart in predicting the vertical dispersion rates and patterns, for both hydrogen and methane. The observed slower mixture evolution in the 2-D simulation may be due to the known tendency
of large-scale features to grow more rapidly in a 3-D domain than in a 2-D one. Moreover, a fuel concentration front can be considered to be similar to a moving object in the flow and it is known that 2-D models tend to overestimate the drag force on submerged objects confining the RTI and forcing them to move slower. It is also interesting to note in Fig. 14 that the helium first-arrival-boundary is fairly close to that of hydrogen. This supports the view that helium may be adequately used instead of hydrogen in safety experiments to reduce the difficulties and hazards as has been emphasized by others [15]. The observed velocity of the fuel-first-arrival-boundary obtained from the 3-D simulation along with 3-D deviation is
Fig. 8 – Development of the corresponding 3-D flammable zones.
international journal of hydrogen energy 34 (2009) 6568–6579
6575
Fig. 9 – Development of the corresponding 2-D flammable zones.
plotted in Fig. 15 against the analytical-experimental results given in reference [18]. The experimental value of the limit mixture velocity appears to be sufficiently close to the 3-D constant velocity growth pattern predicted. This confirms the accuracy of the numerical solution for solving this problem. Unlike the fuel-first-arrival-boundary, i.e. the 1% fuel molar fraction front there is little experimental and numerical
data on the lean flammable boundary, i.e. 5% fuel molar front for hydrogen in air, the boundary which is more important in terms of safety issues. Accordingly, focus is made on the dynamics of this specific concentration surface. The lean flammable boundary showed a similar behaviour to that of the fuel-first-arrival-boundary (1% fuel molar front) and had a vertical velocity which was lower than the velocity of the
Fig. 10 – Corresponding 3-D simulation results of hydrogen dispersing into air, (a) Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5, respectively, Gr0 [ 1.7 3 104, Re0 [ 330, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
6576
international journal of hydrogen energy 34 (2009) 6568–6579
Fig. 11 – 2-D Simulation results of hydrogen dispersing into air, Molar fraction distribution, (1), (2) and (3) represent 0.01, 0.1, 0.5 fuel molar fractions, respectively, Gr0 [ 1.7 3 104, Re0 [ 330, d [ 40 mm, L [ 50 mm, and H [ 20 mm.
first-arrival-boundary by 5–10%. This velocity, as stated in Section 2.2 does not remain constant but after some time subsides towards zero. However, for the consideration of the early stages of the potential spread of the fire hazard, this constant value of the lean flammable front needs to be considered.
Figs. 16 and 17 display the 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of gravitational acceleration and the molecular weight of the fuel, respectively. Fig. 16 was produced using 6 simulation points at M ¼ 2(hydrogen), 4(helium), 16(methane), 8, 22 and 27. M ¼ 8, 22 and 27 was considered to have been the result of
Fig. 12 – Development of the corresponding 3-D flammable zones.
international journal of hydrogen energy 34 (2009) 6568–6579
6577
Fig. 13 – Development of the corresponding 2-D flammable zones.
homogenous mixing of proper portions of hydrogen and air and values of the corresponding fluid and flow properties such as diffusion coefficient, viscosity, etc. were obtained accordingly. Thus, the 5% molar fraction front represents the lean limit for only methane and not for other gases considered here. Fig. 17 was produced using 4 simulation points at g ¼ 1, 5, 10, and 20 m/s2. The 5% molar fraction limit velocity increases as the gravitational acceleration is increased or the molecular weight is decreased (density difference is increased) but the variation is not linear. Different formulations were tried to correlate the 5% molar fraction front velocity of the gas with density difference to the density of the air, Dr=r and the gravitational acceleration to the reference gravity ( gref ¼ 10 m/s2 in atmospheric conditions), g=gref , for the cases considered. It was
28
Elevation(cm)
24
3D 2D
observed that the 3-D and 2-D 5% molar fraction front velocities can be best represented by the following formulations that are shown to agree well with the calculated results in Figs. 16 and 17.
V3D ¼ a3D
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !0:4814 0:4960 Dr g Dr g za3D r gref r gref
i:e: V3D za3D
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dr g r gref
(9)
H2
He
20 16 12 CH4
8 4 0 0.0
0.5
1.0
1.5
2.0
Time(s) Fig. 14 – Elevation of the fuel-first-arrival-boundary as the 1% fuel molar concentration as a function of time.
Fig. 15 – Velocity of the fuel-first-arrival-boundary taken as the 1% fuel concentration as predicted by the 3-D simulation relative to the analytical-experimental correlation given in reference [18], as a function of molecular weight of the fuel.
6578
international journal of hydrogen energy 34 (2009) 6568–6579
Equ. 9
Plume velocity (cm/s)
25 3D
20 15 2D
10
Equ. 10
5 0
0
5
10
15
20
25
30
M(kg/kmol) Fig. 16 – 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of molecular weight of the fuel, g [ 10 m/s2.
V2D ¼ a2D
!0:4093 0:4316 Dr g r gref
(10)
Where, a3D ¼ 225.9 mm/s and a2D ¼ 182.1 mm/s. Thus for the cases considered, the 3-D 5% molar fraction front velocity appears to vary effectively with pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðDr=rÞ ðg=gref Þ. Fig. 18 displays the percent deviation observed in the 5% molar fraction front velocity as predicted by the 2-D compared pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to that of the 3-D simulation as a function of ðDr=rÞ ðg=gref Þ when the molecular weight of the fuel, M, and the gravitational acceleration, g, are changed. It can be observed that the deviation in the 5% molar fraction front velocity increases as g is increased or M is decreased, but the amount of variation is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi almost the same for equal changes in the ðDr=rÞ ðg=gref Þ. This was further illustrated by considering two separate cases in which both M and g were changed at the same time, as shown in Fig. 18. Again it can be seen that the deviation between the 2-D and 3-D models in the 5% molar fraction front pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi velocity is directly a function of ðDr=rÞ ðg=gref Þ and not g and M individually. Different formulations were tested to correlate the deviation observed in the 5% molar fraction front velocity as predicted by the 2-D compared to that of the 3-D simulation as a function of ðDr=rÞ ðg=gref Þ. For the cases considered Eq. (11)
Plume velocity (cm/s)
25 Equ. 9
3D
20 15 10 2D
5 0
Equ. 10 0
5
10
15
20
25
g(m/s2) Fig. 17 – 5% molar fraction front velocity as predicted by the 2-D and 3-D results as a function of gravitational acceleration, M [ 16.
Fig. 18 – Deviation in the 5% molar fraction front velocity prediction by the 2-D and 3-D simulations.
fits best with the data presented in Fig. 17 which is a near linear correlation. Dr g %Dev ¼ 23:778 r gref
!1:013
Dr g z23:778 r gref
(11)
Accordingly, if 5% deviation is considered tolerable in modelling, then ðDr=rÞ ðg=gref Þ should be around 0.21. Thus for atmospheric conditions, g ¼ 10 m/s2, the evolution of the 5% molar fraction front of gases with M ¼ 23 kg/kmol or higher can be predicted with less than 5% error using the 2-D model considered in this paper. Similarly, this 2-D modelling would be sufficient to predict, the evolution of the 5% molar fraction front of hydrogen, helium and methane in effective gravity fields of g ¼ 2.3 m/s2, g ¼ 2.5 m/s2 and g ¼ 4.8 m/s2, respectively with the same 5% deviation from the 3-D predictions.
4.
Conclusions
The predicted results of the transient spatial concentration fields when a fixed mass of hydrogen is released with a negligible pressure difference into the atmosphere of an open vertical cylindrical enclosure showed that: - The temporal changes of the fuel concentration field and hence the formation and decay of the flammable mixture zones are strongly dependant on the dynamics of the fuel fingering produced. - The boundary of the first arrival of the fuel, when considered to be represented by the 1% fuel molar concentration contours, travel for the initial stages of the dispersion at approximately constant rate. The prediction of such a rate is better made with 3-D modelling. The 2-D axis-symmetric model may under-predict the value by up to 15% for methane in the cases considered and did not show the complex structure of the flow observed with the corresponding 3-D simulations particularly in the mixing zone behind the fuel-first-arrival-boundary. - The 3-D simulations displayed more complex velocity and concentration fields particularly in the mixing zone
international journal of hydrogen energy 34 (2009) 6568–6579
behind the fuel-first-arrival-boundary compared to the 2D ones - The lean flammable boundary (5% fuel molar front for hydrogen in air) also appeared to travel in the vertical direction approximately with a constant velocity which was around 0.20 m/s and 0.16 m/s for hydrogen and methane spreading into air, respectively, at ambient conditions. This velocity was lower than what was observed for the first-arrival-boundary by 5–10% for the cases considered. - The velocity of the 5% molar fraction boundary was correlated and seen to be best presented by: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V3D ¼ 225:9 ðDr=rÞ ðg=gref Þ mm/s. - The deviation in the 5% molar fraction front velocity prediction by the 2-D versus 3-D modelling increases as the gravitational acceleration is increased or the molecular weight is decreased, but the amount of variation is the same for equal changes in ðDr=rÞ ðg=gref Þ. This deviation was correlated according to: Dr g %Dev ¼ 23:78 r gref
!
- The evolution of the 5% molar fraction front of gases with M ¼ 23 kg/kmol or higher at atmospheric conditions can be predicted with less than 5% deviation using the 2-D model. - Helium propagation rates were found to be sufficiently close to those of hydrogen for it to be used adequately instead of hydrogen in experimentation.
Acknowledgements The financial assistance of Natural Sciences and Engineering Research Council of Canada (NSERC) and Canadian Natural Resources Limited (CNRL), are greatly acknowledged.
references
[1] Karim GA. Hydrogen as a spark ignition engine fuel. Int J Hydrogen Energy 2003;28(5):569–77. [2] MacIntyrea I, Tchouvelevb AV, Hayc DR, Wongd J, Grante J, Benardf P. Canadian hydrogen safety program. Int J Hydrogen Energy 2007;32:2134–43. [3] Karim GA, Tsang PW. Flame propagation through atmospheres involving concentration gradients formed by mass transfer phenomena. J Fluid Eng 1975:615–7. [4] Meroney RN. Validation of fluid modeling techniques for assessing hazards of dense gas cloud dispersion. J Hazard Mater 1987;13:377–417.
6579
[5] Ilegbusi OJ, Mat MD. A comparison of prediction and measurement of kinematical mixing of two fluids. Appl Math Model 2000;24:199–233. [6] Alhajraf S, Al-Awadhi L, Al-Fadala S, Al-Khubaizi A, Khan A, Baby S. Real-time response system for the prediction of the atmospheric transport of hazardous materials. J Loss Prev Process Indust 2005;18:520–7. [7] Bunama R, Karim GA. Investigation of the formation and venting of flammable mixtures formed within liquid fuel vessels. J Eng Gas Turbines Power 1999;121:68–73. [8] Schefer RW, Houf WG, Williams TC, Bourne B, Colton J. Characterization of high pressure, under-expanded hydrogen-jet flames. Int J Hydrogen Energy 2007;32(12): 2081–93. [9] Proust CH, Lacome JM, Jamois D, Perrette L. Process of the formation of large unconfined clouds following a massive spillage of liquid hydrogen on the ground. Spain. In: International conference on hydrogen safety; 2007. [10] Wilkening H, Baraldi DCFD. Modeling of accidental hydrogen release from pipelines. Int J Hydrogen Energy 2007;32: 2206–13. [11] Kikukawa S. Consequence analysis and safety verification of hydrogen fueling stations using CFD simulation. Int J Hydrogen Energy 2008;33:1425–34. [12] Liu H, Schreiber W. The effect of ventilation system design on hydrogen dispersion in a sedan. Int J Hydrogen Energy 2008;33:5113–9. [13] Venetsanos AG, Baraldi D, Adams P, Heggem PS, Wilkening HCFD. Modeling of hydrogen release, dispersion and combustion for automotive scenarios. J Loss Prev Process Indust 2008;21:162–84. [14] Cisse P, Karim GA. The rapid formation and dispersion of flammable zones within cylindrical vertical enclosures following the release of a fixed mass of hydrogen and other gaseous fuels into air. Int J Hydrogen Energy 2007;32:630–6. [15] Swain MR, Filoso P, Grilliot ES, Swain MN. Hydrogen leakage into simple geometric enclosures. Int J Hydrogen Energy 2003;28:229–48. [16] Incropera FP, DeWitt DP, Bergman TL, Lavine AS. Fundamentals of heat and mass transfer. New York, NY: John Wiley; 2006. [17] Young YN, Tufo H, Dubey A, Rosner R. On the miscible Rayleigh–Taylor instability. J Fluid Mech 2001;447:377–408. [18] Goncharov VN. Analytical model of nonlinear, single-mode, classical Rayleigh–Taylor instability at arbitrary Atwood numbers. Phys Rev Lett 2002;88:134502/1–134502/4. [19] Fardisi S, Karim GA. Dynamics of the flammable plumes resulting from the convective dispersion of a fixed mass of the buoyant gaseous fuel, methane, into air. J Hazard Mater 2009;169:852–9. [20] Karim GA, Wierzba I, Boon S. Some considerations of the lean flammability limits of mixtures involving hydrogen. Int J Hydrogen Energy 1985;10:117–23. [21] Chen Y, Falconer RA. Advection–diffusion modeling using the modified QUICK scheme. Int J Numer Meth Fluid 1992;15: 1171–96. [22] Barton I. Comparison of SIMPLE and PISO-type algorithms for transient flows. Int J Numer Meth Fluid 1998;26:459–83. [23] Gebhart B. Instability, transition, and turbulence in buoyancy-induced flows. Annu Rev Fluid Mech 1973;5: 213–46.