Modelling the dispersion of flashing jets using CFD

Available online at www.sciencedirect.com

Journal of Hazardous Materials 154 (2008) 1198–1209

Modelling the dispersion of flashing jets using CFD R.K. Calay ∗ , A.E. Holdo AADE, University of Hertfordshire, Hatfield AL10-9AB, UK Received 11 March 2006; received in revised form 9 November 2007; accepted 9 November 2007 Available online 22 November 2007

Abstract Risk assessments related to industrial environments where gas is kept in liquid form under high pressure rely on the results from predictive tools. Computational Fluid Dynamics (CFD) is one such predictive tool and it is currently used for a range of applications. One of the most challenging application areas is the simulation of multiphase flows resulting from a breach or leakage in a pressurised pipeline or a vessel containing liquefied gas. The present paper deals with the modelling of the post-flashing scenario of a jet emanating from a circular orifice. In addition to being based on the equations governing fluid flow, the models used are those related to turbulence, droplet transport, evaporation, break-up and coalescence. Some of these models are semi-empirical and based on the data from applications other than flashing. However, these are the only models that are currently available in commercial codes and that would be used by consulting engineers for the type of modelling discussed above, namely the dispersion of a flashing release. A method for calculating inlet boundary conditions after flashing is also presented and issues related to such calculations are discussed. The results from a number of CFD based studies are compared with available experimental results. The results show that whilst a number of features of the experimental results can be reproduced by the CFD model, there are also a number of important shortcomings. The shortcomings are highlighted and discussed. Finally, an optimum approach to modelling of this type is suggested and methods to overcome modelling difficulties are proposed. © 2007 Elsevier B.V. All rights reserved. Keywords: CFD modelling; Flashing liquid jet; Dispersion

1. Introduction There are several types of industrial situations where fluid leakage could emanate from pipe flanges, valves or as a result of corrosion in transport pipe systems. In situations where the leak contains pressurised liquefied fluids which are highly toxic and flammable, such an accidental release may result in sudden and rapid explosions called flashing. The prediction of the magnitude of such explosions and spreading of flow constituents after the explosion is critical for the design and implementation of safety procedures. The mathematical models and predictive computer codes that may be applied to such releases to predict design parameters for safe storage, transport of pressurised fluids is of special interest for research. The flows resulting from the breach or a leakage are likely to be two-phase, consisting of liquid droplets and vapour with the droplets generated by flash boiling and mechanical fragmentation of the liquid. Thus the mathemat-

∗

Corresponding author. Tel.: +44 1707 281098; fax: +44 1707 285086. E-mail address: [email protected] (R.K. Calay).

0304-3894/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jhazmat.2007.11.053

ical models should incorporate an adequate understanding of the physical processes that are involved in the flashing region and subsequent spreading of a two-phase jet. In the case of spreading of a flashing jet there is continuous mixing between the ambient air and the jet and colliding and evaporation of liquid droplets and such mixing poses modelling challenges. There are various modelling approaches depending upon the different combinations of the constituents of the two-phase flow. Thus, if the constituents of the flashing liquid jet are known at the onset of leak, a suitable Computational Fluid Dynamics (CFD) model can be used to predict the dispersion of the jet. However, to be able to successfully apply CFD techniques, we require boundary conditions (such as velocity of the leaked liquid, vapour mass fraction, droplet size and distribution) at the onset. Mathematical models are not available that can successfully predict the inlet boundary conditions at the instant of leak for the CFD models, from the conditions which are known prior to the leak. Various theoretical approaches may be used to estimate conditions of flashing at the inlet boundary but there is no guidance regarding the applicability of different models to a wide range of storage conditions and different size of leaks that may occur creating

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

Nomenclature cp CD d d32 D g h H k L m M p q r Rt Re Ri T Tsh v vrd V We x x, y, z x* /d*

specific heat (J kg−1 K−1 ) drag coefficient droplet or particle diameter (␮m) Sauter mean diameter (␮m) nozzle diameter (m) gravity acceleration (m s−2 ) heat transfer coefficient (W m−2 K−1 ) enthalpy (J kg−1 ) thermal conductivity (W m−1 K−1 ) characteristic length (m) mass (kg) Mach number pressure (N m−2 ) heat flux (W m−2 ) droplet radii (m) droplets/gas mass ratio Reynolds number Richardson number temperature (K) super-heat degree (K) velocity (m s−1 ) droplets velocity ratio volume (m−3 ) Weber number distance (m) coordinates dimensionless displacement parameter

Greek symbol αt thermal diffusivity (m2 s−1 ) μ dynamic viscosity (Pa s−1 ) ρ density (kg m−3 ) σ surface tension (N m−2 ) Subscripts a air c critical g gas l liquid max maximum s super-heated sat saturation w water 0 initial ∞ ambient

jets with different characteristics. Unfortunately this is an area where high quality data for flashing release are rare because making measurements close to the flashing release is very difficult due to sharp gradients in time and space. Moreover, the very presence of sensors modifies the flow. Laser based techniques can overcome this shortcoming but acquiring data remains very difficult as the distance from the breach is decreased.

1199

