- Email: [email protected]
Vented gaseous deflagrations
Journal of Hazardous Materials A116 (2004) 1–10
Vented gaseous deflagrations: modelling of hinged inertial vent covers V.V. Molkov∗ , A.V. Grigorash, R.M. Eber, D.V. Makarov University of Ulster, FireSERT, Shore Road, Newtownabbey, Co. Antrim, Northern Ireland BT37 0QB, UK Received 2 December 2003; received in revised form 3 August 2004; accepted 3 August 2004 Available online 19 October 2004
Abstract The model of explosion pressure build up in enclosures with inertial vent covers and the CINDY code implementing the model are validated against experiments by H¨ochst and Leuckel (1998) in a 50 m3 vessel with a pair of ceiling-mounted upwards-opening hinged doors in a ‘butterfly’ configuration with surface densities of 73 and 124 kg/m2 under conditions of initially quiescent and turbulent mixtures. The model and the code are further validated against an experiment by Zalosh (1978) in a 33.5 m3 room-like enclosure with a pair of wall-mounted rectangular doors, in a parallel configuration, each hinged at its bottom edge with a surface density of 23.1 kg/m2 and initially quiescent mixture. A formula for the torque acting upon a rotating venting door is derived under conditions of vent cover jet formation. The vent cover jet effect decreases the torque three times compared to an elementary approach valid at the start of vent cover movement. It is demonstrated that, similar to translating vent covers, the vent cover jet effect is crucial for prediction of interdependent vent cover displacement in time and pressure transients. © 2004 Elsevier B.V. All rights reserved. Keywords: Explosion/deflagration; Inertial vent cover; Vent cover jet effect; Venting generated turbulence; Surface density/inertia
1. Introduction The first theories of vented gaseous deflagration dynamics were developed by Yao [1] and Pasman et al. [2] in 1974 with the turbulence factor as a lumped parameter. A detailed theory of vented deflagration in spherical vessels was published by Bradley and Mitcheson [3] in 1978. In 1981, Molkov and Nekrasov [4] suggested a deflagration theory with two lumped parameters, i.e. the turbulence factor (χ) and generalized discharge coefficient (µ). This theory has been shown in the following years to predict the dynamics of gaseous deflagrations reasonably well for both closed and vented enclosures for a wide range of explosion conditions. The history of development of this model is outlined in Molkov [5]. Recently, Razus and Krause [6] published a comparative study of a number of lumped parameter models. They demonstrated that model compared favourably to its analogues in predict∗
Corresponding author. Tel.: +44 28 9036 8731; fax: +44 28 9036 8700. E-mail address: [email protected] (V.V. Molkov).
0304-3894/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jhazmat.2004.08.027
ing explosion overpressures. In 2002, Hirano [7] presented the correlation for turbulence generated during vented deflagrations, published for the first time in [5]; in his invited paper, Williams [8] cited findings in scaling of explosions based on the model [9]. Russo and coworkers [10], in their study of pressure piling in connected vessels, demonstrated that only “the universal correlation by Molkov [11] for vent sizing at initially elevated pressures is quite able to reproduce almost all overpressures observed in the second vessel in 115 experiments”. The accuracy of model predictions has been found to be significantly higher than predictions made using the approach offered in NFPA 68 [12] and [13]. Modifications to the model to account for translating inertial vent covers were presented recently in Molkov et al. [14]. Herein a derivation of the equations governing behaviour of hinged inertial covers will be presented. Discussion of the development of empirical coefficients supporting the hinged inertial vent model is also presented. Finally, the results of the model and the code validation against published experimental data are shown.
2
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Nomenclature a
radius of spherical vessel of equivalent volume (m) A fraction of cross-section area of vent occupied by burnt gas Ajet empirical coefficient b door length (m) cui speed of sound in gas (m/s) Cjet empirical coefficient Ei combustion products expansion coefficient at initial conditions fjet (ϕ) transient factor F area (m2 ) g acceleration gravity, 9.80665 m/s2 H/D ratio of height to diameter H × W × L height, width and length of enclosure Jm moment of inertia of the vent (kg m2 ) l coordinate along door width, 0 ≤ l ≤ L (m) L door width (length of the pivoting side) (m) M molecular mass (kg/kmol) m mass (kg) n relative gas mass inside the vessel p pressure (Pa) rb burnt gas equivalent sphere radius (m) R universal gas constant, 8314.4 J/K/kmol R outflow contribution R# outflow parameter Sui mixture burning velocity at initial conditions (m/s) t time (s) T temperature (K) Tfull full torque (N m) Tgravity torque exerted by the gas pressure on the hinged door (N m) Tpressure torque generated by the force of gravity (N m) u1 flow velocity through the vent cross-section (m/s) uL flow velocity in the changing venting area (m/s) V enclosure volume (m3 ) w inertia (surface density) (kg/m2 ) W transient venting parameter Z auxiliary quantity Greek α ε γ µ ϕ χ π0 π
angular acceleration (radian/s2 ) overall thermokinetic exponent ratio of specific heats generalized discharge coefficient angle of opening of a hinged door turbulence factor pi number, 3.141593. . . dimensionless pressure, p/pi
ρ σ τ ω
density of gas (kg/m3 ) dimensionless density, ρ/ρi dimensionless time, =t Sui /a angular speed (radian/s)
Subscripts and superscripts 1 in the vent cross-section a at the atmospheric pressure level b burnt gases closed at closed door conditions f flame front full full torque gravity torque generated by the force of gravity i initial state j vent number, summation index jet pertaining to the jet effect L at changing venting area new new value at the end of the integration step N nominal, 100% open vent pressure gas pressure force S sphere u unburnt gases v venting, latching
2. Equations of vented deflagration dynamics for multiple vents Derivation of the equations of vented gaseous deflagration dynamics for the case of a single non-inertial venting device was published in [4] and [5]. The derivation relied upon: conservation laws (mass, volume, energy), ideal gas state equation, and standard orifice equations for calculation of mass flow rate for subsonic and sonic regimes of outflow. Ratio the real flame front area Ff (t) at any moment t to the surface area Fs (t) = 4π0 rb2 of a sphere of radius rb (where b, burnt gases; s, sphere and π0 , the number ‘pi’) to which burnt gas inside the enclosure could be gathered at the same moment. This ratio is called a turbulence factor χ(t) = Ff (t)/Fs (t). For the ideal case of laminar spherical flame propagation the turbulence factor is constant during the course of deflagration and equal χ = 1. For real large-scale explosion problems values of χ two orders higher, i.e. up to 100, can be expected [15]. The simultaneous discharges from multiple vents add up, and that results in the following system of governing equations [16]: 2/3
dπ χ(τ)Zπε+1/γu (1 − nu π−1/γu ) − γb W (τ)R = 3π , dτ π1/γu − (γu − γb /γu )nu (1)
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
R = R#u
j [1 − Aj (τ)]µj Fj (τ)
j µj Fj (τ) 1/γu − nu j Aj (τ)µj Fj (τ) # π + Rb , nb j µj Fj (τ)
(2)
dnb 2/3 = 3 χ(τ)πε+1/γu (1 − nu π−1/γu ) dτ j Aj (τ)µj Fj (τ) # − Rb W (τ) , j µj Fj (τ) dnu 2/3 = −3 χ(τ)πε+1/γu (1 − nu π−1/γu ) dτ [1 − A (τ)]µ F (τ) j j j j + R#u W (τ) , j µj Fj (τ)
sonic flow conditions. For subsonic regime, the outflow parameter is equal to 1/2 2γ pa 2/γ pa (γ+1)/γ # (5) πσ − R = γ −1 pi π pi π and for sonic regime, it is equal to γ+1/γ−1 1/2 2 # R = γ , πσ γ +1
(3)
(4)
where π is the dimensionless pressure (=p/pi , where i is the initial state, p is the pressure (Pa), pi is the initial pressure in the vessel (Pa)), τ is the dimensionless time (=t Sui /a, where Sui laminar burning velocity at initial conditions (m/s) and a is the radius of spherical vessel of equivalent volume (m)), ε is the overall thermokinetic exponent, γ u is the adiabatic exponent (ratio of specific heats) for unburnt mixture (where u is the unburnt gases), γ b is the adiabatic exponent (ratio of specific heats) for burnt mixture (where b is the the burnt gases), nu is the relative mass of unburnt mixture inside the vessel (=mu /mi , where m is the mass (kg)), nb is the relative mass of burnt mixture inside the vessel (=mb /mi ), A is the fraction of cross-section area of vent occupied by burnt gas (j is the vent number or summation index), µj is the generalized discharge coefficient for the jth vent, F is the vent area (m2 ), R is the outflow contribution, where R# is the outflow parameters (defined below), and Z is the auxiliary quantity: γu γb − 1 1−γu /γu γb − γu Z = γb Ei − + π , γb γu − 1 γu − 1 where Ei is the combustion products expansion coefficient at initial conditions. W (τ) is the transient venting parameter cui j µj Fj (τ) 1 W (τ) = √ , √ 3 V 2/3 36π0 γu Sui where cui = (γ u RTui /Mui )1/2 is the speed of sound in unburnt gas (m/s) (where R = the universal gas constant (J/K/kmol), =8314.41; Tui is the temperature of unburnt gas at initial conditions (K) and Mui is the molecular mass of unburnt mixture at initial conditions (kg/kmol)), and V is the enclosure volume (m3 ). The outflow parameters R#u and R#b for unburned and burned mixture in Eqs. (1), (3) and (4) arise from the orifice equations and are calculated differently for subsonic and
3
(6)
where σ is the relative density of gases (σ u for unburnt 1/γ 1/γ gases, = ρu /ρi = πu ; σ b for burnt gases, = ρb /ρi = πb ; ρu 3 is the density of unburnt gases (kg/m ); ρb is the density of burnt gases (kg/m3 ); ρi is the initial density of unburnt gases (kg/m3 ), =mi /V) and pa is the atmospheric pressure outside the vessel (Pa). The unburned and burned versions R#u and R#b of R# are obtained from Eqs. (5) and (6) by substituting the unburnt and burnt versions of γ and σ in these formulae, respectively. The condition of transition from subsonic to sonic flow regime is pa 1 + γ γ/γ−1 π≥ (7) pi 2 Again, this is calculated separately for unburned and burned mixture, with γ u and γ b . 3. Modelling of hinged vent covers Eqs. (1), (3) and (4) above depend on the current venting area that changes with time. The character of this change should be specified. This allows vent covers of any type to be ‘plugged in’ the calculations, as long as the value of the current venting area F(t) can be calculated. For each vent cover, in conditions of pressure growing with time, at some moment tvj , when the gas pressure is equal to the pre-set ‘latch release’ pressure pvj (Pa), the release of vent cover ‘j’ occurs, and outflow of gases from the enclosure through the vent ‘j’ begins into the surrounding atmosphere. Depending on the vent cover type, the venting area either immediately becomes equal to the nominal vent area FN (non-inertial vent covers, or rupture membranes) (where N is the nominal area) or increases gradually with time while vent cover moves away by the pressure force. The focus of this paper is on hinged covers. Translating covers are presented in our earlier paper [14]. A hinged ‘door’ or ‘cover’ is an inertial cover modelled as a solid rectangle able to swing about one of its edges, the hinge, fixed on the enclosure, e.g. as shown in Fig. 1. 3.1. Venting area Denote b the length of the hinged side (m), i.e. the length of the door; L the length of the pivoting side (m), i.e. the width of the door. Then the nominal area of the vent opening and the area of the hinged door is FN = bL. It is further assumed that
4
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
the vent cover mass is distributed uniformly over the surface of the cover, with surface density or inertia w (kg/m2 ). Let ϕ be the angle between the vent opening and the hinged door. It is assumed that the current venting area F(ϕ) is the gap area between the edges of the cover and the vent opening. The gap, as shown in Fig. 1, is formed from one rectangular region, based on the door edge opposite to the hinge and two triangular regions, based on the pivoting edges of the door. The venting area is then:
ϕ ϕ F (ϕ) = min FN , 2L sin b + L cos (8) 2 2 This area is zero for a closed vent (ϕ = 0) and is allowed to increase until ϕ = ϕN ; it reaches the maximum value equal to the nominal vent area FN . Also assume that for angles ϕ > ϕN , the venting area stays equal to FN . Furthermore, assume that the door is inelastically arrested at ϕ = 90◦ . 3.2. Pressure distribution for F(ϕ) ≤ FN When the vent is closed, the gas pressure is uniform throughout the door surface, and is equal to p(t), the pressure inside the enclosure (Pa). Furthermore, the gas mass discharge rate is zero. When the vent is open, the picture changes. First, the static pressure of the escaping gases on the door is smaller than the pressure at ‘stagnation’ conditions inside the enclosure. Second, the pressure is not uniform along the door surface any more. The first issue is addressed by using the Bernoulli or mass conservation law for the gas flowing between the enclosure inside the vent cross-section and the current venting area. At the vent cross-section, let p = p1 , and u = u1 (where u1 is the flow velocity through the vent cross-section (m/s)), and flows are low enough to warrant the assumption of incompressibility. The Bernoulli’s equation for these two levels allows p1 to be expressed as: p1 = p(t) −
ρ u21 , 2
Outside the vessel the gas pressure is p = pa , and the average gas velocity in the changing venting area is u = uL . From the Bernoulli’s equation relating the inside and the outside of the vessel, uL can be expressed as: 2(p(t) − pa ) 1/2 uL = (10) ρ The mass conservation law between the vent opening and the outside of the vessel gives an expression for u1 : u1 =
uL F (ϕ) bL
(11)
Substituting (10) in (11) and the results in (9) the pressure p1 in the vent cross-section becomes p1 = p(t) − (p(t) − pa )
F 2 (ϕ) (bL)2
(12)
This pressure depends on both the current explosion pressure and the current angle of the door opening. The second issue is dealt with by assuming that the pressure along the door surface changes as a linear function of the position on the width of the door: p(l, t) = p1 −
p1 − pa l L
(13)
Here l ≤ L is the current position, and p1 is defined by (12). The assumption of a linear pressure distribution along the width of an inertial hinged vent cover is a simplification of a more complex three-dimensional distribution. A simple modelling including a portion of an enclosure, hinged vent cover, and a portion of the surrounding environment, was performed using FluentTM . A constant internal pressure of 1.3 atm and external pressure of 1.0 atm were assumed. The problem was solved using a 2D approach, and the steady pressure distribution along the cover was determined for several cover opening angles. The results are shown in Fig. 2. Ultimately
(9)
where ρ is the gas density constant (kg/m3 ).
Fig. 1. Hinged door.
Fig. 2. Pressure distribution on inertial hinged vent cover as a function of cover opening angle for constant enclosure pressure. Steadystate problem, p/p0 = 1.3, 2D segregated solver, k–ε turbulence model, ♦, , —transients ‘below’ the cover, i.e. at cover surface hit by the escaping gases; ,,䊉—transients ‘above’ the cover.
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
a three-dimensional modelling, using CFD, will be necessary to accurately determine the pressure distribution on the hinged cover. As a first approximation, however, a linear distribution will force the model cover displacement to more closely match experimental experience than a model assuming the enclosure pressure everywhere on the cover.
5
these considerations in mind, it is possible to replace (16) and (17) with a single formula bL2 F 2 (ϕ) [p(t) − pa ] 1 − , 6fjet (ϕ) (bL)2
Tpressure =
(18)
where the transient factor fjet (ϕ) is 3.3. Pressure force torque for F(ϕ) ≤ FN
The torque exerted by the gas pressure on the door is (N m): L Tpressure =
[p(l, t) − pa ] bl dl
(14)
0
The pressure torque is always positive, since the gas flow momentum always works for opening the vent. Using formulae (12) and (13) by substitution in (14), the torque the gas applies to turn the door on its hinge can be written as: bL2 F 2 (ϕ) Tpressure = (15) [p(t) − pa ] 1 − 6 (bL)2 Our assumption of linear change of pressure in (13) results in very easy derivations. However, the true pressure distribution is not linear, the whole vent cover-pushing phenomenon being three-dimensional, non-stationary and dependent on the geometry of the cover and vent opening. Therefore, we decided to settle on the linear distribution augmented by an empirical ‘jet factor’ Cjet as follows: bL2 F 2 (ϕ) Tpressure = (16) [p(t) − pa ] 1 − 6Cjet (bL)2 The Cjet will compensate for the true non-linear, nonstationary and geometry-dependent character of the hinged door movement. Notice that when the door is closed, ϕ = 0, formula (15) gives three times smaller value than the correct torque for the closed door should be: Tpressure,closed =
bL2 [p(t) − pa ] 2
(17)
Therefore, we have to use different formulae for a closed or almost closed door and a sufficiently wide open door. When the door is shut or is opened within some small range of angles, we apply formula (17). Above a certain angle, when F(ϕ) reaches a certain fraction Ajet of the full area bL, the velocity of gases escaping along the door surface becomes significant, and we have to apply formula (16). To ensure the continuity of transition from the ‘pressure’ regime of formula (17) to the ‘jet’ regime of formula (16), we assume this transition to be linear with respect to F(ϕ) and controlled by the ‘jet fraction’ parameter Ajet . We also have to ensure that formula (16) never produces a value greater than formula (17) for the same pressure and cover dimensions. To that end, we restrict the values of the ‘jet factor’ by the condition Cjet > 1/3. With
1 Cjet F (ϕ) fjet (ϕ) = min max , , Cjet 3 bLAjet
(19)
Indeed, when ϕ = 0, we have F(ϕ) = 0, the transient factor fjet (ϕ) equals 1/3, and formula (18) assumes the form of (17). At a certain angle ϕ, when F(ϕ) = Ajet bL, fjet (ϕ) becomes equal to Cjet , and stays at this value for any greater angle of opening. Respectively, formula (18) assumes the form of (16). For all the interim angles, formula (19) produces interim values of fjet (ϕ) between 1/3 and Cjet . The values for the controlling parameters Ajet and Cjet have to be determined empirically. If comparison with an experiment shows that Ajet assumes relatively low values, say less than 10% of the nominal area, and Cjet turns out to be reasonable, say less than two, then formula (18) may be deemed a plausible approximation of the physical reality. 3.4. Balance of torques for F(ϕ) ≤ FN In rotational movement of our hinged door, the balance of torques (moments about the axis of rotation) can be expressed as follows: Tfull = Tpressure + Tgravity
(20)
where Tfull is the full torque, Tpressure is the torque generated by the pressure of escaping gases, formula (18), and Tgravity is the torque generated by the force of Earth’s gravity. The full torque (N m) is Tfull = Jm × α, where α is the angular acceleration of the vent (radian/s2 ) and Jm is the moment of inertia of the vent (kg m2 ). The general formula for Jm (kg m2 ) is Jm = m l2 dm, where dm is the elementary mass (kg) forming the body, and l is the distance of this elementary mass from the axis of rotation (hinge) (m). With the assumption of uniform mass distribution over the vent surface, the result is: dm = wdF , where w is the inertia of vent cover and dF is the area of a surface element. Then Jm = w l2 dF FN
Since the vent is rectangular, this integral can be rewritten as L Jm = wb
l2 dl = 0
wbL3 3
6
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
As the inertia w is expressed as w = m/(bL) where m = mass (kg), the moment of inertia simplifies to Jm =
mL2 3
The full torque is therefore Tfull =
αmL2 3
The torque Tgravity of the gravity force depends on where the vent is mounted and what side of the vent is hinged. If the vent is mounted on a wall, and is hinged at its side edge, such that its swinging motion is in the horizontal plane, then Tgravity = 0. If the vent is mounted on a wall, and is hinged at its top edge, then Tgravity = −
mgL sin(ϕ) , 2
where g is the acceleration gravity (m/s2 ). Tgravity is negative because for a top-hinged wall-mounted vent the gravity works against the opening vent. If the vent is bottom-hinged to a wall, then the difference is only in the sign Tgravity = +
mgL sin(ϕ) , 2
since gravity helps the vent to open in this case. Analogously, if the vent is mounted on the ceiling and opens upwards, then Tgravity = −
mgL cos(ϕ) , 2
and for vents mounted in the floor and opening downwards, Tgravity = +
mgL cos(ϕ) 2
For wall-mounted vents, gathering formula (18) and the above formulae for the torques in the Eq. (20) yields: αmL2 bL2 (p(t) − pa ) F 2 (ϕ) mgL sin(ϕ) = 1− − , 2 3 6fjet (ϕ) 2 (bL) (21) with a little change in notation for g. We assign the ‘direction’ for the gravitational force, such that g > 0 corresponds to the situation when the hinged edge is at the top; g < 0 when the hinged edge is at the bottom; and g = 0 when the hinged edge is at a side. In a similar way, for ceiling/floor-mounted vents, sin ϕ in formula (21) should be changed to cos ϕ, and it is assumed that g > 0 for ceiling-mounted vents that open upwards; g < 0 for floor-mounted vents that open downwards; and g = 0 for gravity compensated for by a spring or a balance. A simple algebra in Eq. (21) yields the following formula for the angular acceleration of the vent. For a wall-mounted
vent, 3.0 α= 2.0
b (p(t) − pa ) F 2 (ϕ) g 1− − sin(ϕ) , m3fjet (ϕ) L (bL)2
and for a ceiling-mounted vent, sin (ϕ) should be changed to cos (ϕ). The angular speed ω (radians/s) and the angle of opening ϕ (radians), are t ω(t) =
t α(τ)dτ
and
ϕ(t) =
0
ω(τ)dτ,
respectively.
0
Within each step of numerical integration of the governing equations, the angular acceleration is assumed constant. Therefore, uniformly accelerated motion formulae are used to obtain the new values ωnew and ϕnew of the angular speed and the angle of opening at the end of the integration step from the values at the beginning of this step: ωnew = ω + α dt and ϕnew = ϕ + ω dt + α dt2 /2.0. At the end of each integration step the following assignments take place: ω = ωnew ; ϕ = ϕnew . 3.5. Torque when F(ϕ) > FN The applicability of formula (18) is limited in the angle ϕ of the door opening, in that the formula will work only under the assumption that the current venting area F(ϕ) is less than the nominal area FN . At a certain angle ϕN such that F(ϕN ) = FN the pressure at the vent cross-section is equal to the atmospheric pressure, and the gas flow through the vent is unrestricted, as though there had been no vent cover at all. The CINDY code assumes, when this point is reached, that the vent cover displacement continues; yet the energy imparted to the cover through momentum earlier in the deflagration is affecting cover motion. Since the primary interest in the current research was in the influence of vent cover inertia while the cover position and displacement could influence transient enclosure pressures, detailed analyses of cover displacement once the vent was found to be 100% opened have been neglected. Further modelling would be required to assure an accurate representation of cover behaviour.
4. Comparison with experiments 4.1. Values of the empirical coefficients The empirical coefficients Cjet and Ajet in formula (19) were determined through matching of CINDY simulations of H¨ochst and Leuckel [17] experiments 3-B and 3-D (‘experiment 3-B’ or ‘3-B’ denotes experimental results plotted in Fig. 3b of their paper; ‘experiment 3-D’ or ‘3-D’ denotes experimental results plotted in Fig. 3d of their paper). Ajet for hinged covers represents a threshold below which jetting
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
7
Fig. 3. (a) Cjet —H¨ochst and Leuckel, experiment 3-B: opening angle, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (b) Cjet —H¨ochst and Leuckel, experiment 3-B: pressure, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (c) Cjet —H¨ochst and Leuckel, experiment 3-B: χ, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (d) Ajet —H¨ochst and Leuckel, experiment 3-B: opening angle, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (e) Ajet —H¨ochst and Leuckel, experiment 3-B: pressure, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (f) Ajet —H¨ochst and Leuckel, experiment 3-B: χ, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
flows of gases escaping through the opening vent are assumed to not yet be established. The variation of the vent cover displacement and the enclosure pressure for varying Cjet , assuming constant values of Ajet , are shown in Fig. 3a and b for displacement and pressure, respectively (in these and following figures simulated curves are shown until the moment of full burnout of the mixture inside enclosure). Once the values for the coefficients were selected, the values of χ and µ to backfit the simulations to the experimental data were found. Fig. 3c shows the backfitted values of χ that were found (results were obtained with a constant discharge coefficient µ = 1.2). As the figures show, increasing Cjet decreases the amount of force applied to the door—so the door takes longer to open. A slower opening door ‘generates’ less enclosure turbulence—thus the reductions in turbulence with increasing Cjet , as shown in Fig. 3c. The best fit seems to come with a value of Cjet of 1.4. A similar process, holding Cjet constant, was used to explore the impact of varying Ajet . Fig. 3d and e compares
computed and experimental displacements and pressures, respectively, for different values of Ajet , while Fig. 3f shows the resultant turbulence levels corresponding to the selected values of Ajet . For Ajet values at or below 0.01, the pressure is overestimated, yet the cover is not moving fast enough. Thus, values of Ajet at or below 0.01 should not be used. However, values of Ajet from 0.05 to 0.10 could be successfully used, although the turbulence rises as Ajet rises (Fig. 3f). As the figures show, increasing Ajet (in effect the amount of time the door is exposed to enclosure rather than door jetting conditions), ‘increases’ enclosure turbulence—the cover opens faster. The best fit seems to be achieved at an Ajet of 0.05. Note that given the data shown in Fig. 3, an estimation of the sensitivity of the enclosure turbulence to the Ajet and Cjet may be estimated. An average shift in χ as Ajet or Cjet varies may be determined by observing the shift in χ in the backfitted curves, and comparing the shift that occurs for various values of the coefficients. As a result, it was calculated that
8
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Fig. 4. (a). H¨ochst and Leuckel, experiments 3-B, 3-D: opening angle, µ = 1.2 ( –vent starts to open, 䊉–vent 100% open). (b). H¨ochst and Leuckel, experiments 3-B, 3-D: pressure, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
when Ajet shifts ±10%, χ shifts ±1.5%. Similarly, when Cjet shifts ±10%, χ shifts in opposite direction ±6%. From this, it would seem that the enclosure turbulence is more sensitive to the amount of force applied to the vent cover than to the amount of time that force is applied. The discrepancies between the experimental and simulation results in pressure dynamics, especially those in the peak areas, could be explained by heat losses to enclosure walls, which were not modelled in this study. From explosion safety engineering point of view it is acceptable as calculated pressure peaks are conservative relative to experimental data. 4.2. Validation 1: H¨ochst and Leuckel’s experiments The model has been validated against H¨ochst and Leuckel’s experiments [17]. For translation panels, a detailed description of that validation was presented in [14]. Their apparatus consisted of a 50 m3 silo of reinforced concrete with H/D = 4. The vent covers for these experiments were a pair of hinged vent covers arranged in a ‘butterfly’ configuration on the top surface of the silo. The edge of each cover most remote from the hinges was padded so that when the covers opened to 90◦ , the impact of the covers on each other was minimised. Experiment 3-B was a quiescent mixture of 10.7% methane–air, with a total venting area (for the two covers) F = 1.91 m2 , an inertia w = 124 kg/m2 , and torque of 532 N m. Experiment 3-D was a turbulent mixture (mixed by a fan within the enclosure for the purpose) of 10.6% methane–air, with a total venting area (for the two covers) F = 1.91 m2 , an inertia w = 73 kg/m2 , and torque of 314 N m. The mixtures were ignited by an electric match with energy of 75 J, located 3.5 m from the floor at the silo centre line. Sui was 0.38 m/s; Ei was 7.4; γ u was 1.39; γ b was 1.25;cui was 353 m/s; F was 2.45 m2 ; ε was 0.3; and the molecular mass (M) was 27.24 kg/kmol. Fig. 4a compares the calculated opening angle transients with the experimental results, while Fig. 4b compares the calculated enclosure pressures with the experimental results. A good match of both displacement and pressure has been achieved, using the Cjet and Ajet coefficient values of 1.4 and 0.05, respectively. As expected, the pre-mixed enclosure
gases in experiment 3-D results in faster opening of the cover, and higher and earlier pressures than those experienced by the quiescent mixture in experiment 3-B. These results are well matched by the CINDY computations. Fig. 5 shows how χ varied as the calculations progressed, and the best fit values of χ for experiments 3-B and 3-D. The enclosure gases were initially turbulent in experiment 3-D—thus the higher χ was needed to backfit the pressure and displacement are not surprising. In the CINDY code, χ is implemented as piecewise-linear; over some ranges of backfitted values, many small increments in χ can be replaced by a few longer segments with little or no detrimental change to the backfit for either pressure or displacement. Data from future experiments will be required to further determine if a two- or three- step χ curve established by simple rules could predict deflagration dynamics. 4.3. Validation 2: Zalosh’s experiment The model has also been validated against an experiment of Zalosh [18]. His apparatus consisted of a rectangular concrete bunker with interior dimensions H × W × L = 3.1 m × 2.0 m × 5.4 m and a volume of 33.5 m3 . The only relief was a pair of wall-mounted blowoff panels, arranged vertically in parallel, hinged on the bottom edges, opening outwards and downwards, arrested by cables in the 90◦ opened position. Each door was sized to be 1.29 m2 , for a total venting area of 2.58 m2 . An initially quiescent near stoichiometric 10.0% mixture was ignited by a 12 J spark
Fig. 5. H¨ochst and Leuckel, experiments 3-B, 3-D: χ, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Fig. 6. Zalosh: pressure, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
9
• The obtained empirical coefficients Cjet and Ajet assume plausible values of 1.4 and 0.05. Deviation of the empirical parameters Cjet and Ajet from values derived herein results in only small variances in the backfitted values of the turbulence factor χ. • Observation of gradually changing turbulence factor χ over the processed experiments suggests that χ may be capable of being represented by a simple curve. The processing of additional experimental pressure and displacement test data would be useful to assess how predictable the gradual change of χ may be.
Acknowledgements The authors wish to thank Dr. Siegfried H¨ochst for providing additional information about the structure of the vent panel arresting mechanisms used in the experiments he performed with Dr. Wolfgang Leuckel. Funding by the Engineering and Physical Sciences Research Council (EPSRC) through Grant #GR/R53333/01 is gratefully acknowledged. Fig. 7. Zalosh: χ, µ = 1.2. ( —vent starts to open, 䊉—vent 100% open).
located in the centre of the bunker. In the calculations, the burning velocity, Su , was 0.38 m/s; the combustion products expansion coefficient at initial conditions, Ei , was 7.4; the unburnt mixture adiabatic exponent, γ u , was 1.39; the burnt gas adiabatic exponent, γ b , was 1.25; the speed of sound, cui , was 353 m/s; the overall thermokinetic exponent, which gives the dependence of burning velocity on pressure at adiabatic compression conditions, ε, was 0.3; and the molecular mass, M, was 27.56 kg/kmol. Coefficient values of Cjet = 1.4 and Ajet = 0.05, determined by comparisons with experiments by H¨ochst and Leuckel [17], were used. Fig. 6 compares the calculated pressure transients with the experimental results for hinged covers. No displacement data were published in Zalosh [18], so there is only a pressure comparison. The comparison to the experimentally observed pressure is reasonable. Fig. 7 shows how χ varied as the deflagration progressed, and the best fit values of χ for the experiment. The pressure rise at the end of the experiment seems to be due to turbulisation at the end of the process and a significant increase in χ is required to match the pressure.
5. Conclusions • Modelling of vented deflagrations with inertial venting devices has been performed for the case of multiple hinged doors. • The developed model and the CINDY code have been validated against large-scale experiments of H¨ochst and Leuckel [17] and Zalosh [18]. Good matches in both explosion overpressure transients and cover displacement transients have been achieved.
References [1] C. Yao, Explosion venting of low-strength equipment and structures, Loss Prev. 8 (1974) 1–9. [2] H.L. Pasman, Th.M. Groothuizen, H. de Gooijer, Design of pressure relief vents, in: C.H. Buschman (Ed.), Loss Prevention and Safety Promotion in the Process Industries, New York, 1974, pp. 185–189. [3] D. Bradley, A. Mitcheson, The venting of gaseous explosions in spherical vessels, Combust. Flame 32 (1978) 221–236, 237255. [4] V.V. Molkov, V.P. Nekrasov, Dynamics of gas combustion in a constant volume in the presence of exhaust, Combust. Explosion Shock Waves 17 (4) (1981) 363–369. [5] V.V. Molkov, Theoretical generalization of international experimental data on vented explosion dynamics, in: V.V. Molkov (Ed.), Proceedings of the First International Seminar on Fire and Explosion Hazards (17–21 July 1995, Moscow, Russia), All-Russian Research Institute for Fire Protection—Russian Association for Fire Safety Science, 1995, pp. 166–181. [6] D.M. Razus, U. Krause, Comparison of empirical and semi-empirical calculation methods for venting of gas explosions, Fire Saf. J. 36 (2001) 1–23. [7] T. Hirano, Combustion Science for Safety, in: J.H. Chen, M.D. Colket (Eds.), Proceedings of 29th International Symposium on Combustion, July 21–25, 2002, Sapporo, Japan), The Combustion Institute, Pittsburgh, 2002, pp. 167–180. [8] F. Williams, Approaches to the scaling of fires and explosions through mechanistic considerations, in: D. Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 445–456. [9] V.V. Molkov, A.Ya. Korolchenko, S.V. Alexandrov, Venting of deflagrations in buildings and equipment: universal correlation, in: Y. Hasemi (Ed.), Proceedings of the Fifth International Symposium on Fire Safety Science, Melbourne, Australia, March 3–7, 1997, pp. 1249–1260. [10] A. Di Benedetto, E. Salzano, G. Russo, Analysis of pressure piling by means of empirical correlation for vented vessels, in: D.
10
[11]
[12]
[13] [14]
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10 Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 413–424. V.V. Molkov, Unified correlations for vent sizing of enclosures against gaseous deflagrations at atmospheric and elevated pressures, J. Loss Prev. Process Industries 14 (2001) 567–574. V.V. Molkov, Explosion safety engineering: NFPA 68 and improved vent sizing technology, in: Interflam 1999, Proceedings of the Eighth International Fire Science Conference, June 29–July 1, 1999, Edinburgh Conference Centre, Scotland, pp. 1129–1134. NFPA 68, Guide for Venting of Deflagrations, 2002 Edition, National Fire Protection Association, USA, 2002. V.V. Molkov, A.V. Grigorash, F. Tamanini, R. Dobashi, R.M. Eber, Vented gaseous deflagrations: modelling of translating iner-
[15] [16]
[17] [18]
tial vent covers, J. Loss Prev. Process Industries 16 (5) (2003) 395– 402. V.V. Molkov, Venting of Gaseous Deflagrations, in: DSc Dissertation, VNIIPO, Moscow, Russia, 1996 (in Russian). A. Grigorash, R. Eber, V. Molkov, A theoretical model of vented gaseous deflagrations in enclosures with inertial vent covers, in: D. Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 445–456. S. H¨ochst, W. Leuckel, On the effect of venting large vessels with mass inert panels, J. Loss Prev. Process Industries 11 (1998) 89–97. R.G. Zalosh, Gas explosion tests in room-size vented enclosures, in: Proceedings of AIChE Loss Prevention Symposium, 1978, pp. 98–110.
Vented gaseous deflagrations: modelling of hinged inertial vent covers V.V. Molkov∗ , A.V. Grigorash, R.M. Eber, D.V. Makarov University of Ulster, FireSERT, Shore Road, Newtownabbey, Co. Antrim, Northern Ireland BT37 0QB, UK Received 2 December 2003; received in revised form 3 August 2004; accepted 3 August 2004 Available online 19 October 2004
Abstract The model of explosion pressure build up in enclosures with inertial vent covers and the CINDY code implementing the model are validated against experiments by H¨ochst and Leuckel (1998) in a 50 m3 vessel with a pair of ceiling-mounted upwards-opening hinged doors in a ‘butterfly’ configuration with surface densities of 73 and 124 kg/m2 under conditions of initially quiescent and turbulent mixtures. The model and the code are further validated against an experiment by Zalosh (1978) in a 33.5 m3 room-like enclosure with a pair of wall-mounted rectangular doors, in a parallel configuration, each hinged at its bottom edge with a surface density of 23.1 kg/m2 and initially quiescent mixture. A formula for the torque acting upon a rotating venting door is derived under conditions of vent cover jet formation. The vent cover jet effect decreases the torque three times compared to an elementary approach valid at the start of vent cover movement. It is demonstrated that, similar to translating vent covers, the vent cover jet effect is crucial for prediction of interdependent vent cover displacement in time and pressure transients. © 2004 Elsevier B.V. All rights reserved. Keywords: Explosion/deflagration; Inertial vent cover; Vent cover jet effect; Venting generated turbulence; Surface density/inertia
1. Introduction The first theories of vented gaseous deflagration dynamics were developed by Yao [1] and Pasman et al. [2] in 1974 with the turbulence factor as a lumped parameter. A detailed theory of vented deflagration in spherical vessels was published by Bradley and Mitcheson [3] in 1978. In 1981, Molkov and Nekrasov [4] suggested a deflagration theory with two lumped parameters, i.e. the turbulence factor (χ) and generalized discharge coefficient (µ). This theory has been shown in the following years to predict the dynamics of gaseous deflagrations reasonably well for both closed and vented enclosures for a wide range of explosion conditions. The history of development of this model is outlined in Molkov [5]. Recently, Razus and Krause [6] published a comparative study of a number of lumped parameter models. They demonstrated that model compared favourably to its analogues in predict∗
Corresponding author. Tel.: +44 28 9036 8731; fax: +44 28 9036 8700. E-mail address: [email protected] (V.V. Molkov).
0304-3894/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jhazmat.2004.08.027
ing explosion overpressures. In 2002, Hirano [7] presented the correlation for turbulence generated during vented deflagrations, published for the first time in [5]; in his invited paper, Williams [8] cited findings in scaling of explosions based on the model [9]. Russo and coworkers [10], in their study of pressure piling in connected vessels, demonstrated that only “the universal correlation by Molkov [11] for vent sizing at initially elevated pressures is quite able to reproduce almost all overpressures observed in the second vessel in 115 experiments”. The accuracy of model predictions has been found to be significantly higher than predictions made using the approach offered in NFPA 68 [12] and [13]. Modifications to the model to account for translating inertial vent covers were presented recently in Molkov et al. [14]. Herein a derivation of the equations governing behaviour of hinged inertial covers will be presented. Discussion of the development of empirical coefficients supporting the hinged inertial vent model is also presented. Finally, the results of the model and the code validation against published experimental data are shown.
2
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Nomenclature a
radius of spherical vessel of equivalent volume (m) A fraction of cross-section area of vent occupied by burnt gas Ajet empirical coefficient b door length (m) cui speed of sound in gas (m/s) Cjet empirical coefficient Ei combustion products expansion coefficient at initial conditions fjet (ϕ) transient factor F area (m2 ) g acceleration gravity, 9.80665 m/s2 H/D ratio of height to diameter H × W × L height, width and length of enclosure Jm moment of inertia of the vent (kg m2 ) l coordinate along door width, 0 ≤ l ≤ L (m) L door width (length of the pivoting side) (m) M molecular mass (kg/kmol) m mass (kg) n relative gas mass inside the vessel p pressure (Pa) rb burnt gas equivalent sphere radius (m) R universal gas constant, 8314.4 J/K/kmol R outflow contribution R# outflow parameter Sui mixture burning velocity at initial conditions (m/s) t time (s) T temperature (K) Tfull full torque (N m) Tgravity torque exerted by the gas pressure on the hinged door (N m) Tpressure torque generated by the force of gravity (N m) u1 flow velocity through the vent cross-section (m/s) uL flow velocity in the changing venting area (m/s) V enclosure volume (m3 ) w inertia (surface density) (kg/m2 ) W transient venting parameter Z auxiliary quantity Greek α ε γ µ ϕ χ π0 π
angular acceleration (radian/s2 ) overall thermokinetic exponent ratio of specific heats generalized discharge coefficient angle of opening of a hinged door turbulence factor pi number, 3.141593. . . dimensionless pressure, p/pi
ρ σ τ ω
density of gas (kg/m3 ) dimensionless density, ρ/ρi dimensionless time, =t Sui /a angular speed (radian/s)
Subscripts and superscripts 1 in the vent cross-section a at the atmospheric pressure level b burnt gases closed at closed door conditions f flame front full full torque gravity torque generated by the force of gravity i initial state j vent number, summation index jet pertaining to the jet effect L at changing venting area new new value at the end of the integration step N nominal, 100% open vent pressure gas pressure force S sphere u unburnt gases v venting, latching
2. Equations of vented deflagration dynamics for multiple vents Derivation of the equations of vented gaseous deflagration dynamics for the case of a single non-inertial venting device was published in [4] and [5]. The derivation relied upon: conservation laws (mass, volume, energy), ideal gas state equation, and standard orifice equations for calculation of mass flow rate for subsonic and sonic regimes of outflow. Ratio the real flame front area Ff (t) at any moment t to the surface area Fs (t) = 4π0 rb2 of a sphere of radius rb (where b, burnt gases; s, sphere and π0 , the number ‘pi’) to which burnt gas inside the enclosure could be gathered at the same moment. This ratio is called a turbulence factor χ(t) = Ff (t)/Fs (t). For the ideal case of laminar spherical flame propagation the turbulence factor is constant during the course of deflagration and equal χ = 1. For real large-scale explosion problems values of χ two orders higher, i.e. up to 100, can be expected [15]. The simultaneous discharges from multiple vents add up, and that results in the following system of governing equations [16]: 2/3
dπ χ(τ)Zπε+1/γu (1 − nu π−1/γu ) − γb W (τ)R = 3π , dτ π1/γu − (γu − γb /γu )nu (1)
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
R = R#u
j [1 − Aj (τ)]µj Fj (τ)
j µj Fj (τ) 1/γu − nu j Aj (τ)µj Fj (τ) # π + Rb , nb j µj Fj (τ)
(2)
dnb 2/3 = 3 χ(τ)πε+1/γu (1 − nu π−1/γu ) dτ j Aj (τ)µj Fj (τ) # − Rb W (τ) , j µj Fj (τ) dnu 2/3 = −3 χ(τ)πε+1/γu (1 − nu π−1/γu ) dτ [1 − A (τ)]µ F (τ) j j j j + R#u W (τ) , j µj Fj (τ)
sonic flow conditions. For subsonic regime, the outflow parameter is equal to 1/2 2γ pa 2/γ pa (γ+1)/γ # (5) πσ − R = γ −1 pi π pi π and for sonic regime, it is equal to γ+1/γ−1 1/2 2 # R = γ , πσ γ +1
(3)
(4)
where π is the dimensionless pressure (=p/pi , where i is the initial state, p is the pressure (Pa), pi is the initial pressure in the vessel (Pa)), τ is the dimensionless time (=t Sui /a, where Sui laminar burning velocity at initial conditions (m/s) and a is the radius of spherical vessel of equivalent volume (m)), ε is the overall thermokinetic exponent, γ u is the adiabatic exponent (ratio of specific heats) for unburnt mixture (where u is the unburnt gases), γ b is the adiabatic exponent (ratio of specific heats) for burnt mixture (where b is the the burnt gases), nu is the relative mass of unburnt mixture inside the vessel (=mu /mi , where m is the mass (kg)), nb is the relative mass of burnt mixture inside the vessel (=mb /mi ), A is the fraction of cross-section area of vent occupied by burnt gas (j is the vent number or summation index), µj is the generalized discharge coefficient for the jth vent, F is the vent area (m2 ), R is the outflow contribution, where R# is the outflow parameters (defined below), and Z is the auxiliary quantity: γu γb − 1 1−γu /γu γb − γu Z = γb Ei − + π , γb γu − 1 γu − 1 where Ei is the combustion products expansion coefficient at initial conditions. W (τ) is the transient venting parameter cui j µj Fj (τ) 1 W (τ) = √ , √ 3 V 2/3 36π0 γu Sui where cui = (γ u RTui /Mui )1/2 is the speed of sound in unburnt gas (m/s) (where R = the universal gas constant (J/K/kmol), =8314.41; Tui is the temperature of unburnt gas at initial conditions (K) and Mui is the molecular mass of unburnt mixture at initial conditions (kg/kmol)), and V is the enclosure volume (m3 ). The outflow parameters R#u and R#b for unburned and burned mixture in Eqs. (1), (3) and (4) arise from the orifice equations and are calculated differently for subsonic and
3
(6)
where σ is the relative density of gases (σ u for unburnt 1/γ 1/γ gases, = ρu /ρi = πu ; σ b for burnt gases, = ρb /ρi = πb ; ρu 3 is the density of unburnt gases (kg/m ); ρb is the density of burnt gases (kg/m3 ); ρi is the initial density of unburnt gases (kg/m3 ), =mi /V) and pa is the atmospheric pressure outside the vessel (Pa). The unburned and burned versions R#u and R#b of R# are obtained from Eqs. (5) and (6) by substituting the unburnt and burnt versions of γ and σ in these formulae, respectively. The condition of transition from subsonic to sonic flow regime is pa 1 + γ γ/γ−1 π≥ (7) pi 2 Again, this is calculated separately for unburned and burned mixture, with γ u and γ b . 3. Modelling of hinged vent covers Eqs. (1), (3) and (4) above depend on the current venting area that changes with time. The character of this change should be specified. This allows vent covers of any type to be ‘plugged in’ the calculations, as long as the value of the current venting area F(t) can be calculated. For each vent cover, in conditions of pressure growing with time, at some moment tvj , when the gas pressure is equal to the pre-set ‘latch release’ pressure pvj (Pa), the release of vent cover ‘j’ occurs, and outflow of gases from the enclosure through the vent ‘j’ begins into the surrounding atmosphere. Depending on the vent cover type, the venting area either immediately becomes equal to the nominal vent area FN (non-inertial vent covers, or rupture membranes) (where N is the nominal area) or increases gradually with time while vent cover moves away by the pressure force. The focus of this paper is on hinged covers. Translating covers are presented in our earlier paper [14]. A hinged ‘door’ or ‘cover’ is an inertial cover modelled as a solid rectangle able to swing about one of its edges, the hinge, fixed on the enclosure, e.g. as shown in Fig. 1. 3.1. Venting area Denote b the length of the hinged side (m), i.e. the length of the door; L the length of the pivoting side (m), i.e. the width of the door. Then the nominal area of the vent opening and the area of the hinged door is FN = bL. It is further assumed that
4
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
the vent cover mass is distributed uniformly over the surface of the cover, with surface density or inertia w (kg/m2 ). Let ϕ be the angle between the vent opening and the hinged door. It is assumed that the current venting area F(ϕ) is the gap area between the edges of the cover and the vent opening. The gap, as shown in Fig. 1, is formed from one rectangular region, based on the door edge opposite to the hinge and two triangular regions, based on the pivoting edges of the door. The venting area is then:
ϕ ϕ F (ϕ) = min FN , 2L sin b + L cos (8) 2 2 This area is zero for a closed vent (ϕ = 0) and is allowed to increase until ϕ = ϕN ; it reaches the maximum value equal to the nominal vent area FN . Also assume that for angles ϕ > ϕN , the venting area stays equal to FN . Furthermore, assume that the door is inelastically arrested at ϕ = 90◦ . 3.2. Pressure distribution for F(ϕ) ≤ FN When the vent is closed, the gas pressure is uniform throughout the door surface, and is equal to p(t), the pressure inside the enclosure (Pa). Furthermore, the gas mass discharge rate is zero. When the vent is open, the picture changes. First, the static pressure of the escaping gases on the door is smaller than the pressure at ‘stagnation’ conditions inside the enclosure. Second, the pressure is not uniform along the door surface any more. The first issue is addressed by using the Bernoulli or mass conservation law for the gas flowing between the enclosure inside the vent cross-section and the current venting area. At the vent cross-section, let p = p1 , and u = u1 (where u1 is the flow velocity through the vent cross-section (m/s)), and flows are low enough to warrant the assumption of incompressibility. The Bernoulli’s equation for these two levels allows p1 to be expressed as: p1 = p(t) −
ρ u21 , 2
Outside the vessel the gas pressure is p = pa , and the average gas velocity in the changing venting area is u = uL . From the Bernoulli’s equation relating the inside and the outside of the vessel, uL can be expressed as: 2(p(t) − pa ) 1/2 uL = (10) ρ The mass conservation law between the vent opening and the outside of the vessel gives an expression for u1 : u1 =
uL F (ϕ) bL
(11)
Substituting (10) in (11) and the results in (9) the pressure p1 in the vent cross-section becomes p1 = p(t) − (p(t) − pa )
F 2 (ϕ) (bL)2
(12)
This pressure depends on both the current explosion pressure and the current angle of the door opening. The second issue is dealt with by assuming that the pressure along the door surface changes as a linear function of the position on the width of the door: p(l, t) = p1 −
p1 − pa l L
(13)
Here l ≤ L is the current position, and p1 is defined by (12). The assumption of a linear pressure distribution along the width of an inertial hinged vent cover is a simplification of a more complex three-dimensional distribution. A simple modelling including a portion of an enclosure, hinged vent cover, and a portion of the surrounding environment, was performed using FluentTM . A constant internal pressure of 1.3 atm and external pressure of 1.0 atm were assumed. The problem was solved using a 2D approach, and the steady pressure distribution along the cover was determined for several cover opening angles. The results are shown in Fig. 2. Ultimately
(9)
where ρ is the gas density constant (kg/m3 ).
Fig. 1. Hinged door.
Fig. 2. Pressure distribution on inertial hinged vent cover as a function of cover opening angle for constant enclosure pressure. Steadystate problem, p/p0 = 1.3, 2D segregated solver, k–ε turbulence model, ♦, , —transients ‘below’ the cover, i.e. at cover surface hit by the escaping gases; ,,䊉—transients ‘above’ the cover.
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
a three-dimensional modelling, using CFD, will be necessary to accurately determine the pressure distribution on the hinged cover. As a first approximation, however, a linear distribution will force the model cover displacement to more closely match experimental experience than a model assuming the enclosure pressure everywhere on the cover.
5
these considerations in mind, it is possible to replace (16) and (17) with a single formula bL2 F 2 (ϕ) [p(t) − pa ] 1 − , 6fjet (ϕ) (bL)2
Tpressure =
(18)
where the transient factor fjet (ϕ) is 3.3. Pressure force torque for F(ϕ) ≤ FN
The torque exerted by the gas pressure on the door is (N m): L Tpressure =
[p(l, t) − pa ] bl dl
(14)
0
The pressure torque is always positive, since the gas flow momentum always works for opening the vent. Using formulae (12) and (13) by substitution in (14), the torque the gas applies to turn the door on its hinge can be written as: bL2 F 2 (ϕ) Tpressure = (15) [p(t) − pa ] 1 − 6 (bL)2 Our assumption of linear change of pressure in (13) results in very easy derivations. However, the true pressure distribution is not linear, the whole vent cover-pushing phenomenon being three-dimensional, non-stationary and dependent on the geometry of the cover and vent opening. Therefore, we decided to settle on the linear distribution augmented by an empirical ‘jet factor’ Cjet as follows: bL2 F 2 (ϕ) Tpressure = (16) [p(t) − pa ] 1 − 6Cjet (bL)2 The Cjet will compensate for the true non-linear, nonstationary and geometry-dependent character of the hinged door movement. Notice that when the door is closed, ϕ = 0, formula (15) gives three times smaller value than the correct torque for the closed door should be: Tpressure,closed =
bL2 [p(t) − pa ] 2
(17)
Therefore, we have to use different formulae for a closed or almost closed door and a sufficiently wide open door. When the door is shut or is opened within some small range of angles, we apply formula (17). Above a certain angle, when F(ϕ) reaches a certain fraction Ajet of the full area bL, the velocity of gases escaping along the door surface becomes significant, and we have to apply formula (16). To ensure the continuity of transition from the ‘pressure’ regime of formula (17) to the ‘jet’ regime of formula (16), we assume this transition to be linear with respect to F(ϕ) and controlled by the ‘jet fraction’ parameter Ajet . We also have to ensure that formula (16) never produces a value greater than formula (17) for the same pressure and cover dimensions. To that end, we restrict the values of the ‘jet factor’ by the condition Cjet > 1/3. With
1 Cjet F (ϕ) fjet (ϕ) = min max , , Cjet 3 bLAjet
(19)
Indeed, when ϕ = 0, we have F(ϕ) = 0, the transient factor fjet (ϕ) equals 1/3, and formula (18) assumes the form of (17). At a certain angle ϕ, when F(ϕ) = Ajet bL, fjet (ϕ) becomes equal to Cjet , and stays at this value for any greater angle of opening. Respectively, formula (18) assumes the form of (16). For all the interim angles, formula (19) produces interim values of fjet (ϕ) between 1/3 and Cjet . The values for the controlling parameters Ajet and Cjet have to be determined empirically. If comparison with an experiment shows that Ajet assumes relatively low values, say less than 10% of the nominal area, and Cjet turns out to be reasonable, say less than two, then formula (18) may be deemed a plausible approximation of the physical reality. 3.4. Balance of torques for F(ϕ) ≤ FN In rotational movement of our hinged door, the balance of torques (moments about the axis of rotation) can be expressed as follows: Tfull = Tpressure + Tgravity
(20)
where Tfull is the full torque, Tpressure is the torque generated by the pressure of escaping gases, formula (18), and Tgravity is the torque generated by the force of Earth’s gravity. The full torque (N m) is Tfull = Jm × α, where α is the angular acceleration of the vent (radian/s2 ) and Jm is the moment of inertia of the vent (kg m2 ). The general formula for Jm (kg m2 ) is Jm = m l2 dm, where dm is the elementary mass (kg) forming the body, and l is the distance of this elementary mass from the axis of rotation (hinge) (m). With the assumption of uniform mass distribution over the vent surface, the result is: dm = wdF , where w is the inertia of vent cover and dF is the area of a surface element. Then Jm = w l2 dF FN
Since the vent is rectangular, this integral can be rewritten as L Jm = wb
l2 dl = 0
wbL3 3
6
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
As the inertia w is expressed as w = m/(bL) where m = mass (kg), the moment of inertia simplifies to Jm =
mL2 3
The full torque is therefore Tfull =
αmL2 3
The torque Tgravity of the gravity force depends on where the vent is mounted and what side of the vent is hinged. If the vent is mounted on a wall, and is hinged at its side edge, such that its swinging motion is in the horizontal plane, then Tgravity = 0. If the vent is mounted on a wall, and is hinged at its top edge, then Tgravity = −
mgL sin(ϕ) , 2
where g is the acceleration gravity (m/s2 ). Tgravity is negative because for a top-hinged wall-mounted vent the gravity works against the opening vent. If the vent is bottom-hinged to a wall, then the difference is only in the sign Tgravity = +
mgL sin(ϕ) , 2
since gravity helps the vent to open in this case. Analogously, if the vent is mounted on the ceiling and opens upwards, then Tgravity = −
mgL cos(ϕ) , 2
and for vents mounted in the floor and opening downwards, Tgravity = +
mgL cos(ϕ) 2
For wall-mounted vents, gathering formula (18) and the above formulae for the torques in the Eq. (20) yields: αmL2 bL2 (p(t) − pa ) F 2 (ϕ) mgL sin(ϕ) = 1− − , 2 3 6fjet (ϕ) 2 (bL) (21) with a little change in notation for g. We assign the ‘direction’ for the gravitational force, such that g > 0 corresponds to the situation when the hinged edge is at the top; g < 0 when the hinged edge is at the bottom; and g = 0 when the hinged edge is at a side. In a similar way, for ceiling/floor-mounted vents, sin ϕ in formula (21) should be changed to cos ϕ, and it is assumed that g > 0 for ceiling-mounted vents that open upwards; g < 0 for floor-mounted vents that open downwards; and g = 0 for gravity compensated for by a spring or a balance. A simple algebra in Eq. (21) yields the following formula for the angular acceleration of the vent. For a wall-mounted
vent, 3.0 α= 2.0
b (p(t) − pa ) F 2 (ϕ) g 1− − sin(ϕ) , m3fjet (ϕ) L (bL)2
and for a ceiling-mounted vent, sin (ϕ) should be changed to cos (ϕ). The angular speed ω (radians/s) and the angle of opening ϕ (radians), are t ω(t) =
t α(τ)dτ
and
ϕ(t) =
0
ω(τ)dτ,
respectively.
0
Within each step of numerical integration of the governing equations, the angular acceleration is assumed constant. Therefore, uniformly accelerated motion formulae are used to obtain the new values ωnew and ϕnew of the angular speed and the angle of opening at the end of the integration step from the values at the beginning of this step: ωnew = ω + α dt and ϕnew = ϕ + ω dt + α dt2 /2.0. At the end of each integration step the following assignments take place: ω = ωnew ; ϕ = ϕnew . 3.5. Torque when F(ϕ) > FN The applicability of formula (18) is limited in the angle ϕ of the door opening, in that the formula will work only under the assumption that the current venting area F(ϕ) is less than the nominal area FN . At a certain angle ϕN such that F(ϕN ) = FN the pressure at the vent cross-section is equal to the atmospheric pressure, and the gas flow through the vent is unrestricted, as though there had been no vent cover at all. The CINDY code assumes, when this point is reached, that the vent cover displacement continues; yet the energy imparted to the cover through momentum earlier in the deflagration is affecting cover motion. Since the primary interest in the current research was in the influence of vent cover inertia while the cover position and displacement could influence transient enclosure pressures, detailed analyses of cover displacement once the vent was found to be 100% opened have been neglected. Further modelling would be required to assure an accurate representation of cover behaviour.
4. Comparison with experiments 4.1. Values of the empirical coefficients The empirical coefficients Cjet and Ajet in formula (19) were determined through matching of CINDY simulations of H¨ochst and Leuckel [17] experiments 3-B and 3-D (‘experiment 3-B’ or ‘3-B’ denotes experimental results plotted in Fig. 3b of their paper; ‘experiment 3-D’ or ‘3-D’ denotes experimental results plotted in Fig. 3d of their paper). Ajet for hinged covers represents a threshold below which jetting
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
7
Fig. 3. (a) Cjet —H¨ochst and Leuckel, experiment 3-B: opening angle, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (b) Cjet —H¨ochst and Leuckel, experiment 3-B: pressure, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (c) Cjet —H¨ochst and Leuckel, experiment 3-B: χ, Ajet = 0.05, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (d) Ajet —H¨ochst and Leuckel, experiment 3-B: opening angle, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (e) Ajet —H¨ochst and Leuckel, experiment 3-B: pressure, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open). (f) Ajet —H¨ochst and Leuckel, experiment 3-B: χ, Cjet = 1.4, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
flows of gases escaping through the opening vent are assumed to not yet be established. The variation of the vent cover displacement and the enclosure pressure for varying Cjet , assuming constant values of Ajet , are shown in Fig. 3a and b for displacement and pressure, respectively (in these and following figures simulated curves are shown until the moment of full burnout of the mixture inside enclosure). Once the values for the coefficients were selected, the values of χ and µ to backfit the simulations to the experimental data were found. Fig. 3c shows the backfitted values of χ that were found (results were obtained with a constant discharge coefficient µ = 1.2). As the figures show, increasing Cjet decreases the amount of force applied to the door—so the door takes longer to open. A slower opening door ‘generates’ less enclosure turbulence—thus the reductions in turbulence with increasing Cjet , as shown in Fig. 3c. The best fit seems to come with a value of Cjet of 1.4. A similar process, holding Cjet constant, was used to explore the impact of varying Ajet . Fig. 3d and e compares
computed and experimental displacements and pressures, respectively, for different values of Ajet , while Fig. 3f shows the resultant turbulence levels corresponding to the selected values of Ajet . For Ajet values at or below 0.01, the pressure is overestimated, yet the cover is not moving fast enough. Thus, values of Ajet at or below 0.01 should not be used. However, values of Ajet from 0.05 to 0.10 could be successfully used, although the turbulence rises as Ajet rises (Fig. 3f). As the figures show, increasing Ajet (in effect the amount of time the door is exposed to enclosure rather than door jetting conditions), ‘increases’ enclosure turbulence—the cover opens faster. The best fit seems to be achieved at an Ajet of 0.05. Note that given the data shown in Fig. 3, an estimation of the sensitivity of the enclosure turbulence to the Ajet and Cjet may be estimated. An average shift in χ as Ajet or Cjet varies may be determined by observing the shift in χ in the backfitted curves, and comparing the shift that occurs for various values of the coefficients. As a result, it was calculated that
8
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Fig. 4. (a). H¨ochst and Leuckel, experiments 3-B, 3-D: opening angle, µ = 1.2 ( –vent starts to open, 䊉–vent 100% open). (b). H¨ochst and Leuckel, experiments 3-B, 3-D: pressure, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
when Ajet shifts ±10%, χ shifts ±1.5%. Similarly, when Cjet shifts ±10%, χ shifts in opposite direction ±6%. From this, it would seem that the enclosure turbulence is more sensitive to the amount of force applied to the vent cover than to the amount of time that force is applied. The discrepancies between the experimental and simulation results in pressure dynamics, especially those in the peak areas, could be explained by heat losses to enclosure walls, which were not modelled in this study. From explosion safety engineering point of view it is acceptable as calculated pressure peaks are conservative relative to experimental data. 4.2. Validation 1: H¨ochst and Leuckel’s experiments The model has been validated against H¨ochst and Leuckel’s experiments [17]. For translation panels, a detailed description of that validation was presented in [14]. Their apparatus consisted of a 50 m3 silo of reinforced concrete with H/D = 4. The vent covers for these experiments were a pair of hinged vent covers arranged in a ‘butterfly’ configuration on the top surface of the silo. The edge of each cover most remote from the hinges was padded so that when the covers opened to 90◦ , the impact of the covers on each other was minimised. Experiment 3-B was a quiescent mixture of 10.7% methane–air, with a total venting area (for the two covers) F = 1.91 m2 , an inertia w = 124 kg/m2 , and torque of 532 N m. Experiment 3-D was a turbulent mixture (mixed by a fan within the enclosure for the purpose) of 10.6% methane–air, with a total venting area (for the two covers) F = 1.91 m2 , an inertia w = 73 kg/m2 , and torque of 314 N m. The mixtures were ignited by an electric match with energy of 75 J, located 3.5 m from the floor at the silo centre line. Sui was 0.38 m/s; Ei was 7.4; γ u was 1.39; γ b was 1.25;cui was 353 m/s; F was 2.45 m2 ; ε was 0.3; and the molecular mass (M) was 27.24 kg/kmol. Fig. 4a compares the calculated opening angle transients with the experimental results, while Fig. 4b compares the calculated enclosure pressures with the experimental results. A good match of both displacement and pressure has been achieved, using the Cjet and Ajet coefficient values of 1.4 and 0.05, respectively. As expected, the pre-mixed enclosure
gases in experiment 3-D results in faster opening of the cover, and higher and earlier pressures than those experienced by the quiescent mixture in experiment 3-B. These results are well matched by the CINDY computations. Fig. 5 shows how χ varied as the calculations progressed, and the best fit values of χ for experiments 3-B and 3-D. The enclosure gases were initially turbulent in experiment 3-D—thus the higher χ was needed to backfit the pressure and displacement are not surprising. In the CINDY code, χ is implemented as piecewise-linear; over some ranges of backfitted values, many small increments in χ can be replaced by a few longer segments with little or no detrimental change to the backfit for either pressure or displacement. Data from future experiments will be required to further determine if a two- or three- step χ curve established by simple rules could predict deflagration dynamics. 4.3. Validation 2: Zalosh’s experiment The model has also been validated against an experiment of Zalosh [18]. His apparatus consisted of a rectangular concrete bunker with interior dimensions H × W × L = 3.1 m × 2.0 m × 5.4 m and a volume of 33.5 m3 . The only relief was a pair of wall-mounted blowoff panels, arranged vertically in parallel, hinged on the bottom edges, opening outwards and downwards, arrested by cables in the 90◦ opened position. Each door was sized to be 1.29 m2 , for a total venting area of 2.58 m2 . An initially quiescent near stoichiometric 10.0% mixture was ignited by a 12 J spark
Fig. 5. H¨ochst and Leuckel, experiments 3-B, 3-D: χ, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10
Fig. 6. Zalosh: pressure, µ = 1.2 ( —vent starts to open, 䊉—vent 100% open).
9
• The obtained empirical coefficients Cjet and Ajet assume plausible values of 1.4 and 0.05. Deviation of the empirical parameters Cjet and Ajet from values derived herein results in only small variances in the backfitted values of the turbulence factor χ. • Observation of gradually changing turbulence factor χ over the processed experiments suggests that χ may be capable of being represented by a simple curve. The processing of additional experimental pressure and displacement test data would be useful to assess how predictable the gradual change of χ may be.
Acknowledgements The authors wish to thank Dr. Siegfried H¨ochst for providing additional information about the structure of the vent panel arresting mechanisms used in the experiments he performed with Dr. Wolfgang Leuckel. Funding by the Engineering and Physical Sciences Research Council (EPSRC) through Grant #GR/R53333/01 is gratefully acknowledged. Fig. 7. Zalosh: χ, µ = 1.2. ( —vent starts to open, 䊉—vent 100% open).
located in the centre of the bunker. In the calculations, the burning velocity, Su , was 0.38 m/s; the combustion products expansion coefficient at initial conditions, Ei , was 7.4; the unburnt mixture adiabatic exponent, γ u , was 1.39; the burnt gas adiabatic exponent, γ b , was 1.25; the speed of sound, cui , was 353 m/s; the overall thermokinetic exponent, which gives the dependence of burning velocity on pressure at adiabatic compression conditions, ε, was 0.3; and the molecular mass, M, was 27.56 kg/kmol. Coefficient values of Cjet = 1.4 and Ajet = 0.05, determined by comparisons with experiments by H¨ochst and Leuckel [17], were used. Fig. 6 compares the calculated pressure transients with the experimental results for hinged covers. No displacement data were published in Zalosh [18], so there is only a pressure comparison. The comparison to the experimentally observed pressure is reasonable. Fig. 7 shows how χ varied as the deflagration progressed, and the best fit values of χ for the experiment. The pressure rise at the end of the experiment seems to be due to turbulisation at the end of the process and a significant increase in χ is required to match the pressure.
5. Conclusions • Modelling of vented deflagrations with inertial venting devices has been performed for the case of multiple hinged doors. • The developed model and the CINDY code have been validated against large-scale experiments of H¨ochst and Leuckel [17] and Zalosh [18]. Good matches in both explosion overpressure transients and cover displacement transients have been achieved.
References [1] C. Yao, Explosion venting of low-strength equipment and structures, Loss Prev. 8 (1974) 1–9. [2] H.L. Pasman, Th.M. Groothuizen, H. de Gooijer, Design of pressure relief vents, in: C.H. Buschman (Ed.), Loss Prevention and Safety Promotion in the Process Industries, New York, 1974, pp. 185–189. [3] D. Bradley, A. Mitcheson, The venting of gaseous explosions in spherical vessels, Combust. Flame 32 (1978) 221–236, 237255. [4] V.V. Molkov, V.P. Nekrasov, Dynamics of gas combustion in a constant volume in the presence of exhaust, Combust. Explosion Shock Waves 17 (4) (1981) 363–369. [5] V.V. Molkov, Theoretical generalization of international experimental data on vented explosion dynamics, in: V.V. Molkov (Ed.), Proceedings of the First International Seminar on Fire and Explosion Hazards (17–21 July 1995, Moscow, Russia), All-Russian Research Institute for Fire Protection—Russian Association for Fire Safety Science, 1995, pp. 166–181. [6] D.M. Razus, U. Krause, Comparison of empirical and semi-empirical calculation methods for venting of gas explosions, Fire Saf. J. 36 (2001) 1–23. [7] T. Hirano, Combustion Science for Safety, in: J.H. Chen, M.D. Colket (Eds.), Proceedings of 29th International Symposium on Combustion, July 21–25, 2002, Sapporo, Japan), The Combustion Institute, Pittsburgh, 2002, pp. 167–180. [8] F. Williams, Approaches to the scaling of fires and explosions through mechanistic considerations, in: D. Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 445–456. [9] V.V. Molkov, A.Ya. Korolchenko, S.V. Alexandrov, Venting of deflagrations in buildings and equipment: universal correlation, in: Y. Hasemi (Ed.), Proceedings of the Fifth International Symposium on Fire Safety Science, Melbourne, Australia, March 3–7, 1997, pp. 1249–1260. [10] A. Di Benedetto, E. Salzano, G. Russo, Analysis of pressure piling by means of empirical correlation for vented vessels, in: D.
10
[11]
[12]
[13] [14]
V.V. Molkov et al. / Journal of Hazardous Materials A116 (2004) 1–10 Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 413–424. V.V. Molkov, Unified correlations for vent sizing of enclosures against gaseous deflagrations at atmospheric and elevated pressures, J. Loss Prev. Process Industries 14 (2001) 567–574. V.V. Molkov, Explosion safety engineering: NFPA 68 and improved vent sizing technology, in: Interflam 1999, Proceedings of the Eighth International Fire Science Conference, June 29–July 1, 1999, Edinburgh Conference Centre, Scotland, pp. 1129–1134. NFPA 68, Guide for Venting of Deflagrations, 2002 Edition, National Fire Protection Association, USA, 2002. V.V. Molkov, A.V. Grigorash, F. Tamanini, R. Dobashi, R.M. Eber, Vented gaseous deflagrations: modelling of translating iner-
[15] [16]
[17] [18]
tial vent covers, J. Loss Prev. Process Industries 16 (5) (2003) 395– 402. V.V. Molkov, Venting of Gaseous Deflagrations, in: DSc Dissertation, VNIIPO, Moscow, Russia, 1996 (in Russian). A. Grigorash, R. Eber, V. Molkov, A theoretical model of vented gaseous deflagrations in enclosures with inertial vent covers, in: D. Bradley, D. Drysdale, V. Molkov (Eds.), Proceedings of the Fourth International Seminar on Fire and Explosion Hazards, Londonderry, Northern Ireland, UK, September 8–12, 2003, pp. 445–456. S. H¨ochst, W. Leuckel, On the effect of venting large vessels with mass inert panels, J. Loss Prev. Process Industries 11 (1998) 89–97. R.G. Zalosh, Gas explosion tests in room-size vented enclosures, in: Proceedings of AIChE Loss Prevention Symposium, 1978, pp. 98–110.