- Email: [email protected]
Development of heavy vapour clouds in very low wind speeds
Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Development of heavy vapour clouds in very low wind speeds Graham Atkinson Health and Safety Executive, Harpur Hill, Buxton, UK
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 February 2017 Received in revised form 7 April 2017 Accepted 17 April 2017 Available online 21 April 2017
The majority of very large vapour cloud explosions occur in such low wind speeds that vapour flow is gravitydriven. The rate of stripping (or detrainment) of heavy gas by the wind is so low that almost all of the gas released remains in a gas “blanket”, which typically slumps in all directions around the source. Such gravitydriven flammable clouds may extend many hundreds of metres from the source in all directions without significant dilution: their development cannot be assessed with conventional windy dispersion models. Previous work on the basic problem of gravity-driven vapour transport on a flat open surface is reviewed. Existing methods suffer from a number of deficiencies including arbitrary specification of current height and velocity and the neglect of surface friction. A new method of analysis is presented. Close to the source the flow can be calculated from the initial conditions at the source by integration of momentum and continuity equations. Entrainment of fresh air and friction increase the stability of the flow (Richardson number) and the rate of dilution declines towards zero. When the Richardson number approaches unity the current becomes critical and further out from the source the flow is under downstream control. The new method is used to calculate the total dilution of a vapour cloud before entrainment of fresh air is completely suppressed. This determines whether clouds are flammable at long-range. The effects of surface roughness, bund height, vapour density etc. are described. Predictions are compared with data from incidents. Crown Copyright © 2017 Published by Elsevier Ltd. This is an open access article under the Open Government License (OGL) (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/ 3/).
Keywords: Vapour cloud Gas cloud Axisymmetric Gravity-driven Richardson number Vapour cloud explosion
1. Introduction A recent review of very large, high-loss vapour cloud explosions (Atkinson and Hall, 2015; Atkinson et al., 2016) showed that more than 70% of these incidents occurred in such low wind speeds that vapour flow was gravity-driven. In these cases the rate of detrainment of heavy gas by the wind was so low that almost all of the gas released remained in a gas “blanket”, which typically slumped in all directions around the source. The extent of spread is a function of ground elevation rather than the direction of any weak, overlying wind. If the site is flat, the cloud spreads out in all directions equally. Such low winds are quite rare in most parts of the world e typically 5e10%. The high frequency of recorded incidents in very low winds is a consequence of the large cloud size associated with these weather conditions and consequent high probability of ignition. Fig. 1 shows the extent of the cloud in two examples of this type
E-mail address: [email protected].
of vapour transport: in each case overfilling of gasoline tanks in very still conditions led to the generation of gravity-driven clouds that (in one case) extended up to 500 m from the source. The locations of the vapour sources are marked with yellow dots. Development of clouds like this is an important part of the risk profile at many flammable sites and may dominate the total off-site risk. In the UK for example, land-use planning decisions near gasoline storage depots are made on the basis of gravity-driven vapour flow scenarios (HSE SPC/TECH/GEN 43 and HSE, 2007). In this case the extent and concentration of fuel of the cloud are simply and directly based on experience in the Buncefield incident (Buncefield Major Incident Investigation Board, 2007).
2. Existing methods of assessment The scale and destructive power of this kind of incident means there is a pressing need for methods to assess the vapour transport that can account for site-specific factors such as vapour production rate, surface roughness and topography.
http://dx.doi.org/10.1016/j.jlp.2017.04.011 0950-4230/Crown Copyright © 2017 Published by Elsevier Ltd. This is an open access article under the Open Government License (OGL) (http://www.nationalarchives.gov.uk/ doc/open-government-licence/version/3/).
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
163
Fig. 1. Gravityedriven vapour currents: Left, San Juan (CSB, 2015). Right, Buncefield (Buncefield Major Incident Investigation Board, 2007).
A simple approach is described in FABIG Technical Note 12 (FABIG, 2013) which provides a means of calculating the source term (vapour volume production rate) for tank overfill incidents in calm conditions. Although the volume and concentration of a cloud can be predicted with some confidence using this method, the cloud depth and extent cannot. On the basis of experience at a number of incidents the FABIG TN12 method recommends that assessments for fairly flat sites should be made on the assumption that the cloud spreads symmetrically and has a constant depth of 2 m. This cannot be generally true and there are records of gravity driven clouds with depths significantly above and below this level. The process of accumulation of a gas blanket on flat ground is modelled in DEGADIS (Havens and Spicer, 1988). In principle this model can be applied to purely gravity-driven vapour transport, but it has almost always been applied to intermediate wind speeds where there is initially significant accumulation in a gas blanket which then grows to a size where detrainment balances the input of gas from the source. In DEGADIS it is assumed that the gas blanket has a uniform height H (at any given time) but this height varies as a function of time depending on the volume production rate and the progress of the front of the gravity current. The front velocity Uf is modelled by assuming the Richardson number at the front, Rif, is constant:
Rif ¼
close to the source at 300 s; this is much larger than the observed velocity e which was <1 m/s. However, the maximum extent of the cloud spread into flat areas was reasonably well predicted by DEGADIS. This is the key output when the code is used to determine the lateral size of a gas blanket that will act as a secondary source in a windy dispersion analysis. The problem of axisymmetric spreading of a gravity current has been treated theoretically using a shallow water type model (Slim and Huppert, 2011). As in DEGADIS, their analysis used a Richardson number condition at the front, but with Rif ¼ 1.2. If entrainment of ambient fluid and friction are ignored, this work showed that the gas flows outward with approximately constant speed and with a cloud depth that decreases rapidly and continuously until close to the front, where there is a jump in height and rapid reduction in speed. These inviscid, non-entraining solutions are exact, but do not correspond in a useful way to observations of the structure of vapour clouds during incidents. For example, it is predicted that, as time goes on, an ever-higher proportion of the total volume of gas released is concentrated in an area just behind the front. This is not consistent with direct observations of these clouds or with the thermal and pressure damage caused during explosions. Slim and Huppert go on to consider the effect of entrainment on cloud development, using the assumption that entrainment and
g0H ¼ 0:76 Uf2
where g0 is the relative density of the gas. g0 ¼ g,rcloudr ramb amb Applying this frontal equation together with conservation of volume of released gas allows the radius and height of the cloud to be calculated. Some typical results for a case where g' ¼ 0.5 and the gas volume production rate Q ¼ 200 m3/s are shown in Fig. 2. The predicted cloud height falls steadily with time, while the velocity close to the source is high and increases with time. These cloud densities and volume production rates roughly correspond to conditions in the Buncefield incident e for which CCTV records provide data on the development of the cloud (Gant, 2006). The area around the tank did not correspond completely to the idealised flat, open environment assumed in the model but, even allowing for this, there are some significant discrepancies between the model predictions and the observations. The model predicts a gradual decline in the cloud height with time. This is certainly not what was observed as, on the contrary, the height slowly increased. The predicted height close to the source is 0.4 m at t ¼ 300 s which is less than the observed height of >1 m. The model predicts a gas velocity of around 3 m/s
Fig. 2. Cloud height and velocity calculated by the DEGADIS method.
164
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
dilution of the current continues until a point immediately behind the front. This greatly improves the potential usefulness of the solutions e reducing the thinning of the current and accumulation of gas behind the front. However, their analysis neglects the effect of surface friction and the predicted vapour flows generally have relatively low Richardson numbers - below the limit at which significant entrainment dies away. Without including surface friction, the analysis does not predict the most important feature of a gravity current, which is the potential for gas to travel for distances of many hundreds of metres without significant dilution. CFD offers the potential to account for the details of site topography e which is particularly important in gravity-driven vapour transport. A large number of studies have been carried out involving modelling of large-scale tests in windy conditions (e.g. Gant et al., 2014). However, there appear to be no large-scale experimental data sets corresponding to continuous releases in near calm conditions when vapour transport was dominated by gravity. A few studies have been carried out (Gant, 2006; Venart and Rogers, 2011) in which CFD was used to analyse the vapour cloud development at Buncefield and other incidents in calm conditions. Comparisons in these gravity-driven cases have focussed on the progress of the vapour front because this is generally the only quantitative data that is available. Unfortunately, movement of the front is rather insensitive to the extent to which cloud structure and especially dilution are successfully modelled. Modelling of gravity-driven flows over distances of many hundreds of metres is likely to be a significant challenge. Entrainment rates should be extremely low compared with windy dispersion analyses and the effect of any numerical diffusion needs to be closely considered. The first challenge is to demonstrate that a CFD code can in fact reproduce the longrange transport of vapour without significant dilution that the incident data indicate. Additional large-scale tests, which could provide an opportunity for code validation, would be very valuable.
3. New analysis A new approach to calculation of gravity-driven dispersion in low wind/calm conditions is summarised in Fig. 3. In the discussion the local Richardson number is defined as
Ri ¼
g0H U2
air is suppressed and secondly for Ri > 1 the speed of interfacial gravity waves exceeds the flow speed and such waves can propagate upstream. Initially vapour flow is supercritical (Ri < 1): this could be flow over a bund wall or from slumping of a heavy, upwards facing jet. The development of the supercritical gravity current is controlled by the initial conditions, entrainment and friction. The method of calculation for this section of the flow simply involves forward integration of the momentum and continuity equations from the conditions established at the source. The rate of entrainment and frictional stress can be assumed to be simple functions of the mean flow velocity. The effect of entrainment and friction is to slow and deepen the gravity current and eventually Ri rises to around 1. The flow becomes critical e any further increase in Ri would allow propagation of gravity waves upstream. These disturbances would rapidly reduce friction and any residual entrainment, which would tend to reverse the increase in Ri. In this way the flow progresses under downstream control with Ri remaining close to 1. Near the front, the flow is time dependent and again Ri is of order 1. 3.1. Supercritical flow close to the source The basic mechanics of the supercritical current can be investigated by assuming constant values of velocity and concentration through the depth of the current e this is a similar approach to that taken by (for example) Slim and Huppert (2011) or Johnson and Hogg (2013). Well behind the front, the flow depth and velocity are constant (in time) which leads to equations (1) and (2) - Byron Bird et al. (2007):
Continuity
d ðUHrÞ ¼ EUr dr
where E is the ratio of entrainment velocity to flow speed
Momentum
Fig. 3. Steady flow regimes (behind the front).
d 2 d 1 U Hr þ r g0H 2 ¼ u2* r dr dr 2
(2)
u* is the shear velocity (i.e. the surface shear stress is ru2* )
Eliminating
where U and H are the speed and depth of the vapour current at a given location. The physical significance of this parameter is twofold: firstly, at high values of Ri mixing between the vapour current and overlying
(1)
dU dH gives ¼ dr dr
Hr þ k2 þ E 2 Ri 2 1 Ri
(3)
where k ¼ uU*
Also
dRi Ri Ri Ri,H Ri ¼ þ3 þE 1þ þ k2 dr r Hð1 RiÞ r 2
(4)
These equations are easily numerically integrated (e.g. with a spreadsheet) from a set of initial conditions at the edge of the source area. The dilution is calculated from the H and Ri using the conservation of gas g0UHr ¼ const and the definition of Ri (i.e. Ri ¼ g0H=U 2 ). To begin with additional air is entrained but, as the Richardson number increases, the rate of dilution decreases and eventually approaches zero. At this point the gas concentration in the wider cloud is fixed. As the value of Ri approaches unity and the equations (3) and (4) become singular. The physical interpretation of this is covered in Section 6, where it is argued that further development of the current is under downstream control, with a Richardson number close to unity. The low Ri limit of the equation for H is
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
dH H ¼ þ k2 þ 2E dr r The maximum value of k is around 0.1 and the low Ri limit of E is around 0.08. This means that (for low Ri) changes in current depth are dominated by entrainment; friction has little effect. At large r and with k2 ≪ E the solution of this equation is:
H ¼ Er
165
measurements of entrainment by Lofquist, 1960. Examination of the original data suggest that the small volume entrainment rates measured (nominally at high values of Ri) can actually be accounted for by localised entrainment near the start of the current using the simple detrainment formula proposed by Briggs et al., 1990 for a deep fluid. Once a stably stratified interface is established the entrainment rates are probably extremely low and driven only by diffusion. For this reason, the entrainment correlation shown in (5) is used in the supercritical analysis.
which is the familiar result for a conventional radial wall jet (with minimal density variations).
3.3. Evaluating the friction coefficient k ¼ uU*
3.2. Evaluating E at higher Richardson numbers
The value of k the friction coefficient can be estimated from the roughness length for the surface zrough and cloud depth H. The roughness length is about 10% of the height of obstacles on the surface.
The value of the entrainment coefficient E is a function of Ri. Various relationships have been proposed for this function. Two examples are shown in Fig. 4.
U 1 H ¼ ln u* k zrough
! where k ¼ 0:4 is von Karman0 s constant
Slim and Huppert (2011) or Turner (1986)
0:08 0:1Ri ;0 E ¼ max 1 þ 5Ri
(7) (5)
Johnson and Hogg (2013)
E¼
0:073 1 þ 27Ri
(6)
These proposals vary significantly from one another, but in both cases the rate of entrainment at Ri ¼ 1 is very low or zero and friction may dominate any reduction in the thinning of the current. It is worth noting that expressions like (6) showing E proportional to 1/Ri for large Ri are used in many dispersion codes e.g. DEGADIS, PHAST (DNV, 2013). This appears to be largely based on
Some values for H ¼ 1 m and H ¼ 2 m are shown below in Table 1 for a range of roughness lengths that would be appropriate for short (mown) grass (zrough ¼ 0.003) to moderately overgrown grass (zrough ¼ 0.03). 3.4. Variation of cloud depth and velocity in a steady critical flow (Ri~1) The element of the DEGADIS formulation that leads to discrepancies with observations is the arbitrary assumption that the cloud height is constant with radius. For circular symmetric spreading this leads to unrealistically high velocities close to the source. It is of interest to examine the consequence of different assumptions about the variation of cloud velocity and height with radius. We assume initially that only areas of the cloud quite close to the front are changing in time. This means that the quantity r,UðrÞ,HðrÞ is constant with radius through most of the cloud. The appropriateness of this assumption has to be reviewed at the end of the analysis. If UðrÞ is:
UðrÞ ¼ U0
R0 r
m (8)
Then HðrÞ must vary as:
HðrÞ ¼ H0
1m R0 r
(9)
R0 is the fixed radius at which the flow becomes critical. At large times, R0 will be small compared with the frontal radius. The variation of cloud height with radius is shown schematically in Fig. 5.
Table 1 Values of friction coefficient k ¼ u* as a function of roughness length and cloud U depth. Roughness length (m)
Fig. 4. Entrainment coefficient as a function of Ri.
0.003 0.01 0.03
Cloud depth 1m
2m
k ¼ 0.069 k ¼ 0.087 k ¼ 0.114
k ¼ 0.062 k ¼ 0.075 k ¼ 0.095
166
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Fig. 5. Schematic of the upper surface of a cloud - the depth H decreases with radius.
If the volume flow rate is Q, then the initial height H0 and velocity U0 (both at radius R0) are related as
Rifr 3m1
Q ¼ 2p,R0 ,H0 ,U0
For large times it is possible to express the cloud height and velocity in the simple, time-independent forms (1) and (2) only if
(10)
Conservation of volume implies Rmax Z
UðrÞ ¼ U0
HðrÞ,2p,r,dr ¼ Qt
(11)
R0
Substituting HðrÞ ¼ H0
Gives H0 ¼
R0 r
ðm þ 1ÞQt Rmþ1 2pR1m max 0
1m
Rmax > > R0
(12)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ramb g HðRmax ; tÞ
ramb
dRmax ðtÞ ¼ CE dt
1 3
and HðrÞ ¼ H0
R0 r
2 3
(19)
3 1 4 4 12 g0 4 1 3 ,CE , ,Q 4 ,t 4 3 2p
(20)
The expressions (19) allow the depth and velocity of the steady part of the current to be calculated from the radius and conditions when the flow becomes critical. 4. Results of analysis
ðwith CE 1Þ
(13)
Values for the empirical constant CE vary from 0.91 (Slim and Huppert) to 1.15 (Havens and Spicer)
dRmax ðtÞ ¼ CE dt
R0 r
In this case the Richardson number of the flow does not change with radius. The front edge of the current advances as
Rmax ¼
A constant Ri condition at the front is appropriate.
dRmax ðtÞ ¼ CE dt
(18)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ramb R0 1m g H0 ðtÞ ramb Rmax sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g0ðm þ 1Þ,Q ,t
Examples illustrating the form of solutions to (3) and (4) for some example inputs are shown below. The Slim/Turner form of E(Ri) is used - Equation (5). Example 1: Flow over a bund e roughly corresponding to the Buncefield case
k (friction)
(14)
2p,R2max
0.08
Initial values Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
35
209
1
0.5
0.55
This can be integrated to give to give (at large t)
Rmax ¼
rffiffiffi 1 1 g0 4 1 1 3 4 ðm þ 1Þ4 ,CE2 , ,Q 4 ,t 4 3 2p
(15)
3
The time variation of Rmax (ft 4 ) is independent of m and in line with experimental findings. This is why DEGADIS gives the correct functional dependence of the front location with time even though the assumptions about cloud depth are not correct. In general H0 and U0 vary with time if ms13
H0 ft
13m 4
U0 ft
3m1 4
(16) (17)
The analysis was based on the assumption that cloud height close to the source does not vary with time e so that r,UðrÞ,HðrÞ is independent of radius. This is only consistent for m ¼ 13 The Richardson number Ri ¼ g0H varies with radius if ms13 U2
The Richardson number increases to 1 at a range of about 115 m (Fig. 6). Overall the concentration in the current reduces to around 77% of its initial value before entrainment finishes. This flow could extend for very large distances with minimal entrainment (molecular diffusion is small and has been ignored). It is presumed that once the Richardson number reaches 1 it remains at this value and the flow varies according to (19). Example 2: Jet of LPG (vertically upwards) In still conditions a vertical jet of LPG will eventually flow back to ground level around the source. In this case the source term radius may be much less than the previous example and the initial vapour density difference (g0) higher. A wide range of conditions may occur at the point where the gravity current returns to ground level e depending on the outlet pressure, degree of re-entrainment etc. The example below uses typical values for low pressure jet, rather than analysing a particular source.
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
k (friction)
0.08
Initial values
k (friction)
Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
5
209
1.8
2
0.26
The Richardson number increases to 1 at a range of about 140 m (Fig. 7). Overall the concentration in the current reduces to around 32% of its initial value before entrainment finishes. After this the flow could extend for very large distances with minimal entrainment. In this case the cloud would remain flammable at very long range. Example 3: Upward jet of LPG e ground friction ignored This hypothetical case illustrates the importance of including frictional stress
0
167
Initial values Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
5
209
1.8
2
0.26
The Richardson number never reaches Ri ¼ 1 (Fig. 8). After 400 m the concentration in the current has reduced to around 8% of its initial value and entrainment is continuing. Equation (4) for the rate of increase of Ri with radius is reproduced below
dRi Ri Ri Ri:H Ri ¼ þ3 þE 1þ þ k2 dr r Hð1 RiÞ r 2
Fig. 6. a: Variation of gravity current height and Richardson number for Example 1. b: Variation of gravity current velocity and concentration for Example 1.
168
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Fig. 7. a: Variation of gravity current height and Richardson number for Example 2. b: Variation of gravity current velocity and concentration for Example 2.
It is clear that if k ¼ 0 and the entrainment constant E goes to zero at Ri > 0.8, then Ri cannot exceed 0.8. In this (hypothetical) case there would have to be a slight thickening of the current (a jump) somewhere further out (but behind the front) to match the frontal condition Ri~1. Similar conclusions apply whatever the details of the source; for example if ground friction is neglected for the Buncefield case (bund radius ¼ 35 m) entrainment also continues at long range. 5. Calculation the final gas concentration The most important feature of these results is the dilution of the current before a non-entraining flow is established since this
(combined with the source conditions) determines whether or not the current will be flammable over a very large area. The next sections describe how this dilution depends on site specific circumstances. 5.1. Effect of surface roughness The rougher the surface, the more rapidly the flow slows and deepens. The total extent of dilution before entrainment is completely suppressed is consequently a function of the friction coefficient. Table 2 shows the variation in dilution for the bunded (Buncefield) case at various levels of friction coefficient covering the range of surface roughness 0.003 me0.03 m (see Table 1).
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
169
Fig. 8. a: Variation of gravity current height and Richardson number for Example 3. b: Variation of gravity current velocity and concentration for Example 3.
5.2. Dependence of near-source dilution on bund or vapour barrier height
current for various bund heights (Buncefield-type release: g 0 ¼ 0.5, Q ¼ 200 m3/s, bund radius ¼ 35 m).
If heavy vapour spills over a barrier of significant height the initial speed of the current is increased and this produces additional entrainment. Table 3 shows results for final concentration in the
5.3. Effect of vapour density
Table 2 Final dilution of vapour currents. Friction coefficient k ¼ 0.05 0.075 0.1 0.125
u* U
Final concentration of current at long range (relative to initial concentration) 53% 74% 87% 95%
Lighter vapour currents (lower values of g0 ) such as LNG spills, will undergo larger dilutions before a stable current is formed. Table 4 shows results for final concentration in the current for various initial values of g0 (Q ¼ 200 m3/s, bund radius ¼ 35 m, Ris ¼ 0.55). 5.4. Effect of source (bund) radius If the bund radius is small the current depth decreases more rapidly with radius and hence the Richardson number tends to be lower and the rate of entrainment higher. Table 5 shows the effect
170
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Table 3 Effect of bund height on final dilution of vapour currents. Fall over bund wall (m)
Initial velocity (m/s)
Initial vapour current height (m)
Initial Richardson number Ris
Final dilution factor
0.65 0.9 1.6 3.6 8.1
0.8 0.95 1.26 1.9 2.85
1.19 1 0.75 0.5 0.33
0.93 0.55 0.24 0.069 0.02
0.80 0.77 0.66 0.46 0.29
Table 4 Effect of initial density on dilution. Reduced gravity g0 g0 ¼ g,rcloudr ramb
Final concentration of current at long range (relative to initial concentration)
0.1 0.25 0.5 0.85
0.46 0.65 0.77 0.83
amb
of bund radius on final dilution for flows with a similar initial Richardson number e note in this case the volume flow varies in proportion to the radius. This means that a small volume flow in a small bund dilutes more quickly than a proportionally larger release in a larger bund. 6. Mechanics of a far-field constant Ri number flow The simple time-independent forms for current depth and velocity in equation (19) are independent of the surface and roughness and any residual entrainment. This means they cannot be generally consistent with the lumped continuity and momentum equations (3) and (4) e which include entrainment and friction constants respectively. Instead such a constant Ri number flow would be controlled by the downstream boundary condition (at the front). This apparent mismatch may be explained by considering the effects of a disturbance of the flow that results in a sub-critical flow (Ri > 1) at a given location. Such increase in Ri would corresponds to an increase in layer depth or a reduction in flow speed. For a subcritical flow (Ri > 1) a time dependent change of this sort would produce a gravity wave that propagated back upstream. This wave would lead to both deepening and deceleration of the flow upstream of the original disturbance and associated variations in static pressure across the whole of the layer. Friction and entrainment depend critically on the flow speed in a relatively thin boundary layer close to the lower and upper surfaces of the current respectively. These low speed flows are particularly sensitive to adverse pressure gradients: so a backward propagating gravity wave would produce relatively large proportional reductions in flow close to boundaries and consequently reductions in the frictional stress. In this way any increase in Ri above 1 would rapidly lead to a substantial reduction in upstream frictional stress and this would
Table 5 Effect of bund radius on dilution. Bund radius (m)
Initial volume flow (m3/s)
Initial height (m)
Initial Richardson number Ris
Final dilution factor
3 5 10 20 40
18 30 60 120 240
1 1 1 1 1
0.55 0.55 0.55 0.55 0.55
0.39 0.45 0.55 0.66 0.8
have the effect of strongly reducing Ri in the area where the original disturbance occurred. The effects of temporary disturbances taking Ri above unity would therefore tend to die away with time. A disturbance that took Ri below unity in a particular location would leave a supercritical flow. Friction or entrainment would then tend to increase the Ri number towards 1 and such a disturbance would also tend to die away downstream. The lumped momentum equation (4) is appropriate where there is a relatively fixed relationship between the surface friction and the velocity in the main current. This is a reasonable approximation for low Ri flows (Turner, 1973, 1986) but is not appropriate in critical or sub-critical flows where propagation of gravity waves upstream can have a relatively large effect on boundary layer flow and surface friction. Several authors have provided helpful descriptions of the flow structure within established non-entraining two-dimensional gravity currents (e.g. Thorpe, 1993; Hallworth et al., 1996). The main flow of undiluted heavy fluid is separated from ambient air by a mixing layer extending to about 30% of the total depth so the gradient Richardson number in this layer is around 0.3. Briggs et al. (1990) use this finding to derive important estimates of the total amount of entrainment required to establish the stable flow. For an axisymmetric current it is, at first sight, hard to see how such a structure could be maintained over a long distance. Geometric factors lead to a thinning of the current and this would require an adjustment of the depth of the mixing layer (and potentially additional entrainment) unless the speed of the current simultaneously reduced to maintain stability. Only friction can drive such a reduction in speed and the shear needs to produce just the right degree of slowing of the flow whatever the nature of the surface. In fact the stable flow structure can be maintained if the Richardson number of the layer is constant and the arguments above describe how this may occur in a way that does not depend on the roughness of the surface. 7. Comparison with incident data 7.1. Buncefield, U.K. (2005) Gant reports CCTV data on the early progress of the vapour cloud across a flat, open are to the west of the source (Gant, 2006). Six minutes after over-topping the bund, vapour was observed in a thin layer up to about 200 m from the source (Fig. 9). Fig. 10 shows a comparison between CFD, observations and equation (20) with CE ¼ 0.91 and CE ¼ 1.15, corresponding to the range of uncertainty in this parameter. The correspondence is within the uncertainly associated with the data and equation (20) probably provides a reasonable estimate of cloud spread in very flat, unobstructed conditions. The observed cloud depth is not easy to judge in Fig. 9, but it appears to roughly correspond to the results in Fig. 6a, i.e. a depth of around 0.5 m at a range of 200 m. Later on, the progress of the vapour cloud at Buncefield was strongly constrained by gradients around the site and the blocking effect of large buildings. Fig. 11 shows the correspondence between the final cloud shape and the lie of the land to the south and west.
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
171
Fig. 9. Before and after arrival of the low lying vapour current in an area around 200 m from the source.
7.2. San Juan, Puerto Rico (2009) The vapour cloud (burned area) at San Juan extended north across flat ground by around 480 m in 1800 s. Equation (20) predicts a distance of travel of between 590 m and 660 m for values of CE between 0.91 and 1.15. At this stage the edge of the vapour current would be predicted to be only about 300 mm deep. Fig. 12 shows an aerial view of the burned area at the north edge of the cloud. The vapour travel was in a wetland area covered by dense vegetation, comprising plants significantly higher than the cloud. For this type of groundcover, drag would occur over the full depth of the vapour current e rather than being confined to the boundaries. The mechanism suggested in Section 6, which limits the Richardson number around unity, no longer applies and higher values of Ri could be expected. The vegetation might consequently act as a kind of porous barrier to vapour flow. This could cause increases in depth upstream and slow the progress of the front. 8. Discussion This work provides a quantitative explanation for the observation that in calm conditions, vapour clouds can travel very large distances on flat ground without showing high levels of dilution. If there is constraint on the gravity-driven flow in the far-field e.g.
Fig. 10. A comparison between observed and predicted cloud spread across a flat, open part of the Buncefield site.
Fig. 11. Correspondence between cloud geometry and topography. (Colour scale shows terrain height in metres. Dashed line is vapour cloud). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. Burned vegetation at the northern edge of the San Juan vapour cloud.
172
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
upward slopes or walls or perhaps dense high vegetation, then the flow will back-up, producing a deeper layer, again without significant dilution. Many sites will not conform to the ideal of flat, unobstructed flow analysed here. In this case, equation (20) will over-predict the likely reach of a cloud and under-predict its depth. The recommendations in FABIG Technical Note 12 may be a better general guide to cloud size in this case. CFD provides a means of allowing for site topography and this could be a useful means of assessing the area that might be affected by a cloud. This prediction is reasonably insensitive to the level of entrainment. CFD solutions that indicate rapid dilution of the vapour current below the flammable limit should be examined carefully for the effects of numerical diffusion. Disclaimer This paper and the work it describes were funded by the UK Health and Safety Executive (HSE). Its contents, including any opinions and/or conclusions expressed, are those of the author alone and do not necessarily reflect HSE policy. Nomenclature
k m Q r Rmax R0 Ri Rif
an empirical frontal constant CE ¼ Ri1/2 f entrainment coefficient (the ratio of entrainment velocity to flow speed) acceleration due to gravity (9.81 m/s2) reduced density g0 ¼ g,rcloudr ramb amb height of the gravity current height of the gravity current when the flow becomes critical friction coefficient k ¼ uU* index in power law representations of velocity or height volume flow of heavy gas radius from centre of source radius of the vapour current front radius at which the flow becomes critical Richardson number Ri ¼ g0H U2 Richardson number at the front Rif ¼ g0H U2
Ris t Uf U u* U0
Richardson number at the outer edge of the source time speed of progress of the front of the vapour current flow speed in the gravity current (a function of radius) shear velocity (the surface shear stress is ru2* ) vapour current speed when the flow becomes critical
CE E g g0 H H0
f
zrough
roughness length of the surface
ramb density of ambient air k von Karman's constant r, rcloud density of vapour current
References Atkinson, G., Hall, J., 2015. A Review of Large Vapour Cloud Incidents. HSL Report MH15/80. http://primis.phmsa.dot.gov/meetings/MtgHome.mtg?mtg¼111. Atkinson, G., Cowpe, E., Halliday, J. and Painter, D., (2016) A Review of Very Large Vapour Cloud Explosions. 19th Mary Kay O'Connor International Process Safety Symposium, College Station TX October 2016. Briggs, G.A., Thompson, R.S., Snyder, W.H., 1990. Dense gas removal from a valley by crosswinds. J. Haz. Mat. 24, 1e38. Buncefield Major Incident Investigation Board, 2007. The Buncefield Incidente11th December 2005ethe Final Report of the Major Incident Investigation Board, vol. 1, ISBN 978-07176-6270-8. http://www.hse.gov.uk/comah/buncefield/miibfinal-volume1.pdf. Byron Bird, R., Stewart, W.E., Lightfoot, E.N., 2007. Transport Phenomena, second ed. John Wiley and Sons, New York. CSB, 2015. Caribbean Petroleum Tank Terminal Explosion and Multiple Tank Fires (October 23rd 2009). US. Chemical Safety and Hazard Investigation Board. Final Investigation Report. DNV, 2013. Phast Version 7 Software. DNV Software, London, UK. Available from: http://www.dnv.com/software. Accessed 19 August 2013. FABIG, 2013. FABIG Technical Note 12-Vapour Cloud Development in Over-filling Incidents. April 2013. Available from: http://fabig.com/video-publications/ TechnicalGuidance#. Gant, S., 2006. Buncefield Investigation: Dispersion of the Vapour Cloud. HSL Report CM/06/13. Gant, S.E., Narasimhamurthy, V.D., Skjold, T., Jamois, D., Proust, C., 2014. Evaluation of multi-phase atmospheric dispersion models for application to carbon capture and storage. J. Loss Prev. Process Industries 32 (November 2014), 286e298. Hallworth, M.A., Huppert, H.E., Phillips, J.C., Sparks, R.S.J., 1996. Entrainment into two-dimensional and axisymmetric turbulent gravity currents. J. Fluid Mech. 308, 289e311. Havens, J., Spicer, T., 1988. A Dispersion Model for Elevated Dense Gas Jet Chemical Releases. EPA Report EPA-450/4-88-006a. HSE SPC/TECH/GEN 43 Land use planning advice around large scale petrol storage sites http://www.hse.gov.uk/foi/internalops/hid_circs/technical_general/spc_ tech_gen_43/. HSE 2007 http://www.hse.gov.uk/consult/condocs/cd211.htm. Johnson, C.G., Hogg, A.J., 2013. Entraining gravity currents. J.Fluid Mech. 731, p471e508. Lofquist, K., 1960. Flow and stress near an interface between stratified liquids. Phys. Fluids 3, 158e175. Slim, A.C., Huppert, H.E., 2011. Axisymmetric, constantly supplied gravity currents at high Reynolds number. J. Fluid Mech. 675, 540e551. Thorpe, S.A., 1973. Experiments on instability and turbulence in a stratified shear flow. J. Fluid Mech. 61, 731e751. Turner, J., 1973. Buoyancy Effects in Fluids. Cambridge University Press. Turner, J., 1986. Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows. J. Fluid Mech. 173, 431e447. Venart, J.E.S., Rogers, R.J., 2011. The Buncefield Explosion: Vapour Cloud Dispersion and Other Observation. Hazards XXII, IChemE Symposium Series No. 156.
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Development of heavy vapour clouds in very low wind speeds Graham Atkinson Health and Safety Executive, Harpur Hill, Buxton, UK
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 February 2017 Received in revised form 7 April 2017 Accepted 17 April 2017 Available online 21 April 2017
The majority of very large vapour cloud explosions occur in such low wind speeds that vapour flow is gravitydriven. The rate of stripping (or detrainment) of heavy gas by the wind is so low that almost all of the gas released remains in a gas “blanket”, which typically slumps in all directions around the source. Such gravitydriven flammable clouds may extend many hundreds of metres from the source in all directions without significant dilution: their development cannot be assessed with conventional windy dispersion models. Previous work on the basic problem of gravity-driven vapour transport on a flat open surface is reviewed. Existing methods suffer from a number of deficiencies including arbitrary specification of current height and velocity and the neglect of surface friction. A new method of analysis is presented. Close to the source the flow can be calculated from the initial conditions at the source by integration of momentum and continuity equations. Entrainment of fresh air and friction increase the stability of the flow (Richardson number) and the rate of dilution declines towards zero. When the Richardson number approaches unity the current becomes critical and further out from the source the flow is under downstream control. The new method is used to calculate the total dilution of a vapour cloud before entrainment of fresh air is completely suppressed. This determines whether clouds are flammable at long-range. The effects of surface roughness, bund height, vapour density etc. are described. Predictions are compared with data from incidents. Crown Copyright © 2017 Published by Elsevier Ltd. This is an open access article under the Open Government License (OGL) (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/ 3/).
Keywords: Vapour cloud Gas cloud Axisymmetric Gravity-driven Richardson number Vapour cloud explosion
1. Introduction A recent review of very large, high-loss vapour cloud explosions (Atkinson and Hall, 2015; Atkinson et al., 2016) showed that more than 70% of these incidents occurred in such low wind speeds that vapour flow was gravity-driven. In these cases the rate of detrainment of heavy gas by the wind was so low that almost all of the gas released remained in a gas “blanket”, which typically slumped in all directions around the source. The extent of spread is a function of ground elevation rather than the direction of any weak, overlying wind. If the site is flat, the cloud spreads out in all directions equally. Such low winds are quite rare in most parts of the world e typically 5e10%. The high frequency of recorded incidents in very low winds is a consequence of the large cloud size associated with these weather conditions and consequent high probability of ignition. Fig. 1 shows the extent of the cloud in two examples of this type
E-mail address: [email protected].
of vapour transport: in each case overfilling of gasoline tanks in very still conditions led to the generation of gravity-driven clouds that (in one case) extended up to 500 m from the source. The locations of the vapour sources are marked with yellow dots. Development of clouds like this is an important part of the risk profile at many flammable sites and may dominate the total off-site risk. In the UK for example, land-use planning decisions near gasoline storage depots are made on the basis of gravity-driven vapour flow scenarios (HSE SPC/TECH/GEN 43 and HSE, 2007). In this case the extent and concentration of fuel of the cloud are simply and directly based on experience in the Buncefield incident (Buncefield Major Incident Investigation Board, 2007).
2. Existing methods of assessment The scale and destructive power of this kind of incident means there is a pressing need for methods to assess the vapour transport that can account for site-specific factors such as vapour production rate, surface roughness and topography.
http://dx.doi.org/10.1016/j.jlp.2017.04.011 0950-4230/Crown Copyright © 2017 Published by Elsevier Ltd. This is an open access article under the Open Government License (OGL) (http://www.nationalarchives.gov.uk/ doc/open-government-licence/version/3/).
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
163
Fig. 1. Gravityedriven vapour currents: Left, San Juan (CSB, 2015). Right, Buncefield (Buncefield Major Incident Investigation Board, 2007).
A simple approach is described in FABIG Technical Note 12 (FABIG, 2013) which provides a means of calculating the source term (vapour volume production rate) for tank overfill incidents in calm conditions. Although the volume and concentration of a cloud can be predicted with some confidence using this method, the cloud depth and extent cannot. On the basis of experience at a number of incidents the FABIG TN12 method recommends that assessments for fairly flat sites should be made on the assumption that the cloud spreads symmetrically and has a constant depth of 2 m. This cannot be generally true and there are records of gravity driven clouds with depths significantly above and below this level. The process of accumulation of a gas blanket on flat ground is modelled in DEGADIS (Havens and Spicer, 1988). In principle this model can be applied to purely gravity-driven vapour transport, but it has almost always been applied to intermediate wind speeds where there is initially significant accumulation in a gas blanket which then grows to a size where detrainment balances the input of gas from the source. In DEGADIS it is assumed that the gas blanket has a uniform height H (at any given time) but this height varies as a function of time depending on the volume production rate and the progress of the front of the gravity current. The front velocity Uf is modelled by assuming the Richardson number at the front, Rif, is constant:
Rif ¼
close to the source at 300 s; this is much larger than the observed velocity e which was <1 m/s. However, the maximum extent of the cloud spread into flat areas was reasonably well predicted by DEGADIS. This is the key output when the code is used to determine the lateral size of a gas blanket that will act as a secondary source in a windy dispersion analysis. The problem of axisymmetric spreading of a gravity current has been treated theoretically using a shallow water type model (Slim and Huppert, 2011). As in DEGADIS, their analysis used a Richardson number condition at the front, but with Rif ¼ 1.2. If entrainment of ambient fluid and friction are ignored, this work showed that the gas flows outward with approximately constant speed and with a cloud depth that decreases rapidly and continuously until close to the front, where there is a jump in height and rapid reduction in speed. These inviscid, non-entraining solutions are exact, but do not correspond in a useful way to observations of the structure of vapour clouds during incidents. For example, it is predicted that, as time goes on, an ever-higher proportion of the total volume of gas released is concentrated in an area just behind the front. This is not consistent with direct observations of these clouds or with the thermal and pressure damage caused during explosions. Slim and Huppert go on to consider the effect of entrainment on cloud development, using the assumption that entrainment and
g0H ¼ 0:76 Uf2
where g0 is the relative density of the gas. g0 ¼ g,rcloudr ramb amb Applying this frontal equation together with conservation of volume of released gas allows the radius and height of the cloud to be calculated. Some typical results for a case where g' ¼ 0.5 and the gas volume production rate Q ¼ 200 m3/s are shown in Fig. 2. The predicted cloud height falls steadily with time, while the velocity close to the source is high and increases with time. These cloud densities and volume production rates roughly correspond to conditions in the Buncefield incident e for which CCTV records provide data on the development of the cloud (Gant, 2006). The area around the tank did not correspond completely to the idealised flat, open environment assumed in the model but, even allowing for this, there are some significant discrepancies between the model predictions and the observations. The model predicts a gradual decline in the cloud height with time. This is certainly not what was observed as, on the contrary, the height slowly increased. The predicted height close to the source is 0.4 m at t ¼ 300 s which is less than the observed height of >1 m. The model predicts a gas velocity of around 3 m/s
Fig. 2. Cloud height and velocity calculated by the DEGADIS method.
164
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
dilution of the current continues until a point immediately behind the front. This greatly improves the potential usefulness of the solutions e reducing the thinning of the current and accumulation of gas behind the front. However, their analysis neglects the effect of surface friction and the predicted vapour flows generally have relatively low Richardson numbers - below the limit at which significant entrainment dies away. Without including surface friction, the analysis does not predict the most important feature of a gravity current, which is the potential for gas to travel for distances of many hundreds of metres without significant dilution. CFD offers the potential to account for the details of site topography e which is particularly important in gravity-driven vapour transport. A large number of studies have been carried out involving modelling of large-scale tests in windy conditions (e.g. Gant et al., 2014). However, there appear to be no large-scale experimental data sets corresponding to continuous releases in near calm conditions when vapour transport was dominated by gravity. A few studies have been carried out (Gant, 2006; Venart and Rogers, 2011) in which CFD was used to analyse the vapour cloud development at Buncefield and other incidents in calm conditions. Comparisons in these gravity-driven cases have focussed on the progress of the vapour front because this is generally the only quantitative data that is available. Unfortunately, movement of the front is rather insensitive to the extent to which cloud structure and especially dilution are successfully modelled. Modelling of gravity-driven flows over distances of many hundreds of metres is likely to be a significant challenge. Entrainment rates should be extremely low compared with windy dispersion analyses and the effect of any numerical diffusion needs to be closely considered. The first challenge is to demonstrate that a CFD code can in fact reproduce the longrange transport of vapour without significant dilution that the incident data indicate. Additional large-scale tests, which could provide an opportunity for code validation, would be very valuable.
3. New analysis A new approach to calculation of gravity-driven dispersion in low wind/calm conditions is summarised in Fig. 3. In the discussion the local Richardson number is defined as
Ri ¼
g0H U2
air is suppressed and secondly for Ri > 1 the speed of interfacial gravity waves exceeds the flow speed and such waves can propagate upstream. Initially vapour flow is supercritical (Ri < 1): this could be flow over a bund wall or from slumping of a heavy, upwards facing jet. The development of the supercritical gravity current is controlled by the initial conditions, entrainment and friction. The method of calculation for this section of the flow simply involves forward integration of the momentum and continuity equations from the conditions established at the source. The rate of entrainment and frictional stress can be assumed to be simple functions of the mean flow velocity. The effect of entrainment and friction is to slow and deepen the gravity current and eventually Ri rises to around 1. The flow becomes critical e any further increase in Ri would allow propagation of gravity waves upstream. These disturbances would rapidly reduce friction and any residual entrainment, which would tend to reverse the increase in Ri. In this way the flow progresses under downstream control with Ri remaining close to 1. Near the front, the flow is time dependent and again Ri is of order 1. 3.1. Supercritical flow close to the source The basic mechanics of the supercritical current can be investigated by assuming constant values of velocity and concentration through the depth of the current e this is a similar approach to that taken by (for example) Slim and Huppert (2011) or Johnson and Hogg (2013). Well behind the front, the flow depth and velocity are constant (in time) which leads to equations (1) and (2) - Byron Bird et al. (2007):
Continuity
d ðUHrÞ ¼ EUr dr
where E is the ratio of entrainment velocity to flow speed
Momentum
Fig. 3. Steady flow regimes (behind the front).
d 2 d 1 U Hr þ r g0H 2 ¼ u2* r dr dr 2
(2)
u* is the shear velocity (i.e. the surface shear stress is ru2* )
Eliminating
where U and H are the speed and depth of the vapour current at a given location. The physical significance of this parameter is twofold: firstly, at high values of Ri mixing between the vapour current and overlying
(1)
dU dH gives ¼ dr dr
Hr þ k2 þ E 2 Ri 2 1 Ri
(3)
where k ¼ uU*
Also
dRi Ri Ri Ri,H Ri ¼ þ3 þE 1þ þ k2 dr r Hð1 RiÞ r 2
(4)
These equations are easily numerically integrated (e.g. with a spreadsheet) from a set of initial conditions at the edge of the source area. The dilution is calculated from the H and Ri using the conservation of gas g0UHr ¼ const and the definition of Ri (i.e. Ri ¼ g0H=U 2 ). To begin with additional air is entrained but, as the Richardson number increases, the rate of dilution decreases and eventually approaches zero. At this point the gas concentration in the wider cloud is fixed. As the value of Ri approaches unity and the equations (3) and (4) become singular. The physical interpretation of this is covered in Section 6, where it is argued that further development of the current is under downstream control, with a Richardson number close to unity. The low Ri limit of the equation for H is
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
dH H ¼ þ k2 þ 2E dr r The maximum value of k is around 0.1 and the low Ri limit of E is around 0.08. This means that (for low Ri) changes in current depth are dominated by entrainment; friction has little effect. At large r and with k2 ≪ E the solution of this equation is:
H ¼ Er
165
measurements of entrainment by Lofquist, 1960. Examination of the original data suggest that the small volume entrainment rates measured (nominally at high values of Ri) can actually be accounted for by localised entrainment near the start of the current using the simple detrainment formula proposed by Briggs et al., 1990 for a deep fluid. Once a stably stratified interface is established the entrainment rates are probably extremely low and driven only by diffusion. For this reason, the entrainment correlation shown in (5) is used in the supercritical analysis.
which is the familiar result for a conventional radial wall jet (with minimal density variations).
3.3. Evaluating the friction coefficient k ¼ uU*
3.2. Evaluating E at higher Richardson numbers
The value of k the friction coefficient can be estimated from the roughness length for the surface zrough and cloud depth H. The roughness length is about 10% of the height of obstacles on the surface.
The value of the entrainment coefficient E is a function of Ri. Various relationships have been proposed for this function. Two examples are shown in Fig. 4.
U 1 H ¼ ln u* k zrough
! where k ¼ 0:4 is von Karman0 s constant
Slim and Huppert (2011) or Turner (1986)
0:08 0:1Ri ;0 E ¼ max 1 þ 5Ri
(7) (5)
Johnson and Hogg (2013)
E¼
0:073 1 þ 27Ri
(6)
These proposals vary significantly from one another, but in both cases the rate of entrainment at Ri ¼ 1 is very low or zero and friction may dominate any reduction in the thinning of the current. It is worth noting that expressions like (6) showing E proportional to 1/Ri for large Ri are used in many dispersion codes e.g. DEGADIS, PHAST (DNV, 2013). This appears to be largely based on
Some values for H ¼ 1 m and H ¼ 2 m are shown below in Table 1 for a range of roughness lengths that would be appropriate for short (mown) grass (zrough ¼ 0.003) to moderately overgrown grass (zrough ¼ 0.03). 3.4. Variation of cloud depth and velocity in a steady critical flow (Ri~1) The element of the DEGADIS formulation that leads to discrepancies with observations is the arbitrary assumption that the cloud height is constant with radius. For circular symmetric spreading this leads to unrealistically high velocities close to the source. It is of interest to examine the consequence of different assumptions about the variation of cloud velocity and height with radius. We assume initially that only areas of the cloud quite close to the front are changing in time. This means that the quantity r,UðrÞ,HðrÞ is constant with radius through most of the cloud. The appropriateness of this assumption has to be reviewed at the end of the analysis. If UðrÞ is:
UðrÞ ¼ U0
R0 r
m (8)
Then HðrÞ must vary as:
HðrÞ ¼ H0
1m R0 r
(9)
R0 is the fixed radius at which the flow becomes critical. At large times, R0 will be small compared with the frontal radius. The variation of cloud height with radius is shown schematically in Fig. 5.
Table 1 Values of friction coefficient k ¼ u* as a function of roughness length and cloud U depth. Roughness length (m)
Fig. 4. Entrainment coefficient as a function of Ri.
0.003 0.01 0.03
Cloud depth 1m
2m
k ¼ 0.069 k ¼ 0.087 k ¼ 0.114
k ¼ 0.062 k ¼ 0.075 k ¼ 0.095
166
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Fig. 5. Schematic of the upper surface of a cloud - the depth H decreases with radius.
If the volume flow rate is Q, then the initial height H0 and velocity U0 (both at radius R0) are related as
Rifr 3m1
Q ¼ 2p,R0 ,H0 ,U0
For large times it is possible to express the cloud height and velocity in the simple, time-independent forms (1) and (2) only if
(10)
Conservation of volume implies Rmax Z
UðrÞ ¼ U0
HðrÞ,2p,r,dr ¼ Qt
(11)
R0
Substituting HðrÞ ¼ H0
Gives H0 ¼
R0 r
ðm þ 1ÞQt Rmþ1 2pR1m max 0
1m
Rmax > > R0
(12)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ramb g HðRmax ; tÞ
ramb
dRmax ðtÞ ¼ CE dt
1 3
and HðrÞ ¼ H0
R0 r
2 3
(19)
3 1 4 4 12 g0 4 1 3 ,CE , ,Q 4 ,t 4 3 2p
(20)
The expressions (19) allow the depth and velocity of the steady part of the current to be calculated from the radius and conditions when the flow becomes critical. 4. Results of analysis
ðwith CE 1Þ
(13)
Values for the empirical constant CE vary from 0.91 (Slim and Huppert) to 1.15 (Havens and Spicer)
dRmax ðtÞ ¼ CE dt
R0 r
In this case the Richardson number of the flow does not change with radius. The front edge of the current advances as
Rmax ¼
A constant Ri condition at the front is appropriate.
dRmax ðtÞ ¼ CE dt
(18)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ramb R0 1m g H0 ðtÞ ramb Rmax sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g0ðm þ 1Þ,Q ,t
Examples illustrating the form of solutions to (3) and (4) for some example inputs are shown below. The Slim/Turner form of E(Ri) is used - Equation (5). Example 1: Flow over a bund e roughly corresponding to the Buncefield case
k (friction)
(14)
2p,R2max
0.08
Initial values Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
35
209
1
0.5
0.55
This can be integrated to give to give (at large t)
Rmax ¼
rffiffiffi 1 1 g0 4 1 1 3 4 ðm þ 1Þ4 ,CE2 , ,Q 4 ,t 4 3 2p
(15)
3
The time variation of Rmax (ft 4 ) is independent of m and in line with experimental findings. This is why DEGADIS gives the correct functional dependence of the front location with time even though the assumptions about cloud depth are not correct. In general H0 and U0 vary with time if ms13
H0 ft
13m 4
U0 ft
3m1 4
(16) (17)
The analysis was based on the assumption that cloud height close to the source does not vary with time e so that r,UðrÞ,HðrÞ is independent of radius. This is only consistent for m ¼ 13 The Richardson number Ri ¼ g0H varies with radius if ms13 U2
The Richardson number increases to 1 at a range of about 115 m (Fig. 6). Overall the concentration in the current reduces to around 77% of its initial value before entrainment finishes. This flow could extend for very large distances with minimal entrainment (molecular diffusion is small and has been ignored). It is presumed that once the Richardson number reaches 1 it remains at this value and the flow varies according to (19). Example 2: Jet of LPG (vertically upwards) In still conditions a vertical jet of LPG will eventually flow back to ground level around the source. In this case the source term radius may be much less than the previous example and the initial vapour density difference (g0) higher. A wide range of conditions may occur at the point where the gravity current returns to ground level e depending on the outlet pressure, degree of re-entrainment etc. The example below uses typical values for low pressure jet, rather than analysing a particular source.
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
k (friction)
0.08
Initial values
k (friction)
Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
5
209
1.8
2
0.26
The Richardson number increases to 1 at a range of about 140 m (Fig. 7). Overall the concentration in the current reduces to around 32% of its initial value before entrainment finishes. After this the flow could extend for very large distances with minimal entrainment. In this case the cloud would remain flammable at very long range. Example 3: Upward jet of LPG e ground friction ignored This hypothetical case illustrates the importance of including frictional stress
0
167
Initial values Radius (m)
Vol. flow (m3/s)
Depth (m)
g0 (m/s2)
Ris
5
209
1.8
2
0.26
The Richardson number never reaches Ri ¼ 1 (Fig. 8). After 400 m the concentration in the current has reduced to around 8% of its initial value and entrainment is continuing. Equation (4) for the rate of increase of Ri with radius is reproduced below
dRi Ri Ri Ri:H Ri ¼ þ3 þE 1þ þ k2 dr r Hð1 RiÞ r 2
Fig. 6. a: Variation of gravity current height and Richardson number for Example 1. b: Variation of gravity current velocity and concentration for Example 1.
168
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Fig. 7. a: Variation of gravity current height and Richardson number for Example 2. b: Variation of gravity current velocity and concentration for Example 2.
It is clear that if k ¼ 0 and the entrainment constant E goes to zero at Ri > 0.8, then Ri cannot exceed 0.8. In this (hypothetical) case there would have to be a slight thickening of the current (a jump) somewhere further out (but behind the front) to match the frontal condition Ri~1. Similar conclusions apply whatever the details of the source; for example if ground friction is neglected for the Buncefield case (bund radius ¼ 35 m) entrainment also continues at long range. 5. Calculation the final gas concentration The most important feature of these results is the dilution of the current before a non-entraining flow is established since this
(combined with the source conditions) determines whether or not the current will be flammable over a very large area. The next sections describe how this dilution depends on site specific circumstances. 5.1. Effect of surface roughness The rougher the surface, the more rapidly the flow slows and deepens. The total extent of dilution before entrainment is completely suppressed is consequently a function of the friction coefficient. Table 2 shows the variation in dilution for the bunded (Buncefield) case at various levels of friction coefficient covering the range of surface roughness 0.003 me0.03 m (see Table 1).
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
169
Fig. 8. a: Variation of gravity current height and Richardson number for Example 3. b: Variation of gravity current velocity and concentration for Example 3.
5.2. Dependence of near-source dilution on bund or vapour barrier height
current for various bund heights (Buncefield-type release: g 0 ¼ 0.5, Q ¼ 200 m3/s, bund radius ¼ 35 m).
If heavy vapour spills over a barrier of significant height the initial speed of the current is increased and this produces additional entrainment. Table 3 shows results for final concentration in the
5.3. Effect of vapour density
Table 2 Final dilution of vapour currents. Friction coefficient k ¼ 0.05 0.075 0.1 0.125
u* U
Final concentration of current at long range (relative to initial concentration) 53% 74% 87% 95%
Lighter vapour currents (lower values of g0 ) such as LNG spills, will undergo larger dilutions before a stable current is formed. Table 4 shows results for final concentration in the current for various initial values of g0 (Q ¼ 200 m3/s, bund radius ¼ 35 m, Ris ¼ 0.55). 5.4. Effect of source (bund) radius If the bund radius is small the current depth decreases more rapidly with radius and hence the Richardson number tends to be lower and the rate of entrainment higher. Table 5 shows the effect
170
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
Table 3 Effect of bund height on final dilution of vapour currents. Fall over bund wall (m)
Initial velocity (m/s)
Initial vapour current height (m)
Initial Richardson number Ris
Final dilution factor
0.65 0.9 1.6 3.6 8.1
0.8 0.95 1.26 1.9 2.85
1.19 1 0.75 0.5 0.33
0.93 0.55 0.24 0.069 0.02
0.80 0.77 0.66 0.46 0.29
Table 4 Effect of initial density on dilution. Reduced gravity g0 g0 ¼ g,rcloudr ramb
Final concentration of current at long range (relative to initial concentration)
0.1 0.25 0.5 0.85
0.46 0.65 0.77 0.83
amb
of bund radius on final dilution for flows with a similar initial Richardson number e note in this case the volume flow varies in proportion to the radius. This means that a small volume flow in a small bund dilutes more quickly than a proportionally larger release in a larger bund. 6. Mechanics of a far-field constant Ri number flow The simple time-independent forms for current depth and velocity in equation (19) are independent of the surface and roughness and any residual entrainment. This means they cannot be generally consistent with the lumped continuity and momentum equations (3) and (4) e which include entrainment and friction constants respectively. Instead such a constant Ri number flow would be controlled by the downstream boundary condition (at the front). This apparent mismatch may be explained by considering the effects of a disturbance of the flow that results in a sub-critical flow (Ri > 1) at a given location. Such increase in Ri would corresponds to an increase in layer depth or a reduction in flow speed. For a subcritical flow (Ri > 1) a time dependent change of this sort would produce a gravity wave that propagated back upstream. This wave would lead to both deepening and deceleration of the flow upstream of the original disturbance and associated variations in static pressure across the whole of the layer. Friction and entrainment depend critically on the flow speed in a relatively thin boundary layer close to the lower and upper surfaces of the current respectively. These low speed flows are particularly sensitive to adverse pressure gradients: so a backward propagating gravity wave would produce relatively large proportional reductions in flow close to boundaries and consequently reductions in the frictional stress. In this way any increase in Ri above 1 would rapidly lead to a substantial reduction in upstream frictional stress and this would
Table 5 Effect of bund radius on dilution. Bund radius (m)
Initial volume flow (m3/s)
Initial height (m)
Initial Richardson number Ris
Final dilution factor
3 5 10 20 40
18 30 60 120 240
1 1 1 1 1
0.55 0.55 0.55 0.55 0.55
0.39 0.45 0.55 0.66 0.8
have the effect of strongly reducing Ri in the area where the original disturbance occurred. The effects of temporary disturbances taking Ri above unity would therefore tend to die away with time. A disturbance that took Ri below unity in a particular location would leave a supercritical flow. Friction or entrainment would then tend to increase the Ri number towards 1 and such a disturbance would also tend to die away downstream. The lumped momentum equation (4) is appropriate where there is a relatively fixed relationship between the surface friction and the velocity in the main current. This is a reasonable approximation for low Ri flows (Turner, 1973, 1986) but is not appropriate in critical or sub-critical flows where propagation of gravity waves upstream can have a relatively large effect on boundary layer flow and surface friction. Several authors have provided helpful descriptions of the flow structure within established non-entraining two-dimensional gravity currents (e.g. Thorpe, 1993; Hallworth et al., 1996). The main flow of undiluted heavy fluid is separated from ambient air by a mixing layer extending to about 30% of the total depth so the gradient Richardson number in this layer is around 0.3. Briggs et al. (1990) use this finding to derive important estimates of the total amount of entrainment required to establish the stable flow. For an axisymmetric current it is, at first sight, hard to see how such a structure could be maintained over a long distance. Geometric factors lead to a thinning of the current and this would require an adjustment of the depth of the mixing layer (and potentially additional entrainment) unless the speed of the current simultaneously reduced to maintain stability. Only friction can drive such a reduction in speed and the shear needs to produce just the right degree of slowing of the flow whatever the nature of the surface. In fact the stable flow structure can be maintained if the Richardson number of the layer is constant and the arguments above describe how this may occur in a way that does not depend on the roughness of the surface. 7. Comparison with incident data 7.1. Buncefield, U.K. (2005) Gant reports CCTV data on the early progress of the vapour cloud across a flat, open are to the west of the source (Gant, 2006). Six minutes after over-topping the bund, vapour was observed in a thin layer up to about 200 m from the source (Fig. 9). Fig. 10 shows a comparison between CFD, observations and equation (20) with CE ¼ 0.91 and CE ¼ 1.15, corresponding to the range of uncertainty in this parameter. The correspondence is within the uncertainly associated with the data and equation (20) probably provides a reasonable estimate of cloud spread in very flat, unobstructed conditions. The observed cloud depth is not easy to judge in Fig. 9, but it appears to roughly correspond to the results in Fig. 6a, i.e. a depth of around 0.5 m at a range of 200 m. Later on, the progress of the vapour cloud at Buncefield was strongly constrained by gradients around the site and the blocking effect of large buildings. Fig. 11 shows the correspondence between the final cloud shape and the lie of the land to the south and west.
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
171
Fig. 9. Before and after arrival of the low lying vapour current in an area around 200 m from the source.
7.2. San Juan, Puerto Rico (2009) The vapour cloud (burned area) at San Juan extended north across flat ground by around 480 m in 1800 s. Equation (20) predicts a distance of travel of between 590 m and 660 m for values of CE between 0.91 and 1.15. At this stage the edge of the vapour current would be predicted to be only about 300 mm deep. Fig. 12 shows an aerial view of the burned area at the north edge of the cloud. The vapour travel was in a wetland area covered by dense vegetation, comprising plants significantly higher than the cloud. For this type of groundcover, drag would occur over the full depth of the vapour current e rather than being confined to the boundaries. The mechanism suggested in Section 6, which limits the Richardson number around unity, no longer applies and higher values of Ri could be expected. The vegetation might consequently act as a kind of porous barrier to vapour flow. This could cause increases in depth upstream and slow the progress of the front. 8. Discussion This work provides a quantitative explanation for the observation that in calm conditions, vapour clouds can travel very large distances on flat ground without showing high levels of dilution. If there is constraint on the gravity-driven flow in the far-field e.g.
Fig. 10. A comparison between observed and predicted cloud spread across a flat, open part of the Buncefield site.
Fig. 11. Correspondence between cloud geometry and topography. (Colour scale shows terrain height in metres. Dashed line is vapour cloud). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. Burned vegetation at the northern edge of the San Juan vapour cloud.
172
G. Atkinson / Journal of Loss Prevention in the Process Industries 48 (2017) 162e172
upward slopes or walls or perhaps dense high vegetation, then the flow will back-up, producing a deeper layer, again without significant dilution. Many sites will not conform to the ideal of flat, unobstructed flow analysed here. In this case, equation (20) will over-predict the likely reach of a cloud and under-predict its depth. The recommendations in FABIG Technical Note 12 may be a better general guide to cloud size in this case. CFD provides a means of allowing for site topography and this could be a useful means of assessing the area that might be affected by a cloud. This prediction is reasonably insensitive to the level of entrainment. CFD solutions that indicate rapid dilution of the vapour current below the flammable limit should be examined carefully for the effects of numerical diffusion. Disclaimer This paper and the work it describes were funded by the UK Health and Safety Executive (HSE). Its contents, including any opinions and/or conclusions expressed, are those of the author alone and do not necessarily reflect HSE policy. Nomenclature
k m Q r Rmax R0 Ri Rif
an empirical frontal constant CE ¼ Ri1/2 f entrainment coefficient (the ratio of entrainment velocity to flow speed) acceleration due to gravity (9.81 m/s2) reduced density g0 ¼ g,rcloudr ramb amb height of the gravity current height of the gravity current when the flow becomes critical friction coefficient k ¼ uU* index in power law representations of velocity or height volume flow of heavy gas radius from centre of source radius of the vapour current front radius at which the flow becomes critical Richardson number Ri ¼ g0H U2 Richardson number at the front Rif ¼ g0H U2
Ris t Uf U u* U0
Richardson number at the outer edge of the source time speed of progress of the front of the vapour current flow speed in the gravity current (a function of radius) shear velocity (the surface shear stress is ru2* ) vapour current speed when the flow becomes critical
CE E g g0 H H0
f
zrough
roughness length of the surface
ramb density of ambient air k von Karman's constant r, rcloud density of vapour current
References Atkinson, G., Hall, J., 2015. A Review of Large Vapour Cloud Incidents. HSL Report MH15/80. http://primis.phmsa.dot.gov/meetings/MtgHome.mtg?mtg¼111. Atkinson, G., Cowpe, E., Halliday, J. and Painter, D., (2016) A Review of Very Large Vapour Cloud Explosions. 19th Mary Kay O'Connor International Process Safety Symposium, College Station TX October 2016. Briggs, G.A., Thompson, R.S., Snyder, W.H., 1990. Dense gas removal from a valley by crosswinds. J. Haz. Mat. 24, 1e38. Buncefield Major Incident Investigation Board, 2007. The Buncefield Incidente11th December 2005ethe Final Report of the Major Incident Investigation Board, vol. 1, ISBN 978-07176-6270-8. http://www.hse.gov.uk/comah/buncefield/miibfinal-volume1.pdf. Byron Bird, R., Stewart, W.E., Lightfoot, E.N., 2007. Transport Phenomena, second ed. John Wiley and Sons, New York. CSB, 2015. Caribbean Petroleum Tank Terminal Explosion and Multiple Tank Fires (October 23rd 2009). US. Chemical Safety and Hazard Investigation Board. Final Investigation Report. DNV, 2013. Phast Version 7 Software. DNV Software, London, UK. Available from: http://www.dnv.com/software. Accessed 19 August 2013. FABIG, 2013. FABIG Technical Note 12-Vapour Cloud Development in Over-filling Incidents. April 2013. Available from: http://fabig.com/video-publications/ TechnicalGuidance#. Gant, S., 2006. Buncefield Investigation: Dispersion of the Vapour Cloud. HSL Report CM/06/13. Gant, S.E., Narasimhamurthy, V.D., Skjold, T., Jamois, D., Proust, C., 2014. Evaluation of multi-phase atmospheric dispersion models for application to carbon capture and storage. J. Loss Prev. Process Industries 32 (November 2014), 286e298. Hallworth, M.A., Huppert, H.E., Phillips, J.C., Sparks, R.S.J., 1996. Entrainment into two-dimensional and axisymmetric turbulent gravity currents. J. Fluid Mech. 308, 289e311. Havens, J., Spicer, T., 1988. A Dispersion Model for Elevated Dense Gas Jet Chemical Releases. EPA Report EPA-450/4-88-006a. HSE SPC/TECH/GEN 43 Land use planning advice around large scale petrol storage sites http://www.hse.gov.uk/foi/internalops/hid_circs/technical_general/spc_ tech_gen_43/. HSE 2007 http://www.hse.gov.uk/consult/condocs/cd211.htm. Johnson, C.G., Hogg, A.J., 2013. Entraining gravity currents. J.Fluid Mech. 731, p471e508. Lofquist, K., 1960. Flow and stress near an interface between stratified liquids. Phys. Fluids 3, 158e175. Slim, A.C., Huppert, H.E., 2011. Axisymmetric, constantly supplied gravity currents at high Reynolds number. J. Fluid Mech. 675, 540e551. Thorpe, S.A., 1973. Experiments on instability and turbulence in a stratified shear flow. J. Fluid Mech. 61, 731e751. Turner, J., 1973. Buoyancy Effects in Fluids. Cambridge University Press. Turner, J., 1986. Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows. J. Fluid Mech. 173, 431e447. Venart, J.E.S., Rogers, R.J., 2011. The Buncefield Explosion: Vapour Cloud Dispersion and Other Observation. Hazards XXII, IChemE Symposium Series No. 156.