The aim of the present study to review critically of the existing techniques for simulating accidental atmospheric ‘flashing’ jet releases and subsequent dispersion using commercially available CFD tools. In particular the models for flash atomisation and droplet dynamics downstream are assessed for their applicability to simulate mechanisms in dispersing flashing jet and proposals are made on how best these techniques can be used in realistic situations. Experimental data from the small and large-scale experiments using different fluids including propane, butane and R134A [FLIE; EU funded project] are used to enhance the understanding of the flashing process and to evaluate the validity of assumptions. The focus is to develop guidelines for using these tools more realistically to assess safety issues associated with accidental leaks of flashing liquids in industrial situations. 2. Theoretical background Flashing is violent break-up of the pressurised liquid into small droplets which occurs when a super-heated liquid comes in contact with ambient conditions (atmospheric temperature and pressure) at a point of a leak. The combination of hydrodynamic instabilities and thermal non-equilibrium conditions lead to flashing. The meta-stable liquid can only come to its equilibrium condition by releasing its super-heat through evaporation which consists of boiling and vapourising of the droplets as they disperse in ambient air and provide an explosive characteristic to the process. Fig. 1 shows a schematic representation of a flashing liquid jet issuing from a circular nozzle. The flashing location depends upon the geometry of the leak, initial conditions and the liquid properties. In Fig. 1 flashing is shown occurring outside the nozzle but flashing may begin inside the pipe or vessel. It has been observed experimentally that the ‘expansion’ region is made up of large droplets and liquid ligaments moving with increasing velocity. The velocity starts to decrease in the entrainment region due to mixing with ambient air as jet propagates. The axial temperature keeps decreasing well below the boiling temperature at ambient pressure further downstream and reaches its minimum as droplets evaporate.

Fig. 1. Schematic showing the expansion and entrainment region of the flashing jet and typical temperature and velocity variation along the jet axis.

1200

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

Beyond this point the temperature rises to ambient value. The droplet sizes and velocity mean values also decrease due to evaporation in the entrainment region. The radial variation of droplet size, temperature and velocity at different locations along the axis of the jet shows approximately Gaussian distribution. The foregoing flow behaviour shows that the flashing region is a complex problem and the transition between initial superheated liquid state and the following two-phase jet is not well understood. However, varying combinations of boiling and evaporation processes are present alongside severe mechanical break-up of droplets and turbulent mixing between the jet constituents and between these and the ambient air. Actual flashing phenomenon is usually characterised by the boiling temperature and the degree of super-heat (Tsh ), which is difference between actual liquid temperature and the saturation temperature at storage pressure: Tsh = Tl − Tsat |p

(1)

The super-heat level achieved in a fluid depends upon operating conditions and fluid properties as there is inter-dependence between these. It is important to know this dependence in order to estimate safe operational levels for industrial situations. 3. Modelling concepts As mentioned a flashing liquid jet consists of expansion and entrainment regions and the flow mechanisms are quite different within these regions. Models are available for two-phase flows which may be employed to model the dispersion of flashing jet. However, existing models do not involve all physical mechanisms that are present in the expansion region, i.e. where flashing takes place. It can be argued that it is not possible to include models to take all real physical mechanisms into account. But, it is necessary to achieve a complete understanding of what exactly happens in flashing processes in order to develop an appropriate model which is capable of dealing with the most important physical aspects and simple to be used in numerical simulations. Few computer models are available for simulating flashing jet dispersion (TRAUMA and ECJET by AEA Technology) that also calculate the behaviour of a flashing jet at point of release and Consequence Modelling Package PHAST by DnV, which incorporates a flash model [1] that calculates the depressurisation (external expansion; flashing) from the exit pressure to the ambient pressure. In the following section a review of methods that are used for the expansion and entrainment regions is detailed. 3.1. Expansion Most studies employ homogeneous single-phase and twophase models that assume both phases to move at the same velocity and homogeneous equilibrium methods allow for an equivalent single-fluid calculation. In contrast the constituents in a flashing jet are essentially two-phase and there are studies that provide a useful insight into the flow mechanisms in the expansion region. Fauske and Epstein [2] reported some prac-

tical calculation guidelines for two-phase flashing flows where the discharge rate is determined by three factors: the effect of sub-cooling, the effect of vapour–liquid equilibrium and nonequilibrium effect. An empirical factor in these equations is the discharge coefficient, Cd , which is less than 1.0 to account for viscosity and turbulence losses. The discharge coefficient is considered as a function of the Reynolds number and the diameter of the aperture. Wildgen and Straub [3] studied the flashing of super-heated water jets and reported that three types of boiling mechanics can affect the flashing process, this being particle boiling, surface nucleation and wall boiling. The incidence of these mechanisms depends on the pressure and temperature of the fluid before the leak. The particle boiling only depends on the pressure ratio of the jet and is initiated by activated nuclei; bubbles can grow to a size five times the diameter of the jet [3]. If there is no dissolved gas in the liquid, boiling activation is described by a pressure relation, known as the boiler state, which takes the state before the orifice into account: p0 p0 − psat = p psat − p∞

(2)

Dunbar et al. [4] modelled the flash evaporation of super-heated liquid assuming that the liquid flashed instantaneously into vapour. They applied Bernoulli’s equation to obtain mass flow rates for the sub-critical discharge and assumed the ideal gas relation to be valid for critical discharge or for choked flow. A simple adiabatic heat balance calculation was used to calculate the mass fraction of the liquid that flashes off. Vandroux-Koenig and Berthoud [5] assumed isentropic flow and homogenous equilibrium at the breach and calculated the maximum jet velocity (200 m s−1 ) based on the available energy from the liquid. However, this assumption yielded a velocity four times higher than the actual measured velocity. Aamir and Watkins [6] investigated the flashing problem in two parts. The first involved the construction of a thermodynamic model based on a quasi-state, separated flow analysis, and the second involved the study of the spray behaviour. They validated their result at 95, 500 and 1028 mm from the nozzle using the experimental data from Hervieu and Veneau [7] and Allen [8,9]. The problem considered was the flash evaporation of super-heated liquid and a simple adiabatic heat balance calculation could be used to calculate the mass fraction of the liquid that flashed off. Hallet [10] proposed continuous thermodynamic techniques in which the mixture composition is represented as a probability density function. However, the approach is only valid when the number of mixture components is large enough to approximate the species concentrations as varying continuously with the distribution variable rather than in discrete steps. Consequently, it is not suitable for a mixture with only a few components. Application of this technique to droplet evaporation allows the description of evolution of vapour composition in space and time as vapourisation progresses. The model involves fully transient calculations with property variations, and hence it is already numerically too complex and resource-intensive for including spray evaporation or combustion calculations.

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

It must be noted that the calculated velocity assuming isentropic choked conditions may be an order of magnitude higher than by applying Bernoulli equation. Information on droplet size and droplet size distribution is not obtained from these models. 3.2. Entrainment Fluid flow models based on conservation laws (i.e. Navier– Stokes equations) can also be used with multi phase fluid flow. Conservation equations are required for each phase and therefore various closure models are needed to account for the heat, mass and momentum transfer between phases, nucleation and droplet break-up model. At present these closure models consist of empirical coefficients based on local flow conditions. Some simplifications to the treatment of each phase can also be made depending on the dominating mechanisms. For spray modelling and droplet transport two models known as Locally Homogeneous Flow (LHF) and Separated Flow (SF) are available. In the LHF models the gas and liquid are assumed to be in dynamic and thermodynamic equilibrium, i.e. at each point in the flow both phases have the same velocity and temperature. In SF models the effects of finite rates of transport between the phases are considered. Commercial CFD packages incorporate these modelling concepts and therefore can be successfully applied to predict the dispersion of the multiphase jet following flashing provided the information on the temperature and velocity of liquid and vapour phases of the jet, droplet size and size distribution at the inlet boundary is available. The most common approach available in CFD codes is to solve the N–S equations for the continuum phase, i.e. carrier fluid and then solve for the trajectories of dispersed phase (gas bubbles or liquid droplets). The breaking of droplets or atomisation due to aerodynamic interaction processes exists in this region. These mechanisms are affected also by coalescence and turbulence. Mechanical break-up models exist in the literature such as classic break-up models Taylor Analogy Break-up (TAB), Reitz and Diwakar (RD) [11–13].

1201

4.1. Experimental data Available published data are now reviewed so that operating upstream conditions prior to accidental release may be correlated with the post-flashing conditions. Correlation is difficult because of the range of experimental set-ups and configurations, variations in nozzle size and shape and measurement regions, container conditions, scale factors, nozzle geometry, measurement techniques and type of fluid. Another difficulty with data is that often no reference is given on scaling criteria. In many instances, safety reasons prevented the use of the actual fluid and this resulted in the density ratio and Richardson number of the model fluid not being simultaneously realized. Change of fluid always changes buoyancy effects of the original fluid. The data have been correlated in terms of temperature, velocity, and droplet sizes variations along axial and radial direction. Relevant references are [8,9,11,12] and a comprehensive review on experimental work on flashing liquids can be found in the HSE report [14]. Some qualitative information of the flashing phenomenon has been gathered from visualisation data [14–16]. The minimum temperature distance (MTD) is an important parameter in characterising flashing. As the droplets break-up and expand the temperature drops due to energy absorbed in expansion and vapourisation until it reaches its minimum when there is a rapid rise in temperature that continues downstream though the rate of increase diminishes as it approaches the ambient temperature. Most data suggest that MTD for any release is between x/D = 150 and 170 but it must be noted that the MTD was measured at a distance more than 170 nozzle diameters with orifice sizes of 2 mm and smaller [14]. The degree of super-heat of the liquid and storage pressure influence the MTD but such information has not always been given [8,9]. As a result, it has not been possible to make a general correlation about the location of MTD. Finally, the full-scale release behaviour can also be affected by ambient atmospheric conditions, e.g. wind and humidity. We assess the accuracy of different theoretical models for estimating the input parameters required by the CFD models.

4. Flashing parameters CFD simulations start post-flashing at the exit plane of the orifice where liquid jet is broken into rapidly growing and exploding bubbles due to rapid evaporation occurring within the liquid. In most cases only the container’s pressure and temperature is known and so mass flow rate, velocity, turbulent characteristics, and phase volume fraction at the jet outlet need to be calculated. The phase and the flow rate of the discharge are dependent on the release scenario. Liquefied gas can be in any thermodynamic state as depicted by the enthalpy–temperature curve. Flow rate is also affected by the geometry of the orifice. Therefore in principle the laws of thermodynamics can be applied to predict the parameters in flashing region. However, the choice of theoretical model depends upon the dominating mechanisms in the flashing region which are not known. The proposed correlations in the literature for flash atomisation are also not clearly specified which results in a range of values for inlet boundary conditions based on the assumptions made.

4.1.1. Flashing fraction or vapour/liquid fraction The inlet boundary condition for two-phase requires the knowledge of mass fraction of liquid phase and vapour phase. Beginning of disintegration (flashing) of liquid into droplets depends on size of the leak and super-heat of the fluid. It can be assumed that jet undergo isentropic or isenthalpic expansion. Further adiabatic expansion can be assumed if the point of thermodynamic equilibrium at atmospheric pressure is close to the leak so that air entrainment is considered negligible. For isenthalpic depressurisation the mass fraction of vapour x, can be related to the specific heat, stagnation pressure and enthalpy of the fluid as follows: x=

cp Ts H

(3)

The influence of geometry is not taken into account in Eq. (3). Whereas, it has been observed that distance of disintegration of jet is inversely proportional to the size of nozzle [14]. But

1202

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

enough data were not available to correlate the influence of the geometry of the leak with the fraction vapour at the jet exit. 4.1.2. Exit velocity Several assumptions can be made to estimate exit velocity. When the flow is choked, the gas phase expands to ambient pressure within a downstream distance of about two orifice diameters [1]. This causes a strong acceleration of the two-phase flow and usually an increase of the jet diameter. Some experimental data suggest that the flow speed can be increased by a factor as high as ten in this region [1], and this has important consequences on the downstream dispersion. Therefore isentropic or isenthalpic expansions may be assumed and generally no entrainment can be assumed in this region, which is justified by the observed strong radial expansion. Maximum velocity at the exit for isentropic choked flow at the nozzle would be sonic velocity (i.e. Mach number, M = 1) and can be given by;  γRTsat vmax = (4) M Whereas maximum velocity can be calculated based on the maximum kinetic energy that can be obtained during liquid cooling or expansion:  cp vmax = Ts (5) Tsat Discharge of pure (i.e. non-flashing) liquids through a sharp edged orifice can be described by the classical work of Bernoulli. Assuming a pressurised liquid phase right at the leak, the exit velocity can be considered equivalent to the stagnation pressure:  2P v= (6) ρ The emerging jet is two-phase and it is not known if the droplets and gas velocities are equal. For homogeneous equilibrium assumption these velocities are assumed to be equal. Some CFD packages require as input the droplets/gas mass ratio, Rt , which may be obtained from the vapour fraction defined as follows:  droplets mass flow rate  Rt = (7) gas mass flow rate at the nozzle Vrd =

Vdroplets Vgas

(8)

openings, a two-phase jet develops immediately downstream of the exit, without any liquid core flowing out of the orifice. This indicates that for larger diameter orifices fragmentation due to flashing begins inside the vessel. Atomisation depends strongly upon the degree of super-heat and exit velocity which is dependent on initial pressure. However, existing break-up models are largely for mechanical break-up and typically estimate the Sauter mean diameter. Typically the maximum stable droplet size is predicted by a critical Weber number, We, which represents the ratio of inertia over surface tension forces: dmax =

σ Wemax v2 ρ g

(9)

Another empirical relationship between the dimensionless SMD and the jet Weber number and Reynolds number for large orifice, low-pressure releases is given by: SMD 0.64 = 82.23CD We−0.07 Re−0.5 d0

(10)

Weber number can be correlated with expansion energy instead of jet velocity to include flash atomisation into account. In any case the Weber number value has to be assumed which means calculated droplet diameter depends on the assumed Weber number value. Following empirical relationship to calculate mean droplet diameter after flashing depressurisation as a function of jet velocity and fluid properties [9]:  σ 1 d32 = 0.585 (11) vjet ρl Two established functions namely the Rosin–Rammler and log-normal distribution may be used to represent sprays from flashing releases which influence the prediction of downstream rainout. 5. Numerical modelling The test cases from large-scale and small-scale experiments using propane and R134A [FLIE project] have been numerically modelled. Since the dispersion of the flashing jet is modelled, it is reasonable to assume that flashing has begun at the leak and the droplets are already created. 5.1. Computational model

4.1.3. Atomisation Initial droplet size and size distribution is required as input in the CFD models for predicting droplet evaporation and/or rainout from the leak. The diameter of the droplet is expressed as mean diameter (MD) which is arithmetic mean diameter, volume median diameter (VMD) and mass median diameter (MMD), which is the volume median diameter drop size when measured in terms of volume (or mass). Other definitions that are used to express droplet size are Sauter mean diameter (SMD) expressed as d32 . Various data suggest that the droplet concentration is directly proportional to the nozzle size and the initial pressure. For large

Commercially available CFD software Star-CD was used for simulations presented in this paper. The Eulerian–Lagrangian approach in conjunction with implemented sub-models for evaporation, mechanical break-up and coalescence of the droplets was used for flow predictions of dispersion of the flashing liquid jet [13]. The computational domain and boundary conditions are shown in Fig. 2. The domain size was chosen as a function of the jet exit diameter being 600D and 100D in the axial and radial directions. Co-flow was only used for Case 1. The discretisation algorithm is based on the finite volume method. In these simula-

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

1203

Test cases were run for different boundary conditions for velocity, droplet size and size distribution from Table 1 to investigate how inlet parameters influence the prediction of the dispersion of the jet downstream. Comparison of the results with available experimental data was made to develop methodology for choosing models appropriately to estimate boundary conditions for similar cases.

Fig. 2. Schematic of computation domain.

tions same velocity was used at the inlet for the droplets and the continuous phase. Turbulence was modelled with k–ε model. 5.2. Boundary conditions With no guidance available to calculate inlet flow variables for the given upstream storage conditions, the inlet boundary conditions were estimated using the different theoretical models discussed earlier. Table 1 lists the results for propane stored at a pressure of 800 kPa and super-heat of 65 ◦ C prior to flashing. The calculated number of droplet is not directly put in; instead a number of computational droplets are introduced at the inlet to make up the total mass fraction flow rate of the liquid. There are no models for estimating droplet size due to thermodynamic atomisation so this size was estimated using models for mechanical atomisation and a droplet size class along with liquid phase mass flow rate is specified at the inlet boundary. Rosin–Rammler distribution function was used for the droplet size distribution   D n R = exp (12) Dm

5.2.1. Computational inlet The size of the jet exit influences the onset of flashing, which may start at the physical leak point or some distance away downstream. In the latter case liquid may expand by several diameters away from the jet exit before flashing starts. In the computational model it is assumed that flashing begins and liquid jet starts to disintegrate into droplets at the model inlet which may be several diameters larger in size than the actual jet exit. Statistical analysis of published experimental velocity data was used to model the initial expansion of the liquid jet. Velocity data at four axial positions 0–3, 3–75, 75–125 and 125–257 diameters from the jet exit were considered. “Edge data” represents the extreme data points in the radial direction at four axial locations. Edge 5% represents the radial distances at the four locations where measured velocity is 5% of the maximum axial velocity. Similarly Edge 1% represents the distances where velocity is 1% of the maximum velocity. The results are given in Fig. 3. Thus the diameter of the computational inlet can be calculated as θ Dmodel = 2nDactual tan + Dactual 2

(13)

where n is empirical index taken as 0–3 and θ as 55–60◦ . 5.3. Test cases 5.3.1. Case 1 This case represents the leakage of pressurised propane at 800 kPa and super-heat of ∼65 K to atmospheric conditions through a 2 mm diameter orifice. The annular region surround-

Table 1 Comparison of inlet parameters calculated for propane at 800 kPa and 65 ◦ C Parameter

Assumptions

Vapour mass fraction, x

Adiabatic expansion

Velocity (m s−1 )

1. Due to stagnation pressure 2. Due to maximum kinetic energy during liquid cooling 3. For isentropic choked flow

Turbulence intensity (%) Turbulence length scale

Default Same as model inlet diameter

Droplet diameter (␮m)

Atomisation

Equation

Vmax =

0.17

c p



Vmax =

52 Ts

176

γRTsat M

249

Tsat

10 σWemax V 2 ρg

Maximum, dmax Mean d32

Calculated value

cp Ts x= H  2P V = ρ

1 0.585 vjet

10–500



σ ρl

46

1204

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209 Table 3 Test Case 2 parameters Variable

Model inlet (m s−1 )

Inlet velocity Temperature, Tboil (K) Flow rate (kg s−1 ) Vapour fraction Droplet size (␮m) Turbulence intensity (%) Length scale (m)

34 246 0.025 0.20 200, 100, 50, 10 10 0.009

Fig. 3. Spreading of jet.

ing the jet inlet is the co-flow region. The boundary conditions for the co-flow region are specified at ambient conditions to flow smoothly around the jet in the direction of jet flow. Observations confirmed that the flashing started in region and continued some distance away from the orifice. Inlet temperature was taken as the boiling temperature at atmospheric pressure. Thus the model inlet diameter can be assumed at anywhere between 0 and 3D away from the real orifice assuming the initial expansion angle as 55–60◦ . In this case inlet diameter can be taken as 2–9 mm. Simulations were performed for models of three diameter values, 2, 4.5 and 9 mm. Inlet boundary conditions at the computational nozzle are calculated applying different assumptions discussed in the previous section. A series of simulations were made with these inlet values for comparison with the experimental data. Table 2 shows the inlet parameters chosen for the Case 1. 5.3.2. Case 2 This case is based on from the experiments performed with R134A with nozzles of diameter size 1 and 2 mm. These tests investigated the near field flashing phenomenon. The super-heat and back-pressure for these experiments ranged between 40 and 50 ◦ C and 700 and 900 kPa. Flow visualisation showed that fluid comes out as liquid and flashing starts about 3D away from the nozzle. Thus model inlet size was calculated at 3D away from the actual nozzle. Table 3 shows the inlet conditions used for Case 2. 5.3.3. Mesh sensitivity tests To check the grid independence, simulation of the free jet with no droplets was run with three different mesh resolutions: 1.8 × 105 (Mesh 1), 3.6 × 105 (Mesh 2) and 8 × 105 (Mesh 3).

Fig. 4. Axial distribution of non-dimensional velocity for three mesh sizes.

For this work the inlet velocity was assumed to be equivalent to stagnation pressure and calculated as 52 m s−1 . The axial velocity and temperature profiles were compared for the three mesh sizes. Fig. 4 shows the axial velocity profiles for three meshes and Mesh 2 (3.6 × 105 cells) was found to provide results which were not significantly different to those using the finer mesh (Mesh 3), consequently Mesh 2 was used. 6. Results Parametric analysis was performed for a range of inlet conditions to investigate the sensitivity to inlet conditions. Selection of results and comparison with experimental conditions is presented to highlight the difficulties in choosing correct inlet conditions. 6.1. Case 1 Numerical results were obtained for velocity, temperature and concentration profiles. The velocity data were not available

Table 2 Test Case 1 parameters Variable

Model inlet

Inlet velocity (m s−1 ) Temperature, Tboil (K) Flow rate (liquid propane) (kg s−1 ) Vapour fraction Maximum droplet size (␮m) Turbulence intensity (%) Length scale (m)

52 ∼233 0.04 0.17 100, 50, 10 10 0.005

Co-flow (air) 176 ∼233 0.04 0.17 100, 50, 10 10 0.005

1 298 – – – 1 0.1

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

Fig. 5. Velocity variations in the axial direction for two different inlet velocity values along the jet axis.

in the flashing region. Fig. 5 shows the velocity profile for the two values of inlet velocity. The higher value was calculated assuming the jet exit velocity to be based on maximum kinetic energy from cooling of liquid at flashing. Whereas the lower value was based on the assumption that leaked fluid comes out at a velocity based on the storage pressure. Both velocity profiles show a similar decreasing trend. Experimental data in this case are available for distance between 0.04 and 0.4 m from the orifice and is plotted on the same graph. The predicted results in this case only indicate that the measured value is close to the inlet velocity which was calculated from the stagnation pressure of the liquid. This assumption is applicable when the surface area of jet issuing from the leak point is small relatively to the volume of liquid contained within the jet prior to flashing. In most cases for small leaks flashing initially starts at the corners of the jet whereas the core of the jet propagates as a liquid and this assumption is realistic for these cases. The initial or inlet velocity value also has significant effect on temperature and concentration predictions downstream. Therefore predicted temperature profiles may also be compared using two values of inlet velocity in order to ensure the correct assumption for estimating velocity at the model inlet. Initial increase in velocity in the expansion region was not modelled in this case. Also no experimental data were available for comparison. Fig. 6 shows the axial temperature profile for the two simulations. The profile using the higher inlet velocity is very different

Fig. 6. Temperature profiles inn the axial direction for two inlet velocity values.

1205

to that using the lower value. The later profile exhibits a more realistic variation and shows a decrease due to energy absorbed by the droplets during evaporation. Further downstream (>1.5 m) the temperature rises to the ambient value. Initially, the rise is very sharp due to the evaporation. Using the inlet velocity of 176 m s−1 , the evaporation rate of droplets is not great, so no drop in temperature is seen initially and a slow decrease due to evaporation becomes apparent later. The initial rapid temperature increase does not correspond to reality. Experimental data for temperature were not available for the comparison but the predicted temperature profile for a lower inlet velocity is similar to other experimental observations. The predicted minimum temperature distance corresponds to the limit of droplet existence [8,9]. The temperature distribution is also sensitive to the prescribed initial droplet size and distribution at the inlet boundary. Fig. 7a and b shows results for simulations with two different droplet sizes and different exponents (n) in the Rosin–Rammler distribution function. The predictions are as expected for different size droplets. The influence of different exponent in the Rosin–Rammler distribution function was less significant for both droplet sizes.

Fig. 7. Axial temperature profile for initial droplet distribution based on Rosin–Rammler function with n = 2, 4, 6, 8 and Ui = 52 m s−1 ; (a) dm = 50 ␮m and (b) dm = 100 ␮m.

1206

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

Fig. 8. Droplet size variations in the axial direction of the jet.

The temperature profile for an initial droplet size of 50 ␮m is closer to experimental observations. In the present simulation MTD was predicted at 200–250 diameters away from the nozzle (i.e. x = 0.4 m for n = 8 and x = 0.5 m for n = 2). Compared to the experimental observation of 150–170D given in Section 4, the MTD is over predicted. There are other factors, i.e. atmospheric conditions and droplet break-up and transport that influence the position of MTD. The CFD model does not take all these factors into account. However, the uncertainties in the imposed inlet boundary conditions, i.e. initial droplet size, exit velocity and vapour mass fraction have a more significant effect on the accuracy of the predictions compared to these factors. In the test cases droplet sizes 10, 50 and 100 ␮m were used at the inlet. Fig. 8 shows the predicted droplet size distribution for the initial droplet size of 100 ␮m. The CFD model considers the breaking and evaporation of the droplets as jet propagates and mixes with the ambient air thereby causing droplet size decrease in the direction of the jet. However, the droplet size distribution depends on the inlet droplet size that is imposed as a boundary condition and theoretical models are not available to predict inlet droplet size due to flashing occurring at the inlet. It is suggested that for CFD models a maximum droplet size at the inlet based on mechanical atomisation alone may be used. The empirical model (Eq. (11)) for calculating SMD d32 at the model inlet gives a value of 43 ␮m which shows good agreement to the measured droplet size of 50 ␮m at the same position. However, Eq. (11) is empirical and since its range of validity is not quoted, it is not advisable to use this model in all cases.

Fig. 9. Measured count and mass distributions for the droplet diameter at x = 14D and 18D (Ref. FLIE 2004, VKI data).

Fig. 10. Axial velocity profiles for two droplet sizes at the inlet.

Fig. 10 shows the comparison between the CFD predicted axial velocity profile for two droplet sizes (100 and 200 ␮m) at the inlet with the same velocity for both cases. There is little difference between the axial velocity profiles. However, there is a clear difference in the temperature profiles for the two cases as seen in Fig. 11. The trend of the experimental data is different to the predicted results for both droplet sizes at the inlet. Particularly for the droplets of size

6.2. Case 2 In this case experimental measurements are available in the near field of the jet allowing comparison with theoretically estimated inlet conditions (Table 3). Comparisons are made at distance 14 mm (x = 14D) and 18 mm (x = 18D), where the droplet size distribution and velocity data are available. The measured droplet sizes vary from 50 to 750 ␮m at x = 14D with more than 50% of droplets lying between 50 and 200 ␮m (Fig. 9). The droplets of size 200 ␮m also represent 90% of the mass flow rate.

Fig. 11. Predicted axial temperature profile for two droplet sizes at the inlet.

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

1207

Fig. 12. Velocity–count distribution at distance x = 187D.

100 ␮m the predicted minimum temperature is much lower than the measured temperature and reaches its minimum at a distance of 0.2–0.25 m from the model inlet. The temperature does not drop to the same minimum level for the 200 ␮m size droplets. On the other hand experimental data show a more gentle temperature drop, which continues to decrease over a longer distance of 0.5 m. However, eventually as the jet mixes with the ambient air, its temperature begins to slowly increase to acquire the ambient temperature. Since measurements were taken only up to 0.55 m downstream of the nozzle, the minimum temperature distance was not reached in this case, which seemed to be more than 170D of the actual nozzle (2 mm). However, since the flashing starts few diameters away from the actual nozzle, we may consider the diameter of the model inlet (9 mm) as a reference for MTD. Experimental measurements are available for the velocity– count distribution at the distance of 187 and 507 diameters from the (real) nozzle. The CFD model predicts velocity of the continuous phase. The velocity of dispersed phase, i.e. droplets may be expressed statistically. Fig. 12 shows the comparison between experimental and predicted the velocity–count distribution of droplets at a distance of 187D, i.e. 187 mm from the nozzle. The predicted velocity for the continuous phase at this location is ∼24 m s−1 , which is in excellent agreement with the measured velocity of 24 m s−1 at this distance. Whilst the droplet velocity (count) as predicted from CFD model shows good qualitative similarity to the measured distribution, but the predicted distribution exhibits much narrower range (20–30 m s−1 ) than the measured range (2–24 m s−1 ). It is observed from Fig. 12 that more than 80% of the droplets are moving at 20 m s−1 velocity whereas the predicted droplet velocity for the same count is significantly higher ∼28 m s−1 . It seems that the CFD model does not predict the slowing down of the droplets in a manner similar to reality and there could be several reasons for this. In the CFD model the same velocity was used at the inlet for the droplet and the continuous phase velocity as in Homogeneous Model. Further downstream momentum transfer between the droplet and continuous phase is modelled using drag models which often require empirical constants that are case dependent. Thus, insufficient momentum transfer between the droplet and the continuous phase, insufficient mixing or entrainment of the jet or drag model not

Fig. 13. Comparison of predicted droplet size distribution at x/D = 187 and 507.

affecting the droplet momentum may yield droplet velocities different to the experimental measurements. Also taking the difference in predicted temperature into account, it seems that evaporation rate or heat and mass transfer is not sufficiently well modelled. A similar difference is also evident in Fig. 13 at the two locations, i.e. x = 187D and 507D. The measured droplet count shows an exponential decrease in the percent count with respect to the droplet size. The measurements show a large count of smaller size droplets (<100 ␮m) and show almost an even distribution of droplets of sizes 100–300 ␮m at both locations. As the droplets move downstream the size gets smaller due to break-up and evaporation. The CFD results do not show a similar pattern of decrease in the size. The size distributions at both locations show an essentially uniform distribution of droplet sizes. There is not much change in the distribution profile, which is still more uniform than measured profile. It must be pointed out that the experimental measurements were the point measurements whereas the predicted values represent the average of all droplet size counts at the cross-sectional plane at 187D and 507D. This results in a more uniform droplet size distribution profile so it is not possible to make a conclusive comparison between the droplet size distribution data. 7. Discussion This study has focussed on using commercially available CFD tools to predict the dispersion of flashing jets due to release of liquefied gases in real situations. The general trend in the CFD results used are as expected because these transport models are well validated and verified for fluid flow. However, the accuracy of predictions depends upon the accuracy of the boundary conditions used in the models. The inlet boundary condition of the CFD model is an area where uncertainties exist regarding what happens during flashing. As has been seen the use of various assumptions leads to the use of different theoretical models that yield different values for inlet variables. With no guidance for the choice of correct models for different leak scenarios, in this paper we have investigated

1208

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

the use of various assumptions to estimate inlet conditions; then CFD predictions downstream of the inlet are compared with the experimental data in order to establish the range of applicability of the assumptions. However, more data from a range of scenarios are required so that further correlations can be made with the CFD tools.

Further investigations are needed to extend knowledge considering the influence of various parameters for modelling of droplets-laden flows. However, to take all these factors into consideration might result in expensive simulations and the choice should be made for the models that represent important, real situations with an accuracy dictated by the application and economics.

7.1. Fluid flow model and droplet transport 7.3. Inlet parameters/boundary conditions The existing approach has reproduced trends and main features of the observed behaviour, although the implemented sub-models have been derived from physical phenomena found in different industrial fields. Some of the models involve empirical constants and are only valid for specific range of parameters and the predictions may show uncertainties if used beyond the range of applicability. Both experimental results show that temperature drops gently to a minimum value and then remains constant for a longer distance when compared to the predicted temperature profiles which shows a sharper decay followed by a rise. The evaporation models available in CFD software assume saturated temperature of the liquid at atmospheric pressure as the minimum temperature reached by the droplet of certain size and further cooling is not taken into account. Evaporation of the droplets also depends upon the droplet size which in turn is affected by the break-up models used in the simulations. Droplet transport models should take account of these mechanisms. The results predicted droplet size spectrum at different crosssections in the axial direction did not show good agreement with the experimental data. Considering the similar discrepancy observed in the temperature profile, it seems that the CFD model predicts evaporation rate and rate of breaking up of larger droplets into smaller droplets different to that found experimentally. This discrepancy is reflected in the predicted droplet velocity. 7.2. Turbulence The turbulent characteristics of the flashing jet affect mixing and entrainment with the surrounding air in the same manner as for single-phase jet. The effect of dilution with air enhances evaporation thereby cooling the entrained air. The correct modelling of turbulence, or the effects of turbulence, is critical in predicting the dispersion characteristics of the flashing jet. The existing models allow the influence of turbulent velocity fluctuations on droplet movement but modification of turbulence due to the effect of the presence of a secondary phase (i.e. droplets, bubbles or particles) is not considered, because the k–ε model is developed on the basis of transport of single fluid. It is assumed that 10 ␮m droplets do not modify the flow because of the low value of droplet’s Stokes’ number [5] therefore generally this model is used as the default model in the CFD tools. Modifications are available where the moving secondary phase is considered as the moving obstacles in the flow, but the effects generated by mass exchange at the interface of droplets and main flow due to the relative order of the size of droplets compared to the order of turbulent eddies is not known.

The accurate estimation of flow conditions at the flashing plane significantly influences the predictions critical for the risk and safety implications of a leak scenario. In this study vapour mass fraction was estimated assuming adiabatic expansion to the atmosphere taking only super-heat into account. It has shown that this assumption provides a good estimate for vapour fractions for an orifice diameter of 1–20 mm, a pressure of 700–1000 kPa and a super-heat of ∼40–65 ◦ C. Existing correlations to predict the droplet size do not consider orifice length-to-diameter ratio L/d0 , surface roughness, downstream distance and the influence of liquid properties. 8. Conclusions In the present study comparison between experimental observations and the CFD results suggest that at present initial spray conditions are critical. Thus, when using CFD modelling in risk assessment studies the boundary conditions and sub-models used should be chosen in such a way that the worst possible scenarios are captured and risk associated variables are overpredicted. Following guidelines are thus proposed: • Inlet conditions can be estimated from the storage conditions, i.e. the super-heat of the liquid and pressure which are the only properties known prior to the leak. For a breach size of 1–20 mm diameter, a pressure of 700–1000 kPa and a superheat of ∼40–65 ◦ C the following assumptions can safely be applied to estimate inlet conditions for CFD dispersion model. ◦ The flashing normally would start closer to the breach as the size of the breach increases. The computational inlet diameter can be fixed at the distance of 3D from the point of leak for the size diameter D < 3 mm. ◦ For breach size of 4 mm < D < 20 mm, the computational inlet diameter can be estimated at the real leak point. For bigger breaches flashing may start even inside the vessel. ◦ Vapour mass fraction can be calculated using adiabatic expansion to the atmospheric conditions. ◦ Inlet velocity can be calculated based on liquid stagnation pressure and same velocity can be imposed for both phases at the inlet. • When both phases have the same initial velocities large droplets at inlet spread over much larger distances than the smaller droplets. Existing empirical correlations for estimating droplet size are based only on mechanical atomisation. Few simulations can be sensibly made with different droplet

R.K. Calay, A.E. Holdo / Journal of Hazardous Materials 154 (2008) 1198–1209

sizes based on a range of Weber number to predict the worst case scenario. • The droplet size distribution at the inlet does not affect the flow predictions. More experimental work in order to formulate and validate a modelling strategy for a correct choice of the existing submodels is thus required. Further numerical simulations together with experiments should be undertaken to modify existing theoretical models and also to develop better correlations to estimate the model inlet parameters. Acknowledgements The work was under taken within the scope of collaborative project FLIE (2004) funded by the EU under FP5 [EU Contract no. EVG1-CT-2000-00025]. The experimental work referred to was carried out at INERIS, France and VKI, Belgium. Support from Dr. J.B. Davis is also acknowledged in the preparation of the manuscript. References [1] C.J. Wheatley, A theoretical study of NH3 concentration in moist air arising from accidental releases of liquefied NH3 , using the computer code TRAUMA, Health and Safety Executive, SRD/HSE/R 393, 1987. [2] H.K. Fauske, M. Epstein, Source term considerations in connection with chemical accidents and vapor cloud modeling, in: Proceedings of the International Conference on Vapor Cloud Modeling, Cambridge, MA, November 2–4, AIChE, New York, 1987, p. 251.

1209

[3] A. Wildgen, J. Straub, The boiling mechanism in super-heated free jets, Int. J. Multiphase Flow 15 (2) (1989) 193–207. [4] C.A. Dunbar, A.P. Watkins, J.F. Miller, Theoretical investigation of the spray from a pressurized metered-dose inhaler, Atomization Sprays 7 (1997) 417–436. [5] S. Vandroux-Koenig, G. Berthoud, Modelling of a two phase momentum jet close to the breach, in the containment of liquefied gas, J. Loss Prevent. Process Ind. 10 (1997) 17–29. [6] M.A. Aamir, A.P. Watkins, Numerical analysis of depressurisation of highly pressurised liquid propane, Int. J. Heat Fluid Flow 21 (2000) 420–431. [7] H. Hervieu, T. Veneau, Experimental determination of the droplet size and velocity distributions at the exit of the bottom discharge pipe of a liquefied propane storage tank during a sudden, J. Loss Prevent. Process Ind. 9 (6) (1996) 413–425. [8] J.T. Allen, Laser-based measurements in two-phase flashing propane jets. Part one: velocity profile, J. Loss Prevent. Process Ind. 11 (1998) 291– 297. [9] J.T. Allen, Laser-based measurements in two-phase flashing propane jets. Part two: droplet size distribution, J. Loss Prevent. Process Ind. 11 (1998) 299–306. [10] W.L.H. Hallett, A simple model for the vaporization of droplet with a large number of components, Combust. Flame 50 (121) (2000) 334–344. [11] R. Brown, L. York, Sprays formed by flashing jets, Am. Inst. Chem. Eng. 8 (2) (1962) 149–153. [12] R.D. Reitz, Modelling atomisation processes in high-pressure vaporizing sprays, Atomisation Spray Technol. 3 (1987) 309–337. [13] Star-CD Manual CD-Adapco. [14] H.M. Witlox, P.J. Bowen, Flashing Liquid Jets and Two-phase Dispersion—A review, Health and Safety Executive and DetNorske Veritas, HSE Books, 2003. ISBN0 0717622509. [15] D. Yildiz, P. Rambaud, J. van Beek, Break-up droplet size and velocity characterisation of a two-phase flashing R134A jet, in: 5th International Conference on Multiphase Flow, Yokohama, Japan, 30 May–4 June, 2004. [16] G.N. Abramowich, Theory of Turbulent Jets, MIT Press, 1963.